圖姆-庫克演算法

圖姆-庫克演算法（Toom–Cook），有時也被稱為Toom-3演算法，由安德魯·圖姆命名，他提出了這種演算法的基本原理，而史蒂芬·庫克則最先用簡潔的形式描述並改進了這種演算法，將其作為大整數的乘法算法。

圖姆-庫克演算法的原理是：對於給定的兩個大整數$a$和$b$，將$a$和$b$分成$k$個較小的部分，每個部分的長度為$l$，並對這些部分執行運算。隨著$k$的增長，可以組合許多乘法子運算，從而降低演算法的整體複雜度，然後再次使用圖姆-庫克演算法遞歸計算乘法子運算，依此類推。Toom-3和圖姆-庫克兩個術語有時會被錯誤的混用，但事實上Toom-3只是圖姆-庫克演算法在$k = 3$時的特例。

Toom-3將9次乘法降低至僅需5次，使其在$\Theta(n^{\log(5)/\log(3)}) \approx \Theta(n^{1.46})$的時間裡運行。通常，Toom-$k$的時間複雜度為$\Theta(c(k)n^e)$，其中$e = \log(2k - 1) / \log(k)$。$n^e$是在乘法子運算上花費的時間，$c$則是花費在對小常數進行的加法和乘法運算上的時間。著名的Karatsuba演算法實際上是圖姆-庫克演算法的特例，在Karatsuba演算法中，原始乘數被拆分成兩個較小的數，而原本的4次乘法運算縮減為3次，使之在$\Theta(n^{\log(3)/\log(2)}) \approx \Theta(n^{1.58})$的時間內完成運算。Toom-1等價於普通的長乘法，具有$\Theta(n^2)$的複雜度。

儘管可以通過增加$k$來使指數$e$任意接近1，但函數$c$增長速度非常快。混合級別圖姆-庫克演算法的增長率直到2005年仍然是一個廣為研究的開放性問題。根據高德納所描述演算法的一種實現，其複雜度可降低至$\Theta(n 2^{\sqrt{2 \log n}} \log n)$。

由於工作時的開銷，當乘數包括較小的數時，圖姆-庫克算法會比長乘法更慢，因此它適用於中等規模的乘法。對於更大規模的數據，則有漸進更快的史恩哈格·施特拉森演算法（複雜度為$\Theta(n \log n \log \log n)$）。

這一演算法由安德魯·圖姆1963年首次描述，並在史蒂芬·庫克1966年的博士學位論文中得到漸進等效的改進。

## 細節

本節將討論對於任意給定 $k$ 值， Toom-$k$究竟是如何運作的，這是馬可·波德拉托對圖姆-庫克多項式乘法的簡化描述。這個演算法包括五個主要步驟：

- 拆分
- 求值
- 點乘
- 插值
- 重組

在典型的大整數實現中，每個整數都表示為$b$進位的數字序列（$b$通常取較大的數）。在此示例中，$b = 10000$，因此每個數字序列對應一組十進位數字（在實踐中，$b$通常取$2$的冪）。設要相乘的兩個大整數$m$、 $n$分別是：

$m$｜=｜$12$｜$3456$｜$7890$｜$1234$｜$5678$｜$9012$

$n$｜=｜$9$｜$8765$｜$4321$｜$9876$｜$5432$｜$1098$

這對乘數實際上比圖姆-庫克演算法通常要處理的數據小很多，在此使用學校里學習的普通乘法可能會更快，但這個示例仍有助於說明圖姆-庫克演算法的工作原理。

### 拆分

第一步是選擇基數$B = b^i$，使得兩個數字$m$和$n$可以分成$k$段大小不超過$B$的數字（例如在Toom-3演算法中，拆分段數應至多為3）。$i$常常根據如下公式求得：

$$i = \max\left\{\left\lfloor\frac{\left\lfloor\log_b m\right\rfloor}{k}\right\rfloor, \left\lfloor\frac{\left\lfloor\log_b n\right\rfloor}{k}\right\rfloor\right\} + 1.$$

