This is not the end.

0%

Berlekamp-Massey算法

Posted on

前言

如果现在你有一个序列，你知道他们满足$m$阶的

递推关系，而你需要通过这个序列中的数解出这个序列的递推式。这就是BM算法解决的问题。

形式化地来说，有一个序列$a$，你需要找出一组$\{ r_1,r_2,r_3\cdots r_m \} $，满足$\forall i>m ,a_i = \sum_{1\le j\le m} a_{i-j}r_j$。

过程

BM算法利用了增量法。我们先求出一个对于前$i$项成立的$R_c=\{r_1,r_2,\cdots r_m \}$，然后把第$i+1$项带入。如果此时递推关系仍然成立，我们就继续去看$i+2$项；否则，我们就通过一系列调整，得到一个$R_{c+1}$，对于前$i+1$项都成立，然后我们再去看后面的项。

记我们得到的第$c$个递推关系式为$R_c$，记第$c$个递推关系式第一次不成立的位置为$fail_c$。初始的时候递推关系为空。

现在我们解决如何调整$R$使得递推关系对第$i$项成立。我们记$delta_i = a_i - \sum_{j=1}^m r_j a_{i-j}$。考虑另构造一个递推关系$R’$满足得到的序列在前$i-1$项取值为$0$，在第$i$项取值为$1$，这样一个对前$i$项成立的就可以是$R_c + delta_{i}\cdot R’$。

考虑利用我们之前的修正。对于$R_{d}$这个递推关系式，它满足在前$fail_d-1$项的取值与$a$相同，在第$fail_d$这一项的取值是$delta_{fail_d}$。那么我们可以构造一个递推关系$R’=\{ 0,0,0,\cdots 0,1,-R_d\} $，这个递推关系的前$i-fail_d-1$项是$0$，然后的一项是$1$，接下来是整个$R_d$取反。这个递推关系的长度是$i-fail_d + |R_d|$的。我们构造出来的这个$R’$，对于前$i-1$个位置中的任意一个位置$x$，取值是$a_{x-(i-fail_d)}- \sum_j r_j a_{x-(i-fail_d)-j} = delta_{x-(i-fail_d)}$。由于$x<i$，所以$x-(i-fail_d)<fail_d$，所以这个递推式在$x$处的取值是$0$。考虑这个递推式在第$i$项的取值：$a_{i-(i-fail_d)} -\sum_j r_j a_{i-(i-fail_d)} = delta_{fail_d}$。我们只需要把这个递推式的每一项的系数再除以$delta_{fail_d}$就可以得到一个在第$i$项取值为$1$，在其他项取值为$0$的递推关系式。

概括一下我们调整的方法：

  • 在$R_1,R_2,\cdots R_{c-1}$中选出一个$R_d$。然后构造出$R’ = \{ 0,0,0\cdots 1,-R_{d}\}$。
  • 将$R_{c+1}$调整为$R_c + {delta_i\over delta_{fail_d}} R’$。
  • 调整后，递推关系的长度是$\max \{ |R_c|, i-fail_d+|R_d|\}$。

这里关于$d$的选择，有些资料里面是直接选择$c-1$，也有的地方说应该选择$i-fail_d +|R_d|$，我也没有搞清楚哪种写法是正确的；很多地方都说BM算法解出的是最短线性递推式，然而我也没有弄懂它为什么是最短的。等以后弄懂了再回来补充吧。

时间复杂度

如果要解出一个$m$阶的递推式，我们至少需要$2m$项。我们每一次调整的时候，我们需要修改$m$个位置。所以总复杂度大约是$O(m^2)$的。

模板

bzoj1494生成树计数
一种做法就是暴力行列式求解打表 + 解线性递推式

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
namespace BM {
	vector<int> b[N];
	int delta[N],fail[N];
	void solve(int *a,int n,vector<int> &c) {
		int cnt=0;
		for(int i=1;i<=n;++i) {
			delta[i]=a[i];
			for(int j=0;j<b[cnt].size();++j)
				Dec(delta[i],b[cnt][j]*(ll)a[i-j-1]%mod);
			if(delta[i]==0) continue;
			fail[cnt++]=i;
			if(cnt==1) { b[cnt].resize(i); continue; }
			int id=0,len=1e9;
			for(int k=1;k<cnt-1;++k)
				if(i-fail[k]+(int)b[k].size()<len)
					id=k,len=i-fail[k]+(int)b[k].size();
			int tmp=delta[i]*(ll)Pow(delta[fail[id]],mod-2)%mod;
			b[cnt].resize(i-fail[id]-1); b[cnt].PB(tmp);
			for(int j=0;j<b[id].size();++j)
				b[cnt].PB(b[id][j]*(ll)(mod-tmp)%mod);
			if(b[cnt].size()<b[cnt-1].size()) b[cnt].resize(b[cnt-1].size());
			for(int j=0;j<b[cnt-1].size();++j)
				Add(b[cnt][j],b[cnt-1][j]);
		}
		c=b[cnt];
	}
}

Upd 2020.5.7

参考资料：zzq2019集训队论文

对于有限数列 $\{a_0,a_1,\cdots a_{n-1} \}$，设前缀 $\{a_0,a_1,\cdots a_i\}$ 的最短递推式为 $r^{(i)} = \{r^{(i)}_1,r^{(i)}_2,\cdots \mid a_j = \sum_{k=1}^{|r^{(i)}|} r^{(i)}_k a_{j-k}, j\le i\}$，设 $l_i = |r^{(i)}|$。显然有 $l_{i-1} \le l_i$。

lemma 1.1

如果 $r^{(i-1)}$ 不是 $\{a_0,a_1,\cdots a_{i-1}\}$ 的最短递推式，那么 $l_i \ge \max\{l_{i-1},i+1-l_{i-1}\}$。

证明：假设 $l_i \le i-l_{i-1}$，$p = r^{(i-1)}, q=r^{(i)}$则

这一步是因为 $i-j \ge i-l_{i-1} \ge l_i$。

这一步是因为 $i-k \ge i-l_i \ge l_{i-1}$。

则 $r^{(i-1)}$ 对 $\{a_0,a_1,\cdots a_i\}$ 也成立，与已知矛盾。

这个引理告诉了我们最短递推式长度的下界，而这个下界是可以用 Berlekamp-Massey 算法达到的。

Berlekamp-Massey 算法的过程：

将 $i$ 从 $1$ 枚举到 $n$：

  • 如果 $r^{(i-1)}$ 对 $\{a_0,a_1,\cdots a_i\}$ 仍然成立，令 $r^{(i)} = r^{(i-1)}$ 即可
  • 否则，用上一次 $r^{(p-1)}$ 对 $\{a_0,a_1,\cdots a_{p}\}$ 不再成立的 $r^{(p-1)}$ 来调整 $r^{(i-1)}$ 以得到 $r^{(i)}$

归纳证明其达到了下界：

  • 假设 $l_1,l_2,\cdots l_{i-1}$ 均已达到了下界，而 $r^{(i)} \neq r^{(i-1)}$
  • 由算法流程可知和归纳假设可知，$l_{i-1} = l_p, l_i = i-p+l_{p-1}, l_p = p+1 - l_{p-1}$
  • 则 $l_i = i - p + l_{p-1} = i + 1 - l_{p} = i+1-l_{i-1}$