Berlekamp-Massey 算法（求数列的最短递推式）

用于求数列的最短递推式。

本文参考自 https://www.cnblogs.com/jz-597/p/14983564.html。

增量法，设 \(R_i\) 表示第 \(i\) 个历史递推式，当前为 \(R_{cnt}\)。

设 \(\Delta\) 表示真实的 \(a_i\) 与用 \(R_{cnt}\) 求出的 \(a_i\) 的差值，即 \(\Delta=a_i-\sum_{j=1}^{|R_{cnt}|}a_{i-j}R_{cnt,j}\)。

• 若 \(\Delta=0\)，则无事发生，跳过。
• 若 \(\Delta\neq 0\)，那么我们需要调整这个递推式使得它对 \(a_i\) 也成立。

接下来看第二种情况怎么处理。

有一个想法是找到另一个递推式，使得这个递推式在算 \(\cdots,a_{i-1}\) 时取 \(0\)，在算 \(a_i\) 时取 \(\Delta\)。

我们可以找到一个历史递推式 \(R_{lst}\)，设其失配位置为 \(fail_{lst}\)，失配时的差值为 \(\Delta_{lst}\)（显然 \(\Delta_{lst}\neq 0\)）。

设 \(d=\dfrac{\Delta}{\Delta_{lst}}\)，那么我们可以找到这么一个递推式 \(R'=\{0,0,\cdots,0,d,-d\cdot R_{lst,1},-d\cdot R_{lst,2},\cdots\}\)，其中前面是 \(i-fail_{lst}-1\) 个 \(0\)。

可以发现，使用这个递推式得到的 \(a_i'=\sum_{j=1}^{|R'|}R'_{j}a_{i-j}=d(a_{fail_{lst}}-\sum_{j=1}^{|R_{lst}|}R_{lst,j}a_{fail_{lst}-j})=d\cdot \Delta_{lst}=\Delta\)。

而对于任意 \((i-k)\in [|R'|+1,i-1]\)，\(a_{i-k}'=\sum_{j=1}^{|R'|}R'_j{a_{(i-k)-j}}=d(a_{fail_{lst}-k}-\sum_{j=1}^{|R_{lst}|}R_{lst,j}a_{(fail_{lst}-k)-j})=0\)。

那么我们令 \(R_{cnt+1}\gets R_{cnt}+R'\) 即可。

实质上是利用了以前的信息 \(R_{lst}\)：注意到失配前的差值都是 \(0\)，所以我们用差值来定义递推式。

那么我们怎么找到最短的 \(R'\) 呢？由于 \(|R'|=i+\big(|R_{lst}|-fail_{lst}\big)\)，所以我们直接找 \(|R_{lst}|-fail_{lst}\) 最小的就好了。

至于为什么求出来的递推式一定是最短的不会证。

【XSY3403】求数列的最短递推式 代码：

#include<bits/stdc++.h>

#define N 10010

using namespace std;

namespace modular
{
    const int mod=1000000007;
    inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
    inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
    inline int mul(int x,int y){return 1ll*x*y%mod;}
    inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
}using namespace modular;

inline int poww(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1) ans=mul(ans,a);
        a=mul(a,a);
        b>>=1;
    }
    return ans;
}

inline int read()
{
    int x=0,f=1;
    char ch=getchar();
    while(ch<'0'||ch>'9')
    {
        if(ch=='-') f=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9')
    {
        x=(x<<1)+(x<<3)+(ch^'0');
        ch=getchar();
    }
    return x*f;
}

int n,a[N],delta[N],fail[N];

vector<int>r[N];

int main()
{
    n=read();
    for(int i=1;i<=n;i++) a[i]=read();
    int cnt=0,lst=cnt;
    for(int i=1;i<=n;i++)
    {
        delta[cnt]=0;
        for(int j=0;j<r[cnt].size();j++)
            Add(delta[cnt],mul(r[cnt][j],a[i-j-1]));
        delta[cnt]=dec(a[i],delta[cnt]);
        if(!delta[cnt]) continue;
        fail[cnt]=i;
        if(!cnt)
        {
            r[++cnt].resize(i);
            continue;
        }
        r[cnt+1]=r[cnt];
        r[cnt+1].resize(max((int)r[cnt+1].size(),i+(int)r[lst].size()-fail[lst]));
        int d=mul(delta[cnt],poww(delta[lst],mod-2));
        Add(r[cnt+1][i-fail[lst]-1],d);
        for(int j=0;j<r[lst].size();j++)
            Add(r[cnt+1][i-fail[lst]+j],mul(dec(0,d),r[lst][j]));
        if((int)r[cnt].size()-fail[cnt]<(int)r[lst].size()-fail[lst]) lst=cnt;
        cnt++;
    }
    printf("%d\n1000000006 ",(int)r[cnt].size());
    for(int i=0;i<r[cnt].size();i++)
        cout<<r[cnt][i]<<" \0"[i==r[cnt].size()-1];
    return 0;
}