一种新的线性递推计算方法

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

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

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

新线性递推算法的主要思想是考虑将线性递推转化为生成函数的形式 $\frac{P(x)}{Q(x)}$。注意我们容易在 $\mathsf{M}(n)$ 的时间内通过数列的前 $d$ 项得到 $P(x)$，只需计算 $\left[({\sum }_{i<d}{f}_i{x}^i)Q(x)\right]\mathrm{m}\mathrm{o}\mathrm{d}{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)},U(x^2)=Q(x)Q(-x)$$

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

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

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

$\to \frac{4}{3}\mathsf{M}(d)$

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

因此我们减少了 $\frac{1}{3}$ 的 FFT 运算量，变为 $4\mathsf{E}(2d)=\frac{4}{3}\mathsf{M}(d)$。

$\to \mathsf{M}(d)$

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

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

$\to \frac{2}{3}\mathsf{M}(d)$

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

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

$$\begin{array}{rl}{\hat{A}}_{2k}^{(2n)}& ={\hat{A}}_k\\ {\hat{A}}_{2k+1}^{(2n)}& =\sum _i({\omega }_{2n}^i{A}_i){\omega }_n^{ik}\end{array}$$

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

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

---

你可能认为这是改了一个问题所以变容易了，其实并不，我们可以将原问题：计算 $x^N\mathrm{m}\mathrm{o}\mathrm{d}rev(Q)$ 的时间也降低常数！

我们考虑计算 $\frac{1}{Q(x)}$ 的第 $N,N+1,\dots ,N+d-1$ 项，由 $\frac{1}{Q(x)}=\frac{Q(-x)}{Q(x)Q(-x)}$ 可知，我们只需算出 $\frac{1}{Q(x)Q(-x)}$ 的连续 $2d$ 项，也就是 $U(x^2)=Q(x)Q(-x)$ 后 $\frac{1}{U(x)}$ 的连续 $d$ 项。

其中，$Q\to U$ 这个过程每轮还是 $\frac{1}{3}\mathsf{M}(d)$ 的，这个过程已经算出了 $Q(-x)$ 的 DFT，而我们要和 $\frac{1}{U(x^2)}$ 相乘，注意到这里 $2d$ 长的 DFT 足够，因为后面的部分会溢出掉，所以这个 loop 的总共消耗就是恰好的 $\mathsf{M}(d)$。即总耗时是 $(\log N+O(1))\mathsf{M}(d)$。

接下来注意到从第 $N$ 项截取的 $\frac{1}{Q(x)}$ 记为 $P(x)$，那么 $rev(P(x)Q(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^d)$ 既得。

---

上面那一小段是原论文的见解。但是事实上我们之间将一开始的算法转置就是多项式取模了。