题解 P4723 【【模板】常系数齐次线性递推】

· · 题解

这是一篇对于新线性递推论文的解读，你可以在这里找到原文。

一个众所周知的线性递推计算方法是考虑的通过快速幂计算多项式取模 x^N \bmod Q(x)，这被称为 Fiduccia 算法（1985）。设做多项式乘法的时间为 \mathsf M(n)，在一般环上的计算需要约 \mathsf (3\log N + O(1))\mathsf M(d) 的运算量，即使在可以做 FFT 的运算下也只能做到 \mathsf (2\log N + O(1))\mathsf M(d) 的运算量。

而新的算法可以做到一般环上 (2\log N + O(1))\mathsf M(d) 的运算量，和在可以做 FFT 的运算下做到 (\frac 23 \log N + O(1))\mathsf M(d) 的运算量。也就是说在我们通常有 NTT 模数的情况下，同样的 FFT 实现可以让我们的常数大约是原先方法的 \frac {\mathbf 1}{\mathbf3}。

新线性递推算法的主要思想是考虑将线性递推转化为生成函数的形式 \frac{P(x)}{Q(x)}。注意我们容易在 \mathsf M(n) 的时间内通过数列的前 d 项得到 P(x)，只需计算 \left[(\sum_{i<d} f_i x^i) Q(x)\right] \bmod x^d。接下来我们考虑这样一件事：

\frac {P(x)}{Q(x)} = \frac{P(x)Q(-x)}{Q(x)Q(-x)}

这发生了什么呢？我们注意记 V(x)=Q(x)Q(-x)，可以发现 V(x)=V(-x)，说明 V(x) 只有偶次项有值，因此我们就得到了分解

\frac{P(x)}{Q(x)} = \frac {E(x^2)}{U(x^2)} + x\frac {O(x^2)}{U(x^2)}, \quad U(x^2)=Q(x)Q(-x)

因为这分别填满了二进制的 0 和 1 位，所以我们只需递归到一侧即可。

不难发现朴素实现就是上下都做多项式乘法，每一个二进制位约消耗 2\mathsf M(d) 的时间。

接下来我们看看可以如何优化常数。我们记 DFT 的时间为 \mathsf E(2n)，可得 \mathsf M(n)=3\mathsf E(2n)。

\rightarrow \frac 43 \mathsf M(d)

我们注意当知道 Q(x) 的 DFT 时，可以直接得到 Q(-x) 的 DFT。记 A_k 的 DFT 数组为 \widehat A_k，读者不难自行验证 B_k=(-1)^kA_k 的 DFT 数组对应有 \widehat B_k = \widehat A_{k \oplus \frac n2}。

因此我们减少了 \frac 13 的 FFT 运算量，变为 4\mathsf E(2d)=\frac 43 \mathsf M(d)。

\rightarrow \mathsf M(d)

注意到 Q(x)Q(-x) 最后只剩下偶数项，而 P(x)Q(-x) 我们要么只要奇数项要么只要偶数项，假设我们要提取 A_k 的所有偶数项，那么设 B_k = \frac{A_k + (-1)^k A_k}{2}，可得 \widehat B_k = \frac{\widehat A_k + \widehat A_{k \oplus \frac n2}}2，然后对前 n/2 项做长为 n 的 IDFT 即可，奇数项是类似的。

因此我们将其中的两个 \mathsf E(2d) 换成了 \mathsf E(d)，运算量变为 2\mathsf E(2d) + 2\mathsf E(d) = \mathsf M(d)。

\rightarrow \frac 23\mathsf M(d)

我们考虑不做上一步的 IDFT，始终以 DFT 的形式维护 P,Q。那么我们是否能够更小常数地从 A 的长为 n 的 DFT 推出长为 2n 的 DFT 呢？

事实上我们可以在 2\mathsf E(n) 内完成。只需注意到以下两件事：

\begin{aligned} \widehat A^{(2n)}_{2k} &= \widehat A_k\\ \widehat A^{(2n)}_{2k+1} &= \sum_i (\omega_{2n}^{i}A_i) \omega_n^{ik} \end{aligned}

因此我们先 IDFT 得到 A，然后乘以 \omega_{2n}^i 再做 DFT 就可以了。

现在整个迭代流程只有对 P,Q 的 DFT 数组倍长用到了 FFT，运算量为 4\mathsf E(d) = \frac 23 \mathsf M(d)。