1.25M(n) 的多项式求逆
\(\newcommand{\me}{\mathrm{e}}\newcommand{\bbF}{\mathbb F}\newcommand{\calF}{\mathcal F}\newcommand{\sfE}{\mathsf E}\newcommand{\sfM}{\mathsf M}\)
已知 \(f\in R[[x]]\) 的前 \(n\) 项, 欲求 \(g = 1/f\) 的前 \(n\) 项.
分块原理
这一方法首先由 [1] 引入.
传统的 Newton 迭代递归形式为
因此递归过程中, 主定理使得迭代的常数翻倍, 为 \(T(n) = (2c + o(1))\sfM(n)\). 而在接下来描述的一类分块方法中, 我们有参数 \(r = r(n)\), 取 \(m\) 为大于等于 \(n/r\) 的最小的 \(2\) 的幂, 然后基于其进行计算. 递归形式为
例如取 \(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\), 有
那么先计算出每个 \(\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\), 我们有
设 \(\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\) 遍历完后, 我们对所有计算量进行统计:
- 计算所有的 \(\calF(f_{[i]})\) 和 \(\calF(g_{[i]})\), 总共 \(2\sfE(2n)\).
- 每个 \(\psi\), 总共 \(\sfE(2n)\).
- 计算出 \(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\) 次项. 那么有
那么在计算前 \(n\) 次项时, 计算量无外乎还是 \(10\sfE(n)\). 接下来分析后面部分 \(g\) 的计算.
- 还需要对剩余的 \(f\) 部分也计算 DFT, 消耗 \(2\sfE(2n)\).
- 根据 \(f,g_0\) 计算出 \(\delta\), 注意这里 \(f,g_0\) 的 DFT 前面已经算了, 消耗 \(2\sfE(2n)\).
- 计算 \(\delta_1\) 的 DFT, 消耗 \(\sfE(2n)\).
- 避开直接计算 \(\delta_2\), 转而计算 \(\delta_1^2 - \delta_2\), 消耗 \(\sfE(2n)\), 再计算
其 DFT, 消耗 \(\sfE(2n)\). (注意这里如果计算了 \(\delta_2\) 就和原来的做法常数一样了.) - 计算 \(g\) 的后面 \(2n\) 次项, 消耗 \(2\sfE(2n)\).
共计 \(T(3n) = (\frac{13}3 + o(1))\sfM(n)\), 因此 \(T(n) = (1\frac49 + o(1))\sfM(n)\).
注记. 如果采取二阶 Newton 迭代, 读者自行计算不难发现有 \(T(2n)=(18+o(1))\sfE(n)\), 也即 \(T(n)=(1\frac12 +o(1))\sfM(n)\). 比三阶 Newton 迭代效率略差.
改进分块
这里介绍 [2] 引入的方法. 首先的想法是半在线地计算 \(f\cdot g\cdot h\) 不一定要做 \(2m\) 的 DFT, 而是可以改为 \(3m\), 这样就可以在已经有 \(f,g,h\) 的 DFT 的情况下, 在 \((3+o(1))\sfE(n) + O(rn)\) 的时间内算出 \(f\cdot g\cdot h\).
由二阶 Newton 迭代 \(g = g(1-gf)\) 将精度翻倍, 我们只需计算整个 \(f\) 的 \(3m\) 长 DFT, \(g\) 的前半部分的 \(3m\) 长 DFT, 就能算出整个 \(g\).
整个过程中, \(\calF(f)\) 消耗 \(3\sfE(n)\), \(\calF(g)\) 消耗 \(1.5\sfE(n)\), \(g\) 消耗 \(3\sfE(n)\), 总共消耗 \(7.5\sfE(n)\). 也即 \((1.25 + 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.
[2] I. S. Sergeev. “Fast algorithms for elementary operations on complex power series”. In: 20.1 (2010), pp. 25–60. doi: doi:10.1515/dma.2010.002. url: https://doi.org/10.1515/dma.2010.002.