我們的示例將演繹Toom-3演算法的運算過程，因此確定$B = b^2 = 10^8$，接著把$m$和$n$拆分為3段，即$m_i$和$n_i$，則有：

$$\begin{align} m_2 & {} = 123456 \\ m_1 & {} = 78901234 \\ m_0 & {} = 56789012 \\ n_2 & {} = 98765 \\ n_1 & {} = 43219876 \\ n_0 & {} = 54321098 \end{align}$$

然後，我們把這些數作為$(k - 1)$階多項式$p$和$q$的係數，with the property that $p(B) = m$ and $q(B) = n$：

$$p(x) = m_2x^2 + m_1x + m_0 = 123456x^2 + 78901234x + 56789012 \,$$

$$q(x) = n_2x^2 + n_1x + n_0 = 98765x^2 + 43219876x + 54321098 \,$$

定義這些多項式的目的在於：如果計算出它們的乘積$r(x) = p(x)q(x)$，我們的答案就會是$r(B) = m \times n$。

如果乘數位數不同，對於$m$、$n$分別取不同的$k$值十分有用，我們將其稱為$k_m$和$k_n$。例如，演算法「Toom-2.5」是指$k_m = 3$且$k_n = 2$時的圖姆-庫克演算法。這時$B = b^i$中的$i$通常被確定為：

$$i = \max\left\{\left\lfloor\frac{\left\lceil\log_b m\right\rceil}{k_m}\right\rfloor, \left\lfloor\frac{\left\lceil\log_b n\right\rceil}{k_n}\right\rfloor\right\}.$$

### 求值

圖姆-庫克演算法包含一種常用的方法，來計算多項式$p(x)$、$q(x)$的乘積。注意，次數為$d$的多項式可以通過$d + 1$個空間中的點確定（例如一次多項式是一條直線，它由兩個點確定）。這個方法是在各個點上求值$p(\cdot)$和$q(\cdot)$，然後把這些點相乘以獲得多項式乘積上的點，最後進行插值以找到其係數。

由於$\deg(pq) = \deg(p) + \deg(q)$，我們將需要$\deg(p) + \deg(q) + 1 = k_m + k_n - 1$個點來確定最終結果$d$。在Toom-3的情況下，$d = 5$。無論選擇什麼點，該演算法都可以工作（有一些小例外，請參閱插值中的矩陣可逆性約束），但為了簡化演算法，最好選擇較小的整數值，例如$0$、$1$、$-1$和$-2$。

無窮大是一個常被使用的不尋常點，其記作$\infty$或$1/0$。求多項式$p$在無窮大時的值，實際上意味著令$p(x) / x^{\deg p}$的上限為$x$ 且趨向無窮大。因此，$p(\infty)$總是其高階係數的值（$m_2$是上文中的係數）。

在我們的Toom-3示例中，我們將使用點$0$、 $1$、 $-1$、 $-2$和$\infty$，這些選擇簡化了求值，如下式子：

$$\begin{array}{lrlrl} p(0) & = & m_0 + m_1(0) + m_2(0)^2 & = & m_0 \\ p(1) & = & m_0 + m_1(1) + m_2(1)^2 & = & m_0 + m_1 + m_2 \\ p(-1) & = & m_0 + m_1(-1) + m_2(-1)^2 & = & m_0 - m_1 + m_2 \\ p(-2) & = & m_0 + m_1(-2) + m_2(-2)^2 & = & m_0 - 2m_1 + 4m_2 \\ p(\infty) & = & m_2 & & \end{array}$$

對於$q$也是如此。在示例中，我們得到的值是：

$p(0)$｜=｜$m_0$｜=｜$56789012$｜=｜$56789012$

$p(1)$｜=｜$m_0 + m_1 + m_2$｜=｜$56789012 + 78901234 + 123456$｜=｜$135813702$

$p(-1)$｜=｜$m_0 - m_1 + m_2$｜=｜$56789012 - 78901234 + 123456$｜=｜$-21988766$

