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) = (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)\).

注記. 如果采取二階 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.