Toom-Cook 大整數乘法
演算法引入
在 Karatsuba 分治乘法 這篇文章中，我介紹了 Karatsuba 分治乘法。通過將兩個數分成兩段，它的時間複雜度可以達到 \(T(n) = O(n^{\log_23})=O(n^{1.585})\)。這篇文章將推廣 Karatsuba 演算法，進一步討論分治乘法，介紹時間複雜度更低的 Toom-Cook 演算法。其實 Toom-Cook 演算法不是一個單一的演算法，它是一個解決分治高精度乘法問題的一個思想，基於這個思想我們可以給出無數種不同的演算法，而它們的思想和原理大同小異。下面的文章主要會介紹 Toom-Cook 3 Way，最後進行歸納。
寫在前面：這篇文章介紹的演算法對於演算法競賽、實際工作不會有非常大的幫助，文章主要供讀者擴充套件思維；但是如果讀者想要深入學習高精度演算法或者實現一個比較快速的高精度整數類的話，這篇文章應該會提供比較大的幫助。
注意：
- 建議先閱讀 Karatsuba 分治乘法 這篇文章。
- 下文中所有 \(m_1\)，\(m_2\) 指兩個數的資料規模，\(n=\max(m_1,m_2)\)
演算法介紹
在 Karatsuba 演算法中，我們選擇了計算 \((U_0+U_1)(V_0+V_1)\) 減少一次遞迴乘法。為什麼我們要這樣算？這個思路是怎麼想出來的？其實它是可以“湊”出來的。然而，如果我們把一個數分成三段甚至更多，“湊”就變得極為困難，所以我們要找一種“不用湊的方法”。
Toom-Cook 演算法就解決了這個問題。它可以分為大概五個部分：劃分(split)、求值(evaluation)、逐點相乘(pointwise multiply)、插值(interpolate)、重組(recomposition)
劃分
顧名思義，就是把數分成幾段。在這裡，我們將兩個數 U, V 都分成三段，也就是說
\[k = \max(m_1, m_2)/3 \]\[U=U_0+U_1\cdot \beta^k+U_2\cdot \beta^{2k} \]\[V=V_0+V_1\cdot \beta^k+V_2\cdot \beta^{2k} \]其中 \(m_1\)，\(m_2\) 是資料規模。
注意，我們可以設兩個多項式
\[p(x)=U_0+U_1x+U_2x^2 \]\[q(x)=V_0+V_1x+V_2x^2 \]當 \(x=\beta^k\) 時，\(p(x)=U\)，\(q(x)=V\)，我們只需要計算出 \(r(x)=p(x)\cdot q(x)\)
求出 W0 ~ W4 我們就可以得到答案。
求值
這是一個比較關鍵的步驟。我們可以任選 5 個（下面會說明為什麼是 5 個）不同的 \(x\) 帶入 \(p(x)\)，\(q(x)\)。為了簡化計算，我們應該選擇比較小的 \(x\) 帶入。在這裡，我選擇 0, 1, -1, 2, -2 這幾個數帶入。
\[p(0)=U_0 \]\[p(1)=U_0+U_1+U_2 \]\[p(-1)=U_0-U_1+U_2 \]\[p(2)=U_0+2U_1+4U_2 \]\[p(-2)=U_0-2U_1+4U_2 \]\[q(0)=V_0 \]\[q(1)=V_0+V_1+V_2 \]\[q(-1)=V_0-V_1+V_2 \]\[q(2)=V_0+2V_1+4V_2 \]\[q(-2)=V_0-2V_1+4V_2 \]我們發現了一個特別重要的事：計算過程中可能出現負數。具體的實現必須要考慮這一點，並實現有符號整數的加減法。
為了簡化計算，我們可以選擇 \(\infty\) 這個特殊的求值點。你可能會疑問，\(x=\infty\) 時兩個多項式的值不應該都是 \(\infty\) 嗎？但其實，我們求的是
\[{\lim_{x\to +\infty}}\frac{p(x)}{x^2} \]和
\[{\lim_{x\to +\infty}}\frac{q(x)}{x^2} \]它們的結果分別是 \(U_2\) 和 \(V_2\)。我們直接記作 \(p(\infty)\) 和 \(q(\infty)\)。
它們相乘的結果：
\[r(\infty)={\lim_{x\to +\infty}}\frac{r(x)}{x^4}=W_4=p(\infty)\cdot q(\infty)=U_2\cdot V_2 \]亂七八糟說了這麼多，其實就是想表達：\(W_4\) 可以直接由 \(U_2,V_2\) 確定，這樣表示主要是想和求值點的表達統一一下。
我們可以用 \(\infty\) 這個求值點代替上面 -2 的求值點，以減少一些計算。
逐點相乘
我們知道了 \(p(0), ..., p(\infty)\)，知道了 \(q(0), ..., q(\infty)\) ，那麼我們只要把它們依次相乘，就可以得到 \(r(0), ..., r(\infty)\)，而這個逐點相乘的操作是遞迴進行的。
插值
這是另一個比較關鍵的步驟。我們知道了 \(r(0), ..., r(\infty)\) ，我們就可以求 \(W_0, ..., W_4\) 了。具體如何實現？答案是：解方程！
。我們設 \(r(0),r(1),r(-1),r(2),r(\infty)=a,b,c,d,e\)，我們可以得到如下方程組：
\[\begin{align} r(0)=W_0=a\\ r(1)=W_0+W_1+W_2+W_3+W_4=b\\ r(-1)=W_0-W_1+W_2-W_3+W_4=c\\ r(2)=W_0+2W_1+4W_2+8W_3+16W_4=d\\ r(\infty)=W_4=e \end{align} \]直接使用 Wolfram|Alpha, Mathematica 等軟體解就可以了。結果是：
\[\begin{align} & W_0=a\\ & W_1=\frac{-3a+6b-2c-d+12e}{6}\\ & W_2=\frac{-2a+b+c-2e}{2}\\ & W_3=\frac{3a-3b-c+d-12e}{6}\\ & W_4=e \end{align} \]至此 \(W_0\) ~ \(W_4\) 就全部被求出來了。
這裡我們還回答了上面的問題：為什麼要選 5 個值帶入 \(p(x)\)，\(q(x)\)？因為確定一個四次多項式需要 5 個值，多了會導致無意義的多餘計算，少了則無法確定全部係數。
我們也可以換一個理解：設一個矩陣
\[\begin{bmatrix} 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{bmatrix} \]計算出它的逆矩陣，乘以向量 \(\begin{bmatrix}a,b,c,d,e\end{bmatrix}\) 即可。
重組
令 \(x=\beta^k\) 帶入 \(r(x)\) 即可。這一步是 O(n) 的，不再贅述。
時間複雜度分析
這個演算法的求值、插值、重組部分的操作都是 O(n) 的，劃分部分可以做到 O(1) （直接調整指標） 或者 O(n) （把內容複製一遍）。我們可以得到遞迴式：
\[T(n)=5T(n/3)+O(n) \]其中 n 是兩個資料規模中較大的一個。解得：
\[T(n)=O(n^{\log_35})=O(n^{1.465}) \]比 Karatsuba 演算法的 \(O(n^{1.585})\) 略快一些。
更快速的解決方案
我們可以發現：這個演算法的求值、插值部分的操作雖然都是 O(n) 的，但是操作次數很多，常數很大，而我們發現 \(p(1)\) 和 \(p(-1)\) 有重複部分，所以理論上來講我們可以通過一些技巧減少一些操作次數。
[1] 這個網頁展示了幾個不同的 Toom-Cook 演算法最優的求值和插值順序。下面是 bt-3wap-z.gp 這個頁面的內容。只要令 \(x=\beta^k\)，\(y=1\) 就可以實現演算法。需要注意的是，這裡的 <<1, >>1 是二進位制左移右移的意思，作者這麼寫主要是考慮到目前很多主流的大整數運算庫都是二進位制高精度整數。
\\ (C) 2007-2011 Marco Bodrato <http://marco.bodrato.it/>
\\ This code is released under GNU-GPL 3.0+ licence.
U = U2*x^2 + U1*x*y + U0*y^2
V = V2*x^2 + V1*x*y + V0*y^2
\\ P(x,y): P0=(0,1); P1=(2,1); P2=(1,1); P3=(-1,1); P4=(1,0)
\\ Evaluation[1]: 5*2 add/sub, 2 shift; 5mul (n)
W0 = U0 + U2 ; W4 = V0 + V2
W3 = W0 - U1 ; W2 = W4 - V1
W0 = W0 + U1 ; W4 = W4 + V1
W1 = W3 * W2 ; W2 = W0 * W4
W0 =(W0 + U2)<<1 - U0; W4 =(W4 + V2)<<1 - V0
W3 = W0 * W4
W0 = U0 * V0 ; W4 = U2 * V2
\\ Interpolation[1,2]: 8 add/sub, 3 shift, 1 Sdiv (2n)
W3 =(W3 - W1)/3
W1 =(W2 - W1)>>1
W2 = W2 - W0
W3 =(W3 - W2)>>1 - W4<<1
W2 = W2 - W1
\\ Early recomposition[3] (save 0.5 sub):
W3 = W4*x + W3*y
W1 = W2*x + W1*y
W1 = W1 - W3
\\ Final recomposition:
W = W3*x^3+ W1*x*y^2 + W0*y^4
W == U*V
\\ References: http://bodrato.it/papers/#Toom-Cook
\\ [1] Marco BODRATO, "Towards Optimal Toom-Cook Multiplication
\\ for Univariate and Multivariate Polynomials in Characteristic
\\ 2 and 0"; "WAIFI'07 proceedings" (C.Carlet and B.Sunar, eds.)
\\ LNCS#4547, Springer, Madrid, Spain, June 2007, pp. 116-133.
\\ [2] M. BODRATO, A. ZANONI, "Integer and Polynomial
\\ Multiplication: Towards Optimal Toom-Cook Matrices";
\\ "Proceedings of the ISSAC 2007 conference", pp. 17-24
\\ ACM press, Waterloo, Ontario, Canada, July 29-August 1, 2007
\\ [3] Marco BODRATO, "High degree Toom'n'half for balanced and
\\ unbalanced multiplication"; "Proceedings of the 20th IEEE symposium
\\ on Computer Arithmetic", pp 15-22, Tuebingen, Germany, 2011
這篇文章就不提供 C++ 原始碼了，一方面是這個頁面給出的虛擬碼對於演算法的具體實現已經非常詳細且清晰了，C++ 原始碼不會起到很大的幫助；另一方面，如果真的要實現它的話，直接按照虛擬碼照貓畫虎就可以了，不需要刻意背這些步驟。
Toom-Cook 演算法的變種
我們在上述內容介紹的演算法將兩個數各分成三段，但其實我們理解了 Toom-Cook 演算法的思路後就可以推廣，將兩個數各分成四段、五段甚至更多。如果演算法將兩個數各分成 k 段，它的時間複雜度是 \(O(n^{\log_k{2k-1}})=O(n^{\ln(2k-1)/{\ln{k}}})\)。但是分的段數越多，求值和插值兩個過程就越繁瑣，演算法常數就越大。另外，分的段數越多，時間複雜度的提升也越少。比如 Toom-Cook 4 Way 的複雜度是 \(O(n^{1.404})\)，相比於 Toom-Cook 3 Way 提升較小。很多高精度運算庫（如Java BigInteger 和 libtommath）只實現了二重迴圈乘法，Karatsuba 演算法和 Toom-Cook 3 Way 演算法。gmplib 還實現了 toom-4, toom-6 和 toom-8 及各種不平衡乘法的演算法（下面馬上提到）。
上述演算法我們都是將兩個數分成相同的段數。我們還可以將兩個數分成不同的段數，比如說把兩個數中較大的數分成 4 段，較小的數分成 2 段。這就是不平衡的分治乘法 (unbalanced multiplcation)。當兩個數的資料規模相差較大時，這種分段的方法可以提高一點效率。[1] 這個網頁有 Toom-Cook-4x2 的具體步驟。另外值得注意的是，Toom-Cook-3x3 和 Toom-Cook-4x2 的插值部分完全相同（前提是使用的求值點相同），這是因為求值點相同，要解的方程相同，插值部分自然完全相同。實現中可以把插值部分單獨寫成一個函式以減少程式碼量。
參考文獻
[1] Optimal Toom-Cook Polynormial Multiplication, Marco Bodrato.
[2] Toom-Cook multiplication, Wikipedia.