当前位置：首页 > news >正文

Berlekamp–Massey 算法

news 2026/6/2 1:07:22

Berlekamp–Massey 算法

模板link
简单来说就是给你一个序列，求它的最短线性递推式，复杂度 $O(n^2)$。
为了说清楚这个算法，我用文字来演示计算一个例子的过程，让读者了解这个算法的基本思想，然后介绍如何用代码实现。

算法流程

先给一个例子：设 $a = (1, 3, 7, 15, 31，63)$。用 $r$ 来表示我们当前的这个递推式。
从左到右遍历每一个，对于每个数，我们进行以下操作：

算出根据当前递推式，这个位置理论上的值 $sum = \sum_{j = 1}^{k}r_j \times a_{i-j}$ （$k$ 为递推式长度）。
用实际值 $a_i$ 减去理论值，得到当前位置的差值 $delta$。
若 $delta$ 不为 $0$，尝试修改递推式使其继续成立，同时记录下 $ a_i - \sum_{j = 1}^{k}r_j \times a_{i-j} = delta$ 这个差值。

好，现在我们开始计算。首先遍历到 $a_1 = 1$，由于这是第一个数我们什么都不知道，就在递推式里放一个 $0$ 好了。于是当前 $r = (0)$。同时记录 $a_1 = 1$。

现在遍历到 $a_2 = 3$，根据递推式，算出理论值 $sum = 0$，发现 $delta = 3$。如何修改递推式呢？这时就要用到BM算法的核心思想：利用之前的差值来修改递推式。

之前我们记录了 $a_1 = 1$，所以可以利用它得到 $a_2 = a_1 \times 3$，所以修改递推式为 $r = (3)$，表示 $a_n = a_{n-1}\times 3$ （$n \geq $ 1）。同时我们记录下这次的差值 $a_2 = 3$。

现在我们遍历到 $a_3 = 7$，根据递推式，算出理论值 $sum = 3 \times a_2 = 9$，发现 $delta = 2$，即 $a_3 = a_2 \times 3 - 2$。同样利用之前的差值 $a_2 = 3$，修改为 $a_3 = a_2 \times 3 - a_2 \times \frac{2}{3} = a_2 \times \frac{7}{3}$，所以修改递推式为 $r = (\frac{7}{3})$，表示 $a_n = a_{n-1}\times \frac{7}{3}$ （$n \geq 2$）(注意 $n \geq 2$ 时才成立，所以第一项是不符合的，相当于第一二项是初始值)。同时记录 $a_3 - a_2 \times 3 = -2$ 这个差值。

之前有多个错误，用哪个才能保证递推式最短呢？这点会在后面代码的部分讲。

现在我们遍历到 $a_4 = 15$，同样的方法算出 $delta = -\frac{4}{3}$，即 $a_4 = a_3 \times \frac{7}{3} - \frac{4}{3}$。同样之前用记录的 $a_3 - a_2 \times 3 = -2$ 来修改，即 $a_4 = a_3 \times \frac{7}{3} + \frac{2}{3} \times (a_3 - a_2 \times 3) = a_3 \times 3 - a_2 \times 2$。所以修改递推式为 $r = (3, -2)$，表示 $a_n = a_{n-1} \times 3 + a_{n-2} \times (-2)$（$n\geq2$）。

遍历后面的数，发现都有 $delta = 0$，所以递推式符合要求，最终的递推式就是 $r=(3,-2)$。

代码实现

我们用 $r$ 记录递推式，$cnt$ 为递推式数量，$r_{cnt}$ 即为当前的递推式，用 $fail$ 记录递推式第一次失效的地方，用 $delta_i$ 记录 $a_i$ 与其理论值的差值。（请注意 $r$ 和 $fail$ 的下标表示的是第几个递推式，而 $delta$ 的下标是原数列的第几项）。

遍历到 $i$ 时，先像上面说的一样计算当前的 $delta_i$，如果 $delta_i=0$，说明递推式依然成立，直接continue即可；否则，则记录 $fail_{cnt} = i$。

  int sum = 0;for(int j = 0; j < r[cnt].size(); j++){sum = (sum + a[i - j - 1] * r[cnt][j] % mod) % mod;}delta[i] = a[i] - sum;if(!delta[i]) continue;fail[cnt] = i;

