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]);
    }
};
分类: 文章

3 条评论

Remmina · 2019年1月11日 11:56 下午

勘误:

  • $deg(H) > deg(G)$的情况确实存在,比如 1 1 2 3 5 8 13 21 34 55 1 1 ...,会在第 12 个数字处取到(这个估计您已经知道了)
  • 代码所得到的不一定是最优解(其实就是上面那个问题)(相比您也知道了)

我就是提醒您一下以防您忘了修正(大雾)

    boshi · 2019年1月12日 7:42 上午

    已经修正,并通过多种数据初步验证了代码无误。

【题解】 [NOI2007]生成树计数 矩阵树定理+高斯消元+Berlekamp-Massey(BM)+矩阵快速幂 BZOJ – – K-XZY · 2019年1月12日 12:13 上午

[…] 关于 BM 的理解我认为 Boshi 大佬写的挺好的:https://www.k-xzy.xyz/archives/7500 […]

发表回复

Avatar placeholder

您的电子邮箱地址不会被公开。 必填项已用 * 标注