Berlekamp–Massey 算法
Berlekamp–Massey 算法
简单来说就是给你一个序列，求它的最短线性递推式，复杂度 \(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 tmp = \frac{delta_i}{delta_{fail_{id}}}\)，两边同乘 \(tmp\) 就有：
在递推式 \(r\) 中，第 \(t\) 项表示的是序列中第 \(i - t\) 项的系数，所以 \(fail_{id}\) 的系数就在递推式的第 \(i-fail_{id}\) 位。
所以我们构造的 \(r'\) 应该长这样：
新的递推式就是 \(r_{cnt+1}=r_{cnt}+r'\)。
代码实现时，我们可以先加上上一个递推式，然后直接修改 \(i - fail_{id}\) 后面的部分。
另外可以发现，这个新的递推式的长度就是 \(i-fail_{id}+len(r_{id})\,\,\, (len\text{表示递推式长度})\)，所以每次找最短的来修改就能求出最短递推式了。
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;
}
最后给出 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\) 使得：
也就是大项 \(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){
if(m < n) return a[m + 1] % mod;
if(len == 0) return 0;
if(len == 1) return a[1] * qpow(p[1], m, mod) % mod;
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){
if(m < n) return a[m + 1] % mod;
if(len == 0) return 0;
if(len == 1) return a[1] * qpow(p[1], m, mod) % mod;
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;
}
浙公网安备 33010602011771号