跳转至

# Berlekamp–Massey 算法

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

### 定义

定义一个数列 {𝑎 0…𝑎 𝑛−1}{a 0…a n−1} 的递推式为满足下式的序列 {𝑟 0…𝑟 𝑚}{r 0…r m}：

∑𝑚 𝑗=0 𝑟 𝑗 𝑎 𝑖−𝑗 =0,∀𝑖 ≥𝑚∑j=0 m r j a i−j=0,∀i≥m

其中 𝑟 0 =1 r 0=1．𝑚 m 称为该递推式的 阶数．

数列 {𝑎 𝑖}{a i} 的最短递推式即为阶数最小的递推式．

### 做法

与上面定义的稍有不同，这里定义一个新的递推系数 {𝑓 0…𝑓 𝑚−1}{f 0…f m−1}，满足：

𝑎 𝑖 =∑𝑚−1 𝑗=0 𝑓 𝑗 𝑎 𝑖−𝑗−1,∀𝑖 ≥𝑚 a i=∑j=0 m−1 f j a i−j−1,∀i≥m

容易看出 𝑓 𝑖 = −𝑟 𝑖+1 f i=−r i+1，并且阶数 𝑚 m 与之前的定义是相同的．

我们可以增量地求递推式，按顺序考虑 {𝑎 𝑖}{a i} 的每一位，并在递推结果出现错误时对递推系数 {𝑓 𝑖}{f i} 进行调整．方便起见，以下将前 𝑖 i 位的最短递推式记为 𝐹 𝑖 ={𝑓 𝑖,𝑗}F i={f i,j}．

显然初始时有 𝐹 0 ={}F 0={}．假设递推系数 𝐹 𝑖−1 F i−1 对数列 {𝑎 𝑖}{a i} 的前 𝑖 −1 i−1 项均成立，这时对第 𝑖 i 项就有两种情况：

1. 递推系数对 𝑎 𝑖 a i 也成立，这时不需要进行任何调整，直接令 𝐹 𝑖 =𝐹 𝑖−1 F i=F i−1 即可．
2. 递推系数对 𝑎 𝑖 a i 不成立，这时需要对 𝐹 𝑖−1 F i−1 进行调整，得到新的 𝐹 𝑖 F i．

设 Δ 𝑖 =𝑎 𝑖 −∑𝑚 𝑗=0 𝑓 𝑖−1,𝑗 𝑎 𝑖−𝑗−1 Δ i=a i−∑j=0 m f i−1,j a i−j−1，即 𝑎 𝑖 a i 与 𝐹 𝑖−1 F i−1 的递推结果的差值．

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

否则设上一次对递推系数进行修改时，已考虑的 {𝑎 𝑖}{a i} 的项数为 𝑘 k．如果存在一个序列 𝐺 ={𝑔 0…𝑔 𝑚′−1}G={g 0…g m′−1}，满足：

∑𝑚′−1 𝑗=0 𝑔 𝑗 𝑎 𝑖′−𝑗−1 =0,∀𝑖′ ∈[𝑚′,𝑖)∑j=0 m′−1 g j a i′−j−1=0,∀i′∈[m′,i)

并且 ∑𝑚′−1 𝑗=0 𝑔 𝑗 𝑎 𝑖−𝑗−1 =Δ 𝑖∑j=0 m′−1 g j a i−j−1=Δ i，那么不难发现将 𝐹 𝑘 F k 与 𝐺 G 按位分别相加之后即可得到一个合法的递推系数 𝐹 𝑖 F i．

考虑如何构造 𝐺 G．一种可行的构造方案是令

𝐺 ={0,0,…,0,Δ 𝑖 Δ 𝑘, −Δ 𝑖 Δ 𝑘 𝐹 𝑘−1}G={0,0,…,0,Δ i Δ k,−Δ i Δ k F k−1}

其中前面一共有 𝑖 −𝑘 −1 i−k−1 个 0 0，且最后的 −Δ 𝑖 Δ 𝑘 𝐹 𝑘−1−Δ i Δ k F k−1 表示将 𝐹 𝑘−1 F k−1 每项乘以 −Δ 𝑖 Δ 𝑘−Δ i Δ k 后接在序列后面．

