淺談多項式除法 Part 1: 快速演算法

計算兩個不超過 $n$ 次的複係數一元多項式 $f(x)$ 和 $g(x)$ 的乘積，可以藉由 FFT (Fast Fourier Transform, 快速傅立葉變換) 在 $O(n\log n)$ 時間內求得。更一般的情形，設 $R$ 為一個環且 $f(x), g(x) \in R[x]$，則 $f(x)g(x)$ 可以用 Karatsuba 演算法在 $O(n^{\log_23})$ 時間內求得，也比 naïve 演算法的 $O(n^2)$ 快。假設多項式乘法的時間複雜度是 $O(M(n))$，這篇將介紹如何在 $O(M(n))$ 時間內計算 $f(x)\div g(x)$ 的商式和餘式。

先定義好我們的問題：設 $R$ 為一個有 $1$ 的環，$f(x), g(x) \in R[x]$ 滿足

• $\deg f = n$
• $g(x) \neq 0, \deg g = m \leq n$
• $g$ 的 $x^m$ 項係數為 $1$

我們要找出 $q(x), r(x) \in R[x]$ 使得

• $f(x) = q(x)g(x) + r(x)$
• $\deg r < m$

不難從長除法證明 $q$ 和 $r$ 存在且唯一。接著，引入反轉算子 $\mathscr{R}: R[x] \to R[x]$，定義對於所有的 $h(x) \in R[x]$，

$$(\mathscr{R}h)(x) = \begin{cases}0,&\text{if }h(x) = 0,\\x^ph(\frac{1}{x}),&\text{if }\deg h = p.\end{cases}$$ 事實上，若 $h(x) = \alpha_0 + \sum_{i=1}^p\alpha_ix^i \neq 0$，則 $(\mathscr{R}h)(x) = \alpha_p+\sum_{i=1}^p\alpha_{p-i}x^i$。這樣一來，就有

$$(\mathscr{R}f)(x) = \begin{cases}(\mathscr{R}q)(x)\cdot(\mathscr{R}g)(x),&\text{if }r(x) = 0,\\(\mathscr{R}q)(x)\cdot(\mathscr{R}g)(x)+x^{n-\deg r}(\mathscr{R}r)(x),&\text{otherwise}.\end{cases}$$ 注意當 $r(x) \neq 0$ 時，有 $n-\deg r \geq n-m+1$。所以我們總是有 $$(\mathscr{R}f)(x) \operatorname{mod}x^{n-m+1} = ((\mathscr{R}q)(x)\cdot(\mathscr{R}g)(x)) \operatorname{mod}x^{n-m+1}$$ 如果我們能找到 $(\mathscr{R}g)(x)^{-1} \operatorname{mod}x^{n-m+1}$，就有 $$(\mathscr{R}q)(x) \operatorname{mod}x^{n-m+1} = ((\mathscr{R}f)(x)\cdot(\mathscr{R}g)(x)^{-1}) \operatorname{mod}x^{n-m+1}$$ 又因為 $\deg q = n-m$，事實上我們有 $(\mathscr{R}q)(x) \operatorname{mod}x^{n-m+1} = (\mathscr{R}q)(x)$，因此上式可以改成 $$(\mathscr{R}q)(x) = ((\mathscr{R}f)(x)\cdot(\mathscr{R}g)(x)^{-1}) \operatorname{mod}x^{n-m+1}$$ 一旦求得了 $(\mathscr{R}g)(x)^{-1} \operatorname{mod}x^{n-m+1}$，就能在 $O(M(n))$ 時間內計算出 $q(x)$ 以及 $r(x) = f(x) - q(x)g(x)$。接下來我們把目標放在：給定 $h \in R[x]$ 和 $k \in \mathbb{Z}_{\geq 1}$，其中 $h$ 的常數項為 $1$，求出 $h^{-1}\operatorname{mod}x^k$。

不難發現 $h^{-1}\operatorname{mod}x = 1$。令 $l \in \mathbb{Z}_{\geq 1}$ 並假設 $s := h^{-1} \operatorname{mod} x^l$ 為已知，我們介紹一個在 $O(M(l))$ 時間內求出 $h^{-1}\operatorname{mod}x^{2l}$ 的演算法。

令 $h(x) = h_0 + h_1x^l + (\text{higher terms})$，其中 $h_0, h_1 \in R[x]$ 且 $\deg h_0, \deg h_1 \leq l-1$。我們想找某個 $a, b \in R[x]$，其中 $\deg a, \deg b \leq l-1$，使得 $$\begin{equation}\label{eq1}(a+bx^l)h\operatorname{mod}x^{2l} = (a+bx^l)(h_0+h_1x^l)\operatorname{mod}x^{2l} = (ah_0 + (ah_1+bh_0)x^l)\operatorname{mod}x^{2l} = 1\end{equation}$$ 注意 $$(a+bx^l)h\operatorname{mod}x^l = ((a+bx^l)h\operatorname{mod}x^{2l})\operatorname{mod}x^l = ah_0\operatorname{mod}x^l = ah\operatorname{mod}x^l=1$$ 所以 $a = s$。接著，令 $ah_0 = sh_0 = 1 + cx^l$，其中 $c \in R[x]$，則 $(\ref{eq1})$可以改寫成 $$(1 + (sh_1+bh_0+c)x^l)\operatorname{mod}x^{2l} = 1$$ 因此 $sh_1+bh_0+c=0\operatorname{mod}x^l$，亦即 $b = -s(sh_1+c)\operatorname{mod}x^l$，此時的 $a+bx^l$ 就是 $h^{-1}\operatorname{mod}x^{2l}$。因為 $s, h_0, h_1, c$ 的次數都不超過 $l-1$，由 $h^{-1}\operatorname{mod}x^l$ 計算 $h^{-1}\operatorname{mod}x^{2l}$ 的時間複雜度是 $O(M(l))$。

最後，假設計算 $h^{-1}\operatorname{mod}x^k$ 的總時間為 $T(k)$，利用這個倍增法，可知 $T(k) = T(\lceil\frac{k}{2}\rceil)+O(M(k)) = O(M(k))$，也就是計算 $q(x)$ 和 $r(x)$ 只需要 $O(M(n))$ 時間，和乘法一樣。

我是在這裡學到這個演算法的。至於實作部分，有興趣的話可以參考我放在github上的整係數多項式模板。