Berlekamp-Massey
阿新 • • 發佈:2021-07-01
\(BM\)演算法(\(Berlekamp-Massey\))
給定一個長度為\(n\)的序列,在\(O(n^2)\)的時間內求出序列的最短遞推式(前提是要能通過這\(n\)項求出至少一個遞推式)
\(\forall i>m,a_i=\sum_{j=1}^mf_ja_{i-j}\)
考慮增量法,設當前遞推式為\(f_{cnt}\),check到了\(a_i\)
設\(dlt_i=a_i-\sum_{j=1}^{m_{cnt}}f_{cnt,j}a_{i-j}\)
- \(dlt_i=0\),說明\(f_{cnt}\)合法
- 否則考慮在\(f_{cnt}\)的基礎上,加一些東西,構造出一個新的\(f\)
設\(fail_i\)表示第\(i\)個遞推式失效的位置,對於之前求出的一個遞推式\(f_{prv}\)有(設\(x=fail_{prv}\)):
\(dlt_{x}=a_x-\sum_{j=1}^{m_{prv}}f_{prv,j}a_{x-j}\)
把下標代換一下可得到(設\(i-x=k\))
\(dlt_x=a_{i-k}-\sum_{j=1}^{m_{prv}}f_{prv,j}a_{i-k-j}\)
於是可以求出一個\(f'\)使得\(\sum_{j=1}^{m'}f'_ja_{i-j}=dlt_x\)
但是我們想加上的東西是\(dlt_i\),所以係數乘一下就行
於是就構造出來\(f_{cnt+1}=f_{cnt}+\frac{dlt_i}{dlt_x}f'\)
對應一下下標關係就是加上一個\((0\cdots0,k,-kf_{prv,1},-kf_{prv,2}\cdots)\),前面是\(i-fail_{prv}-1\)個\(0\)
想求出最短的遞推式,就取\(i-fail_{prv}+len(prv)\)最小的就行了
for(register int i=1; i<=n; i++){ // for(register int j=0; j<dp[cnt].size(); j++) printf("%d ",dp[cnt][j]);puts(""); if(!cnt){ if(a[i]){ dlt[i] = a[i]; fail[cnt++] = i; dp[cnt].resize(i,0); } continue; } int &sum = dlt[i], num = dp[cnt].size(); sum = a[i]; fail[cnt] = i; for(register int j=0; j<num; j++) sum = dec(sum,1ll*a[i-j-1]*dp[cnt][j]%MOD); if(!sum) continue; int prv = cnt-1, mnlen = i-fail[prv]+(int)dp[prv].size(); for(register int j=0; j<cnt; j++) if(i-fail[j]+(int)dp[j].size()<mnlen) mnlen = i-fail[j]+(int)dp[j].size(), prv = j; int k = 1ll*dlt[i]*ksm(dlt[fail[prv]])%MOD; dp[cnt+1] = dp[cnt]; cnt++; while(dp[cnt].size()<mnlen) dp[cnt].push_back(0); iinc(dp[cnt][i-fail[prv]-1],k); num = dp[prv].size(); for(register int j=0; j<num; j++) ddec(dp[cnt][i-fail[prv]+j],1ll*k*dp[prv][j]%MOD); }