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< ileq 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 ight} ]

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

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

#include<bits/stdc++.h>
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; 
} 
#define ll long long
#define FOR(i,a,b) for(int i=a;i<=b;i++) 
typedef long double ld;
const int N = 3021;
const int mod = 1e9+7;
int n,pn=0;
ll R[N],Rn[N];
int rl,nl;
ll a[N],dlt[N];
ll tp[N],tl;
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(){
	n=read();
	FOR(i,1,n) scanf("%lld",&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;mi=0;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);
			if(tl-i<mi) mi=tl-i;nl=tl;FOR(j,1,tl) Rn[j]=tp[j];id=i; //更新Rn 
		}
	}
	while(R[rl]==0) rl--;
	printf("%d
",rl);
	FOR(i,1,rl){
		printf("%ld ",R[i]);
	}
	return 0;
}