Berlekamp–Massey 算法

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

定义

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

其中。称为该递推式的阶数。

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

做法

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

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

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

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

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

设，即与的递推结果的差值。

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

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

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

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

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

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

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

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

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

参考实现

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 算法就可以了。如果是一个阶方阵，则显然最小多项式的次数不超过。

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

假设有个非零项，则复杂度为。

求稀疏矩阵行列式

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

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

设为阶方阵，且有个非零项，则复杂度为。

求稀疏矩阵的秩

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

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

设有个非零项，且，则复杂度为。

解稀疏方程组

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

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

（证明略）

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

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

参考实现

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;
}

例题

最后更新: 2024年10月9日
创建日期: 2022年3月18日