算法引入

在 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} \]

求出 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.