Berlekamp–Massey 算法

线性递推式

对于数列 \{a_0,a_1 \dots\} ，存在有限数列 \{r_0,r_1\dots r_{m-1}\} ，使得 p\ge m-1,\sum\limits_{i=0}^{m-1}r_ia_{p-i}=0 ，那么我们称数列 r 为数列 a 的线性递归式，若 r_0=1 ，则称 r 为线性递推式，称其阶数为 m-1 。

容易发现，如果设 r 的生成函数为 R ，设 a 的生成函数为 A ，设 S=AR ，则 S 的次数不超过 m-2 。

Berlekamp-Massey 算法

已知数列 \{a_0,a_1\dots a_{n-1}\} ，Berlekamp-Massey 算法会对其每个前缀 \{a_0,a_1\dots a_i\} 求出这一前缀的最短线性递推式 r^{(i)} ，我们设 |r^{(i)}|-1=l_i ，即 l_i 为线性递推式 r^{(i)} 的阶，显然有 l_{i-1}\le l_i 。我们令 r^{(-1)}=\{1\} 。

算法流程

引理：若 r^{(i-1)} 不是 \{a_0,a_1,a_2\dots a_i\} 的线性递推式，那么有 l_i\ge i+1-l_{i-1} 。

证明：反证法，设 l_i\le i-l_{i-1} ，设 r^{(i-1)}=\{p_0,p_1\dots p_{l_{i-1}}\},r^{(i)}=\{q_0,q_1\dots q_{l_i}\} 。

我们有 l_i\le j\le i,a_j=-\sum\limits_{t=1}^{l_i}q_ta_{j-t} 。

由于 l_i\le i-l_{i-1} ，我们可以做出如下变换。

\begin{aligned} -\sum\limits_{j=1}^{l_{i-1}}p_ja_{i-j}&=\sum\limits_{j=1}^{l_{i-1}}p_j\sum\limits_{k=1}^{l_i}q_ka_{i-j-k}\\ &=\sum\limits_{k=1}^{l_i}q_k\sum\limits_{j=1}^{l_{i-1}}p_ja_{i-j-k}\\ &=-\sum\limits_{k=1}^{l_i}q_ka_{i-k}\\ &=a_i \end{aligned}

所以 r^{(i-1)} 是 \{a_0,a_1,a_2\dots a_i\} 的线性递推式，矛盾。

于是，我们有 l_i\ge \max(i+1-l_{i-1},l_{i-1}) 。

下面给出一种构造使得 l_i=\max(i+1-l_{i-1},l_{i-1}) ，显然这样的 r^{(i)} 即为 \{a_0,a_1,a_2\dots a_i\} 的最短线性递推式。

如果 r^{(i-1)} 是 \{a_0,a_1,a_2\dots a_i\} 的线性递推式， r^{(i)}=r^{(i-1)} 。

否则，我们设 a 的生成函数为 A ， r^{(i)} 的生成函数为 R^{(i)} ，设 S^{(i)}\equiv AR^{(i)}\pmod{x^{i+1}} ，则有 S^{(i)} 的次数不超过 l_i-1 ，即 R^{(i)} 的次数 -1 。

考虑由 R^{(i-1)} 推出 R^{(i)} ，由于 r^{(i-1)} 不是 \{a_0,a_1,a_2\dots a_i\} 的线性递推式，我们有 AR^{(i-1)}\equiv S^{(i-1)}+dx^i\pmod{x^{i+1}} 。

若 0\le j< i,l_j=0 ，我们令 r^{(i)}=\{1,0\dots 0\} ，此处有 i+1 个 0 ，显然满足 l_i=i+1-l_{i-1} 。

否则我们考虑上一次 l_p>l_{p-1} 的情形，设当时的情况为 AR^{(p-1)}\equiv S^{(p-1)}+cx^p\pmod{x^{p+1}} ，给两侧分别乘上 x^{i-p}dc^{-1} ，那么有 x^{i-p}dc^{-1}AR^{(p-1)}\equiv x^{i-p}dc^{-1}S^{(p-1)}+dx^i\pmod {x^{i+1}} 。

两式相减可以得到 A(R^{(i-1)}-x^{i-p}dc^{-1}R^{(p-1)})\equiv S^{(i-1)}-x^{i-p}dc^{-1}S^{(p-1)}\pmod {x^{i+1}} 。

显然 S^{(i-1)}-x^{i-p}dc^{-1}S^{(p-1)} 的次数不超过 R^{(i-1)}-x^{i-p}dc^{-1}R^{(p-1)} 的次数 -1 。

所以我们令 R^{(i)}=R^{(i-1)}-x^{i-p}dc^{-1}R^{(p-1)} 即可。

