图姆-库克算法

图姆-库克算法（英语：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 (n2^{\sqrt {2\log n}}\log n)$。

由于工作时的开销，当乘数包括较小的数时，图姆-库克算法会比长乘法更慢，因此它适用于中等规模的乘法。对于更大规模的数据，则有渐进更快的史恩哈格·施特拉森算法（复杂度为$\Theta (n\log n\log \log n)$）。

这一算法由安德鲁·图姆1963年首次描述，并在斯蒂芬·库克1966年的博士学位论文中得到渐进等效的改进。

细节

本节将讨论对于任意给定 $k$ 值， Toom-$k$究竟是如何运作的，这是马可·波德拉托对图姆-库克多项式乘法的简化描述。这个算法包括五个主要步骤：

1. 拆分
2. 求值
3. 点乘
4. 插值
5. 重组

在典型的大整数实现中，每个整数都表示为$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{aligned}m_{2}&{}=123456\\m_{1}&{}=78901234\\m_{0}&{}=56789012\\n_{2}&{}=98765\\n_{1}&{}=43219876\\n_{0}&{}=54321098\end{aligned}}$

然后，我们把这些数作为$(k-1)$阶多项式$p$和$q$的系数，with the property that $p(B)=m$ and $q(B)=n$：

$p(x)=m_{2}x^{2}+m_{1}x+m_{0}=123456x^{2}+78901234x+56789012\,$

$q(x)=n_{2}x^{2}+n_{1}x+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{aligned}\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)\\&{}=\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{aligned}}$

此矩阵的构造与求值步骤中的矩阵相同，不过它是$d\times d$的。我们可以用高斯消元法来求出方程的解，但这样非常昂贵。根据以下事实：只要求值点的选择合适，这个矩阵就是可逆的。因此我们有：

${\begin{aligned}\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}1&0&0&0&0\\{\tfrac {1}{2}}&{\tfrac {1}{3}}&-1&{\tfrac {1}{6}}&-2\\-1&{\tfrac {1}{2}}&{\tfrac {1}{2}}&0&-1\\-{\tfrac {1}{2}}&{\tfrac {1}{6}}&{\tfrac {1}{2}}&-{\tfrac {1}{6}}&2\\0&0&0&0&1\end{matrix}}\right)\left({\begin{matrix}r(0)\\r(1)\\r(-1)\\r(-2)\\r(\infty )\end{matrix}}\right).\end{aligned}}$

接下来即要求得该矩阵的向量积。尽管矩阵中包含分数，但所得的系数却是整数——因此所有这些都可以在整数算术中完成，仅仅是与小常数进行加减乘除。图姆-库克设计时面临的一个困难挑战就是找到有效的操作顺序来计算该乘积。下面是波德拉托为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\\&+&6740415721237444x\\&+&3422416581971852x^{2}\\&+&13128433387466x^{3}\\&+&12193131840x^{4}\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&{\tfrac {1}{2}}&-{\tfrac {1}{2}}&-1\\-1&{\tfrac {1}{2}}&{\tfrac {1}{2}}&0\\0&0&0&1\end{matrix}}\right).$

另见

• 乘法算法
• 任意精度整数
• 长乘法
• Karatsuba 算法
• 施特拉森算法
• 富尔算法

引用

• D. Knuth. The Art of Computer Programming, Volume 2. Third Edition, Addison-Wesley, 1997. Section 4.3.3.A： Digital methods, pg.294.
• R. Crandall & C. Pomerance. Prime Numbers – A Computational Perspective. Second Edition, Springer, 2005. Section 9.5.1： Karatsuba and Toom–Cook methods, pg.473.
• M. Bodrato. Toward Optimal Toom-Cook （页面存档备份，存于互联网档案馆）. In WAIFI'07, Springer, 2007.

外部链接

• http://gmplib.org/manual/Toom-3_002dWay-Multiplication.html（页面存档备份，存于互联网档案馆）