# Berlekamp-Massey算法简单介绍

Berlekamp-Massey算法是用来求一个数列的递推式的。

int n,pn=0,fail[2333];
typedef vector<ld> vld;
vld ps[2333]; ld x[2333],delta[2333];
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++) scanf("%lf",x+i);
int best=0;
for(int i=1;i<=n;i++)
{
ld dt=-x[i];
for(int j=0;j<ps[pn].size();j++)
dt+=x[i-j-1]*ps[pn][j];
delta[i]=dt;
if(fabs(dt)<=1e-7) continue;
fail[pn]=i; if(!pn) {ps[++pn].resize(i); continue;}
vld&ls=ps[best]; ld k=-dt/delta[fail[best]];
vld cur; cur.resize(i-fail[best]-1); //trailing 0
cur.pb(-k); for(int j=0;j<ls.size();j++) cur.pb(ls[j]*k);
if(cur.size()<ps[pn].size()) cur.resize(ps[pn].size());
for(int j=0;j<ps[pn].size();j++) cur[j]+=ps[pn][j];
if(i-fail[best]+(int)ps[best].size()>=ps[pn].size()) best=pn;
ps[++pn]=cur;
}
for(unsigned g=0;g<ps[pn].size();g++)
cout<<ps[pn][g]<<" "; puts("");
}

const int MOD=1e9+7;
ll qp(ll a,ll b)
{
ll x=1; a%=MOD;
while(b)
{
if(b&1) x=x*a%MOD;
a=a*a%MOD; b>>=1;
}
return x;
}
namespace linear_seq {
inline vector<int> BM(vector<int> x)
{
vector<int> ls,cur;
int pn=0,lf,ld;
for(int i=0;i<int(x.size());++i)
{
ll t=-x[i]%MOD;
for(int j=0;j<int(cur.size());++j)
t=(t+x[i-j-1]*(ll)cur[j])%MOD;
if(!t) continue;
if(!cur.size())
{cur.resize(i+1); lf=i; ld=t; continue;}
ll k=-t*qp(ld,MOD-2)%MOD;
vector<int> c(i-lf-1); c.pb(-k);
for(int j=0;j<int(ls.size());++j) c.pb(ls[j]*k%MOD);
if(c.size()<cur.size()) c.resize(cur.size());
for(int j=0;j<int(cur.size());++j)
c[j]=(c[j]+cur[j])%MOD;
if(i-lf+(int)ls.size()>=(int)cur.size())
ls=cur,lf=i,ld=t;
cur=c;
}
vector<int>&o=cur;
for(int i=0;i<int(o.size());++i)
o[i]=(o[i]%MOD+MOD)%MOD;
return o;
}
int N; ll a[SZ],h[SZ],t_[SZ],s[SZ],t[SZ];
inline void mull(ll*p,ll*q)
{
for(int i=0;i<N+N;++i) t_[i]=0;
for(int i=0;i<N;++i) if(p[i])
for(int j=0;j<N;++j)
t_[i+j]=(t_[i+j]+p[i]*q[j])%MOD;
for(int i=N+N-1;i>=N;--i) if(t_[i])
for(int j=N-1;~j;--j)
t_[i-j-1]=(t_[i-j-1]+t_[i]*h[j])%MOD;
for(int i=0;i<N;++i) p[i]=t_[i];
}
inline ll calc(ll K)
{
for(int i=N;~i;--i) s[i]=t[i]=0;
s[0]=1; if(N!=1) t[1]=1; else t[0]=h[0];
for(;K;mull(t,t),K>>=1) if(K&1) mull(s,t); ll su=0;
for(int i=0;i<N;++i) su=(su+s[i]*a[i])%MOD;
return (su%MOD+MOD)%MOD;
}
inline int gao(vector<int> x,ll n)
{
if(n<int(x.size())) return x[n];
vector<int> v=BM(x); N=v.size(); if(!N) return 0;
for(int i=0;i<N;++i) h[i]=v[i],a[i]=x[i];
return calc(n);
}
}

posted @ 2017-05-19 11:10  fjzzq2002  阅读(15095)  评论(2编辑  收藏  举报