草稿。。。

自翻译了 David Harvey 的两篇论文（即 多项式牛顿迭代的分块优化） 之后，我就一直对 inv 9E 的实现好奇，它是 3 阶牛顿迭代，比其他 2 阶迭代复杂度低一截，在分块思路上也继续有效。

因此我就自然而然的产生了一个想法：是否存在 4 阶牛顿迭代呢？

经过一番思索，还真成功的推导出了 $\frac{28}{3}$ E 的逆和 $\frac{31}{3}$ E 的倒数平方根算法。

下面让我介绍这个算法。记号与上一篇博客相同，$\mathsf{E}(n)$ 是一次长 $n$ 的 DFT/IDFT 运算。

倒数

假设我们已经获得了获得了 $f$ 的逆的前 $n$ 项，那么设：

$$g_0 = f^{-1} \bmod {x^n}$$

想扩展到 $4n$ 项，可以假设 $g$ 有如下形式

$$g = g_0(1 + c_0 x^n + c_1x^{2n} + c_2x^{3n}) + O(x^{4n})$$

为了求解多项式 $c_0,c_1,c_2$，我们可以计算 $f g_0$ 的残差：

$$fg_0 = 1 + a_0 x^n + a_1 x^{2n} + a_2 x^{3n} + O(x^{4n})$$

可得方程：

$$fg = (1 + a_0 x^n + a_1 x^{2n} + a_2 x^{3n})(1 + c_0 x^n + c_1x^{2n} + c_2x^{3n}) + O(x^{4n})= 1 + O(x^{4n})$$

展开，对比系数解得

$$\begin{cases} a_0 + c_0 = 0 \\ a_1+c_1+ a_0c_0 = 0 \\ a_2+a_1c_0+a_0c_1+c_2 = 0 \end{cases} \Rightarrow \begin{cases} c_0 = -a_0 \\ c_1 = a_0^2 - a_1 \\ c_2 = -a_2 + 2a_0a_1 - a_0^3 \end{cases}$$

需要注意的是这里 $a_i$ 是长为 $n$ 的多项式块，乘法会向下一块进位，所以其实还需要再取余。记

$$a_0^2 = q_0 + q_1 x^n$$

带入得到，这就是我们要计算的式子了。

$$\begin{aligned} \overline{c_0} &= -a_0,\\ \overline{c_1} &= q_0 - a_1,\\ \overline{c_2} &= q_1 - a_2 + a_0(2a_1 - q_0) \bmod x^n. \end{aligned}$$

计算过程也比较清晰：

0. 已知 $g_0 \bmod x^n$。
1. 计算 $fg_0 \bmod x^{4n}$，，取出 $a_0,a_1,a_2$，开销 3E(4n)。
2. 计算 $a_0^2 = q_0 + q_1x^n$，取出 $q_0, q_1$，开销 2E(2n)。
3. 计算 $a_0(2a_1 - q_0) \bmod x^n$，开销 2E(2n)，这里复用上一步的 $\mathcal{F}_{2n}(a_0)$。
4. 此时我们可以得到 $c = \overline{c_0} + \overline{c_1} x^n + \overline{c_2}x^{2n}$。
5. 计算 $g_0 c$，开销 2E(4n)，这里复用了第一步的 $\mathcal{F}_{4n}(g_0)$。

总而言之，复杂度为：

$$T(n) = 3 \mathsf{E}(4n) + 2 \mathsf{E}(2n) + 2 \mathsf{E}(2n) + 2 \mathsf{E}(4n) = 28 \mathsf{E}(n)$$

因为每次 $n \to 4n$，所以复杂度为

$$T(4n) = T (n) + 28 \mathsf{E}(n), \quad T(n) = \frac{28}{3}\mathsf{E}(n) \approx 9.33 \mathsf{E}(n)$$

Taylor 展开

其实这个方程可以更快的求解，注意到：

$$C = cx^n = f g_0 - 1 = a_0x^n +a_1x^{2n}+a_2x^{3n} \pmod {x^{4n}}$$

然后考虑Taylor展开：

$$\begin{aligned} (1+C)^{-1} &= 1 - C + C^2 - C^3 + O(x^{4n}) \\ &= 1 - a_0x^n + (a_0^2-a_1)x^{2n} + (2a_0a_1 - a_0^3-a_2) x^{3n} + O(x^{4n}) \end{aligned}$$

这样就可以直截了当地得到。

倒数平方根

道理相同，假设我们已经获得了获得了 $f$ 的倒数平方根的前 $n$ 项，那么设：

$$g_0 = f^{-1/2} \bmod {x^n}$$

同样设

$$C = cx^n = f g_0^2 - 1 = a_0x^n +a_1x^{2n}+a_2x^{3n} \pmod {x^{4n}}$$

Taylor 展开得到

