Berlekamp-Massey 算法学习笔记

在李巨的指引下学了一发.

简称 $BM$ 算法,可以以 $O(n^2)$ 的时间复杂度求解一个长度为 $n$ 的数列的最短线性递推式.

对于数列 $\lbrace a_1,a_2,a_3,\dots,a_n \rbrace$ ,称数列 $\lbrace R_1,R_2,R_3,\dots,R_m \rbrace$ 为其线性递推式,当且仅当下式成立:

$$
\forall m+1\le i\le n,a_i=\sum_{j=1}^m R_j\cdot a_{i-j}
$$

若它是所有合法线性递推式中长度 $m$ 最小的线性递推式,则称它为数列 $\lbrace a_1,a_2,a_3,\dots,a_n \rbrace$ 的最短线性递推式.

若未特殊说明,接下来的递推式,线性递推式均指 最短线性递推式 ,

  • $BM$ 算法的主要思路是从前往后依次考虑每个数,修改当前递推式使其合法,并且满足最短.
  • 尝试由 $\lbrace a_1,a_2,a_3,\dots,a_{i-1} \rbrace$ 的递推式 $\lbrace r_1,r_2,r_3,\dots,r_m \rbrace$ 得出 $\lbrace a_1,a_2,a_3,\dots,a_i \rbrace$ 的递推式.
  • 记递推式被更改的次数为 $cnt$ ,第 $i$ 次修改后得到的递推式为 $R_i$ ,规定 $R_0$ 为空.记递推式 $\lbrace r_1,r_2,r_3,\dots,r_m \rbrace$ 为当前递推式,即 $R_{cnt}$ .
  • 记 $delta_i=a_i-\sum_{j=1}^m r_{j}\cdot a_{i-j}$ ,即用当前递推式算得的值与真实值之间的误差.若 $delta_i=0$ ,则不需要对递推式进行修改,直接考虑下一个元素 $a_{i+1}$ .
  • 否则说明 $R_{cnt}$ 在位置 $i$ 出错了,记 $R_{cnt}$ 第一次出错的位置为 $fail_i$ ,则 $fail_{cnt}=i$ .考虑对 $R_{cnt}$ 进行修改,得到 $R_{cnt+1}$ ,并使得它在位置 $i$ 也成立.
  • 若当前 $cnt=0$ ,说明 $a_i$ 是第一个非零元素,直接将 $R_1$ 置为 $\lbrace 0,0,0,\dots0 \rbrace$ ( $i$ 个 $0$ ) 即可.
  • 否则, $cnt>0$ ,说明之前已经修改过递推式,即存在 $R_k$, 记录 $p=fail_k$ ,尝试在 $R_k$ 的基础上修改,在 $a_i$ 的位置上递推出一个 $-delta_{p}$ 每个位置乘上 $\frac {delta_i} {delta_p}$ ,每个位置再加上原来的 $R_{cnt}$ 就得到合法的 $R_{cnt+1}$ .
  • 将 $R_k$ 的元素全部变成它的相反数,再在前面补上一个 $1$ , $-delta_p$ 就到 $p+1$ 位置上来了.再在前面补 $i-p-1$ 个 $0$ , $-delta_p$ 就到位置 $i$ 上来了.
  • 于是得到 $R_{cnt+1}=\frac {delta_i} {delta_p}*\lbrace 0,0,0,\dots,1,-R_{k_1},-R_{k_2},\dots,-R_{k_M} \rbrace + R_{cnt}$ . $0$ 有 $i-p-1$ 个.
  • 为了保证得到的递推式长度最短,我们需要选取恰当的 $k$ .容易看出,得到的 $R_{cnt+1}$ 的长度为 $\max(i-p+M,m)$ , $M$ 为 $R_k$ 的长度, $m$ 为 $R_{cnt}$ 的长度.于是记录 $M-p$ 最短的递推式作为 $R_k$ .
  • 最坏情况需要更新 $n$ 次递推式,时间复杂度 $O(n^2)$ .

用 $BM$ 得到的最短递推式长度最好要明显小于 $n$ 的一半,否则需要再打些表.

为什么?因为若长度为 $\frac n 2$,可以看做 $\frac n 2$ 个变量列出 $\frac n 2$ 个方程,总能找到解. 所以一个随机数列解出的最短递推式长度就是 $n$ 的一半左右,长度在 $\frac n 2$ 左右说明原数列很可能并没有一定的规律,即,大概率对之后的数据不适用.

另,因为计算中涉及 $\frac {delta_i} {delta_p}$ ,所以 $BM$ 在实数域内求解可能有一定的精度误差.若在模质数意义下,则不用考虑.

板子. $vector$ 下标从 $0$ 开始,要注意处理.

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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
int out=0,fh=1;
char jp=getchar();
while ((jp>'9'||jp<'0')&&jp!='-')
jp=getchar();
if (jp=='-')
fh=-1,jp=getchar();
while (jp>='0'&&jp<='9')
out=out*10+jp-'0',jp=getchar();
return out*fh;
}
const int P=1e9+7;
const int inf=P+1;
inline int add(int a,int b)
{
return (a+b)%P;
}
inline int sub(int a,int b)
{
return (a-b+P)%P;
}
inline int mul(int a,int b)
{
return 1LL * a * b % P;
}
int fpow(int a,int b)
{
int res=1;
while(b)
{
if(b&1)
res=mul(res,a);
a=mul(a,a);
b>>=1;
}
return res;
}
int inv(int x)
{
return fpow(x,P-2);
}
#define len(A) A.size()
typedef vector<int> poly;
poly BerlekampMassey(poly a)
{
poly R_k,R;
int n=a.size();
int p=-inf,Delta;
int fail,delta;
for(int i=0;i<n;++i)
{
delta=a[i];
for(int j=0;j<len(R);++j)
delta=sub(delta,mul(R[j],a[i-(j+1)]));
if(!delta)
continue;
fail=i;
if(p==-inf)
{
R.resize(i+1);
p=fail;
Delta=delta;
continue;
}
poly r;
int tmp=mul(delta,inv(Delta));
for(int j=1;j<=i-p-1;++j)
r.push_back(0);
r.push_back(tmp);
for(int j=0;j<len(R_k);++j)
r.push_back(sub(0,mul(tmp,R_k[j])));
if(len(r)<len(R))
r.resize(len(R));
for(int j=0;j<len(R);++j)
r[j]=add(r[j],R[j]);
if(len(R)+p<len(R_k)+fail)
{
R_k=R;
Delta=delta;
p=fail;
}
R=r;
}
return R;
}
int main()
{
int n=read();
poly a;
for(int i=1;i<=n;++i)
a.push_back(read());
poly R=BerlekampMassey(a);
int m=R.size();
printf("%u\n",m);
for(int i=0;i<m;++i)
printf("%d ",R[i]);
puts("");
for(int i=m;i<n;++i)
{
int delta=a[i];
for(int j=0;j<m;++j)
delta=sub(delta,mul(R[j],a[i-(j+1)]));
assert(!delta);
}
return 0;
}