zoukankan      html  css  js  c++  java
  • BM算法线性递推

    学习BM算法正确搜索方式:

    搜索“BM算法线性递推”->随便点开一个博客,得到全名“Berlekamp-Massey算法”->复制搜索。

    其实单纯是记不住全名


    参考资料:

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

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


    概要

    给出个数列,求最短递推式。(即(a_n=sum_{i=1}^mr_ia_{n-i},nge m),下标从0开始)

    (最短递推式不一定唯一。如果知道最短递推式有(m)项,那么写出(2m)项的数列(还是(2m-1)(2m-1))才能写出唯一的最短递推式)

    思想:增量构造,如果矛盾就调整。


    流程

    增量构造。设(R_{i})表示第(i)个历史递推式,设(R_{cnt})表示当前递推式。

    现在把(a_i)带入递推式,得到变化量(Delta =a_i-sum_{j=1}^{|R_{cnt}|} r_ja_{i-j})

    如果(Delta=0),之前递推式仍然有效,继续使用。

    如果(Delta eq 0):如果之前没有其它递推式,令递推式为(i+1)(0)

    否则:

    找之前某个历史递推式(R_{id}),记其失配位置为(fail_{id}),变化量(Delta_{id})

    根据定义得(forall iin[|R_{id}|,fail_{id}-1],a_i-sum_{j=1}^{|R_{id}|}r_ja_{i-j}=0)(i=fail_{id},a_i-sum_{j=1}^{|R_{id}|}r_ja_{i-j}=Delta_{id})

    (mul=frac{Delta}{Delta_{id}})

    构造递推式(R'={0,0,dots,0,mul,-mul*R_{id,1},-mul*R_{id,2},dots})。其中前面是(i-fail_{id})(0)

    结合上面的分析,可以发现这个递推式的贡献为在(a_i)处有(Delta)

    然后(R_{cnt+1}=R_{cnt}+R')。长度为(max(|R_{cnt}|,i-fail_{id}+|R_{id}|))

    具体取哪个(id),当然是取(-fail_{id}+|R_{id}|)最小的。记下这个最小的(id)的相关信息即可。

    时间(O(n^2))


    lj洛谷模板强行将BM算法和远处系数求值套一起。

    害得我去copy了个远处系数求值的板子。

    using namespace std;
    #include<bits/stdc++.h>
    typedef long long ll;
    const int N=10005,mo=998244353;
    ll qpow(ll x,ll y=mo-2){
    	ll r=1;
    	for (;y;y>>=1,x=x*x%mo)
    		if (y&1)
    			r=r*x%mo;
    	return r;
    }
    int n;
    ll a[N];
    bool fir;
    ll mf,f[N],fail,fdelta;
    ll m,r[N];
    namespace chang{
    	const int N=65536;
    	int n,k;
    	int a[N],f[N];
    	int nN,re[N];
    	void setlen(int n){
    		int bit=0;
    		for (nN=1;nN<=n;nN<<=1,++bit);
    		for (int i=1;i<nN;++i)
    			re[i]=re[i>>1]>>1|(i&1)<<bit-1;
    	}
    	void clear(int A[],int n){memset(A,0,sizeof(int)*n);}
    	void copy(int A[],int a[],int n){
    		clear(A,nN);
    		for (int i=0;i<=n;++i)
    			A[i]=a[i];
    	}
    	void dft(int A[],int flag){
    		for (int i=0;i<nN;++i)
    			if (i<re[i])
    				swap(A[i],A[re[i]]);
    		static int wnk[N];
    		for (int i=1;i<nN;i<<=1){
    			ll wn=qpow(3,flag==1?(mo-1)/(2*i):mo-1-(mo-1)/(2*i));
    			wnk[0]=1;
    			for (int k=1;k<i;++k)
    				wnk[k]=wnk[k-1]*wn%mo;
    			for (int j=0;j<nN;j+=i<<1)
    				for (int k=0;k<i;++k){
    					ll x=A[j+k],y=(ll)A[j+k+i]*wnk[k];
    					A[j+k]=(x+y)%mo;
    					A[j+k+i]=(x-y)%mo;
    				}
    		}
    		if (flag==-1)
    			for (int i=0,invn=qpow(nN);i<nN;++i)
    				A[i]=(ll)A[i]*invn%mo;
    		for (int i=0;i<nN;++i)
    			A[i]=(A[i]+mo)%mo;
    	}
    	void multi(int c[],int a[],int b[],int n,int an=-1,int bn=-1){
    		if (an==-1) an=n-1;
    		if (bn==-1) bn=n-1;
    		static int A[N],B[N],C[N];
    		setlen(an+bn);
    		copy(A,a,an),dft(A,1);
    		if (a==b)
    			for (int i=0;i<nN;++i)
    				C[i]=(ll)A[i]*A[i]%mo;
    		else{
    			copy(B,b,bn),dft(B,1);
    			for (int i=0;i<nN;++i)
    				C[i]=(ll)A[i]*B[i]%mo;
    		}
    		dft(C,-1);
    		for (int i=0;i<=min(n-1,an+bn);++i)
    			c[i]=C[i];
    	}
    	void getinv(int c[],int a[],int n,int an=-1){
    		if (an==-1) an=n-1;
    		static int b[N],g[N];
    		int nn=1;for (;nn<n;nn<<=1);
    		clear(b,nn);
    		b[0]=qpow(a[0]);
    		for (int i=1;i<n;i<<=1){
    			clear(g,i<<1);
    			multi(g,b,b,2*i,i-1,i-1);
    			multi(g,g,a,2*i,2*i-1,min(2*i-1,an));
    			for (int j=0;j<2*i;++j)
    				b[j]=(2*b[j]-g[j]+mo)%mo;
    		}
    		for (int i=0;i<n;++i)
    			c[i]=b[i];
    	}
    	void getrev(int A[],int a[],int n){for (int i=0;i<=n;++i) A[i]=a[n-i];}
    	void getdiv(int c[],int a[],int b[],int an,int bn){
    		static int A[N],B[N],C[N];
    		clear(A,an-bn+1),clear(B,an-bn+1);
    		getrev(A,a,an),getrev(B,b,bn);
    		getinv(B,B,an-bn+1,an-bn+1);
    		multi(C,A,B,an-bn+1,an-bn,an-bn);
    		getrev(c,C,an-bn);	
    	}
    	void getmod(int c[],int a[],int b[],int an,int bn){
    		static int D[N];
    		getdiv(D,a,b,an,bn);
    		multi(D,D,b,an+1,an-bn,bn);
    		for (int i=0;i<bn;++i)
    			c[i]=(a[i]-D[i]+mo)%mo;
    	}
    	int g[N],q[N],mx;
    	void dfs(int n){
    		if (n==0){
    			q[0]=1;
    			mx=0;
    			return;
    		}
    		if (n&1){
    			dfs(n-1);
    			for (int i=mx;i>=0;--i)
    				q[i+1]=q[i];
    			q[0]=0;
    			if (mx+1<k)
    				mx++;
    			else{
    				getmod(q,q,g,mx+1,k);
    				mx=k-1;
    			}
    		}
    		else{
    			dfs(n>>1);
    			multi(q,q,q,2*mx+1,mx,mx);
    			if (2*mx<k)
    				mx*=2;
    			else{
    				getmod(q,q,g,2*mx,k);
    				mx=k-1;
    			}
    		}
    	}
    	int main(int _n,int _k,ll _f[],ll _a[]){
    		n=_n,k=_k;
    		for (int i=1;i<=k;++i)
    			f[i]=_f[i];
    		for (int i=0;i<k;++i)
    			a[i]=_a[i],(a[i]+=mo)%=mo;
    		for (int i=0;i<k;++i)
    			g[i]=(mo-f[k-i])%mo;
    		g[k]=1;
    		dfs(n);
    		ll ans=0;
    		for (int i=0;i<k;++i)
    			(ans+=(ll)q[i]*a[i])%=mo;
    		printf("%lld
    ",ans);
    		return 0;
    	}
    }
    int main(){
    	freopen("in.txt","r",stdin);
    	int ask;
    	scanf("%d%d",&n,&ask);
    	for (int i=0;i<n;++i)
    		scanf("%lld",&a[i]);
    	fir=1;
    	for (int i=0;i<n;++i){
    		ll delta=a[i];
    		for (int j=1;j<=m;++j)
    			(delta-=r[j]*a[i-j])%=mo;
    		if (delta==0) continue;
    		if (fir){
    			mf=m,fail=i,fdelta=delta;
    			m=i+1;
    			fir=0;
    			continue;
    		}
    		ll m_=max(m,i-fail+mf);
    		static ll t[N];
    		memset(t,0,sizeof(ll)*(m_+1));
    		ll tmp=delta*qpow(fdelta)%mo;
    		t[i-fail]=tmp;
    		for (int j=1;j<=mf;++j)
    			t[i-fail+j]=-tmp*f[j]%mo;
    		if (m-i<mf-fail){
    			mf=m,fail=i,fdelta=delta;
    			memcpy(f,r,sizeof(ll)*(m+1));
    		}
    		m=m_;
    		for (int j=1;j<=m;++j)
    			r[j]=(r[j]+t[j])%mo;
    	}
    	
    	for (int i=m;i<n;++i){
    		ll sum=a[i];
    		for (int j=1;j<=m;++j)
    			(sum-=r[j]*a[i-j])%=mo;
    		if (sum)
    			printf("%d
    ",i);
    		assert(sum==0);
    	}
    	
    	//printf("%lld
    ",m);
    	for (int i=1;i<=m;++i)
    		printf("%lld ",r[i]=(r[i]+mo)%mo);
    	printf("
    ");
    	chang::main(ask,m,r,a);
    		
    	return 0;
    }
    /*
    Input 1
    7
    1 2 4 9 20 40 90
    
    Output 1
    4 
    0 0 10 0
    
    Input 2 
    18
    2 4 8 16 32 64 128 256 512 2 4 8 16 32 64 128 256 512
    
    Output 2 
    0 0 0 0 0 0 0 0 1
    */
    
  • 相关阅读:
    js replace替换 忽略大小写问题
    Spring security实现国际化问题
    Mac 的mysql5.7没有配置文件,如何解决only_full_group_by 问题
    java设计模式学习
    synchronized的锁问题
    理解java的三种代理模式
    [acm]HDOJ 2059 龟兔赛跑
    [acm]HDOJ 2673 shǎ崽 OrOrOrOrz
    [acm]HDOJ 1200 To and Fro
    [acm]HDOJ 2064 汉诺塔III
  • 原文地址:https://www.cnblogs.com/jz-597/p/14983564.html
Copyright © 2011-2022 走看看