不难验证此时 ∑𝑚′−1 𝑗=0 𝑔 𝑗 𝑎 𝑖−𝑗−1 =Δ 𝑘 Δ 𝑖 Δ 𝑘 =Δ 𝑖∑j=0 m′−1 g j a i−j−1=Δ k Δ i Δ k=Δ i，因此这样构造出的是一个合法的 𝐺 G．将 𝐹 𝑖 F i 赋值为 𝐹 𝑘 F k 与 𝐺 G 逐项相加后的结果即可．

如果要求的是符合最开始定义的递推式 {𝑟 𝑖}{r i}，则将 {𝑓 𝑗}{f j} 全部取相反数后在最开始插入 𝑟 0 =1 r 0=1 即可．

从上述算法流程中可以看出，如果数列的最短递推式的阶数为 𝑚 m，则算法的复杂度为 𝑂(𝑛 𝑚)O(n m)．最坏情况下 𝑚 =𝑂(𝑛)m=O(n)，因此算法的最坏复杂度为 𝑂(𝑛 2)O(n 2)．

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

参考实现

     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 算法求解的是有限项数列的最短递推式．如果待求递推式的序列有无限项，但已知最短递推式的阶数上界，则只需取出序列的前 2 𝑚 2 m 项即可求出整个序列的最短递推式．（证明略）

### 应用

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

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

如果要求向量列 𝒗 𝑖 v i 的最短递推式，设向量的维数为 𝑛 n，我们可以随机一个 𝑛 n 维行向量 𝐮 𝑇 u T，并计算标量序列 {𝒖 𝑇 𝒗 𝑖}{u T v i} 的最短递推式．由 Schwartz–Zippel 引理，二者的最短递推式有至少 1 −𝑛 𝑝 1−n p 的概率相同．

求矩阵列 {𝐴 𝑖}{A i} 的最短递推式也是类似的，设矩阵的大小为 𝑛 ×𝑚 n×m，则只需随机一个 1 ×𝑛 1×n 的行向量 𝐮 𝑇 u T 和一个 𝑚 ×1 m×1 的列向量 𝒗 v，并计算标量序列 {𝒖 𝑇 𝐴 𝑖 𝒗}{u T A i v} 的最短递推式即可．由 Schwartz–Zippel 引理可以类似地得到二者相同的概率至少为 1 −𝑛+𝑚 𝑝 1−n+m p．

#### 优化矩阵快速幂

设 𝒇 𝑖 f i 是一个 𝑛 n 维列向量，并且转移满足 𝒇 𝑖 =𝐴 𝒇 𝑖−1 f i=A f i−1，则可以发现 {𝒇 𝑖}{f i} 是一个不超过 𝑛 n 阶的线性递推向量列．（证明略）

我们可以直接暴力求出 𝒇 0…𝒇 2 𝑛−1 f 0…f 2 n−1，然后用前面提到的做法求出 {𝒇 𝑖}{f i} 的最短递推式，再调用 常系数齐次线性递推 即可．

如果要求的向量是 𝒇 𝑚 f m，则算法的复杂度是 𝑂(𝑛 3 +𝑛 l o g⁡𝑛 l o g⁡𝑚)O(n 3+n log⁡n log⁡m)．如果 𝐴 A 是一个只有 𝑘 k 个非零项的稀疏矩阵，则复杂度可以降为 𝑂(𝑛 𝑘 +𝑛 l o g⁡𝑛 l o g⁡𝑚)O(n k+n log⁡n log⁡m)．但由于算法至少需要 𝑂(𝑛 𝑘)O(n k) 的时间预处理，因此在压力不大的情况下也可以使用 𝑂(𝑛 2 l o g⁡𝑚)O(n 2 log⁡m) 的线性递推算法，复杂度同样是可以接受的．

#### 求矩阵的最小多项式

方阵 𝐴 A 的最小多项式是次数最小的并且满足 𝑓(𝐴) =0 f(A)=0 的多项式 𝑓 f．

