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）

就先更到这里吧。