Berlekamp-Massey算法学习小记

简介

Berlekamp-Massey算法，简称BM算法，可以在$O(N^2)$O(N2)时间内求解一个数列的最短线性递推式。

教程

一篇讲的很详细的博客

Berlekamp-Massey算法

我们采用增量法构造数列$\{a_n\}${an​}最短线性递推式。

假设现在已经得到了数列$a_1,a_2,...,a_{i-1}$a1​,a2​,...,ai−1​的最短线性递推式，且在这之前递推式被修改过$cnt$cnt次，设第$k$k次修改后的递推式为$R_k$Rk​。

如果在加入$a_i$ai​后原来的递推式仍然满足，也就是有$\sum_{k=1}^mr_k*a_{i-k}=a_i$k=1∑m​rk​∗ai−k​=ai​
那么数列$a_1,a_2,...,a_n$a1​,a2​,...,an​的最短线性递推式就是$R_{cnt}$Rcnt​。

不然的话，表明第$cnt$cnt次修改后的递推式在位置$i$i处出错。定义$fail_{k}$failk​表示第$k$k次修改后的递推式在$fail_{k}$failk​处出错，显然有$fail_{cnt}=i$failcnt​=i，同时定义$delta_{cnt}=a_i-\sum_{k=1}^mr_k*a_{i-k}$deltacnt​=ai​−∑k=1m​rk​∗ai−k​

如果$cnt=0$cnt=0，表明$a_i$ai​是数列中第一个非零元素，不难证明当前最短线性递推数列就是$\{0,0,...,0\}${0,0,...,0}也就是$i$i个$0$0。

不然的话，我们可以考虑构造一个新的递推式$R'$R′，若对于任意$m'+1\le j\le i-1$m′+1≤j≤i−1有$\sum_{k=1}^{m'}r'_{k}*a_{j-k}=0$k=1∑m′​rk′​∗aj−k​=0
且满足$\sum_{k=1}^{m'}r'_k*a_{i-k}=delta_{cnt}$k=1∑m′​rk′​∗ai−k​=deltacnt​

那么我们只要把$R_{cnt}$Rcnt​和$R'$R′的对应位分别相加就可以得到一个合法的线性递推式。

考虑通过以下方法构造$R'$R′：

设$mul=\frac{delta_{cnt}}{delta_{cnt-1}}$mul=deltacnt−1​deltacnt​​，使$R'=\{0,0,...,0,mul,-mul*R_{cnt-1}\}$R′={0,0,...,0,mul,−mul∗Rcnt−1​}，也就是最前面有$i-fail_{cnt-1}-1$i−failcnt−1​−1个$0$0，再把数列$\{1,-R_{cnt-1}\}${1,−Rcnt−1​}乘上$mul$mul后接到后面。

为什么这样的$R'$R′是合法的呢？

首先$\sum_{k=1}^{m'}r'_k*a_{i-k}=mul*delta_{cnt-1}=delta_{cnt}$k=1∑m′​rk′​∗ai−k​=mul∗deltacnt−1​=deltacnt​
其次对于所有$m'+1\le j\le i-1$m′+1≤j≤i−1，设$p=j-i+fail_{cnt-1}$p=j−i+failcnt−1​，有$\sum_{k=1}^{m'}r'_k*a_{j-k}=mul*(a_p-\sum_{k=1}^{m_{cnt-1}}r_{cnt-1,k}*a_{p-k})$k=1∑m′​rk′​∗aj−k​=mul∗(ap​−k=1∑mcnt−1​​rcnt−1,k​∗ap−k​)
又因为$R_{cnt-1}$Rcnt−1​对于数列$a_1,...,a_{fail_{cnt-1}}$a1​,...,afailcnt−1​​是一个合法的递推式，故上式等于$0$0。

所以新的递推式就是$R'+R_{cnt}$R′+Rcnt​。

代码

这是模1e9+7意义下的版本
（话说最后求出来线性递推式貌似并不是最短的？）

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>

typedef long long LL;

const int MOD=1000000007;
const int N=3005;

int n,a[N],r[N][N],delta[N],len[N],fail[N];

int ksm(int x,int y)
{
	int ans=1;
	while (y)
	{
		if (y&1) ans=(LL)ans*x%MOD;
		x=(LL)x*x%MOD;y>>=1;
	}
	return ans;
}

int main()
{
	scanf("%d",&n);
	for (int i=1;i<=n;i++) scanf("%d",&a[i]);
	int cnt=0;
	for (int i=1;i<=n;i++)
	{
		int t=a[i];
		for (int j=1;j<=len[cnt];j++) (t+=MOD-(LL)r[cnt][j]*a[i-j]%MOD)%=MOD;
		if (!t) continue;
		delta[cnt]=t;
		fail[cnt]=i;
		if (!cnt)
		{
			len[++cnt]=i;
			continue;
		}
		int mul=(LL)delta[cnt]*ksm(delta[cnt-1],MOD-2)%MOD;
		cnt++;
		len[cnt]=std::max(len[cnt-1],i-fail[cnt-2]+len[cnt-2]);
		for (int j=1;j<=len[cnt-1];j++) r[cnt][j]=r[cnt-1][j];
		(r[cnt][i-fail[cnt-2]]+=mul)%=MOD;
		for (int j=1;j<=len[cnt-2];j++) (r[cnt][j+i-fail[cnt-2]]+=MOD-(LL)r[cnt-2][j]*mul%MOD)%=MOD;
	}
	printf("%d\n",len[cnt]);
	for (int i=1;i<=len[cnt];i++) printf("%d ",r[cnt][i]);
	return 0;
}