[学习笔记] Berlekamp-Massey 算法

都 2202 年了，现代 OIer 早该会会了！参考了 此博客。

引入

Berlekamp-Massey 算法，又称为 BM 算法，其可以在 $O(n^2)$ 时间内求解一个长度为 $n$ 的数列的最短线性递推式。

在当今 OI 界，尚没有很多 BM 算法的应用，但在一些输入的数很少的题目中，BM 能够成为发掘题目性质（找规律）的一大助力，甚至有可能直接解出答案的线性递推式，不失为一种有效的工具。

算法流程

对于数列 $a_{1\cdots n}$，我们称数列 $r_{1\cdots m}$ 为其线性递推式当且仅当 $a_i=\sum _{j=1}^m{r}_j\cdot a_{i-j}$ 对于任意 $m+1\le i\le n$ 均成立。若数列 $a_{1\cdots n}$ 的线性递推式 $r_{1\cdots m}$ 还满足 $m$ 是所有该数列的线性递推式中最小的，则称 $r_{1\cdots m}$ 为数列 $a_{1\cdots n}$ 的最短线性递推式。我们的问题是，在数列 $a$ 已知的情况下，如何求出其最短线性递推式 $r$。

考虑增量法，假设我们已经求出了 $a_{1\cdots i-1}$ 的最短线性递推式 $r_{1\cdots m}$，如何求出 $a_{1\cdots i}$ 的最短线性递推式。

定义 $a_{1\cdots i-1}$ 的最短线性递推式 $r_{1\cdots m}$ 为当前递推式，记递推式被更改的次数为 $c$，第 $i$ 次更改后的递推式为 $R_i$，特别地，定义 $R_0$ 为空，那么当前递推式应当为 $R_c$。

记 ${\mathrm{\Delta }}_i=a_i-\sum _{j=1}^m{r}_j\ast a_{i-j}$，其中 $r_{1\cdots m}$ 为当前递推式，显然若 ${\mathrm{\Delta }}_i=0$，那么当前递推式就是 $a_{1\cdots i}$ 的最短线性递推式。否则，我们认为 $R_c$ 在 $a_i$ 处出错了，定义 ${\mathrm{fail}}_i$ 为 $R_i$ 最早的出错位置，则有 ${\mathrm{fail}}_c=i$。

若 $c=0$，这意味着 $a_i$ 是序列中第一个非零元素，我们可以令 $R_{c+1}=\{0,0,0,...,0\}$，即用 $i$ 个 $0$ 填充线性递推式，此时由于不存在 $j$ 使得 $m+1\le j\le i$，因此 $R_{c+1}$ 显然为 $a_{1\cdots i}$ 的线性递推式，并且由于 $a_i$ 是序列中第一个非零元素，不难证明 $R_{c+1}$ 也是 $a_{1\cdots i}$ 的最短线性递推式。

否则，即 $c>0$，考虑 $R_{c-1}$ 出错的位置 ${\mathrm{fail}}_{c-1}$，记 $\mathrm{mul}=\frac{{\mathrm{\Delta }}_i}{{\mathrm{\Delta }}_{{\mathrm{fail}}_{c-1}}}$。我们希望得到数列 $R^{\prime }=r_{1\cdots m^{\prime }}^{\prime }$，使得 $\sum _{j=1}^{m^{\prime }}{r}_j^{\prime }\cdot a_{k-j}=0$ 对于任意 $m^{\prime }+1\le k\le i-1$ 均成立，并且 $\sum _{j=1}^{m^{\prime }}{r}_j^{\prime }\cdot a_{i-j}={\mathrm{\Delta }}_i$。如果能够找到这样的数列 $R^{\prime }$，那么令 $R_{c+1}=R_c+R^{\prime }$ 即可（其中 $+$ 定义为各位分别相加）。

构造数列 $R^{\prime }$ 如下：$\{0,0,0,...,0,\mathrm{mul},-\mathrm{mul}\cdot R_{c-1}\}$，即填充 $i-{\mathrm{fail}}_{c-1}-1$ 个零，然后将数列 $\{1,-R_{c-1}\}$ 的 $\mathrm{mul}$ 倍放在后面。容易验证其合法性，故令 $R_{c+1}=R_c+R^{\prime }$ 即可。在最坏情况下，我们可能需要对数列进行 $O(n)$ 次修改，因此该算法的时间复杂度为 $O(n^2)$。

代码

以下代码实现了求解给定 $n$ 元数列在实数域上的最短线性递推式。显然，BM 算法只需要数域中每个非零元素均存在乘法逆元即可实现，读者不妨自行实现一下在模质数意义下的 BM 算法。

Code

/*
最黯淡的一个　梦最为炽热
万千孤单焰火　让这虚构灵魂鲜活
至少在这一刻　热爱不问为何
存在为将心声响彻
*/
#include <bits/stdc++.h>
#define pii pair<int, int>
#define mp(x, y) make_pair(x, y)
#define pb push_back
#define eb emplace_back
#define fi first
#define se second
#define int long long
#define mem(x, v) memset(x, v, sizeof(x))
#define mcpy(x, y, n) memcpy(x, y, sizeof(int) * (n))
#define lob lower_bound
#define upb upper_bound
using namespace std;

inline int read() {
	int x = 0, w = 1;char ch = getchar();
	while (ch > '9' || ch < '0') { if (ch == '-')w = -1;ch = getchar(); }
	while (ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar();
	return x * w;
}

const int MN = 2e3 + 5;
const int Mod = 998244353;
const int inf = 1e9;
const double eps = 1e-8;

inline int qPow(int a, int b = Mod - 2, int ret = 1) {
    while (b) {
        if (b & 1) ret = ret * a % Mod;
        a = a * a % Mod, b >>= 1;
    }
    return ret;
}

#define dbg

int N, c, fail[MN];
double val[MN], delta[MN];
vector <double> ans[MN];

signed main(void) {
    N = read();
    for (int i = 1; i <= N; i++) 
        scanf("%lf", &val[i]);
    for (int i = 1; i <= N; i++) {
        double tmp = val[i];
        for (int j = 0; j < ans[c].size(); j++) 
            tmp -= ans[c][j] * val[i - j - 1];
        delta[i] = tmp;
        if (fabs(tmp) <= eps) continue;
        fail[c] = i;
        if (!c) {
            ans[++c].resize(i);
            continue;
        }
        double mul = delta[i] / delta[fail[c - 1]];
        ++c, ans[c].resize(i - fail[c - 2] - 1);
        ans[c].pb(mul);
        for (int j = 0; j < ans[c - 2].size(); j++)
            ans[c].pb(ans[c - 2][j] * -mul);
        if (ans[c].size() < ans[c - 1].size()) ans[c].resize(ans[c - 1].size());
        for (int j = 0; j < ans[c - 1].size(); j++)
            ans[c][j] += ans[c - 1][j];
    }
    for (int i = 0; i < ans[c].size(); i++)
        printf("%.lf ", ans[c][i]);
    return 0;
}