Berlekamp_Massey 算法 (BM算法) 学习笔记

原文链接www.cnblogs.com/zhouzhendong/p/Berlekamp-Massey.html

前言

BM算法用于求解常系数线性递推式。

它可以在 $O(n^2)$ 的时间复杂度内解决问题。

由于许多问题会涉及线性递推，所以 BM 算法将会有不错的应用。

问题模型

给定一个有 $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)$ 的时间复杂度内解决此问题。

BM算法

　　考虑增量法。

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

　　考虑在当前数列末尾加入 $a_i$ 。假设当前递推式长度为 $m$ 。

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

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

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

　　如果 $c = 0$ ，说明 $a_i$ 之前都是 0 ，显然新的递推式由 $i$ 个 $0$ 组成。

　　考虑 $n\neq 0$ 的情况：考虑构造一个递推式 $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'$ 符合要求。

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

　　至此，我们可以在 $O(n^2)$ 的时间复杂度内，求出数列 $a_i$ 的一个较短线性递推式。

　　那么如何求最短的线性递推式呢？

　　只要在对 $id$ 取值时，每次找 $i - Fail[id] + len(R[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;
}

关于模板的测试

到这篇博文写完为止，各大OJ似乎并没有BM算法的模板题。因此这里说明两个测试数据来源：

1. cz_xuyixuan 的博客中的例子、评论中的数据：

Input 1
7
1 2 4 9 20 40 90

Output 1
4 
0 0 10 0

Input 2 
18
2 4 8 16 32 64 128 256 512 2 4 8 16 32 64 128 256 512

Output 2 
0 0 0 0 0 0 0 0 1

2. fjzzq2002 的博客中给出的数据。

相应博客链接见“参考文献”。

参考文献

https://blog.csdn.net/qq_39972971/article/details/80725873

http://www.cnblogs.com/zzqsblog/p/6877339.html

 

转载于:https://www.cnblogs.com/zhouzhendong/p/Berlekamp-Massey.html