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][j]$ 。
如果 $delta_i = 0$ ,那么递推式 $R[c]$ 依然合法,不用修改。
否则,设 $Fail_{c} = i$ 表示递推式 $R[c]$ 第一次失效的位置为 $i$ 。
如果 $c = 0$ ,说明 $a_i$ 之前都是 0 ,显然新的递推式由 $i$ 个 $0$ 组成。
考虑 $c\neq 0$ 的情况:考虑构造一个递推式 $R'$ 使得对于 $|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