跳转至

Berlekamp–Massey 算法

Berlekamp–Massey 算法是一种用于求数列的最短递推式的算法。给定一个长为 的数列,如果它的最短递推式的阶数为 ,则 Berlekamp–Massey 算法能够在 时间内求出数列的每个前缀的最短递推式。最坏情况下 ,因此算法的最坏复杂度为

定义

定义一个数列 的递推式为满足下式的序列

其中 称为该递推式的 阶数

数列 的最短递推式即为阶数最小的递推式。

做法

与上面定义的稍有不同,这里定义一个新的递推系数 ,满足:

容易看出 ,并且阶数 与之前的定义是相同的。

我们可以增量地求递推式,按顺序考虑 的每一位,并在递推结果出现错误时对递推系数 进行调整。方便起见,以下将前 位的最短递推式记为

显然初始时有 。假设递推系数 对数列 的前 项均成立,这时对第 项就有两种情况:

  1. 递推系数对 也成立,这时不需要进行任何调整,直接令 即可。
  2. 递推系数对 不成立,这时需要对 进行调整,得到新的

,即 的递推结果的差值。

如果这是第一次对递推系数进行修改,则说明 是序列中的第一个非零项。这时直接令 即可,显然这是一个合法的最短递推式。

否则设上一次对递推系数进行修改时,已考虑的 的项数为 。如果存在一个序列 ,满足:

并且 ,那么不难发现将 按位分别相加之后即可得到一个合法的递推系数

考虑如何构造 。一种可行的构造方案是令

其中前面一共有 ,且最后的 表示将 每项乘以 后接在序列后面。

不难验证此时 ,因此这样构造出的是一个合法的 。将 赋值为 逐项相加后的结果即可。

如果要求的是符合最开始定义的递推式 ,则将 全部取相反数后在最开始插入 即可。

从上述算法流程中可以看出,如果数列的最短递推式的阶数为 ,则算法的复杂度为 。最坏情况下 ,因此算法的最坏复杂度为

在实现算法时,由于每次调整递推系数时都只需要用到上次调整时的递推系数 ,因此如果只需要求整个数列的最短递推式,可以只存储当前递推系数和上次调整时的递推系数,空间复杂度为

参考实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
vector<int> berlekamp_massey(const vector<int> &a) {
  vector<int> v, last;  // v is the answer, 0-based, p is the module
  int k = -1, delta = 0;

  for (int i = 0; i < (int)a.size(); i++) {
    int tmp = 0;
    for (int j = 0; j < (int)v.size(); j++)
      tmp = (tmp + (long long)a[i - j - 1] * v[j]) % p;

    if (a[i] == tmp) continue;

    if (k < 0) {
      k = i;
      delta = (a[i] - tmp + p) % p;
      v = vector<int>(i + 1);

      continue;
    }

    vector<int> u = v;
    int val = (long long)(a[i] - tmp + p) * power(delta, p - 2) % p;

    if (v.size() < last.size() + i - k) v.resize(last.size() + i - k);

    (v[i - k - 1] += val) %= p;

    for (int j = 0; j < (int)last.size(); j++) {
      v[i - k + j] = (v[i - k + j] - (long long)val * last[j]) % p;
      if (v[i - k + j] < 0) v[i - k + j] += p;
    }

    if ((int)u.size() - i < (int)last.size() - k) {
      last = u;
      k = i;
      delta = a[i] - tmp;
      if (delta < 0) delta += p;
    }
  }

  for (auto &x : v) x = (p - x) % p;
  v.insert(v.begin(), 1);

  return v;  // $\forall i, \sum_{j = 0} ^ m a_{i - j} v_j = 0$
}

朴素的 Berlekamp–Massey 算法求解的是有限项数列的最短递推式。如果待求递推式的序列有无限项,但已知最短递推式的阶数上界,则只需取出序列的前 项即可求出整个序列的最短递推式。(证明略)

应用

由于 Berlekamp–Massey 算法的数值稳定性比较差,在处理实数问题时一般很少使用。为了叙述方便,以下均假定在某个质数 的剩余系下进行运算。

求向量列或矩阵列的最短递推式