$$\begin{aligned} (1+C)^{-1/2} &= 1 - \frac12 C + \frac38 C^2 - \frac{5}{16}C^3 + O(x^{4n}) \\ &= 1 - \frac12(a_0x^n+a_1x^{2n}+a_2x^{3n}) \\ &\quad + \frac38(a_0^2x^{2n}+2a_0a_1x^{3n}) - \frac{5}{16}a_0^3x^{3n} + O(x^{4n}) . \end{aligned}$$

和倒数一样，$a_i$ 是多项式块，所以 $a_0^2$ 的高半会进到下一块。仍然记

$$a_0^2 = q_0 + q_1x^n.$$

那么修正因子

$$g = g_0(1+c_0x^n+c_1x^{2n}+c_2x^{3n}) + O(x^{4n})$$

中的三个块为

$$\begin{aligned} \overline{c_0} &= -\frac12 a_0,\\ \overline{c_1} &= -\frac12 a_1 + \frac38 q_0,\\ \overline{c_2} &= -\frac12 a_2 + \frac38 q_1 + a_0\left(\frac34a_1-\frac{5}{16}q_0\right) \bmod x^n. \end{aligned}$$

如果能直接算出 $C$ 的前三个块，那这里就结束了。问题在于倒数平方根的残差是 $fg_0^2-1$，而 $g_0^2$ 已经有 $2n$ 项，$f$ 又 $4n$ 项，直接做普通卷积需要 $8n$，因此需要另想循环卷积的法子。

注意到 $f$ 乘 $g_0^2$ 最多只有 6 项，开头一项还是 $1$，因此实际上只有 5 块。因此注意到 $4n$ 的循环卷积可以写成：

$$fg_0^2 \bmod (x^{4n} - 1) = 1 + a_3 + (a_0+a_4)x^n + a_1 x^{2n} + a_2x^{3}$$

接下来需要想办法把 $a_4$ 抠出来。一个办法是再算一次点值

$$\begin{aligned} fg_0^2 \bmod (x^{n} - \xi_8) &= f(\eta x)g_0^2(\eta x) \bmod (x^n - 1) \\ &= 1 + \xi_8 a_0 + \xi_8^2a_1 + \xi_8^3a_2 + \xi_8^4a_3 + \xi_8^5a_4 \end{aligned}$$

其中 $\xi_8$ 是 8 阶单位根，$\xi_8^4 = -1$，而 $\eta = \xi_{8n}$，也叫做 $8n$ 阶单位根。由于 $a_1,a_2,a_3$ 是无干扰的，我们可以计算得到 $a_0$ 了。

这个 trick 和 inv 9E 所使用的 DFT(3n) 的 trick 本质相同，这里变成 DFT(4n) 了。

得到 $a_0,a_1,a_2$ 后和前面的方法基本相同，就不再详述。具体计算流程如下：

0. 已知 $g_0 \bmod x^n$。
1. 计算 $fg_0^2 \bmod (x^{4n}-1)$，开销 3E(4n)。
2. 用长度 $n$ 的 twisted DFT 计算 $x^n=\eta$ 处的 $fg_0^2$，开销 3E(n)。
3. 恢复 $a_0,a_1,a_2$。
4. 计算 $a_0^2=q_0+q_1x^n$，开销 2E(2n)。
5. 计算 $a_0(\frac34a_1-\frac{5}{16}q_0)\bmod x^n$，开销 2E(2n)，这里复用上一步的 $\mathcal{F}_{2n}(a_0)$。
6. 拼出 $c=\overline{c_0}+\overline{c_1}x^n+\overline{c_2}x^{2n}$。
7. 计算 $g_0c$，开销 2E(4n)，这里复用第一步的 $\mathcal{F}_{4n}(g_0)$。

于是单次 $n\to 4n$ 的开销为

$$3\mathsf{E}(4n)+3\mathsf{E}(n)+2\mathsf{E}(2n)+2\mathsf{E}(2n)+2\mathsf{E}(4n)=31\mathsf{E}(n).$$

因此总复杂度为

$$T(4n)=T(n)+31\mathsf{E}(n), \quad T(n)=\frac{31}{3}\mathsf{E}(n)\approx 10.33\mathsf{E}(n).$$

性能比较

取 NTT 模数 $P=998244353$，并使用了 AVX2 来优化。

inv

invsqrt

总结

本文通过尝试 4 阶牛顿迭代，得到了 9.33E 的多项式逆和 10.33E 的多项式倒数平方根算法。其中 inv 并没有在复杂度上取得突破，只是实机性能还不错。invsqrt 似乎并未查找到相关文献，可能这个算子不怎么受人关心，一般论坛上似乎都是使用 12E 的标准牛顿迭代，我这个思路可能是互联网上首次提出。

我还有一个 18E 的 pow 算法的 idea，等我尝试验证后再考虑把本文的两个算法实现交到 OJ 上测一下。