接下来考虑如何修改递推式。如果 $cnt = 0$，说明这是第一个数（或前面都是 $0$），显然将递推式改为 $i$ 个 $0$ 就可以了。

  if(!cnt){r[++cnt].resize(i);delta[i] = a[i];continue;}

否则，我们需要像上文说的一样，用利用之前的差值来凑出 $delta_i$，实现修改。

具体的来说，我们对于每一个 $cnt$ 之前的长度为 $k$ 的递推式 $id$，我们都有：

\[\large a_{fail_{id}}- \underbrace{\sum_{j=1}^{k}r[id][j] \times a_{fail_{id} - j}}_{\text{按照递推式}\, id \,\text{的理论值}}=delta_{fail_{id}} \]

设 $\large tmp = \frac{delta_i}{delta_{fail_{id}}}$，两边同乘 $tmp$ 就有：

\[\large tmp \cdot a_{fail_{id}}+(-tmp) \cdot r[id][1] \cdot a_{fail_{id}-1}+\cdots+(-tmp)\cdot r[id][k] \cdot a_{fail_{id}-k}=delta_i \]

在递推式 $r$ 中，第 $t$ 项表示的是序列中第 $i - t$ 项的系数，所以 $fail_{id}$ 的系数就在递推式的第 $i-fail_{id}$ 位。

所以我们构造的 $r'$ 应该长这样：

\[\large (\underbrace{0,0,\cdots,0}_{(i-fail_{id}-1)\text{个}0},tmp,-tmp\cdot r[id][1],-tmp\cdot r[id][2],\cdots,-tmp\cdot r[id][k]) \]

新的递推式就是 $r[cnt+1]=r[cnt]+r'$。

