Berlekamp-Massey 算法

考虑给定序列 \(a\)，求出序列 \(a\) 的最短线性递推式（在固定域 \(\mathbb{F}\) 下，通常是模意义或者实数域）。

引理：对于一个无限长的递推序列 \(a\)，若其满足线性递推且其最短递推式长度为 \(l\)，则取前 \(2l\) 项跑Berlekamp-Massey即可得到其线性递推式。也就是需要前 \(2l\) 项才能完整体现递推关系。

换言之，对于给出的长为 \(n\) 的序列，若 Berlekamp-Massey 给出的线性递推式长度为 \(t\)，则 \(2t\le n\) 即合法。

问题等价于求出最短的多项式 \(C(x)=1-c_1x-c_2x^2-\dots -c_Lx^{L}\)，使得对于 \(i\ge L\)，满足：

\[\sum_{j=0}^Lc_ja_{i-j}=0 \]

定义：

\[d_i=\sum_{j=0}^L a_{i-j}[x^j]C(x) \]

我们维护几个东西：

1. \(C(x)\) 当前求出的答案。
2. \(L\)：当前的递推长度。 \(\deg C(x)\le L\)。
3. \(B(x)\)：上一次修改 \(L\) 之前的答案，也就是修改 \(L\) 之前的 \(C(x)\)。
4. \(b\) ：\(B(x)\) 被修改的那一次对应的 \(d\)。
5. \(m\)：上次修改的时刻（下标）

算法流程：

初始化 \(C(x)=1,L=0,m=-1,b=1,B(x)=1\)：

枚举 \(i:0\to n\)：

1. 计算 \(d_i\)

2. \(d_i=0\)，说明合法，跳过，否则进入第三步。

3. \(d_i\neq 0\)

   现在失配了，我们想要修正 \(C(x)\) 为新的 \(C’(x)\)，在其满足之前的递推关系的情况下，这个位置成功匹配。

   关键的发现是，\(B(x)\) 是在 \(<m\) 处 \(d\) 全部为 \(0\) 但在 \(=m\) 处 \(d=b\) 的多项式。

   那么可以构造出多项式：

   \[\Delta(x)=x^{i-m}B(x) \]

   性质：\(\Delta(x)\) 在 \(<i\) 的下标的 \(d\) 全部都是 \(0\)，在位置 \(i\) 时 \(d=b\)。

   证明：

   \[\begin{aligned} d^{(\Delta)}_j&=\sum_{k=0}^{\deg \Delta} \Delta_k·a_{j-k}\\ &=\sum_{k=0}^{\deg\Delta}a_{j-k}\times [x^k]x^{i-m}B(x)\\ &=\sum_{k=i-m}^{\deg\Delta }a_{j-k}[x^{k-(i-m)}]B(x)\\ &=d^{(B)}_{j-(i-m)} \end{aligned} \]

   而我们知道 \(d_{j-(i-m)}\)，对于 \(j<i\) 而言，\(j-(i-m)<m\)，而 \(d^{(B)}_{<m}=0\)，\(d^{(B)}_m=b\)。

   证毕。

   因此我们就知道：

   \[C'(x)=C(x)+\frac{-d_i^{(C)}x^{i-m}}{b}B(x) \]

   那么：

   \[\deg C'=\max(\deg C,(i-m)+\deg B)=\max(L,i-m+L_m) \]

   因此更新 \(L'=\max(L,i+1-L)\)。

   当 \(L'\neq L\) 时，更新：

   \[B(x)\leftarrow C(x)\\ C(x)\leftarrow B(x)\\ L\leftarrow i+1-L\\ m\leftarrow i \]

   否则就不更新，毕竟我们要的是长度最小的 \(B(x)\)。

   时间复杂度 \(O(nL+L^2)\)，其中 \(L\) 是最终递推长度。

vector<int> BM(vector<int> f){
    vector<int>C={1},B={1};
    int L=0,m=-1,b=1;
    for(int i=0;i<f.size();++i){
        int d=0;
        for(int j=0;j<C.size()&&j<=i;++j)d+=f[i-j]*C[j]%p;
        d%=p;
        if(d==0)continue;
        vector<int>tmp=C;
        int w=d*power(b,p-2)%p,dt=i-m;
        if(C.size()<B.size()+dt)C.resize(B.size()+dt,0);
        for(int j=0;j<B.size();++j)(C[j+i-m]-=w*B[j]%p)%=p;
        if((L<<1)<=i){
            L=i+1-L;m=i;
            B=tmp;b=d;
        }
    }
    vector<int>res;
    for(int i=1;i<C.size();++i)res.push_back((p-C[i])%p);
    while(!res.empty()&&res.back()==0)res.pop_back();
    return res;
}