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},n\ge 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\neq 0\):如果之前没有其它递推式,令递推式为\(i+1\)个\(0\)。

否则:

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

根据定义得\(\forall i\in[|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\n",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\n",i);
		assert(sum==0);
	}
	
	//printf("%lld\n",m);
	for (int i=1;i<=m;++i)
		printf("%lld ",r[i]=(r[i]+mo)%mo);
	printf("\n");
	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
*/
上一篇:**1279 - Asteroid Rangers (最小生成树)


下一篇:动力系统笔记1