Berlekamp Massey 算法

一些定义

$F,G$：这样的大写字母一般代表一个数列

$deg(f)$：数列 f 的长度，下标从 1 开始

功能

为一个序列，求一个非平凡的齐次线性递推方程。

解释

非平凡的齐次线性递推方程，即找到一种规律，使这个序列满足这个规律。同时这个规律必然是通过前 $n-1$或更短的前缀找到的，然后后面几项依旧满足这个规律。

形式化定义

已知 $F$，求最短的 $G$使：
$$
\forall deg(F)\geq n > deg(G),\ F_n = \sum_{i=1}^{deg(G)}G_iF_{n-i}
$$
也就是说，找到一个序列 $G$，使得通过 $G$和 $F$的前几项，可以 “ 生成” 出 $F$的后面的所有项。这里需要注意两点，第一点，” 生成” 指的是通过上面的递推方式推算出的值与原来的 $F$完全一样。第二点，构造出的 $G$不一定能生成出 $F$的前几项，因为 $F_0,F_{-1}$等项是未定义的。

原理

• 1. 如果存在 $F,G_0,G_1$，且 $n=deg(F)=deg(G_0)=deg(G_1)$， $G_0,G_1$根据 $F$生成的数列的第 $n+1$项分别为 $a,b$，那么 $(G_0+G_1)$根据 $F$生成的数列的第 $n+1$项就是 $a+b$。
• 2. 将 $G$移动一位，第一位补 0，则 $G$根据 $F$生成的数列整体移动一位。
• 3. 如果 $G$通过 $F$生成的第 $[m+1,k]$项与 $F$相同，那么 $G$移动一位，第一项补-1 再根据 $F$生成的数列第 $[m+2,k+1]$项为 0。

以上三个原理的证明不难，在此略过。根据以上几个原理，我们就可以找到一种通过调整 $G$获得正确的递推方程的方法。

步骤

• ①随便给出一个递推数列 $G={1}$，和一个用于调整 $G$的数列 $H={\frac{1}{F[1]}}$。
• ②算法不断用 $G$通过 $F$生成第 $i\in[1,n]$位，生成的数列叫做 $I$
  • $I[i] = F[i]$：数列 $G$暂时通过了检测，将 $H$移动一位，第一位补 $0$，返回②
  • $I[i] \not = F[i]$：数列 $G$没有通过检测，通过将 $G$加上 $H$的倍数，使得 $I[i]=F[i]$，然后 $H$变为 $G$，并移动一位，第一位补 $-1$，返回②
• 最终得到的 $G$，当 $deg(G) < deg(F)/2$时，就是最短的 $G$。其它时候，$G$只是一个能够满足递推关系的序列。

疑问

• 1. 如果在第二步的第二种情况下，$deg(H)>deg(G)$，是否还需要交换 $G,H$?

第二步的第二种情况，$deg(H)>deg(G)$的情况可能存在，如果输入的数列为 $1,2,4,8,1,2,4,8,16,32,64,128,256$，在这种情况下交换 $G,H$会使答案变成 $2,0,0,0,0,0$，但是正确答案为 $2,0,0,0,0$。我还没有构造出其它可以区分这两种写法的数据。如果大家发现了，可以在这里留言。

• 2.BM 算法什么时候可以保证答案正确，且是最短的？

对于我下面贴出的奇葩代码来说，如果将输入数据的前导 0 全部去除，该算法的输出的长度若不大于输入数据的一半，则可保证答案是最短的。(当然，我并不会证明正确性，这份代码的正确性也没有得到过验证，仅仅通过了若干组对拍、题目、典型数列的验证)

实现

namespace Berlekamp_Massey
{
    ll A[MX], F[MX], G[MX], T[MX];
    int lenA, lenF, lenG, lenT;

    ll calc(int p, int len, ll* x)
    {
        ll ret = 0;
        for(int i=1; i<=len&&(p-i); i++)
            ret = (ret + A[p-i]*x[i]) % MOD;
        return ret%MOD;
    }

    void work()
    {
        int beg = 1;
        while(!A[beg] && beg<=lenA) beg++;
        memmove(A+1, A+beg, sizeof(ll)*(lenA-beg+1));
        lenA -= beg-1; 
        G[1] = inv(A[1]), lenG = 1;
        F[1] = 0, lenF = 1;
        for(int i=2; i<=lenA; i++)
        {
            ll delta = (A[i]-calc(i, lenF, F)+MOD)%MOD;
            if(!delta)
            {
                memmove(G+2, G+1, sizeof(ll)*lenG);
                G[1] = 0;
                lenG ++;
            }
            else
            {
                for(int j=1; j<=lenF; j++) T[j] = F[j];
                for(int j=1; j<=lenG; j++) T[j] = (T[j] + delta*G[j]) % MOD;
                lenT = max(lenF, lenG);
                if(lenF <= lenG)
                {
                    memmove(G+2, F+1, sizeof(ll)*lenF);
                    G[1] = MOD-1;
                    lenG = lenF+1;
                    ll now = inv(calc(i+1, lenG, G));
                    for(int j=1; j<=lenG; j++) G[j] = G[j] * now % MOD;
                }
                else
                {
                    memmove(G+2, G+1, sizeof(ll)*lenG);
                    G[1] = 0;
                    lenG++;
                }
                memmove(F+1, T+1, sizeof(ll)*lenT);
                lenF = lenT;
            }
        }
        printf("%d\n", lenF);
        for(int i=1; i<=lenF; i++) printf("%lld ", F[i]);
        putchar('\n');
    }

    void input()
    {
        scanf("%d", &lenA);
        for(int i=1; i<=lenA; i++) scanf("%lld", &A[i]);
    }
};