Berlekamp-Massey算法学习笔记
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\not=0\)
所以是不可能用之前的数推出\(a_i\)的
这种情况下直接把\(f_0,f_1,...f_i\)记为\(0\)
这样\(i\le 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\n",G.size());
for(int i=0;i<G.size();++i)printf("%d ",G[i]);
}
然后遇到一个题就可以Berlekamp-Massey+矩阵快速幂水过
当然多项式取模优化线性递推更好
毕竟Berlekamp-Massey的复杂度是\(O(k^2)\)
而矩阵快速幂的复杂度是\(O(k^3\log n)\)
综合代码效率和实现难度
和\(O(k^2\log 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\n",BerlekampMassey::main(vector<int>(F,F+(sizeof(F))/4),n));
}
