\(\newcommand{\me}{\mathrm{e}}\newcommand{\bbF}{\mathbb F}\newcommand{\calF}{\mathcal F}\newcommand{\sfE}{\mathsf E}\newcommand{\sfM}{\mathsf M}\)当环 \(R\) 中存在 \(n=2^k\) 次单位根 \(\omega_n\) (例如常用的 \(\bbF_{998244353}\)), 我们容易进行 \(O(n\log n)\)已知 \(f\in R[[x]]\) 的前 \(n\) 项, 欲求 \(g = 1/f\) 的前 \(n\) 项, [1] 明确给出了一个 \((1\frac{4}9 + o(1))\sfM(n)\) 的算法, 我们先对其进行介绍, 接着分析其最后提到的可能得到 \((1\frac{1}3 +o(1))\sfM(n)\) 做法的优化方向.

分块原理

传统的 Newton 迭代递归形式为

\[T(n) = (c + o(1))\sfM(n) + T(n/2) \]

因此递归过程中, 主定理使得迭代的常数翻倍, 为 \(T(n) = (2c + o(1))\sfM(n)\). 而在接下来描述的一类分块方法中, 我们有参数 \(r = r(n)\), 取 \(m\) 为大于等于 \(n/r\) 的最小的 \(2\) 的幂, 然后基于其进行计算. 递归形式为

\[T(n) = (c+o(1))\sfM(n) + O(rn) + T(n/r) \]

例如取 \(r = \Theta(\sqrt{\log n})\) 时解得 \(T(n) = (c+o(1))\sfM(n)\).

我们首先介绍一个 \(1\frac23\sfM(n)\) 的简洁做法, 其较为具有一般性.

对于形式幂级数 \(f\), 记 \(f_{[i]}\) 为它的第 \(mi \sim m(i+1)-1\) 次项, 记 \(X=x^m\), 我们如果逐块确定 \(f_{[0]},f_{[1]},\dots\), 可以逐块计算出 \((fg)_{[0]}, (fg)_{[1]},\dots\), 有

\[fg = \sum_{i,j} f_{[i]}g_{[j]}X^{i+j} \]

那么先计算出每个 \(\calF_{2m}(f_{[i]}), \calF_{2m}(g_{[i]})\), 注意到 \((fg)_{[i]}\) 由 \(\sum_{i+j=k} f_{[i]} g_{[j]}\) 和 \(\sum_{i+j=k-1} f_{[i]} g_{[j]}\) 贡献. 借助 DFT 的循环卷积性质可以在 \(O(km)+ \sfE(2m)\) 时间内还原出 \((fg)_{[k]}\).

现在我们已经递归求出 \(g_0 = (f^{-1})_{[0]}\), 由 \(fg = 1\), 我们有

\[f_{[0]}g_{[k]} = - \left[\left(\sum_{i} f_{[i]}X^i\right) \left(\sum_{j<k} g_{[j]}X^j\right)\right]_{[k]} \]

设 \(\psi = \left[\left(\sum_{i} f_{[i]}X^i\right) \left(\sum_{j<k} g_{[j]}X^j\right)\right]_{[k]}\), 我们每轮先计算出 \(\psi\) 后, 只需计算 \(g_{[k]} = -\psi \cdot g_{[0]} \bmod X\).

那么在 \(k\) 遍历完后, 我们对所有计算量进行统计:

1. 计算所有的 \(\calF(f_{[i]})\) 和 \(\calF(g_{[i]})\), 总共 \(2\sfE(2n)\).
2. 每个 \(\psi\), 总共 \(\sfE(2n)\).
3. 计算出 \(g_{[k]}\), 这需要对 \(\psi\) 做 DFT 最后再做回来, 总共 \(2\sfE(2n)\).

综上, 总共消耗 \(10\sfE(n) = 1\frac 23 \sfM(n)\).

三阶 Newton 迭代

如果将 \(g\) 计算到 \(n\) 次项精度, 我们考虑如何将其扩展到 \(3n\) 次精度. 由于 \(fg = 1\), 设 \(g_0\) 为 \(n\) 次精度的解, 有 \(fg_0 = 1 + \delta_1 + \delta_2 + O(x^{3n})\), 其中 \(\delta_1\) 为 \(n\sim 2n-1\) 次项, \(\delta_2\) 为 \(2n\sim 3n-1\) 次项. 那么有

\[f^{-1} = \frac{g_0}{1+\delta_1 + \delta_2} = g_0 (1 - \delta_1 + (\delta_1^2 - \delta_2)) + O(x^{3n}) \]

那么在计算前 \(n\) 次项时, 计算量无外乎还是 \(10\sfE(n)\). 接下来分析后面部分 \(g\) 的计算.

1. 还需要对剩余的 \(f\) 部分也计算 DFT, 消耗 \(2\sfE(2n)\).
2. 根据 \(f,g_0\) 计算出 \(\delta\), 注意这里 \(f,g_0\) 的 DFT 前面已经算了, 消耗 \(2\sfE(2n)\).
3. 计算 \(\delta_1\) 的 DFT, 消耗 \(\sfE(2n)\).
4. 避开直接计算 \(\delta_2\), 转而计算 \(\delta_1^2 - \delta_2\), 消耗 \(\sfE(2n)\), 再计算
   其 DFT, 消耗 \(\sfE(2n)\). (注意这里如果计算了 \(\delta_2\) 就和原来的做法常数一样了.)
5. 计算 \(g\) 的后面 \(2n\) 次项, 消耗 \(2\sfE(2n)\).

共计 \(T(3n) = (\frac{13}3 + o(1))\sfM(n)\), 因此 \(T(n) = (1\frac49 + o(1))\sfM(n)\).

参考文献

[1] David Harvey. “Faster algorithms for the square root and reciprocal of power series”. In: CoRR abs/0910.1926 (2009). arXiv: 0910.1926. url: http://arxiv.org/abs/0910.1926.