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}-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了。

如果有多个满足条件的递推式?那么我们选择最近的一个肯定最优。

我们来举个例子,我有一个数列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;
typedef vector<ld> vld;
vld ps[2333]; int pn=0,fail[SZ];
ld x[SZ],delta[SZ];
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf",x+i);
    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(1); continue;}
        vld&ls=ps[pn-1]; ld k=-dt/delta[fail[pn-1]];
        vld cur; cur.resize(i-fail[pn-1]-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];
        ps[++pn]=cur;
    }
    for(unsigned g=0;g<ps[pn].size();g++)
        cout<<ps[pn][g]<<" "; puts("");
}

什么?你说你的递推式有常数项?这个简单,你遇到第一个非0数的时候直接把它令给常数项就好了。

什么,卡空间?递推式只有最近两项是有用的,可以滚动数组。

由于某种原因,具体应用下一篇再说好了。

posted @ 2017-05-19 11:10 fjzzq2002 阅读(...) 评论(...) 编辑 收藏