Berlekamp-Massey算法学习笔记
声明:博主已退役,这是以前的总结,如有错误望指正,如有问题
不妨看看别人的博客
问题描述
给一个序列({a_0,a_1,...,a_n})
要求最短的序列({f_0,f_1,...,f_m})
使得对于所有(i>m)有
算法流程
主要思想是依次考虑每个(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));
}