Berlekamp-Massey算法簡單介紹

請閱讀本文的同窗們注意：以前這篇博客和所附代碼有點問題，求的不必定是最短遞推式，很是抱歉

看毛爺爺的論文大概斷斷續續看了一個月了，看得不是很懂，google了一波好像很快就看懂了，就先口胡一下這個算法好了。這篇文章的介紹方式大概和毛爺爺的論文不大一致，固然算法的本質是一致的。

參考連接：https://grocid.net/2012/11/22/berlekamp-massey-algorithm-explained/

Berlekamp-Massey算法是用來求一個數列的遞推式的。

我如今有一個序列$x_1,x_2...x_n$，首先咱們先來定義遞推式。序列x的遞推式是個序列$p_1,p_2...p_m$。一個合法遞推式知足對於$i\geq{m+1}$，$x_i=\sum_{j=1}^mp_jx_{i-j}$。一個遞推式的次數定義爲m，注意這裏的每一個p均可覺得0。

咱們定義一個遞推式帶入一個序列的第i項所得的結果爲$\sum_{j=1}^mx_{i-j}p_j-x_i$，若是$i\leq{m}$定義爲0。那麼合法遞推式就是帶進去每一項都得0的遞推式對吧。

好，如今咱們想要求遞推式。怎麼求呢，咱們開始隨便口胡一個遞推式，就{}好了。

咱們考慮第一個數，這樣一個一個往下考慮。咱們看看如今這個多項式滿不知足x[i]，知足就留着。

若是不知足，咱們就要想辦法修復這個遞推多項式對吧。

一種naive的想法就是說我開始有個遞推多項式爲g，它對於前i-1項全成立，到了第i項掛了。咱們來算算這個遞推式帶進第i項，假設是s，那麼咱們要想辦法把帶進去獲得的值減去s。

若是咱們憑空變出一個遞推式f，它知足這個遞推式帶進前i-1項全爲0，帶進第i項是1，那麼咱們只要把g減去sf不就好了嗎。

若是咱們不是第一次遇到這種狀況，那麼咱們以前也曾經有一個遞推式，它知足前j-1項，帶進去剛好在第j項不爲0，不妨假設在第j項爲1，若是不爲1除掉便可。那麼咱們只要在這個遞推式稍微移一移項，前面補若干個0就能夠獲得f了。

若是有多個知足條件的遞推式？那麼咱們選擇最近的一個確定最優。

咱們能夠維護當前最優的一個遞推式，使得它補完0以後的長度最少，取這個就好了。

咱們來舉個例子，我有一個數列1,2,4,10,24,50,116。

咱們一個一個數考慮，首先數列裏啥都沒有，遞推式是{}。

而後看到一個1，一臉懵逼，那麼遞推式隨便寫一個{0}好了。

接下來咱們發現到了2這裏掛了。以前1掛掉的經歷告訴咱們x[1]=1，那麼咱們令f={1}，咱們把遞推式改爲{2}好了。

接下來到了4，咱們發現一切良好。

到了10，咱們發現又掛了，帶進去這項獲得了-2。以前2掛的經歷告訴咱們x[2]=2，那麼咱們令f={0,0.5,0}，把遞推式加上2f={0,1,0}，那麼遞推式變成了{2,1,0}。（注意這裏結尾有一個0，由於以前遞推式裏面是{0}）

到了24，一切良好。

到了50，咱們發現又掛了，帶進去獲得了8。以前10掛的經歷告訴咱們-x[4]+2x[3]=-2，那麼兩邊除以-2，0.5x[4]-x[3]=1，因此能夠令f={0,0.5,-1}，那麼8f={0,4,-8}，因此遞推式要減去8f，變成了{2,-3,8}。

到了116，又掛掉了，帶進去獲得了-8。以前50的狗帶經歷告訴咱們-x[6]+2x[5]+x[4]=8，那麼兩邊除以8獲得-0.125x[6]+0.25x[5]+0.125x[4]=1，那麼令f={-0.125,0.25,0.125,0}，8f={-1,2,1,0}，那麼遞推式要加上8f，變成{1,-1,9,0}（這裏結尾的0是由於原來帶了一個0）。

實際程序寫起來仍是挺簡單的。