$p(-2)$｜=｜$m_0 - 2m_1 + 4m_2$｜=｜$56789012 - 2 \times 78901234 + 4 \times 123456$｜=｜$-100519632$

$p(\infty)$｜=｜$m_2$｜=｜$123456$｜=｜$123456$

$q(0)$｜=｜$n_0$｜=｜$54321098$｜=｜$54321098$

$q(1)$｜=｜$n_0 + n_1 + n_2$｜=｜$54321098 + 43219876 + 98765$｜=｜$97639739$

$q(-1)$｜=｜$n_0 - n_1 + n_2$｜=｜$54321098 - 43219876 + 98765$｜=｜$11199987$

$q(-2)$｜=｜$n_0 - 2n_1 + 4n_2$｜=｜$54321098 - 2 \times 43219876 + 4 \times 98765$｜=｜$-31723594$

$q(\infty)$｜=｜$n_2$｜=｜$98765$｜=｜$98765$ .

如上所示，這些值可以包括負值。

為了下文的闡述，把這個求值過程視作矩陣向量乘法較為有用。其中，矩陣的每一行都包含求值點之一的冪，且向量包含多項式的係數：

$$\left(\begin{matrix}p(0) \\ p(1) \\ p(-1) \\ p(-2) \\ p(\infty)\end{matrix}\right) = \left(\begin{matrix} 0^0 & 0^1 & 0^2 \\ 1^0 & 1^1 & 1^2 \\ (-1)^0 & (-1)^1 & (-1)^2 \\ (-2)^0 & (-2)^1 & (-2)^2 \\ 0 & 0 & 1 \end{matrix}\right) \left(\begin{matrix}m_0 \\ m_1 \\ m_2\end{matrix}\right) = \left(\begin{matrix} 1 & 0 & 0 \\ 1 & 1 & 1 \\ 1 & -1 & 1 \\ 1 & -2 & 4 \\ 0 & 0 & 1 \end{matrix}\right) \left(\begin{matrix}m_0 \\ m_1 \\ m_2\end{matrix}\right).$$

The dimensions of the matrix are $d$ by $k_m$ for $p$ and $d$ by $k_n$ for $q$。除最後一列的 $1$ 以外，無窮大的行總是$0$。

#### 更快的求值

與上述公式相比，多點求值可能會減少基本運算（加、減）的次數，更快獲得需要的結果。波德拉托 為Toom-3給出的序列如下所示，它是在運行示例的第一個操作數（多項式$p$上進行的）：

$p_0$｜$\leftarrow$｜$m_0 + m_2$｜=｜$56789012 + 123456$｜=｜$56912468$

$p(0)$｜=｜$m_0$｜=｜$56789012$｜=｜$56789012$

$p(1)$｜=｜$p_0 + m_1$｜=｜$56912468 + 78901234$｜=｜$135813702$

$p(-1)$｜=｜$p_0 - m_1$｜=｜$56912468 - 78901234$｜=｜$-21988766$

$p(-2)$｜=｜$(p(-1) + m_2) \times 2 -m_0$｜=｜$(-21988766 + 123456 ) \times 2 - 56789012$｜=｜$-100519632$

$p(\infty)$｜=｜$m_2$｜=｜$123456$｜=｜$123456$

此序列需要進行五次加/減運算，比簡單求值少一次，同時節省了在計算$p(-2)$時乘以$4$的開銷。

### 點乘

與對多項式$p(\cdot)$和$q(\cdot)$所進行的乘法不同，將$p(a)$ 和$q(a)$被求出的值相乘僅涉及整數相乘——這是原始問題的較小實例。我們遞歸調用我們的乘法過程來使每對已求值的點相乘。在實踐中，隨著乘數減小，演算法將逐漸過渡為教科書長乘法。令$r$為多項式乘積，我們將得到：

$r(0)$｜=｜$p(0)q(0)$｜=｜$56789012 \times 54321098$｜=｜$3084841486175176$

$r(1)$｜=｜$p(1)q(1)$｜=｜$135813702 \times 97639739$｜=｜$13260814415903778$

