Berlekamp-Massey算法学习笔记

Berlekamp-Massey算法学习笔记

声明:博主已退役,这是以前的总结,如有错误望指正,如有问题不妨看看别人的博客

问题描述

给一个序列({a_0,a_1,...,a_n})

要求最短的序列({f_0,f_1,...,f_m})

使得对于所有(i>m)有

[a_i=sum_{k=0}^{m}f_k*a_{i-1-k} ]

算法流程

主要思想是依次考虑每个(a_i)

不断修改({f})

使得其在保证正确的同时尽量短

一开始({f})为空

对每个(a_i)判断当前递推式是否满足条件

如果满足就直接判断下一个

否则需要修改

如果当前({f})为空

说明(a_i)之前数全是(0)

而(a_i ot=0)

所以是不可能用之前的数推出(a_i)的

这种情况下直接把(f_0,f_1,...f_i)记为(0)

这样(ile m)就不用推了

否则说明之前已经有一个错误的({f})

为了方便记为({F})

我们希望能通过({F})在(i)位推出一个不为(0)的值

然后把这次的错误抵消掉

如果({F})出错的位置是(p)且多了(Delta)

这个时候(a_p)一定不等于(sum_{k=0}^{m}F_k*a_{i-k})

就可以对({F})稍作修改

在(a_i)的位置上递推出一个(-Delta)

具体来说

把({F})全部变为相反数再在最前面补(1)

得到的新的({F})可以满足在(p+1)位置推出一个(-Delta)来

再在({F})最前面补(i-p-1)个(0)

(-Delta)就跑到(i)位置来了

把现在得到的({F})除以(Delta)再乘上这次的差

加上({f})就可以把这次的差抵消掉

因为要求递推式最短

我们希望每次能得到最优的({F})

考虑每次修改时({f})的长度变化

变成了(max(len(f),len(F)+i-p))

所以只要记录(len(F)-p)最短的({F})即可

还有Berlekamp-Massey返回的递推式的长度最好要在输入的一半以内

不然还是再多打点表吧

为什么？

感谢LHJ神仙指教

这样多出来的一半就可以列出至少(len(F))个方程组

确定了递推式的唯一性

代码

#include<bits/stdc++.h>

using namespace std;

#define ll long long

const int p=1e9+7;

inline int add(int a,int b){
    return (a+=b)>=p?a-p:a;
}

inline int sub(int a,int b){
    return (a-=b)<0?a+p:a;
}

inline ll qpow(ll a,ll b){
	ll ans=1;
	for(;b;b>>=1,(a*=a)%=p)if(b&1)(ans*=a)%=p;
	return ans;
}

namespace BerlekampMassey{

	typedef vector<int> poly;

	#define len(A) A.size()

	inline poly BM(const poly& A){
		poly F,F0;
		int d0,p0;
		for(int i=0;i<len(A);++i){
			int d=0;
			for(int j=0;j<len(F);++j){
				d=add(d,(ll)F[j]*A[i-j-1]%p);
			}
			d=sub(d,A[i]);
			if(!d)continue;
			if(!len(F)){
				F.resize(i+1);
				d0=d;
				p0=i;
				continue;
			}
			ll t=qpow(d0,p-2)*d%p;
			poly G(i-p0-1);
			G.push_back(t);
			t=sub(0,t);
			for(int j=0;j<len(F0);++j){
				G.push_back(F0[j]*t%p);
			}
			if(len(G)<len(F))G.resize(len(F));
			for(int j=0;j<len(F);++j){
				G[j]=add(G[j],F[j]);
			}
			if(i-p0+len(F0)>=len(F)){
                //注意这里不要把i移项,vector的size是unsigned类型,减成负的就凉了
				F0=F;
				d0=d;
				p0=i;
			}
			F=G;
		}
		return F;
	}
}
using namespace BerlekampMassey;

int F[]={1,2,4,9,20,40,90};

