Berlekamp-Massey算法

Posted 2020-11-22 yuanquming

(BM)算法可以在(O(n^2))的时间里用来求出一个长度为(n)的数列的最短递推式

用处是在题目中打出小范围的表之后求出递推式并配合CH定理来求出最终的答案

以下无特殊说明时均默认下标从(1)开始，用(|B|)表示数列(B)的长度

算法流程

对于某个长为(n)的数列({a_i})，我们称某个长为(m)的数列({b_i})为其递推式，当且仅当对于任意(m+1leq ileq n)，有(a_i=sumlimits_{j=1}^m b_ja_{i-j})恒成立。称({b_i})为其最短线性递推式当且仅当(m)是所有递推式中最小的。注意，(mgeq n)的情况也是允许的，此时由于不存在(m+1leq ileq n)，所以此时必定是(a)的线性递推式

我们考虑增量法，设当前已经求出了(a_{1,...,i-1})的最短线性递推式，考虑加入(a_i)之后的最短线性递推式

记初始时的递推式为(B_0=varnothing)，假设递推式被修改了(k)次，且第(i)次修改后的递推式为(B_i)，那么当前的线性递推式就是(B_k)

设当前递推式为(B_k)，记(d_i=a_i-sumlimits_{j=1}^m b_ja_{i-j})，如果(d_i=0)，那么当前递推式已经合法，可以退出

否则，我们认为(B_k)在(a_i)处出错，记(fail_i)为(B_i)最早出错的位置，则有(fail_k=i)。我们尝试对(B_k)进行修改，使之变为(B_{k+1})，并其是(a_{1,...,i})的最短递推式

如果(k=0)，那么我们可以令(B_1={0,0,...,0})，即用(i)个(0)填充，那么(B_1)显然是(a_{1,..,i})的线性递推式。并且由于(a_i)是序列中第一个非零元素，所以容易证明(B_1)是最短递推式

否则(k>0)，我们记(id)为(0)到(k-1)中的任意一个元素，并记(coef={d_iover d_{fail_{id}}})

如果我们能构造出一个长度为(q)的序列(B')，使得(sum_{j=1}^{q}b_{j}a_{k-j}=0)对任意(q+1leq ileq k-1)成立，且(sum_{j=1}^{q}b_{j}a_{i-j}=d_i)，那么我们令(B_{k+1}=B_k+B')即可（这里加法定义为按位相加）

对于(B')，可以是({0,0,...,0,coef,-coef imes B_{id}})，即先填上(i-fail_{id}-1)个(0)，然后把({1,-B_{id}})乘上(coef)倍接在后面。手玩一下发现它满足我们的条件，那么令(B_{k+1}=B_k+B')即可

由于我们需要最短线性递推式，而(|B_{k+1}|=max(|B_k|,i-fail_{id}+|B_{id}|))，那么我们取(|B_{id}|-fail_{id})最小的即可，这样我们就可以只记录当前的和最小的这两个线性递推式了

由于最坏情况下会修改(O(n))次，总复杂度为(O(n^2))

洛谷5487，注释掉的是(NTT)优化(CH)，没注释掉的是暴力实现(CH)