$r(-1)$｜=｜$p(-1)q(-1)$｜=｜$-21988766 \times 11199987$｜=｜$-246273893346042$

$r(-2)$｜=｜$p(-2)q(-2)$｜=｜$-100519632 \times -31723594$｜=｜$3188843994597408$

$r(\infty)$｜=｜$p(\infty)q(\infty)$｜=｜$123456 \times 98765$｜=｜$12193131840$

如上所示，這些值也可以是負數。對於足夠大的數值，這裡是最昂貴的、唯一與$m$、$n$大小不成線性關係的步驟。

### 插值

這一步最為複雜。與求值相反：給定多項式乘積$r(\cdot)$上的$d$點，我們需要確定其係數。換句話說，我們要在右側求解其向量的矩陣方程：

$$\begin{align} \left(\begin{matrix}r(0) \\ r(1) \\ r(-1) \\ r(-2) \\ r(\infty)\end{matrix}\right) & {} = \left(\begin{matrix} 0^0 & 0^1 & 0^2 & 0^3 & 0^4 \\ 1^0 & 1^1 & 1^2 & 1^3 & 1^4 \\ (-1)^0 & (-1)^1 & (-1)^2 & (-1)^3 & (-1)^4 \\ (-2)^0 & (-2)^1 & (-2)^2 & (-2)^3 & (-2)^4 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}\right) \left(\begin{matrix}r_0 \\ r_1 \\ r_2 \\ r_3 \\ r_4\end{matrix}\right) \\ <pre> & {} = </pre> \left(\begin{matrix} 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 & 1 \\ 1 & -2 & 4 & -8 & 16 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}\right) \left(\begin{matrix}r_0 \\ r_1 \\ r_2 \\ r_3 \\ r_4\end{matrix}\right). \end{align}$$

此矩陣的構造與求值步驟中的矩陣相同，不過它是$d \times d$的。我們可以用高斯消元法來求出方程的解，但這樣非常昂貴。根據以下事實：只要求值點的選擇合適，這個矩陣就是可逆的。因此我們有：

$$\begin{align} \left(\begin{matrix}r_0 \\ r_1 \\ r_2 \\ r_3 \\ r_4\end{matrix}\right) & {} = \left(\begin{matrix} 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 & 1 \\ 1 & -2 & 4 & -8 & 16 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}\right)^{-1} \left(\begin{matrix}r(0) \\ r(1) \\ r(-1) \\ r(-2) \\ r(\infty)\end{matrix}\right) \\ & {} = \left(\begin{matrix} <pre> 1 & 0 & 0 & 0 & 0 \\ \tfrac12 & \tfrac13 & -1 & \tfrac16 & -2 \\ -1 & \tfrac12 & \tfrac12 & 0 & -1 \\ </pre> -\tfrac12 & \tfrac16 & \tfrac12 & -\tfrac16 & 2 \\ <pre> 0 & 0 & 0 & 0 & 1 </pre> \end{matrix}\right) \left(\begin{matrix}r(0) \\ r(1) \\ r(-1) \\ r(-2) \\ r(\infty)\end{matrix}\right). \end{align}$$

接下來即要求得該矩陣的向量積。儘管矩陣中包含分數，但所得的係數卻是整數——因此所有這些都可以在整數算術中完成，僅僅是與小常數進行加減乘除。圖姆-庫克設計時面臨的一個困難挑戰就是找到有效的操作順序來計算該乘積。下面是波德拉托為Toom-3找到的一組順序，通過上面的示例演示：

$r_0$｜$\leftarrow$｜$r(0)$｜=｜$3084841486175176$

$r_4$｜$\leftarrow$｜$r(\infty)$｜=｜$12193131840$

$r_3$｜$\leftarrow$｜$(r(-2) - r(1))/3$｜=｜$(3188843994597408 - 13260814415903778)/3$

=｜$-3357323473768790$

$r_1$｜$\leftarrow$｜$(r(1) - r(-1))/2$｜=｜$(13260814415903778 - (-246273893346042))/2$

