淺談大數除法 Part 2: 快速演算法

前言

要計算兩個大數的乘積，可以先把這兩個大數轉成兩個多項式，計算它們的乘積，再轉回大數，時間複雜度 $O(n\log n)$；然而，大數除法和多項式除法可說是完全不一樣的問題。以下要介紹如何利用牛頓法把除法問題轉化成乘法問題，並利用倍增 (doubling) 來得到 $O(n\log n)$ 的除法演算法。

本文內容大致上只是把之前我寫的這篇 document 換成中文，但一些證明以及虛擬碼 (pseudocode) 則省略。另，實作部分可以參考我在 github 上的大數模板，讀者也可以在 ZJ b960 測試自己的程式。

動機

以下設 $\xi$ 和 $\eta$ 都是 $A$ 進位正整數，其中 $A \geq 9$，$\xi \geq \eta$，且 $\xi$ 和 $\eta$ 的長度分別為 $m$ 和 $n$，目標求出 $\lfloor\xi/\eta\rfloor$ 的 $A$ 進位表示。因為 $$\xi/\eta = \xi\times 1/\eta$$ 如果我們能把 $1/\eta$ 算得夠精確，就能把除法問題轉換成乘法問題。我們希望能找到某個 $x \leq 1/\eta$ 使得 $$\begin{equation}\label{goal}\lfloor\xi/\eta\rfloor-1 \leq \lfloor\xi x\rfloor \leq \lfloor\xi/\eta\rfloor\end{equation}$$ 不難發現只要 $x$ 滿足 $$\begin{equation}\label{tol1}0 \leq 1/\eta-x \leq A^{-m}\end{equation}$$ 就有 $$\xi/\eta-\xi x = \xi(1/\eta-x) \leq A^m\cdot A^{-m} = 1$$ 也就是 $(\ref{goal})$ 會成立。接下來我們把目標放在找到 $x$ 滿足 $(\ref{tol1})$。

牛頓法

考慮函數 $f(x) := \frac{1}{x} - \eta$。我們找一個 $x_0 \in (0, \frac{1}{\eta}]$，並對於所有的 $k \geq 0$，令 $$x_{k+1}\gets x_k - \frac{f(x_k)}{f'(x_k)} = x_k(2-\eta x_k)$$ 引入 $\epsilon_k := 1-\eta x_k$，則我們有 $$\begin{equation}\label{raweps}\epsilon_{k+1} = 1-\eta x_k(2-\eta x_k) = (1-\eta x_k)^2 = \epsilon_k^2\end{equation}$$ $(\ref{raweps})$ 告訴我們，如果某個 $k_0$ 滿足 $\epsilon_{k_0} \leq A^{-1}$，則對任意 $l \geq 0$ 都有 $\epsilon_{k_0+l} \leq A^{-2^l}$，收斂速度相當快。

捨去低位與倍增

$(\ref{raweps})$ 看起來很棒，可是實際上我們不可能在每次迭代都算出 $x_k$ 的精確值 (不然時間複雜度還是會壞掉)，必須捨去低位並重新估計誤差收斂速率。把 $(\ref{tol1})$ 重寫成 $$\begin{equation}\label{tol2}\frac{A^{n-1}}{\eta}-A^{n-1}x \leq A^{-(m-n+1)}\end{equation}$$ 為了方便，我們 (暫時) 重新設 $\eta \gets \eta/A^{n-1}$ 並設 $x \gets A^{n-1}x$，這樣一來 $(\ref{tol2})$ 就被簡化成 $$\begin{equation}\label{tol3}\frac{1}{\eta}-x \leq A^{-(m-n+1)}\end{equation}$$ 注意這樣一來就有 $1 \leq \eta < A$，因此要滿足 $(\ref{tol3})$，只要找到 $x$ 使得 $$\begin{equation}\label{tol4} 2(1-\eta x)\leq A^{-(m-n+1)} \end{equation}$$ 就行了。

性質一：設 $h$ 為 $\eta$ 的最高兩位 (如當 $A=10$ 且 $\eta = 1.234$ 時，$h = 12$)，取 $x_0 = \lfloor\frac{A^3}{h+1}\rfloor A^{-2}$，則

1. $x_0 \leq \frac{1}{\eta}$
2. $2\epsilon_0 := 2(1-\eta x_0) \leq A^{-1/2}$

性質二：設 $k\geq 1$。假定 $x_{k-1}$ 滿足

1. $x_{k-1} \leq \frac{1}{\eta}$
2. $2\epsilon_{k-1} \leq A^{-2^{k-2}}$

且

1. 令 $y_k$ 為把 $\eta$ 無條件進位至小數第 $2^{k-1}+1$ 位的結果
2. 令 $z_k$ 為把 $y_kx_{k-1}$ 無條件進位至小數第 $2^{k-1}+1$ 位的結果
3. 令 $x_k$ 為把 $x_{k-1}(2-z_k)$ 無條件捨去至小數第 $2^{k-1}+3$ 位的結果

則

1. $x_k \leq \frac{1}{\eta}$
2. $2\epsilon_k \leq A^{-2^{k-1}}$

根據上面兩個性質，我們只要用上面的構造方式，執行迭代直到 $2^{k-1} \geq m-n+1$，就能取 $x \gets x_k$ 以滿足 ($\ref{tol4}$)。

時間複雜度分析

由於 $T(m) = T(\frac{m}{2}) + O(m\log m) = O(m\log m)$，時間複雜度和乘法一樣，但常數大上許多。