//quming
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
const int P=998244353;
inline void upd(R int &x,R int y){(x+=y)>=P?x-=P:0;}
inline int inc(R int x,R int y){return x+y>=P?x+y-P:x+y;}
inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
int ksm(R int x,R int y){
    R int res=1;
    for(;y;y>>=1,x=mul(x,x))(y&1)?res=mul(res,x):0;
    return res;
}
const int N=(1<<16)+5;
//namespace CH{
//  int rt[2][N],r[18][N],inv[18],lg[N],md[N],A[N],B[N],a[N],b[N],lim,d,n,K;
//  void init(){
//      fp(d,1,16){
//          fp(i,1,(1<<d)-1)r[d][i]=(r[d][i>>1]>>1)|((i&1)<<(d-1));
//          inv[d]=ksm(1<<d,P-2),lg[1<<d]=d;
//      }
//      for(R int t=(P-1)>>1,i=1,x,y;i<65536;t>>=1,i<<=1){
//          x=ksm(3,t),y=ksm(332748118,t),rt[0][i]=rt[1][i]=1;
//          fp(k,1,i-1){
//              rt[0][i+k]=mul(rt[0][i+k-1],x);
//              rt[1][i+k]=mul(rt[1][i+k-1],y);
//          }
//      }
//  }
//  void NTT(int *A,int ty){
//      int t;
//      fp(i,0,lim-1)if(i<r[d][i])swap(A[i],A[r[d][i]]);
//      for(R int mid=1;mid<lim;mid<<=1)
//          for(R int j=0;j<lim;j+=(mid<<1))
//              fp(k,0,mid-1){
//                  A[j+k+mid]=inc(A[j+k],P-(t=mul(A[j+k+mid],rt[ty][mid+k])));
//                  upd(A[j+k],t);
//              }
//      if(!ty){
//          t=inv[d];
//          fp(i,0,lim-1)A[i]=mul(A[i],t);
//      }
//  }
//  void Inv(int *a,int *b,int len){
//      if(len==1)return b[0]=ksm(a[0],P-2),void();
//      Inv(a,b,len>>1);
//      static int A[N],B[N];lim=(len<<1),d=lg[lim];
//      fp(i,0,len-1)A[i]=a[i],B[i]=b[i];
//      fp(i,len,lim-1)A[i]=B[i]=0;
//      NTT(A,1),NTT(B,1);
//      fp(i,0,lim-1)A[i]=mul(A[i],mul(B[i],B[i]));
//      NTT(A,0);
//      fp(i,0,len-1)upd(b[i],inc(b[i],P-A[i]));
//      fp(i,len,lim-1)b[i]=0;
//  }
//  void Mod(int *a,int *b,int *c,int n,int m){
//      while(!a[n-1])--n;
//      if(n<m){
//          fp(i,0,n-1)c[i]=a[i];fp(i,n,m-2)c[i]=0;
//          return;
//      }
//      static int A[N],B[N],IB[N],C[N];
//      R int len=1;while(len<=n-m)len<<=1;
//      fp(i,0,n-1)A[i]=a[n-i-1];fp(i,0,m-1)B[i]=b[m-i-1];
//      fp(i,m,len-1)B[i]=0;Inv(B,IB,len);
//      lim=(len<<1),d=lg[lim]; 
//      fp(i,n-m+1,lim-1)A[i]=IB[i]=0;
//      NTT(A,1),NTT(IB,1);
//      fp(i,0,lim-1)A[i]=mul(A[i],IB[i]);
//      NTT(A,0);
//      lim=1;while(lim<n)lim<<=1;d=lg[lim];
//      fp(i,0,n-m)C[i]=A[n-m-i];fp(i,n-m+1,lim-1)C[i]=0;
//      fp(i,0,m-1)B[i]=b[i];fp(i,m,lim-1)B[i]=0;
//      NTT(B,1),NTT(C,1);
//      fp(i,0,lim-1)B[i]=mul(B[i],C[i]);
//      NTT(B,0);
//      fp(i,0,m-2)c[i]=inc(a[i],P-B[i]);
//      fp(i,m-1,lim-1)c[i]=0;
//  }
//  void Mul(int *a,int *b,int *c,int n,int m){
//      static int A[N],B[N];
//      lim=1;while(lim<(n+m))lim<<=1;d=lg[lim];
//      fp(i,0,n-1)A[i]=a[i];fp(i,0,m-1)B[i]=b[i];
//      fp(i,n,lim-1)A[i]=0;fp(i,m,lim-1)B[i]=0;
//      NTT(A,1),NTT(B,1);
//      fp(i,0,lim-1)A[i]=mul(A[i],B[i]);
//      NTT(A,0);
//      fp(i,n+m-1,lim-1)A[i]=0;
//      Mod(A,md,c,n+m-1,K+1);
//  }
//  void ksm(int y){
//      R int sz=2,psz=1;
//      A[1]=1,B[0]=1;
//      for(;y;y>>=1,Mul(A,A,A,sz,sz),sz=sz+sz-1,cmin(sz,K))
//          if(y&1)Mul(A,B,B,sz,psz),psz+=sz-1,cmin(psz,K);
//  }
//  void MAIN(){
//      init();
//      md[K]=1;fp(i,0,K-1)md[i]=inc(0,P-a[K-i]);
//      ksm(n);
//      R int res=0;
//      fp(i,0,K-1)upd(res,mul(B[i],b[i]));
//      printf("%d
",res);
//  }
//  inline void fr(int nn,int kk,int *A,int *B){n=nn,K=kk;fp(i,0,K-1)a[i+1]=A[i+1],b[i]=B[i];}
//}
namespace CH{
    typedef vector<int> poly;
    int a[N],b[N],K,n,res;poly c,d,md;
    poly operator %(R poly a,R poly b){
        R int n=a.size(),m=b.size(),t,iv=inc(0,P-ksm(b[m-1],P-2));
        fd(i,n-1,m-1)if(a[i]){
            t=mul(a[i],iv);
            fp(j,0,m-1)upd(a[i-j],mul(t,b[m-1-j]));
        }
        while(!a.empty()&&!a.back())a.pop_back();
        return a;
    }
    poly operator *(R poly a,R poly b){
        R int n=a.size(),m=b.size();poly c(n+m-1);
        fp(i,0,n-1)fp(j,0,m-1)upd(c[i+j],mul(a[i],b[j]));
        return c%md;
    }
    void MAIN(){
        md.resize(K+1);md[K]=1;fp(i,0,K-1)md[i]=inc(0,P-a[K-i]);
        c.resize(1),d.resize(2),c[0]=1,d[1]=1;
        for(;n;n>>=1,d=d*d)if(n&1)c=c*d;
        res=0;
        fp(i,0,K-1)upd(res,mul(c[i],b[i]));
        printf("%d
",res);
    }
    inline void fr(int nn,int kk,int *A,int *B){n=nn,K=kk;fp(i,0,K-1)a[i+1]=A[i+1],b[i]=B[i];}
}
namespace BM{
    vector<int>now,las,ret;int a[N],del[N],A[N],B[N],fail,mn,sz,n,m;
    void MAIN(){
        scanf("%d%d",&n,&m);
        fp(i,1,n)scanf("%d",&a[i]),upd(a[i],P);
        fail=0;
        fp(i,1,n){
            if(!fail){
                if(a[i])now.resize(i,0),del[i]=a[i],fail=i,mn=-i,sz=i;
                continue;
            }
            R int sum=a[i];
            fp(j,0,sz-1)upd(sum,P-mul(now[j],a[i-j-1]));
            del[i]=sum;
            if(!sum)continue;
            ret=now;
            R int tsz=i+mn,coef=mul(del[i],ksm(del[fail],P-2));
            if(sz<tsz)now.resize(tsz);
            upd(now[i-fail-1],coef);
            fp(j,0,las.size()-1)upd(now[i-fail+j],P-mul(coef,las[j]));
            if(mn>=sz-i)mn=sz-i,las=ret,fail=i;
            cmax(sz,tsz);
        }
        fp(i,0,sz-1)printf("%d%c",now[i]," 
"[i==sz-1]);
        fp(i,0,sz-1)A[i+1]=now[i],B[i]=a[i+1];
        CH::fr(m,sz,A,B);
    }
}
int main(){
//  freopen("testdata.in","r",stdin);
    BM::MAIN(),CH::MAIN();
    return 0;
}

参考文献

https://blog.csdn.net/qq_39972971/article/details/80725873

https://www.cnblogs.com/zhouzhendong/p/Berlekamp-Massey.html