实际上最小多项式就是 {𝐴 𝑖}{A i} 的最小递推式，所以直接调用 Berlekamp–Massey 算法就可以了．如果 𝐴 A 是一个 𝑛 n 阶方阵，则显然最小多项式的次数不超过 𝑛 n．

瓶颈在于求出 𝐴 𝑖 A i，因为如果直接每次做矩阵乘法的话复杂度会达到 𝑂(𝑛 4)O(n 4)．但考虑到求矩阵列的最短递推式时实际上求的是 {𝒖 𝑇 𝐴 𝑖 𝒗}{u T A i v} 的最短递推式，因此我们只要求出 𝐴 𝑖 𝒗 A i v 就行了．

假设 𝐴 A 有 𝑘 k 个非零项，则复杂度为 𝑂(𝑘 𝑛 +𝑛 2)O(k n+n 2)．

#### 求稀疏矩阵行列式

如果能求出方阵 𝐴 A 的特征多项式，则常数项乘上 ( −1)𝑛(−1)n 就是行列式．但是最小多项式不一定就是特征多项式．

实际上如果把 𝐴 A 乘上一个随机对角阵 𝐵 B，则 𝐴 𝐵 A B 的最小多项式有至少 1 −2 𝑛 2−𝑛 𝑝 1−2 n 2−n p 的概率就是特征多项式．最后再除掉 d e t 𝐵 det B 就行了．

设 𝐴 A 为 𝑛 n 阶方阵，且有 𝑘 k 个非零项，则复杂度为 𝑂(𝑘 𝑛 +𝑛 2)O(k n+n 2)．

#### 求稀疏矩阵的秩

设 𝐴 A 是一个 𝑛 ×𝑚 n×m 的矩阵，首先随机一个 𝑛 ×𝑛 n×n 的对角阵 𝑃 P 和一个 𝑚 ×𝑚 m×m 的对角阵 𝑄 Q, 然后计算 𝑄 𝐴 𝑃 𝐴 𝑇 𝑄 Q A P A T Q 的最小多项式即可．

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

设 𝐴 A 有 𝑘 k 个非零项，且 𝑛 ≤𝑚 n≤m，则复杂度为 𝑂(𝑘 𝑛 +𝑛 2)O(k n+n 2)．

#### 解稀疏方程组

问题：已知 𝐴 𝐱 =𝐛 A x=b, 其中 𝐴 A 是一个 𝑛 ×𝑛 n×n 的 满秩 稀疏矩阵，𝐛 b 和 𝐱 x 是 1 ×𝑛 1×n 的列向量．𝐴,𝐛 A,b 已知，需要在低于 𝑛 𝜔 n ω 的复杂度内解出 𝑥 x．

做法：显然 𝐱 =𝐴−1 𝐛 x=A−1 b．如果我们能求出 {𝐴 𝑖 𝐛}{A i b}(𝑖 ≥0 i≥0) 的最小递推式 {𝑟 0…𝑟 𝑚−1}{r 0…r m−1}(𝑚 ≤𝑛 m≤n), 那么就有结论

𝐴−1 𝐛 = −1 𝑟 𝑚−1∑𝑚−2 𝑖=0 𝐴 𝑖 𝐛 𝑟 𝑚−2−𝑖 A−1 b=−1 r m−1∑i=0 m−2 A i b r m−2−i

（证明略）

因为 𝐴 A 是稀疏矩阵，直接按定义递推出 𝐛…𝐴 2 𝑛−1 𝐛 b…A 2 n−1 b 即可．

同样地，设 𝐴 A 中有 𝑘 k 个非零项，则复杂度为 𝑂(𝑘 𝑛 +𝑛 2)O(k n+n 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
    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

> 本页面最近更新：2026/1/7 08:56:54，更新历史
> 发现错误？想一起完善？ 在 GitHub 上编辑此页！
> 本页面贡献者：Tiphereth-A, antileaf, Enter-tainer, AntiLeaf, ZnPdCo
> 本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用