Bostan-Mori 算法与常系数齐次线性递推

这两天比较颓废，没有训练的欲望，因此来写写科技。

前置知识：生成函数、快速傅里叶（数论）变换。

问题的引出

常系数齐次线性递推问题：给定 $f_{1 \sim k}$ ， $a_{0 \sim k - 1}$ ，对于数列 $a$ 的第 $n \geq k$ 项，其满足递推式：

a_{n} = \sum_{i = 1}^{k} f_{i} a_{n - i}

求数列第 $n$ 项的值。

对于这个问题，直接模拟给出了一个 $O (n k)$ 的算法。

将递推式写成矩阵形式，则给出了一个 $O (k^{3} \log n)$ 的算法。适用于 $k$ 较小的情况。

对于 $n$ 很大， $k$ 较大，经典算法是快速幂求 $x^{n}$ 对某个特征多项式取模。（Fiduccia 算法，1985
）该算法的复杂度为 $O (k \log k \log n)$ 。

2020 年，Alin Bostan 和 Ryuhei Mori 提出了一个新算法。该算法拥有与上述算法相同的时间复杂度，并且常数更小，更好理解，所需的前置工具也更少（只需使用快速傅里叶变换）。

引问题 1

给定一个有理分式生成函数 $\frac{F (x)}{G (x)}$ ，求其第 $n$ 项的系数。其中 $F (x), G (x)$ 都是多项式。

也即求： $[x^{n}] \frac{F (x)}{G (x)}$ 。

Bostan-Mori 算法

我们不妨在原式子上下各乘上一个 $G (- x)$ ：

\frac{F (x)}{G (x)} = \frac{F (x) G (- x)}{G (x) G (- x)}

这样做有什么好处呢？我们发现， $G (x) G (- x)$ 是一个偶函数，只有偶次项有值，感觉性质非常好啊！

于是我们令 $V (x^{2}) = G (x) G (- x)$ 。

对于上面的 $F (x) G (- x)$ ，我们不妨把它的奇次项和偶次项分离出来。令：

F (x) G (- x) = U_{0} (x^{2}) + x U_{1} (x^{2})

其中 $U_{0} (x^{2})$ 为偶次项多项式， $U_{1} (x^{2})$ 为奇次项多项式。

我们现在讨论 $n$ 的奇偶性，当 $n$ 为偶数时， $U_{1} (x^{2})$ 对答案没有贡献，因此：

\begin{aligned} [x^{n}] \frac{F (x)}{G (x)} & = [x^{n}] \frac{U_{0} (x^{2})}{V (x^{2})} \\ = [x^{\frac{n}{2}}] \frac{U_{0} (x)}{V (x)} \end{aligned}

同理，当 $n$ 为奇数时：

\begin{aligned} [x^{n}] \frac{F (x)}{G (x)} & = [x^{n}] \frac{x U_{1} (x^{2})}{V (x^{2})} \\ = [x^{\frac{n - 1}{2}}] \frac{U_{1} (x)}{V (x)} \end{aligned}

发现我们通过上述讨论，可以将问题规模缩小到原来的二分之一左右。

那么我们直接递归运用该算法求出更小规模的解即可。

递归边界是 $n = 0$ ，此时易得 $[x^{0}] \frac{F (x)}{G (x)} = \frac{f_{0}}{g_{0}}$ 。
（以上 $f_{0}, g_{0}$ 分别表示 $F (x), G (x)$ 的常数项系数。)

截一张论文里的算法流程：

回到原问题

接下来我们只需要求出数列的生成函数就可以利用 Bostan-Mori 实现数列第 $n$ 项求值了。

实际上这个是简单的，考虑我们一般求生成函数的过程，在程序中模拟即可。

以下 $f$ 为递推系数， $a$ 为原数列， $F (x)$ 为其生成函数。

以 $k = 3$ 为例：

\begin{aligned} F (x) & = a_{0} x^{0} + a_{1} x^{1} + a_{2} x^{2} + a_{3} x^{3} + a_{4} x^{4} + \dots \\ x f_{1} F (x) & = a_{0} f_{1} x^{1} + a_{1} f_{1} x^{2} + a_{2} f_{1} x^{3} + a_{3} f_{1} x^{4} + \dots \\ x^{2} f_{2} F (x) & = a_{0} f_{2} x^{2} + a_{1} f_{2} x^{3} + a_{2} f_{2} x^{4} + \dots \\ x^{3} f_{3} F (x) & = a_{0} f_{3} x^{3} + a_{1} f_{3} x^{4} + \dots \\ (1 - x f_{1} - x^{2} f_{2} - x^{3} f_{3}) F (x) & = a_{0} x^{0} + (a_{1} - a_{0} f_{1}) x^{1} + (a_{2} - a_{1} f_{1} - a_{0} f_{2}) x^{2} + (a_{3} - a_{2} f_{1} - a_{1} f_{2} - a_{0} f_{3}) x^{3} \end{aligned}

不难得出最终生成函数为如下形式：

F (x) = \frac{\sum_{i = 0}^{k} (a_{i} - (\sum_{j = 0}^{i - 1} a_{j} f_{i - j})) x_{i}}{1 - \sum_{i = 1}^{k} x^{i} f_{i}}

然后分子是一个卷积的形式，可以 FFT 算，分母直接得出，然后就可以使用 Bostan-Mori 算法了！

（以下为未经任何常数优化的代码）

ll Bostan_Mori(poly F,poly G,int n){
	while(n){
		poly Q = G; //Q(x)=G(-x)
		for(int i=1;i<Q.size();i+=2)Q[i]=(mod-Q[i])%mod;
		poly V = Q*G; //V(x^2)=G(x)G(-x)
		poly T;
		for(int i=0;i<V.size();i+=2)T.a.push_back(V[i]);
		poly U = F*Q; //U(x^2)=F(x)G(-x)
		poly W;
		for(int i=n%2;i<U.size();i+=2)W.a.push_back(U[i]);
		F=W; F.reduce();
		G=T; G.reduce();
		n/=2;
	}
	if(!F.size() || !G.size())return 0;
	return F[0]*inv(G[0])%mod;
}

参考资料

Bostan, A., & Mori, R. (2020). A simple and fast algorithm for computing the N-th term of a linearly recurrent sequence. arXiv Preprint. arXiv:2008.08822v1 [cs.SC]. https://arxiv.org/pdf/2008.08822
Elegia 题解 P4723 【【模板】常系数齐次线性递推】 https://www.luogu.com.cn/article/oggmoj4k