都 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 \leq i \leq 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\)。
记 \(\Delta_i=a_i−\sum_{j=1}^m r_j∗a_{i−j}\),其中 \(r_{1 \cdots m}\) 为当前递推式,显然若 \(\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 \leq 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{\Delta_i}{\Delta_{\mathrm{fail}_{c-1}}}\)。我们希望得到数列 \(R'=r'_{1 \cdots m'}\),使得 \(\sum_{j=1}^{m'} r'_j \cdot a_{k-j} = 0\) 对于任意 \(m'+1≤k≤i−1\) 均成立,并且 \(\sum_{j=1}^{m'} r'_j \cdot a_{i-j} = \Delta_i\)。如果能够找到这样的数列 \(R'\),那么令 \(R_{c+1}=R_c+R'\) 即可(其中 \(+\) 定义为各位分别相加)。
构造数列 \(R'\) 如下:\(\{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'\) 即可。在最坏情况下,我们可能需要对数列进行 \(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;
}