Berlekamp_Massey 算法 | pinkyHead's blog

问题模型

给定一个有n个元素的数列a，其中第i个元素是 $a_i$ 。

求一个较短或最短的数列b，假设b有m个元素，那么要求满足 $\forall m<i\leq n, a_i = \sum_{j=1}^m a_{i-j} b_j$

要求在 $O(n^2)$ 的时间复杂度内解决此问题

Berlekamp_Massey

假设递推式经过c次更新，第i次更新后递推式为 $R[i]$ 。初始时，定义 $R[0]$ 为空。

假设当前递推式的长度为m，这时候考虑用递推式去算 $a[i]$ 是否成立

设 $delta_i = a_i - \sum_{j=1}^m a_{i-j} R[c][j]$

如果 $delta_i=0$ ，那么 $R[c]$ 依然合法，不用修改

否则，设 $fail_c=i$ 表示递推式 $R[c]$ 第一次失效的位置为i

如果 $c=0$ ，即初始状态。则把 $R[1]$ 设置为i个0,因为对于空串如果成立，则说明i的元素之前都是0

如果 $c\neq 0$ ，则把问题转换成构造出一个递推式 $R'$ ，使得对于 $|R'|+1\leq k < i$ ，有 $\sum_j a_{k-j} R'_j = 0$ ，且 $\sum_j a_{i-j} R'_j = delta_i$

设 $0\leq id < c$ ，设 $tmp = \frac{delta_i}{delta_{fail[id]}}$ ，则可以构造

$R' = \{ 0,0,\cdots, 0, tmp, -tmp\cdot R[id][1],- tmp \cdot R[id][2],\cdots \}$

其中开头有 $i-fail[id]-1$ 个0，tmp之后是-tmp倍的 $R[id]$

令 $R[c+1] = R[c] + R'$ 即可

该构造的正确性证明
首先，当 $k=n$ 时，显然 $R'$ 计算得到的答案为 $tmp*delta_{fail[id]}=delta_i$
其次，当 $|R'|+1\leq k < i$ 时， $R[id]$ 满足使得和式为0。注意当k变小时，相当于与 $R[id]$ 配对的数向左移动了，若 $|R'|+1\leq k ==i-1$ 时成立，则 $|R'|+1\leq kk < i-1$ 时也必然成立。因为每个递推式对于任何下标大于该递推式长度，且小于递推式出错位置的 $a_i$ ，一定成立

故可以在 $O(n^2)$ 内求出数列 $a_i$ 的一个较短线性递推式

若要求最短的线性递推式，则在每次id取值时，每次找 $min{i - fail[id] + len(R[id])}$ 即可。易证明，id取值为 $cnt-1$ 时最优~~但是我目前没有发现有标程直接用这个结论，所以如果怕出错的话还是老老实实挨个找最优的id吧~~

#include <bits/stdc++.h>
#define clr(x) memset(x,0,sizeof (x))
#define For(i,a,b) for (int i=a;i<=b;i++)
#define Fod(i,b,a) for (int i=b;i>=a;i--)
#define pb(x) push_back(x)
#define mp(x,y) make_pair(x,y)
#define fi first
#define se second
#define _SEED_ ('C'+'L'+'Y'+'A'+'K'+'I'+'O'+'I')
#define outval(x) printf(#x" = %d\n",x)
#define outvec(x) printf("vec "#x" = ");for (auto _v : x)printf("%d ",_v);puts("")
#define outtag(x) puts("----------"#x"----------")
#define outarr(a,L,R) printf(#a"[%d...%d] = ",L,R);\
                        For(_v2,L,R)printf("%d ",a[_v2]);puts("");
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef vector <int> vi;
LL read(){
    LL x=0,f=0;
    char ch=getchar();
    while (!isdigit(ch))
        f|=ch=='-',ch=getchar();
    while (isdigit(ch))
        x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return f?-x:x;
}
const int N=0x1233,mod=1e9+7;
void Add(int &x,int y){
    if ((x+=y)>=mod)
        x-=mod;
}
void Del(int &x,int y){
    if ((x-=y)<0)
        x+=mod;
}
int Pow(int x,int y){
    int ans=1;
    for (;y;y>>=1,x=(LL)x*x%mod)
        if (y&1)
            ans=(LL)ans*x%mod;
    return ans;
}
int n,cnt;
int a[N];
int Fail[N],delta[N];
vector <int> R[N];
int main(){
    n=read();
    For(i,1,n)
        a[i]=read();
    R[0].clear();
    cnt=0;
    For(i,1,n){
        if (cnt==0){
            if (a[i]){
                Fail[cnt++]=i;
                delta[i]=a[i];
                R[cnt].resize(0);
                R[cnt].resize(i,0);
            }
            continue;
        }
        int sum=0,m=R[cnt].size();
        delta[i]=a[i];
        Fail[cnt]=i;
        For(j,0,m-1)
            Add(sum,(LL)a[i-j-1]*R[cnt][j]%mod);
        Del(delta[i],sum);
        if (!delta[i])
            continue;
        int id=cnt-1,v=i-Fail[id]+(int)R[id].size();
        For(j,0,cnt-1)
            if (i-Fail[j]+(int)R[j].size()<v)
                id=j,v=i-Fail[j]+(int)R[j].size();
        int tmp=(LL)delta[i]*Pow(delta[Fail[id]],mod-2)%mod;
        R[cnt+1]=R[cnt];
        while (R[cnt+1].size()<v)
            R[cnt+1].pb(0);
        Add(R[cnt+1][i-Fail[id]-1],tmp);
        For(j,0,(int)R[id].size()-1)
            Del(R[cnt+1][i-Fail[id]+j],(LL)tmp*R[id][j]%mod);
        cnt++;
    }
    printf("%d\n",(int)R[cnt].size());
    For(i,0,(int)R[cnt].size()-1)
        printf("%d ",R[cnt][i]);
    puts("");
    return 0;
}