【学习小记】Berlekamp-Massey算法

短命女 2022-10-02 15:41 128阅读 0赞

## Preface ##

BM算法是用来求一个数列的最短线性递推式的。  
形式化的，BM算法能够对于长度为n的有穷数列或者已知其满足线性递推的无穷数列\\(a\\)，找到最短的长度为m的有穷数列\\(c\\)，满足对于所有的\\(i\\geq n\\)，有\\\[a\_i=\\sum\\limits\_\{j=1\}^\{m\}c\_ja\_\{i-j\}\\\]

## Text ##

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

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

假设我们当前构造出来的递推系数C是第\\(cnt\\)版（经过cnt次修正）长度为\\(m\\)，能够满足前\\(a\_0...a\_\{i-1\}\\)项，记做\\(\_\{cnt\}C\\)，初始时\\(\_\{cnt\}C\\)为空,m=0  
记\\(d\_i=a\_i-\\sum\\limits\_\{j=1\}^\{m\}c\_ja\_\{i-j\}\\)

若\\(d\_i=0\\)，那么C符合的很好，不用管它

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

记\\(fail\_\{cnt\}\\)表示\\(\_\{cnt\}C\\)在\\(a\_i\\)处拟合失败。

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

否则我们考虑如何构造，如果能找到一个\\(C'\\)，满足对于\\(m\\leq j\\leq i-1\\)，都有\\(\\sum\\limits\_\{k=1\}^\{m\}c'\_ka\_\{j-k\}=0\\)，且\\(\\sum\\limits\_\{k=1\}^\{m\}c'\_ka\_\{i-k\}=1\\)

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

如何构造呢？我们可以利用之前的C！  
找到某一个\\(k\\in\[0..cnt-1\]\\)  
我们构造设\\(w=\{d\_i\\over d\_\{fail\_k\}\}\\)，构造\\(wC'=\\\{0,0,0,0,...,0,w,-w\*\{\_\{k\}C\}\\\}\\)

其中前面填上了\\(i-fail\_k-1\\)个0，后面相当于是\\(\_kC\\)乘上\\(-w\\)接在了后面。

为什么这是对的？其实很简单，对于\\(a\_i\\)，带进去的算出来的东西相当于是\\\[w\*a\_\{fail\_k\}-w(a\_\{fail\_k\}-d\_\{fail\_k\})=w\*d\_\{fail\_k\}=d\_i\\\]  
而对于\\(m\\leq j\\leq i-1\\)，算出来的是正好是\\(w\*a\_\{j-(i-fail\_k)\}-w\*a\_\{j-(i-fail\_k)\}=0\\)，利用了\\(\_kC\\)在1到\\(fail\_k-1\\)都满足关系式，而在\\(fail\_k\\)相差\\(d\\)的性质。

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

最坏时间复杂度显然是\\(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); } }

转载于:https://www.cnblogs.com/BAJimH/p/10840788.html