BM 学习笔记

两个 BM 哟

1.Bostan-Mori

常系数其次线性递推。

实际上这个算法是用来计算 $[x^n]\frac {F(x)}{G(x)}$ 的。。。

我们考虑一个神奇的多项式：$F(x)F(-x)$。容易发现存在一个多项式 $H(x)$ 满足 $H(x^2)=F(x)F(-x)$。

于是我们对这个分式的分子和分母同时乘上 $G(-x)$，能够得到有 $[x^n]\frac {F(x)G(-x)}{G(x)G(-x)}$。

按照奇偶拆开，假设这玩意儿是 $[x^n]\frac {Q(x^2)+xP(x^2)}{H(x^2)}$ ，分子的两个多项式一个只在偶数次处有值，另一个只在奇数次处有值。

所以我们按照 $n$ 的奇偶性来讨论应该取 $Q(x)$ 还是 $P(x)$。然后得到了规模减小一半的子问题。

最后变成 $[x^0]\frac {F(x)}{G(x)}$，此时计算 $\frac {f_0}{g_0}$ 就行。

那么这个东西如何解决常系数其次线性递推呢？

常系数其次线性递推，说白了就是给你一个递推式和这个数列的前几项，问你第 $n$ 项。

众所周知这个是可以矩快 $O(m^3\log n)$ 的。但是以上算法可以做到 $O(m\log m\log n)$虽然一般情况不需要 NTT，所以正常情况下复杂度是 $O(m^2\log n)$ 的。

首先我们知道斐波那契数列的生成函数是 $F(x)=\frac 1 {1-x-x^2}$，通过这个我们就可以 $O(2^2\log n)$ 计算出斐波那契数列的第 $n$ 项比矩快少三次乘法。

inline int Solve(long long n){
	int f0=1,f1=0,g0=1,g1=mod-1,g2=mod-1;
	for(;n;n>>=1){
		if(n&1)f0=(1ll*g0*f1+1ll*(mod-g1)*f0)%mod,f1=1ll*g2*f1%mod;
		else f0=1ll*g0*f0%mod,f1=(1ll*g2*f0+1ll*(mod-g1)*f1)%mod;
		g1=(2ll*g0*g2+1ll*(mod-g1)*g1)%mod;g0=1ll*g0*g0%mod;g2=1ll*g2*g2%mod;
	}
	return 1ll*f0*pow(g0,mod-2)%mod;
}

其实对于正常递推式都可以写成 $\frac {F(x)}{G(x)}$ 的形式。。。但是不会证。。。

然而我们知道有 $\frac {F(x)}{G(x)} \equiv H(x) \bmod x^{m+1}$，那么我们就可以知道 $F(x)=G(x)H(x)$。

$G(x)$ 实际上就是这个序列的特征多项式，可以通过 $1-P(x)$（$P(x)$ 是递推式）得到。

inline void times(int*f,int*g,const int&n){
	static int sav[M<<1];int i,j;
	for(i=0;i<n;++i)for(j=0;j<n;++j)sav[i+j]=(sav[i+j]+1ull*f[i]*g[j])%mod;
	for(i=0;i<(n<<1);++i)f[i]=sav[i],sav[i]=0;
}
inline int Solve(int*f,int*g,const int&len,long long n){
	static int sav[M<<1];int i;
	for(;n;n>>=1){
		for(i=0;i<len;++i)sav[i]=g[i],sav[i+len]=0;
		for(i=1;i<len;i+=2)sav[i]=mod-sav[i];times(f,sav,len);times(g,sav,len);
		for(i=n&1;i<(n<<1);i+=2)f[i>>1]=f[i];
		for(i=0;i<(n<<1);i+=2)g[i>>1]=g[i],f[(i>>1)+n]=0;
	}
	return 1ll*f[0]*pow(g[0],mod-2)%mod;
}

好处？你可以不用推矩阵和调矩阵的那些细节了，并且更快。

坏处？如果模数是非质数那就寄了。但是老算法是不需要模数为质数的

如果有多个序列之间推来推去？如果你是萌新那还是矩阵快速幂吧，如果你是老鸽那可以大力拆得只剩下一个序列，然后化成分式继续 BM。

CF392C？那题直接 $O(k^2)$ 跑一个暴力求逆不就完了？

2.Berlekamp-Massey

用来猜一个数列的递推式。

如果是数数题然后你状态设错了，但是能够得到正确的答案，然后你又懒得重新设一遍，就可以 BM+BM $O(m^2\log n)$ 大力草过去。

这个东西和拉格朗日插值有点儿像。

这个 BM 的思想是，我把前几项的递推式都搞出来，并且记录下来，如果我不能由上一项的递推式递推过来，那就通过前面的部分做一些奇怪的操作，然后得到正确的递推式。但实际上我们正常情况下记录的都是最短递推式，所以递推式不一定正确，使用这个 BM 之前建议打表尽量多项，尽管是 $O(n^2)$ 的

好，我们假设前面的递推式是 $R_0$ 到 $R_{cnt-1}$。

$delta_i$ 为到第 $i$ 项之前，上一个递推式在 $i$ 处的值减去 $a_i$。

$fail_i$ 的意思是 $R_i$ 在 $fail_i$ 处失效。

如果 $delta_i=0$，那么上一个递推式一定是正确的，直接跳过。

否则我们再构造一个递推式，使得这个递推式在前面 $i-1$ 处的值为 $0$，在 $i$ 处的值为 $delta_i$。

如何得到这个递推式？

如果我前面一个递推式都没有，那么这个递推式一定是一车 $0$ 后面跟了一个 $1$。

否则我们构造一个神秘递推式。随便拉一个递推式过来，假设是第 $x$ 个递推式。

那么我们设 $t = \frac {delta_i} {delta_{fail_x}}$。

我们构造：

$$\{ 0,0,...0,t,-t \times R_x[1],-t \times R_x[2],... \}$$

其中开头有 $i-fail_x-1$ 个 $0$。

于是我们找到 $fail_x+R[x].len$ 最小的 $x$ 就能够得到最短递推式啦。

inline int BM(int*a,const int&n){
	int i,j,id,tmp,cnt=0;
	for(i=0;i^n;++i){
		delta[i]=a[i];
		if(!cnt){
			if(a[i])fail[cnt++]=i,R[cnt].resize(i,0);
			continue;
		}
		delta[fail[cnt]=i]=a[i];
		for(j=0;j^R[cnt].size();++j)delta[i]=Del(delta[i],1ll*a[i-j-1]*R[cnt][j]%mod);
		if(!delta[i])continue;id=cnt-1;
		for(j=0;j^cnt;++j)if(fail[j]-R[j].size()>fail[id]-R[id].size())id=j;
		tmp=1ll*delta[i]*pow(delta[fail[id]])%mod;++cnt;R[cnt]=R[cnt-1];
		while(R[cnt].size()<i-fail[id]+R[id].size())R[cnt].push_back(0);
		R[cnt][i-fail[id]-1]=Add(R[cnt][i-fail[id]-1],tmp);
		for(j=0;j^R[id].size();++j){
			R[cnt][i-fail[id]+j]=Del(R[cnt][j-fail[id]+i],1ll*tmp*R[id][j]%mod);
		}
	}
	for(i=0,id=R[cnt].size();i^id;++i)p[i+1]=R[cnt][i];
	return id;
}

递推式将会存放在 $p$ 中，长度会由 $BM()$ 返回。