BM算法学习笔记

BM算法学习笔记

CSP前就要多学学乱搞（确信）

myy的论文里引进了BM算法，可以\(n^2\)求出齐次线性递推式

学习于大佬的博客：

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

https://www.cnblogs.com/zzqsblog/p/6877339.html

本文没有易证易得，符号非常初等，请放心食用。

我们要求的递推式形如

\[\forall m< i\leq n,a_i=\sum_{j=1}^m a_{i-j}b_j \]

\(b_i\)即为递推式系数。

考虑增量法，设\(R[i]\)表示当前的 \(b\)，对于新的一个要求\(a_i\)，有两种情况：

1. 不用更改，此时\(a_i=\sum_{j=1}^m a_{i-j}R[i][j]\)
2. 需要更改递推式，考虑构造递推式的减量递推式\(F\)

令\(\Delta_i=a_i-\sum_{j=1}^m a_{i-j}R[i][j]\) 为当前\(i\)计算出多的部分，则减量\(F\)需要满足：

1. \[\forall |R'|<k<i,\sum_{j=1}^{|F_j|}a_{k-j}F_j=0 \]

2. \[\sum_{j=1}^{|F_j|}a_{i-j}F_j=\Delta_i \]

即\(F\)需要满足计算\(a_k(k<i)\)的时候都是0，而计算\(a_i\)的时候是\(\Delta_i\)，才能使减去之后\(\Delta'=0\)

如果之前没出现过需要更改的情况，直接搞上\(i\)个0一定是合法的。(其实也只会在\(i=1\)出现)

否则直接用前面某个\(id\)的\(R_{id}\)来构造\(F\)。对于\(R_{id}\)，由于定义，其满足

\[\sum_{j=1}^{|R_{id}|}a_{id-j}R_{id}-a_{id}=\Delta_{id} \]

那么我们把它除以\(\Delta_{id}\)，乘上\(\Delta_i\)，再等效替换出 \(i\)：

\[\sum_{j=1}^{|R_{id}|}a_{i-(i-id)-j}\frac{R_{id}[j]}{\Delta_{id}}-\frac{a_{id}}{\Delta_{id}}=1\\ \sum_{j=1}^{|R_{id}|}a_{i-(i-id)-j}\frac{R_{id}[j]\Delta_i}{\Delta_{id}}-\frac{a_{id}\Delta_i}{\Delta_{id}}=\Delta_i \]

那么原来的第\(j\)位变成了第\(i-id+j\)位，相当于右移了\(i-id\)位，再在前面放上一个\(-\frac{\Delta_i}{\Delta_{id}}\)。而对于 \(k<id\) 都是\(\sum_j a_{k-j}R_{id}[j]-a_k\)等于0的情况，这些情况也还是0，因为乘除对0没有影响。这样我们就构造出了满足要求的\(F\)

那么当前的答案就是\(R_{i}\)减去\(F\)。注意\(R_i\)并不需要加F，因为后面还要用。

总结一下，减量F等于

\[\left\{\overbrace{0,0,\cdots,0}^{i-id-1个0},-\frac{\Delta_i}{\Delta_{id}},\frac{R_{id}[1]\Delta_i}{\Delta_{id}},\frac{R_{id}[2]\Delta_i}{\Delta_{id}},\cdots\right\} \]

那么id选哪个呢？我们要求的是最短递推式（不然整\(n\)个0也行），那么选\(i-id+|R_{id}|\)最小的就行了。每次就比较一下\(-i+|R_i|\)更新\(id\)，具体看代码实现。

zzq博客的例子比较生动形象，我就不再造了。

#include<bits/stdc++.h>
#define FOR(i,a,b) for(register int i=a;i<=b;++i)
#define ROF(i,a,b) for(register int i=a;i>=b;--i)
using namespace std;
int read(){
	int x=0,pos=1;char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) if(ch=='-') pos=0;
	for(;isdigit(ch);ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';
	return pos?x:-x; 
} 
const int N = 1.5e5+200;
const int mod = 998244353;
#define ll long long
namespace POLY{
	int r[N],lim,li,omg[N][2];
	const int nlim=131072;
	inline int ksm(int a,int b){
		int res=1;while(b){
			if(b&1) res=1ll*res*a%mod;
			a=1ll*a*a%mod;b>>=1;
		}return res;
	}
	void init1(){
		omg[0][0]=omg[0][1]=1;
		omg[1][1]=ksm(3,(mod-1)/nlim);omg[1][0]=ksm(omg[1][1],mod-2);
		FOR(i,2,nlim-1) omg[i][1]=1ll*omg[i-1][1]*omg[1][1]%mod,omg[i][0]=1ll*omg[i-1][0]*omg[1][0]%mod;
	}
	inline void init(){
		FOR(i,0,lim-1) r[i]=(r[i>>1]>>1)|((i&1)?(1<<(li-1)):0);
	}
	inline void ntt(int *a,int opt){
		FOR(i,0,lim-1) if(i<r[i]) swap(a[i],a[r[i]]);
		for(int l=2;l<=lim;l*=2){
			int m=l/2;for(int *g=a;g!=a+lim;g+=l){
				FOR(i,0,m-1){
					int t=1ll*omg[nlim/l*i][opt]*g[i+m]%mod;
					g[i+m]=(g[i]-t+mod)%mod;
					g[i]=(g[i]+t)%mod;
				}
			}
		}if(!opt){
			int iv=ksm(lim,mod-2);
			FOR(i,0,lim-1) a[i]=1ll*a[i]*iv%mod;
		}
	}
	static int x[N],y[N];
	void mul(int *a,int *b,int al,int bl){
		lim=1,li=0;while(lim<=al+bl) lim*=2,li++;init();
		FOR(i,0,al) x[i]=a[i];FOR(i,al+1,lim-1) x[i]=0;
		FOR(i,0,bl) y[i]=b[i];FOR(i,bl+1,lim-1) y[i]=0;
		ntt(x,1);ntt(y,1);
		FOR(i,0,lim-1) x[i]=1ll*x[i]*y[i]%mod;
		ntt(x,0);FOR(i,0,lim-1) a[i]=x[i];
	}
	static int FV[N];
	inline void getinv(const int *a,int *b,int n){
		if(n==1) return b[0]=ksm(a[0],mod-2),void();
		getinv(a,b,(n+1)/2);
		lim=1,li=0;while(lim<=n*2) lim*=2,li++;init();
		FOR(i,(n+1)/2,lim-1) b[i]=0;
		FOR(i,0,n-1) FV[i]=a[i];FOR(i,n,lim-1) FV[i]=0;
		ntt(FV,1);ntt(b,1);
		FOR(i,0,lim-1) b[i]=(2ll-1ll*FV[i]*b[i]%mod+mod)%mod*b[i]%mod;
		ntt(b,0);
	}
	static int fr[N],gr[N],tmp[N];
	inline void getmod(int *f,int *g,int n,int m,int *q){
		FOR(i,0,n) fr[n-i]=f[i];
		FOR(i,0,m) gr[m-i]=g[i];
		FOR(i,n-m+1,m) gr[i]=0;
		getinv(gr,tmp,n-m+1);
		FOR(i,0,n-m) gr[i]=tmp[i];
		mul(fr,gr,n,n-m);
		reverse(fr,fr+n-m+1);
		mul(fr,g,n-m,m);
		FOR(i,0,m-1) q[i]=(f[i]-fr[i]+mod)%mod;
	}
	int Mtmp[N];
	inline void mulmod(int *f,int *g,int *h,int n,int *res){
		FOR(i,0,n) Mtmp[i]=g[i];
		mul(Mtmp,f,n,n);
		getmod(Mtmp,h,n*2,n,res);
	}
	int F[N],G[N],H[N];
	inline int polydt(int *c,int *a,int n,int k){
		FOR(i,1,k) F[k-i]=(mod-c[i])%mod; 
		F[k]=1;
		G[1]=1; H[0]=1;
		for(int m=n;m;m>>=1){
			if(m&1){
				mulmod(H,G,F,k,H);
			}
			mulmod(G,G,F,k,G);
		}
		int ans=0;
		FOR(i,0,k-1) ans=(ans+1ll*H[i]*a[i]%mod)%mod;
		return ans;
	} 
}
int n,pn=0;
int R[N],Rn[N];
int rl,nl;
int a[N],dlt[N];
int tp[N],tl;
typedef long double ld;
int ksm(int a,int b){
	int res=1; 
	while(b){
		if(b&1){
			res=1ll*res*a%mod;
		}a=1ll*a*a%mod;b>>=1;
	}
	return res;
}
const ld eps = 1e-7;
int main(){
	POLY::init1();
	n=read();int m=read();
	FOR(i,1,n) scanf("%d",&a[i]);
	int id=0,mi=0;
	FOR(i,1,n){
		ll dlti=(mod-a[i])%mod;
		for(int j=1;j<=rl;j++){
			dlti=(dlti+1ll*a[i-j]*R[j]%mod)%mod;
		}
		dlt[i]=dlti;
		if(dlti==0) continue;
		if(i==1){
			rl=1;R[1]=Rn[1]=0;id=1;nl=rl=1;continue;
		}else{
			FOR(j,1,rl) tp[j]=R[j];tl=rl; //??R 
			int tmp=1ll*dlti*ksm(dlt[id],mod-2)%mod;
			FOR(j,1,nl) R[i-id+j]=(R[i-id+j]-1ll*Rn[j]*tmp%mod+mod)%mod;   //??F?R 
			R[i-id]=(R[i-id]+tmp)%mod;
			rl=max(rl,i-id+nl);
			nl=tl;FOR(j,1,tl) Rn[j]=tp[j];id=i;
			//??Rn 
		}
	}
	while(R[rl]==0) rl--;
	FOR(i,1,rl){
		printf("%d ",R[i]);
	}putchar(10);
	FOR(i,1,rl) a[i-1]=a[i];FOR(i,rl,n) a[i]=0;
	int ans=POLY::polydt(R,a,m,rl);
	printf("%d\n",ans);
	return 0;
}