int main(){
	poly G(F,F+((sizeof F)>>2));
	G=BM(G);
	printf("%d
",G.size());
	for(int i=0;i<G.size();++i)printf("%d ",G[i]);
}

然后遇到一个题就可以Berlekamp-Massey+矩阵快速幂水过

当然多项式取模优化线性递推更好

毕竟Berlekamp-Massey的复杂度是(O(k^2))

而矩阵快速幂的复杂度是(O(k^3log n))

综合代码效率和实现难度

和(O(k^2log n))的多项式取模搭配最佳

虽然我觉得不太可能

但也不排除(k)比较大

必须本地跑出递推式然后上(O(k log k log n))的多项式取模的情况

#include<bits/stdc++.h>

using namespace std;

#define gc c=getchar()
#define r(x) read(x)
#define ll long long

template<typename T>
inline void read(T&x){
    x=0;T k=1;char gc;
    while(!isdigit(c)){if(c=='-')k=-1;gc;}
    while(isdigit(c)){x=x*10+c-'0';gc;}x*=k;
}

const int p=1e9+7;

inline int add(int a,int b){
	a+=b;
	if(a>=p)a-=p;
	return a;
}

inline int sub(int a,int b){
	a-=b;
	if(a<0)a+=p;
	return a;
}

inline ll qpow(ll a,ll b){
	ll ans=1;
	for(;b;b>>=1,(a*=a)%=p)if(b&1)(ans*=a)%=p;
	return ans;
}

namespace BerlekampMassey{

	typedef vector<int> poly;
	
	#define len(A) A.size()
	
	inline poly mul(const poly &A,const poly&B,const poly&F){
		poly ret(len(F)*2-1);
		for(int i=0;i<len(F);++i){
			for(int j=0;j<len(F);++j){
			    ret[i+j]=add(ret[i+j],(ll)A[i]*B[j]%p);
			}
		}
		for(int i=len(F)*2-2;i>=len(F);--i){
			for(int j=0;j<len(F);++j){
			    ret[i-j-1]=add(ret[i-j-1],(ll)ret[i]*F[j]%p);
			}
		}
		ret.resize(len(F));
		return ret;
	}
	
	inline int solve(const poly &A,const poly &F,ll n){
		if(n<len(A))return A[n];
		poly base(len(F)),ans(len(F));
		base[1]=ans[0]=1;
		for(;n;n>>=1){
			if(n&1)ans=mul(ans,base,F);
			base=mul(base,base,F);
		}
		int ret=0;
		for(int i=0;i<len(F);++i)ret=add(ret,(ll)A[i]*ans[i]%p);
		return ret;
	}
	
	inline poly BM(const poly& A){
		poly F,F0;
		int d0,p0;
		for(int i=0;i<len(A);++i){
			int d=0;
			for(int j=0;j<len(F);++j){
				d=add(d,(ll)F[j]*A[i-j-1]%p);
			}
			d=sub(d,A[i]);
			if(!d)continue;
			if(!len(F)){
				F.resize(i+1);
				d0=d;
				p0=i;
				continue;
			}
			ll t=qpow(d0,p-2)*d%p;
			poly G(i-p0-1);
			G.push_back(t);
			t=sub(0,t);
			for(int j=0;j<len(F0);++j){
				G.push_back(F0[j]*t%p);
			}
			if(len(G)<len(F))G.resize(len(F));
			for(int j=0;j<len(F);++j){
				G[j]=add(G[j],F[j]);
			}
			if(i-p0+len(F0)>=len(F)){
				F0=F;
				d0=d;
				p0=i;
			}
			F=G;
		}
		return F;
	}
	
	inline int main(const poly &A,ll n){
		return solve(A,BM(A),n);
	}
}

int F[]={0,1,1,2,3,5,8,13};

int main(){
	//freopen(".in","r",stdin);
	//freopen(".out","w",stdout);
	ll n;r(n);
	printf("%d
",BerlekampMassey::main(vector<int>(F,F+(sizeof(F))/4),n));
}