【学习笔记】Berlekamp-Massey算法

【算法简介】

• Berlekamp-Massey算法，常简称为BM算法，是用来求解一个数列的最短线性递推式的算法。
  
• BM算法可以在\(O(N^2)\)的时间内求解一个长度为\(N\)的数列的最短线性递推式。
  
• 在当今OI竞赛界，尚没有很多BM算法的应用，但在一些输入的数很少的题目中，BM能够成为发掘题目性质的一大助力，甚至有可能直接解出答案的线性递推式，不失为一种不错的工具。

【算法流程】

• 对于数列\(\{a_1,a_2,a_3,...,a_n\}\)，我们称数列\(\{r_1,r_2,r_3,...,r_m\}\)为其线性递推式当且仅当\(a_i=\sum_{j=1}^{m}r_j*a_{i-j}\)对于任意\(m+1≤i≤n\)均成立。
• 若数列\(\{a_1,a_2,a_3,...,a_n\}\)的线性递推式\(\{r_1,r_2,r_3,...,r_m\}\)还满足\(m\)是所有该数列的线性递推式中最小的，则称\(\{r_1,r_2,r_3,...,r_m\}\)为数列\(\{a_1,a_2,a_3,...,a_n\}\)的最短线性递推式。
• 现在考虑我们已经求得了\(\{a_1,a_2,a_3,...,a_{i-1}\}\)的最短线性递推式\(\{r_1,r_2,r_3,...,r_m\}\)，如何求得\(\{a_1,a_2,a_3,...,a_i\}\)的最短线性递推式。
• 定义\(\{a_1,a_2,a_3,...,a_{i-1}\}\)的最短线性递推式\(\{r_1,r_2,r_3,...,r_m\}\)为当前递推式，记递推式被更改的次数为\(cnt\)，第\(i\)次更改后的递推式为\(R_i\)，特别地，定义\(R_0\)为空，那么当前递推式应当为\(R_{cnt}\)。
• 记\(delta_i=a_i-\sum_{j=1}^{m}r_j*a_{i-j}\)，其中\(\{r_1,r_2,r_3,...,r_m\}\)为当前递推式，显然若\(delta_i=0\)，那么当前递推式就是\(\{a_1,a_2,a_3,...,a_i\}\)的最短线性递推式。
• 否则，我们认为\(R_{cnt}\)在\(a_i\)处出错了，定义\(fail_i\)为\(R_i\)最早的出错位置，则有\(fail_{cnt}=i\)。
• 考虑对\(R_{cnt}\)进行修改，使其变为\(R_{cnt+1}\)，并在\(a_i\)处同样成立。
• 若\(cnt=0\)，这意味着\(a_i\)是序列中第一个非零元素，我们可以令\(R_{cnt+1}=\{0,0,0,...,0\}\)，即用\(i\)个0填充线性递推式，此时由于不存在\(j\)使得\(m+1≤j≤i\)，因此\(R_{cnt+1}\)显然为\(\{a_1,a_2,a_3,...,a_i\}\)的线性递推式，并且由于\(a_i\)是序列中第一个非零元素，不难证明\(R_{cnt+1}\)也是\(\{a_1,a_2,a_3,...,a_i\}\)的最短线性递推式。
• 否则，即\(cnt>0\)，考虑\(R_{cnt-1}\)出错的位置\(fail_{cnt-1}\)，记\(mul=\frac{delta_i}{delta_{fail_{cnt-1}}}\)。
• 我们希望得到数列\(R'=\{r'_1,r'_2,r'_3,...,r'_{m'}\}\)，使得\(\sum_{j=1}^{m'}r'_j*a_{k-j}=0\)对于任意\(m'+1≤k≤i-1\)均成立，并且\(\sum_{j=1}^{m'}r'_j*a_{i-j}=delta_i\)。如果能够找到这样的数列\(R'\)，那么令\(R_{cnt+1}=R_{cnt}+R'\)即可（其中\(+\)定义为各位分别相加）。
• 数列\(R'\)可以是下述数列：\(\{0,0,0,...,0,mul,-mul*R_{cnt-1}\}\)即填充\(i-fail_{cnt-1}-1\)个零，然后将数列\(\{1,-R_{cnt-1}\}\)的\(mul\)倍放在后面。此时有\(\sum_{j=1}^{m'}r'_j*a_{i-j}=delta_{fail_{cnt-1}}*mul=delta_i\)，并且\(\sum_{j=1}^{m'}r'_j*a_{k-j}=0\)对于任意\(m'+1≤k≤i-1\)均成立。
• 故令\(R_{cnt+1}=R_{cnt}+R'\)即可。
• 在最坏情况下，我们可能需要对数列进行\(O(N)\)次修改，因此该算法的时间复杂度为\(O(N^2)\)。

【一组实例】

• 以数列\(\{1,2,4,9,20,40,90\}\)为例，我们来具体地理解一下算法流程。
• 初始时，我们有\(R_0=\{\},cnt=0\)。
• \(i=1\)时，将\(a_1=1\)代入递推式，得到\(delta_1=1\)，\(R_0\)在\(i=1\)时出错，记\(fail_0=1\)。由于此时\(cnt=0\)，我们将递推式修改为\(R_1=\{0\}\)。
• \(i=2\)时，将\(a_2=2\)代入递推式，得到\(delta_2=2\)，\(R_1\)在\(i=2\)时出错，记\(fail_1=2\)。此时\(mul=2\)，可以构造得到\(R'=\{2\}\)，递推式被修改为\(R_2=\{2\}\)。
  
• \(i=3\)时，将\(a_3=4\)代入递推式，得到\(delta_3=0\)，\(R_2\)没有出错。
  
• \(i=4\)时，将\(a_4=9\)代入递推式，得到\(delta_4=1\)，\(R_2\)在\(i=4\)时出错，记\(fail_2=4\)。此时\(mul=0.5\)，可以构造得到\(R'=\{0,0.5,0\}\)，递推式被修改为\(R_3=\{2,0.5,0\}\)。
  
• \(i=5\)时，将\(a_5=20\)代入递推式，得到\(delta_5=0\)，\(R_3\)没有出错。
  
• \(i=6\)时，将\(a_6=40\)代入递推式，得到\(delta_6=-4.5\)，\(R_3\)在\(i=6\)时出错，记\(fail_3=6\)。此时\(mul=-4.5\)，可以构造得到\(R'=\{0,-4.5,9\}\)，递推式被修改为\(R_4=\{2,-4,9\}\)。
  
• \(i=7\)时，将\(a_7=90\)代入递推式，得到\(delta_7=9\)，\(R_4\)在\(i=7\)时出错，记\(fail_4=7\)。此时\(mul=-2\)，可以构造得到\(R'=\{-2,4,1,0\}\)，递推式被修改为\(R_5=\{0,0,10,0\}\)。
  
• 因此以数列\(\{1,2,4,9,20,40,90\}\)的递推式即为\(R_5=\{0,0,10,0\}\)。

【代码】

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

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 2005;
const double eps = 1e-8;
template <typename T> void read(T &x) {
	x = 0; int f = 1;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
	for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
	x *= f;
}
int cnt, fail[MAXN];
double val[MAXN], delta[MAXN];
vector <double> ans[MAXN];
int main() {
	int n; read(n);
	for (int i = 1; i <= n; i++)
		scanf("%lf", &val[i]);
	for (int i = 1; i <= n; i++) {
		double tmp = val[i];
		for (unsigned j = 0; j < ans[cnt].size(); j++)
			tmp -= ans[cnt][j] * val[i - j - 1];
		delta[i] = tmp;
		if (fabs(tmp) <= eps) continue;
		fail[cnt] = i;
		if (cnt == 0) {
			ans[++cnt].resize(i);
			continue;
		}
		double mul = delta[i] / delta[fail[cnt - 1]];
		cnt++; ans[cnt].resize(i - fail[cnt - 2] - 1);
		ans[cnt].push_back(mul);
		for (unsigned j = 0; j < ans[cnt - 2].size(); j++)
			ans[cnt].push_back(ans[cnt - 2][j] * -mul);
		if (ans[cnt].size() < ans[cnt - 1].size()) ans[cnt].resize(ans[cnt - 1].size());
		for (unsigned j = 0; j < ans[cnt - 1].size(); j++)
			ans[cnt][j] += ans[cnt - 1][j];
	}
	for (unsigned i = 0; i < ans[cnt].size(); i++)
		cout << ans[cnt][i] << ' ';
	return 0;
}