【学习小记】Berlekamp-Massey算法

Preface

BM算法是用来求一个数列的最短线性递推式的。
形式化的，BM算法能够对于长度为n的有穷数列或者已知其满足线性递推的无穷数列$a$a，找到最短的长度为m的有穷数列$c$c，满足对于所有的$i\geq n$i≥n，有$a_i=\sum\limits_{j=1}^{m}c_ja_{i-j}$ai​=j=1∑m​cj​ai−j​

Text

BM算法的流程十分简洁明了——增量，构造，修正。

方便起见，我们令a的下标从0开始，c的下标从1开始

假设我们当前构造出来的递推系数C是第$cnt$cnt版（经过cnt次修正）长度为$m$m，能够满足前$a_0...a_{i-1}$a0​...ai−1​项，记做$_{cnt}C$cnt​C，初始时$_{cnt}C$cnt​C为空,m=0
记$d_i=a_i-\sum\limits_{j=1}^{m}c_ja_{i-j}$di​=ai​−j=1∑m​cj​ai−j​

若$d_i=0$di​=0，那么C符合的很好，不用管它

否则，我们需要进行一定的修正，$_{cnt}C$cnt​C需要变换到$_{cnt+1}C$cnt+1​C。

记$fail_{cnt}$failcnt​表示$_{cnt}C$cnt​C在$a_i$ai​处拟合失败。

若$cnt=0$cnt=0，说明这是a的第一个非0元素，直接设$m=i+1$m=i+1,在$C$C中填上i+1个0。显然这满足定义式（因为前m项是可以不满足递推式的）。

否则我们考虑如何构造，如果能找到一个$C'$C′，满足对于$m\leq j\leq i-1$m≤j≤i−1，都有$\sum\limits_{k=1}^{m}c'_ka_{j-k}=0$k=1∑m​ck′​aj−k​=0，且$\sum\limits_{k=1}^{m}c'_ka_{i-k}=1$k=1∑m​ck′​ai−k​=1

那么可以构造$_{cnt+1}C=_{cnt}C+d_iC'$cnt+1​C=cnt​C+di​C′，显然这一定满足性质。其中加法为按项数对应加。

如何构造呢？我们可以利用之前的C！
找到某一个$k\in[0..cnt-1]$k∈[0..cnt−1]
我们构造设$w={d_i\over d_{fail_k}}$w=dfailk​​di​​，构造$wC'=\{0,0,0,0,...,0,w,-w*{_{k}C}\}$wC′={0,0,0,0,...,0,w,−w∗k​C}

其中前面填上了$i-fail_k-1$i−failk​−1个0，后面相当于是$_kC$k​C乘上$-w$−w接在了后面。

为什么这是对的？其实很简单，对于$a_i$ai​，带进去的算出来的东西相当于是$w*a_{fail_k}-w(a_{fail_k}-d_{fail_k})=w*d_{fail_k}=d_i$w∗afailk​​−w(afailk​​−dfailk​​)=w∗dfailk​​=di​
而对于$m\leq j\leq i-1$m≤j≤i−1，算出来的是正好是$w*a_{j-(i-fail_k)}-w*a_{j-(i-fail_k)}=0$w∗aj−(i−failk​)​−w∗aj−(i−failk​)​=0，利用了$_kC$k​C在1到$fail_k-1$failk​−1都满足关系式，而在$fail_k$failk​相差$d$d的性质。

此时我们还希望总的长度最短，也就是说$max(m_{cnt},i-fail_k+m_{k})$max(mcnt​,i−failk​+mk​)最短。
我们只需要动态维护最短的$i-fail_k+m_{k}$i−failk​+mk​即可，每次算出$_{cnt+1}$cnt+1​时都与之前的k比较一下谁更短即可，这样贪心可以感受出来是正确的。

最坏时间复杂度显然是$O(nm)$O(nm)的

Code

LL rc[4*N],rp[4*N],le,le1,rw[4*N];
void BM()
{
	le=le1=0;
	memset(rc,0,sizeof(rc));
	memset(rp,0,sizeof(rp));
	int lf=0;LL lv=0;
	fo(i,0,n1)
	{
		LL v=0;
		fo(j,1,le) inc(v,rc[j]*ap[i-j]%mo);
		if(v==ap[i]) continue;
		if(le==0) 
		{
			le=i+1;
			fo(j,1,le) rc[j]=rp[j]=0;
			le1=0,lf=i,lv=(ap[i]-v)%mo;
			continue;
		}
		v=(ap[i]-v+mo)%mo;
		LL mul=v*ksm(lv,mo-2)%mo;
		
		fo(j,1,le) rw[j]=rc[j];
		
		inc(rc[i-lf],mul);
		fo(j,i-lf+1,i-lf+le1) inc(rc[j],(mo-mul*rp[j-(i-lf)]%mo)%mo);
		if(le<i-lf+le1)
		{
			swap(le1,le);
			le=i-lf+le,lf=i,lv=v;
			fo(j,1,le1) rp[j]=rw[j];
		}
		
		v=0;
		fo(j,1,le) inc(v,rc[j]*ap[i-j]%mo);
	}
}