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)\)，再帶入 \(\beta^k\)，就可以得到結果。由上一篇文章得知，帶入的步驟只需要 \(O(n)\) 的復雜度。設

\[r(x)=W_0+W_1\cdot \beta^k+W_2\cdot \beta^{2k}+W_3\cdot \beta^{3k}+W_4\cdot \beta^{4k} \]

求出 \(W_0\) ~ \(W_4\) 我們就可以得到答案。

求值

這是一個比較關鍵的步驟。我們可以任選 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), \dots, p(\infty)\)，知道了 \(q(0), \dots, q(\infty)\) ，那么我們只要把它們依次相乘，就可以得到 \(r(0), \dots, r(\infty)\)，而這個逐點相乘的操作是遞歸進行的。

插值

這是另一個比較關鍵的步驟。我們知道了 \(r(0), \dots, r(\infty)\) ，我們就可以求 \(W_0, \dots, 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-3way-z.gp 這個頁面的內容。這份偽代碼使用 0, 1, -1, 2, \(\infty\) 作為求值點，bt-3way-z-2.gp 將求值點 2 換成了 -2。實現中我們令 \(x=\beta^k\)，\(y=1\)。需要注意的是，這里的 <<1, >>1 是二進制左移右移的意思，作者這么寫主要是考慮到目前很多主流的大整數運算庫都是二進制高精度整數。另外，代碼中出現了 /3。這個操作應該使用高精度除以低精度 \(O(n)\) 地處理，算法保證這個除法運算一定是整除。

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