如果要求向量列 的最短递推式,设向量的维数为 ,我们可以随机一个 维行向量 ,并计算标量序列 的最短递推式。由 Schwartz–Zippel 引理,二者的最短递推式有至少 的概率相同。

求矩阵列 的最短递推式也是类似的,设矩阵的大小为 ,则只需随机一个 的行向量 和一个 的列向量 ,并计算标量序列 的最短递推式即可。由 Schwartz–Zippel 引理可以类似地得到二者相同的概率至少为

优化矩阵快速幂

是一个 维列向量,并且转移满足 ,则可以发现 是一个不超过 阶的线性递推向量列。(证明略)

我们可以直接暴力求出 ,然后用前面提到的做法求出 的最短递推式,再调用 常系数齐次线性递推 即可。

如果要求的向量是 ,则算法的复杂度是 。如果 是一个只有 个非零项的稀疏矩阵,则复杂度可以降为 。但由于算法至少需要 的时间预处理,因此在压力不大的情况下也可以使用 的线性递推算法,复杂度同样是可以接受的。

求矩阵的最小多项式

方阵 的最小多项式是次数最小的并且满足 的多项式

实际上最小多项式就是 的最小递推式,所以直接调用 Berlekamp–Massey 算法就可以了。如果 是一个 阶方阵,则显然最小多项式的次数不超过

瓶颈在于求出 ,因为如果直接每次做矩阵乘法的话复杂度会达到 。但考虑到求矩阵列的最短递推式时实际上求的是 的最短递推式,因此我们只要求出 就行了。

假设 个非零项,则复杂度为

求稀疏矩阵行列式

如果能求出方阵 的特征多项式,则常数项乘上 就是行列式。但是最小多项式不一定就是特征多项式。

实际上如果把 乘上一个随机对角阵 ,则 的最小多项式有至少 的概率就是特征多项式。最后再除掉 就行了。

阶方阵,且有 个非零项,则复杂度为

求稀疏矩阵的秩

是一个 的矩阵,首先随机一个 的对角阵 和一个 的对角阵 , 然后计算 的最小多项式即可。

实际上不用调用矩阵乘法,因为求最小多项式时要用 乘一个向量,所以我们依次把这几个矩阵乘到向量里就行了。答案就是最小多项式除掉所有 因子后剩下的次数。

个非零项,且 ,则复杂度为

解稀疏方程组

问题:已知 , 其中 是一个 满秩 稀疏矩阵, 的列向量。 已知,需要在低于 的复杂度内解出

做法:显然 。如果我们能求出 () 的最小递推式 (), 那么就有结论

(证明略)

因为 是稀疏矩阵,直接按定义递推出 即可。

同样地,设 中有 个非零项,则复杂度为

参考实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
vector<int> solve_sparse_equations(const vector<tuple<int, int, int>> &A,
                                   const vector<int> &b) {
  int n = (int)b.size();  // 0-based

  vector<vector<int>> f({b});

  for (int i = 1; i < 2 * n; i++) {
    vector<int> v(n);
    auto &u = f.back();

    for (auto [x, y, z] : A)  // [x, y, value]
      v[x] = (v[x] + (long long)u[y] * z) % p;

    f.push_back(v);
  }

  vector<int> w(n);
  mt19937 gen;
  for (auto &x : w) x = uniform_int_distribution<int>(1, p - 1)(gen);

  vector<int> a(2 * n);
  for (int i = 0; i < 2 * n; i++)
    for (int j = 0; j < n; j++) a[i] = (a[i] + (long long)f[i][j] * w[j]) % p;

  auto c = berlekamp_massey(a);
  int m = (int)c.size();

  vector<int> ans(n);

  for (int i = 0; i < m - 1; i++)
    for (int j = 0; j < n; j++)
      ans[j] = (ans[j] + (long long)c[m - 2 - i] * f[i][j]) % p;

  int inv = power(p - c[m - 1], p - 2);

  for (int i = 0; i < n; i++) ans[i] = (long long)ans[i] * inv % p;

  return ans;
}

例题

  1. LibreOJ #163. 高斯消元 2
  2. ICPC2021 台北 Gym103443E. Composition with Large Red Plane, Yellow, Black, Gray, and Blue