int n,pn=0,fail[2333];
typedef vector<ld> vld;
vld ps[2333]; ld x[2333],delta[2333];
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf",x+i);
    int best=0;
    for(int i=1;i<=n;i++)
    {
        ld dt=-x[i];
        for(int j=0;j<ps[pn].size();j++)
            dt+=x[i-j-1]*ps[pn][j];
        delta[i]=dt;
        if(fabs(dt)<=1e-7) continue;
        fail[pn]=i; if(!pn) {ps[++pn].resize(i); continue;}
        vld&ls=ps[best]; ld k=-dt/delta[fail[best]];
        vld cur; cur.resize(i-fail[best]-1); //trailing 0
        cur.pb(-k); for(int j=0;j<ls.size();j++) cur.pb(ls[j]*k);
        if(cur.size()<ps[pn].size()) cur.resize(ps[pn].size());
        for(int j=0;j<ps[pn].size();j++) cur[j]+=ps[pn][j];
        if(i-fail[best]+(int)ps[best].size()>=ps[pn].size()) best=pn;
        ps[++pn]=cur;
    }
    for(unsigned g=0;g<ps[pn].size();g++)
        cout<<ps[pn][g]<<" "; puts("");
}

什麼，卡空間？只有兩項是有用的，能夠滾動數組。

因爲某種緣由，具體應用下一篇再說好了。下一篇鴿了

那麼Berlekamp-Massey算法有啥用呢？好比我有一個優秀的dp，f[i][j]表示某種神奇的量，而後對1<=i<=n，j<=m，$f[i][j]=\sum_s f[i-1][s]*val[j][s]$，val是與i無關的量，而後要求f[n][1]之類的。之前咱們用矩陣快速冪，如今咱們就能夠用BM算出f[n][1]的遞推式，而後悶聲m^2logn。

如下附贈一個素質二連板子：

const int MOD=1e9+7;
ll qp(ll a,ll b)
{
    ll x=1; a%=MOD;
    while(b)
    {
        if(b&1) x=x*a%MOD;
        a=a*a%MOD; b>>=1;
    }
    return x;
}
namespace linear_seq {
inline vector<int> BM(vector<int> x)
{
    vector<int> ls,cur;
    int pn=0,lf,ld;
    for(int i=0;i<int(x.size());++i)
    {
        ll t=-x[i]%MOD;
        for(int j=0;j<int(cur.size());++j)
            t=(t+x[i-j-1]*(ll)cur[j])%MOD;
        if(!t) continue;
        if(!cur.size())
        {cur.resize(i+1); lf=i; ld=t; continue;}
        ll k=-t*qp(ld,MOD-2)%MOD;
        vector<int> c(i-lf-1); c.pb(-k);
        for(int j=0;j<int(ls.size());++j) c.pb(ls[j]*k%MOD);
        if(c.size()<cur.size()) c.resize(cur.size());
        for(int j=0;j<int(cur.size());++j)
            c[j]=(c[j]+cur[j])%MOD;
        if(i-lf+(int)ls.size()>=(int)cur.size())
            ls=cur,lf=i,ld=t;
        cur=c;
    }
    vector<int>&o=cur;
    for(int i=0;i<int(o.size());++i)
        o[i]=(o[i]%MOD+MOD)%MOD;
    return o;
}
int N; ll a[SZ],h[SZ],t_[SZ],s[SZ],t[SZ];
inline void mull(ll*p,ll*q)
{
    for(int i=0;i<N+N;++i) t_[i]=0;
    for(int i=0;i<N;++i) if(p[i])
        for(int j=0;j<N;++j)
            t_[i+j]=(t_[i+j]+p[i]*q[j])%MOD;
    for(int i=N+N-1;i>=N;--i) if(t_[i])
        for(int j=N-1;~j;--j)
            t_[i-j-1]=(t_[i-j-1]+t_[i]*h[j])%MOD;
    for(int i=0;i<N;++i) p[i]=t_[i];
}
inline ll calc(ll K)
{
    for(int i=N;~i;--i) s[i]=t[i]=0;
    s[0]=1; if(N!=1) t[1]=1; else t[0]=h[0];
    for(;K;mull(t,t),K>>=1) if(K&1) mull(s,t); ll su=0;
    for(int i=0;i<N;++i) su=(su+s[i]*a[i])%MOD;
    return (su%MOD+MOD)%MOD;
}
inline int gao(vector<int> x,ll n)
{
    if(n<int(x.size())) return x[n];
    vector<int> v=BM(x); N=v.size(); if(!N) return 0;
    for(int i=0;i<N;++i) h[i]=v[i],a[i]=x[i];
    return calc(n);
}
}

你能夠用這個數據來測試你的BM板子：https://files.cnblogs.com/files/zzqsblog/BM-data.zip （mod 1e9+7，最短長度爲712或713）

就先更到這裏吧。