另外可以发现，这个新的递推式的长度就是 $i-fail_{id}+r[id].size()$，所以每次找最短的来修改就能求出最短递推式了。

  int id = cnt - 1, mi = i - fail[id] + r[id].size();for(int j = 0; j < cnt; j++){ //找出最短的if(i - fail[j] + r[j].size() < mi){id = j;mi = i - fail[j] + r[j].size();}}int tmp = delta[i] * qpow(delta[fail[id]], mod - 2, mod) % mod;cnt++;r[cnt] = r[cnt - 1]; //实现和上面的做法有点区别，这里是先加上上一个while(r[cnt].size() < mi) r[cnt].push_back(0); //然后补满长度r[cnt][i - fail[id] - 1] = (r[cnt][i - fail[id] - 1] + tmp) % mod; //最后填充后面的for(int j = 0; j < r[id].size(); j++){ //另外vector是以0开头的，注意下标r[cnt][i - fail[id] + j] = ((r[cnt][i - fail[id] + j] - tmp * r[id][j] % mod) % mod + mod) % mod;}

最后给出BM算法的完整代码：

inline void BM(){for(int i = 1; i <= n; i++){int sum = 0;for(int j = 0; j < r[cnt].size(); j++){sum = (sum + a[i - j - 1] * r[cnt][j] % mod) % mod;}delta[i] = a[i] - sum;if(!delta[i]) continue;fail[cnt] = i;if(!cnt){r[++cnt].resize(i);delta[i] = a[i];continue;}int id = cnt - 1, mi = i - fail[id] + r[id].size();for(int j = 0; j < cnt; j++){if(i - fail[j] + r[j].size() < mi){id = j;mi = i - fail[j] + r[j].size();}}int tmp = delta[i] * qpow(delta[fail[id]], mod - 2, mod) % mod;cnt++;r[cnt] = r[cnt - 1];while(r[cnt].size() < mi) r[cnt].push_back(0);r[cnt][i - fail[id] - 1] = (r[cnt][i - fail[id] - 1] + tmp) % mod;for(int j = 0; j < r[id].size(); j++){r[cnt][i - fail[id] + j] = ((r[cnt][i - fail[id] + j] - tmp * r[id][j] % mod) % mod + mod) % mod;}}len = r[cnt].size();for(int i = 0; i < r[cnt].size(); i++){p[i + 1] = (r[cnt][i] % mod + mod) % mod;cout << p[i + 1] << " ";}cout << '\n';
}

线性递推

以下结论的证明可以直接参考【模板】常系数齐次线性递推的题解。

我们想要找出一组系数 $f$ 使得：

\[\large a_m = \sum_{i=0}^{k-1}f_i\cdot a_i \]

也就是大项 $a_m$ 表示成前 $k$ 个初始项的线性组合。

然后结论就是这组系数是 $x_m$ 在模特征多项式意义下的多项式系数。对于这道题直接暴力快速幂即可。

因为做多项式乘法时两个 $k-1$ 次多项式相乘，最大可能出现 $2k-2$ 项，但是对于大于等于 $k$ 次的项在我们的状态里没有对应位置，所以每一个这样的 $x_i$ 都要用递推式拆成 $x_{i−1},\cdots,x_{i−k}$。

代码：

int solve(int m, int *p){f[0] = 1;g[1] = 1;auto mul = [&](int *a, int *b, int *c){for(int i = 0; i <= 2 * len; i++) tmp[i] = 0;for(int i = 0; i < len; i++){for(int j = 0; j < len; j++){tmp[i + j] = (tmp[i + j] + a[i] * b[j] % mod) % mod;}}for(int i = 2 * len; i >= len; i--){ //把大于等于k次的项用递推式拆开if(tmp[i]){for(int j = len; j >= 0; j--){tmp[i - j] = (tmp[i - j] + tmp[i] * p[j] % mod) % mod;}}}for(int i = 0; i <= 2 * len; i++) c[i] = tmp[i];return 0;}; //学习了大佬的匿名函数写法for(; m; m >>= 1){if(m & 1) mul(f, g, f);mul(g, g, g);}for(int i = 0; i < len; i++){ans = (ans + a[i + 1] * f[i] % mod) % mod;}return ans;
}

可以用NTT优化多项式乘法，不过这个题不需要。

完整代码：

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e4 + 5;
const int mod = 998244353;
int n, m, a[N], p[N], len, ans;
int qpow(int a, int b, int mod){int m = 1;for(; b; b >>= 1){if(b & 1) m = (m * a) % mod;a = (a * a) % mod;}return m;
}
vector<int> r[N];
int cnt, fail[N], delta[N];
inline void BM(){for(int i = 1; i <= n; i++){int sum = 0;for(int j = 0; j < r[cnt].size(); j++){sum = (sum + a[i - j - 1] * r[cnt][j] % mod) % mod;}delta[i] = a[i] - sum;if(!delta[i]) continue;fail[cnt] = i;if(!cnt){r[++cnt].resize(i);delta[i] = a[i];continue;}int id = cnt - 1, mi = i - fail[id] + r[id].size();for(int j = 0; j < cnt; j++){if(i - fail[j] + r[j].size() < mi){id = j;mi = i - fail[j] + r[j].size();}}int tmp = delta[i] * qpow(delta[fail[id]], mod - 2, mod) % mod;cnt++;r[cnt] = r[cnt - 1];while(r[cnt].size() < mi) r[cnt].push_back(0);r[cnt][i - fail[id] - 1] = (r[cnt][i - fail[id] - 1] + tmp) % mod;for(int j = 0; j < r[id].size(); j++){r[cnt][i - fail[id] + j] = ((r[cnt][i - fail[id] + j] - tmp * r[id][j] % mod) % mod + mod) % mod;}}len = r[cnt].size();for(int i = 0; i < r[cnt].size(); i++){p[i + 1] = (r[cnt][i] % mod + mod) % mod;cout << p[i + 1] << " ";}cout << '\n';
}
int f[N], g[N], tmp[N];
int solve(int m, int *p){f[0] = 1;g[1] = 1;auto mul = [&](int *a, int *b, int *c){for(int i = 0; i <= 2 * len; i++) tmp[i] = 0;for(int i = 0; i < len; i++){for(int j = 0; j < len; j++){tmp[i + j] = (tmp[i + j] + a[i] * b[j] % mod) % mod;}}for(int i = 2 * len; i >= len; i--){if(tmp[i]){for(int j = len; j >= 0; j--){tmp[i - j] = (tmp[i - j] + tmp[i] * p[j] % mod) % mod;}}}for(int i = 0; i <= 2 * len; i++) c[i] = tmp[i];return 0;};for(; m; m >>= 1){if(m & 1) mul(f, g, f);mul(g, g, g);}for(int i = 0; i < len; i++){ans = (ans + a[i + 1] * f[i] % mod) % mod;}return ans;
}
signed main(){ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin >> n >> m;for(int i = 1; i <= n; i++) cin >> a[i];BM();cout << solve(m, p);return 0;
}

查看全文

http://www.jsqmd.com/news/619784/