从分式第n项到线性递推——bostan-mori 算法的扩展应用

本文是在洛谷博客、github 博客同时发布的P4723 【模板】常系数齐次线性递推题解。

介绍一个常数小还极其好写的科技：bostan-mori 算法。

与其他线性递推算法不同，利用本方法求解线性递推问题不需要任何矩阵知识，唯一前置知识：多项式乘法

波斯坦-茉莉算法简介

这个东西是求分式第 (n) 项，即 ([x^n]frac {f(x)}{g(x)}) 的，而我们知道分式第 n 项和线性递推式可以很容易地互化的（后文会细说）。于是我们先看如何利用波斯坦-茉莉算法求解分式第 (n) 项。

(egin{aligned}&[x^n]frac{f(x)}{g(x)}\=&[x^n]frac{f(x)g(-x)}{g(x)g(-x)}\=&[x^n]frac{c_{even}(x^2)+xcdot c_{odd}(x^2)}{v(x^2)}\=&[x^{n/2}]frac{c_{even}(x)}{v(x)},N ext{ is even}\&[x^{n/2}]frac{c_{odd}(x)}{v(x)},N ext{ is odd} end{aligned})

然后就可以 (mathcal{O}(klog klog n)) 求出解。

示例代码如下：

int divAt(Poly F,Poly G, ll k){
	int i;
	for(;k;k>>=1){
		Poly R=G;
		// R=G(-x)
		for(i=1;i<R.size();i+=2)R[i]=mod-R[i];
		F*=R,G*=R;
		for(i=k&1;i<F.size();i+=2)F[i/2]=F[i];
		F.resize(i/2);
		for(i=0;i<G.size();i+=2)G[i/2]=G[i];
		G.resize(i/2);
	}
	return F.empty()?0:F[0]*qpow(G[0],mod-2)%mod;
}

应用场景：需要注意这一个算法要求模数是质数。

好，回到线性递推。

从线性递推到分式第 n 项

如果你会两者之间的转换，可以直接跳到第三节。

我们求分式第 (n) 项可以化作线性递推后求解：

设 (frac{f(x)}{g(x)}=h(x))，其中 (deg(f)=m,deg(g)=k)，则 (f=gast h)。

根据多项式乘法，项数 (i) 很大的时候 (f_i=0)， (h) 就是一个递推系数为 (-frac{g_{1cdots k}}{g_0}) 的线性递推式：

(0=f_i=g_0h_i+g_1h_{i-1}+cdots+g_{k}h_{i-k})

(h_i=sum_{j=1}^{k}frac{-g_j}{g_0}h_{i-j})

而 (h) 的前 (0cdots m) 项可以容易地通过 (fast g^{-1}(mod x^k)) 得出，于是可以求解。

(m>k) 的时候可以用 (frac fg) 余数进行计算；(m=k) 的时候求出的是 (h_0cdots h_k)，移一位即可获得答案。

参考实现如下（其中 get_an(F,A,n,k) 为与本题相同的含义）：

int div_at(Poly F,Poly G,ll n){
	int m=F.size(),k=G.size();
	if(m>k){
		Poly P=F/G,R=F-P*G;
		return div_at(R,G,n)+P.at(n);
	}
	Poly h=poly.inv(G,k)*F;
	if(m==k){
		for(int i=0;i<k&&i+1<h.size();i++)h[i]=h[i+1];
		n--;
	}
	h.cut(k-1); // 就是 h %= x^{k-1}
	int invg0=qpow(G[0],mod-2);
	G*=-invg0;
	return get_an(G,h,m,n-1);
} 

从分式第 n 项到线性递推

我们有一个以 (F_{1cdots k}) 为递推系数的线性递推序列 (h)，想要构造 (f,g) 使得 ([x^n]frac fg=h_n)

首先，根据多项式乘法，构造项数 (i) 特别大的时候（(f_i=0)）的递推关系。套用上面的式子：

(h_ig_0=sum_{j=1}^{k}{-g_j}h_{i-j})

令 (g_0=0,g_i=F_i(iin[1,k]))即可。

又：(f=gh(mod x^k))

所以可以求出 (f) 的 ([0,k-1]) 项。

参考实现：

int getAn(Poly F,Poly A,ll n,int k){
	F=Poly{1}-F;
	Poly f=A*F;
	f.cut(k); // 就是 h %= x^k
	return divAt(f,F,n);
} 

其中 divAt 套用波斯坦-茉莉算法即可，非常的方便。

代码

namespace POLY{
	int divAt(Poly F,Poly G, ll k){
		int i;
		for(;k;k>>=1){
			Poly R=G;
			// R=G(-x)
			for(i=1;i<R.size();i+=2)R[i]=mod-R[i];
			F*=R,G*=R;
			for(i=k&1;i<F.size();i+=2)F[i/2]=F[i];
			F.resize(i/2);
			for(i=0;i<G.size();i+=2)G[i/2]=G[i];
			G.resize(i/2);
		}
		return F.empty()?0:F[0]*qpow(G[0],mod-2)%mod;
	}
	int getAn(Poly F,Poly A,ll n,int k){
		F=Poly{1}-F;
		Poly f=A*F;f.cut(k);
		return divAt(f,F,n);
	} 
}
int main(){
	int i,x,n,k;
	scanf("%d%d",&n,&k);
	Poly F(k+1),A(k);
	F[0]=0;
	for(i=1;i<=k;i++){
		scanf("%d",&x);
		F[i]=(x+mod)%mod;
	}
	for(i=0;i<k;i++){
		scanf("%d",&x);
		A[i]=(x+mod)%mod;
	}
	printf("%d
",POLY::getAn(F,A,n,k));
	return 0;
}

完整代码就是封装了一些多项式乘除法，就不贴了，不会多项式的大常数选手写的很垃圾。