【学习小记】Berlekamp-Massey算法

Preface

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

Text

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

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

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

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

否则，我们需要进行一定的修正，cntC_{cnt}CcntC需要变换到cnt+1C_{cnt+1}Ccnt+1C。

记failcnt=ifail_{cnt}=ifailcnt=i表示cntC_{cnt}CcntC在aia_iai处拟合失败。

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

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

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

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

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

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

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

最坏时间复杂度显然是O(nm)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));//rc为当前构造出的递推系数列,le为当前系数列的长度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) //如果是第一次出现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;//计算d_iLL mul=v*ksm(lv,mo-2)%mo;//计算wfo(j,1,le) rw[j]=rc[j];inc(rc[i-lf],mul);//lf为上次存的最短的C'的失配位置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];}}
}