=｜$6753544154624910$

$r_2$｜$\leftarrow$｜$r(-1) - r(0)$｜=｜$-246273893346042 - 3084841486175176$

=｜$-3331115379521218$

$r_3$｜$\leftarrow$｜$(r_2 - r_3)/2 + 2r(\infty)$｜=｜$(-3331115379521218 - (-3357323473768790))/2 + 2 \times 12193131840$

=｜$13128433387466$

$r_2$｜$\leftarrow$｜$r_2 + r_1 - r_4$｜=｜$-3331115379521218 + 6753544154624910 - 12193131840$

=｜$3422416581971852$

$r_1$｜$\leftarrow$｜$r_1 -r_3$｜=｜$6753544154624910 - 13128433387466$

=｜$6740415721237444.$

現在我們知道多項式乘積$r$：

$$\begin{array}{rrr} r(x) = & {} & 3084841486175176 \\ <pre> & + & 6740415721237444x \\ & + & 3422416581971852x^2 \\ & + & 13128433387466x^3 \\ & + & 12193131840x^4 </pre> \end{array}$$

如果我們使用不同的$k_m$、$k_n$或求值點，矩陣和我們的插值將改變。但是它不依賴於輸入，因此可以對任何給定的參數集進行硬編碼。

### 重組

最後，我們將求出$r(B)$的值以獲得最終結果。很顯然，由於$B$是$b$的冪，因此對$B$的冪的乘法同樣也可以應用於所有以$b$為底數的數值。在這個示例中，$b = 10^4$且$B = b^2 = 10^8$。

$3084$｜$8414$｜$8617$｜$5176$

$6740$｜$4157$｜$2123$｜$7444$

$3422$｜$4165$｜$8197$｜$1852$

$13$｜$1284$｜$3338$｜$7466$

$+$｜$121$｜$9313$｜$1840$

$121$｜$9326$｜$3124$｜$6761$｜$1632$｜$4937$｜$6009$｜$5208$｜$5858$｜$8617$｜$5176$

這實際上是$1234567890123456789012$與$987654321987654321098$的乘積。

## $k$在其他取值時的插值矩陣

這裡我們給出了幾種$k_m$和$k_n$取常見較小值的插值矩陣。

### Toom-1

Toom-1（$k_m = k_n = 1$）需要一個求值點，這裡選擇$0$。它退化為長乘法，並且使用恆等矩陣的插值矩陣。

$$\left(\begin{matrix}1\end{matrix}\right)^{-1} = \left(\begin{matrix}1\end{matrix}\right).$$

### Toom-1.5

Toom-1.5（$k_m = 2, k_n = 1$）需要兩個求值點，這裡選擇$0$和$\infty$，且其插值矩陣就是恆等矩陣。

$$\left(\begin{matrix}1 & 0 \\ 0 & 1\end{matrix}\right)^{-1} = \left(\begin{matrix}1 & 0 \\ 0 & 1\end{matrix}\right).$$

這裡亦退化為長乘法：一個因數的兩個係數都乘以另一個因數的兩個係數。

### Toom-2

Toom-2（$k_m = 2, k_n = 2$）需要三個求值點，這裡選擇$0$、$1$和$\infty$。它與 Karatsuba 演算法相同，其插值矩陣為：

$$\left(\begin{matrix} 1 & 0 & 0 \\ 1 & 1 & 1 \\ 0 & 0 & 1 \end{matrix}\right)^{-1} = \left(\begin{matrix} 1 & 0 & 0 \\ -1 & 1 & -1 \\ 0 & 0 & 1 \end{matrix}\right).$$

### Toom-2.5

Toom-2.5（$k_m = 3, k_n = 2$）需要四個求值點，這裡選擇$0$、$1$、$-1$和$\infty$。它的插值矩陣為：

$$\left(\begin{matrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 0 & 0 & 0 & 1 \end{matrix}\right)^{-1} = \left(\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \tfrac12 & -\tfrac12 & -1 \\ -1 & \tfrac12 & \tfrac12 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}\right).$$