由归纳法，我们设 l_p=\max(l_{p-1},p+1-l_{p-1}) ，由于 l_p>l_{p-1} ，所以 l_p=p+1-l_{p-1} ，那么 l_i=\max(l_{i-1},i-p+l_{p-1})=\max(l_{i-1},i+1-l_p)=\max(l_{i-1},i+1-l_{i-1}) 。

时间复杂度 O(n^2) 。

引理：对于一个无限数列，若我们已知其线性递推式阶数不超过 s ，则只需取 \{a_0,a_1\dots a_{2s-1}\} 的最短线性递推式即可。

证明：反证，令 t 为最小的满足 t\ge 2s,r^{(t)}\ne r^{(2s-1)} 的数，由上一个引理，我们就有 l_t\ge t+1-l_{2s-1}\ge t+1-s\ge 2s+1-s\ge s+1 ，矛盾。

扩展：关于线性递推数列的远处单点求值

设数列 \{a_0,a_1\dots \} 的线性递推式为 \{r_0,r_1\dots r_{m-1}\} ，我们现在要求解 a_n 。

我们需要已知 a_0,a_1\dots a_{m-2} 。

设 F(x)=\sum\limits_{i=0}^{m-1}r_ix^{m-1-i} ， G(x)=x^n\bmod F(x)=\sum\limits_{i=0}^{m-2}g_ix^i 。

则答案为 \sum\limits_{i=0}^{m-2}g_ia_i 。

时间复杂度 O(m\log m\log n) 或 O(m^2\log n) ，取决于是否使用 FFT 加速卷积。

示例代码 P5487 【模板】Berlekamp–Massey 算法 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include<iostream>
using std::cin;
using std::cout;
int const mod=998244353;
int pow(int x,int y){
    int res=1;
    while(y){
        if(y&1) res=1ll*res*x%mod;
        x=1ll*x*x%mod,y>>=1;
    }
    return res;
}
int n,m,r[5010],a[10010],len;
int bm(int *a,int n,int *r){
    int len=1,p=0,clen=0,dt=0;
    static int cr[5010];
    r[0]=1;
    for(int i=0;i<n;i++){
        int d=0;
        for(int j=0;j<len;j++) d=(d+1ll*a[i-j]*r[j])%mod;
        if(d==0) continue;
        if(len==1){
            len=i+2,p=i,cr[0]=1,clen=1,dt=d;
            continue;
        }
        int u=1ll*d*pow(dt,mod-2)%mod;
        if(i-p+clen>len){
            static int pr[5010];
            for(int i=0;i<len;i++) pr[i]=r[i];
            for(int j=0;j<clen;j++) pr[i-p+j]=(pr[i-p+j]-1ll*cr[j]*u%mod+mod)%mod;
            for(int i=0;i<len;i++) cr[i]=r[i];
            int tmp=clen;
            clen=len,dt=d;
            len=i-p+tmp;
            p=i;
            for(int i=0;i<len;i++) r[i]=pr[i];
        }else for(int j=0;j<clen;j++) r[i-p+j]=(r[i-p+j]-1ll*cr[j]*u%mod+mod)%mod;
    }
    return len-1;
} 
int getans(int *r,int len,int m,int *a){
    static int x[20010],res[20010];
    x[1]=1,res[0]=1;
    if(len==1) x[0]=1ll*(mod-x[1])*r[1]%mod,x[1]=0;
    int t=2*len-2;
    auto mul=[&t,&len](int *a,int *b,int *p){
        static int tmp[20010];
        for(int i=0;i<=t;i++) tmp[i]=0;
        for(int i=0;i<len;i++) for(int j=0;j<len;j++) tmp[i+j]=(tmp[i+j]+1ll*a[i]*b[j])%mod;
        for(int i=0;i<=t;i++) a[i]=tmp[i];
        for(int i=t;i>=len;i--){
            int u=(mod-a[i])%mod;
            a[i]=0;
            for(int j=1;j<=len;j++) a[i-j]=(a[i-j]+1ll*u*p[j])%mod;
        }
    };
    while(m){
        if(m&1) mul(res,x,r);
        mul(x,x,r);
        m>>=1;
    }
    int ans=0;
    for(int i=0;i<len;i++) ans=(ans+1ll*a[i]*res[i])%mod;
    return ans;
}
int main(){
    std::ios::sync_with_stdio(false),cin.tie(nullptr);
    cin>>n>>m;
    for(int i=0;i<n;i++) cin>>a[i];
    len=bm(a,n,r);
    for(int i=1;i<=len;i++) cout<<(mod-r[i])%mod<<' ';
    cout<<'\n'<<getans(r,len,m,a)<<'\n';
    return 0;
}
build本页面最近更新：，更新历史
edit发现错误？想一起完善？ 在 GitHub 上编辑此页！
people本页面贡献者：AFOI-wiki
copyright本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用

评论