大整数相乘的 Toom-Cook 算法

间完成。对固定宽度整数来说比如 32 位或 64 位整数乘法通常确实可以由少量机器指令完成成本可以近似看成常数。但如果是两个几千位、几万位甚至更大的整数相乘情况就不一样了。此时*背后不再是一条简单的 CPU 乘法指令而是一整套大整数乘法算法。这篇文章要讲的就是大整数乘法中的一个经典算法Toom-Cook 算法。从普通乘法说起当数值远超一个机器字或寄存器能表示的范围时就不能再把一次整数乘法视作 (1)O(1) 操作了。为了简化讨论下面把 n 理解成大整数被拆成的“单元数量”。每个单元可以是一个机器字也可以是一段固定长度的二进制位可以快速相乘。小时候学过的竖式乘法本质上就是普通乘法。如果两个整数各有 n 个单元那么普通乘法大概要做 2n2 级别的单元乘法。对于特别大的整数2n2 会很快变得昂贵。Karatsuba 算法是第一个非常经典的优化。它把一个大整数拆成两半本来需要 4 次子乘法但通过代数变形只需要 3 次子乘法。它的复杂度大约是(log⁡23)≈(1.585)O(nlog2​3)≈O(n1.585)这已经比普通乘法好不少。而 Toom-Cook 可以看成是 Karatsuba 的推广Karatsuba 是拆成 2 块Toom-Cook 可以拆成 3 块、4 块、5 块甚至更多块。不同种类的乘法讨论大数乘法问题时很多不同类型的“乘法”成本并不一样。这个首先要明确。机器字乘法如果两个数都能放进一个机器字比如 32 位或 64 位整数那么一次乘法通常可以看成 (1)O(1)。这是我们平时写固定宽度整数乘法时最接近的情况。大整数乘法如果两个数都有很多个单元那么A * B就不是一次 (1)O(1) 操作而是一个规模随整数长度增长的计算问题。普通竖式乘法、Karatsuba、Toom-Cook 优化的都是这一层的乘法。大整乘小常数比如 212a1​这里的 2 和 4 是固定小常数。它们通常可以通过移位和加法完成。从单个机器字的角度看每一步很便宜但如果整个大整数有 n 个单元扫完整个数组仍然需要 ()O(n) 时间。乘以 β 或 βi后面我们会把大整数写成0122Aa0​a1​βa2​β2这里的 β 是分块基数。乘以 β、2β2 只需要移位、搬移和进位处理而不是普通意义上的大整数乘法。所以当我们说 Toom-Cook 把 9 次乘法减少到 5 次乘法时减少的是第 2 类真正昂贵的大整数子乘法。大整数如何拆成多项式假设我们要计算两个大整数 A 和 B 的乘积。要先把它们拆成若干块。比如把 A 拆成 3 块0122Aa0​a1​βa2​β2这里的 β 可以理解成一个很大的进制基数。例如一个大整数内部用数组存储每个数组元素存一段二进制位那么 β 就对应一块的进位单位。A [a2][a1][a0]假如 10β10这就等价于个十百千万了而实际上基数也可以取大得多的数。同理 0122Bb0​b1​βb2​β2接下来不再把 A 和 B 仅仅看成整数而是看成两个多项式在 xβ 时的取值。也就是()0122A(x)a0​a1​xa2​x2()0122B(x)b0​b1​xb2​x2原来的整数 A 和 B就是 ()A(β) 和 ()B(β)也就是函数 ()()()C(x)A(x)B(x)在 xβ 的值Toom-Cook 的核心三步Toom-Cook 的流程可以概括成三步取值 evaluation逐点相乘 pointwise multiplication插值 interpolation下面以 Toom-3 为例说明。就是把整数拆成 3 块。此时 ()A(x) 和 ()B(x) 都是二次多项式。两个二次多项式相乘结果是一个四次多项式()()()C(x)A(x)B(x)四次多项式有 5 个系数()01223344C(x)c0​c1​xc2​x2c3​x3c4​x4所以只要知道它在 5 个不同点上的值就能恢复出整个多项式。第一步取值Toom-3 常用的取值点是0, 1, -1, 2, infinity这里的infinity不是真的无穷大也不是要带入极限。它只是一个记号表示取最高次项系数。对于()0122A(x)a0​a1​xa2​x2有A(infinity) a2其他几个点也很直观A(0) a0A(1) a0 a1 a2A(-1) a0 - a1 a2A(2) a0 2a1 4a2()B(x) 也做同样的事情B(0) b0B(1) b0 b1 b2B(-1) b0 - b1 b2B(2) b0 2b1 4b2B(infinity) b2注意取值过程里虽然出现了 212a1​、424a2​ 这样的写法但这不是递归意义上的“大整数乘法”。因为 2 和 4 是固定小常数它们通常可以通过移位和加法完成。第二步逐点相乘求完值之后我们得到若干组点值。然后在相同点上分别相乘C(0) A(0) * B(0)C(1) A(1) * B(1)C(-1) A(-1) * B(-1)C(2) A(2) * B(2)C(infinity) A(infinity) * B(infinity)这里一共做了 5 次乘法。需要特别注意的是这里的“乘法”不是 CPU 的单条机器乘法指令而是 5 次规模更小的大整数乘法。每个参与相乘的数长度大约是原来的三分之一再加上一点由求值带来的常数级增长。如果直接拆成 3 块做普通分块乘法本来需要 9 次块乘法a0*b0, a0*b1, a0*b2a1*b0, a1*b1, a1*b2a2*b0, a2*b1, a2*b2Toom-3 的关键是用加减、移位、乘以小常数和插值换掉其中一部分昂贵的大整数乘法把 9 次子乘法压缩成 5 次子乘法。对于大整数来说乘法通常比加减法贵得多。所以用更多线性级别的操作换更少的递归乘法通常是划算的。第三步插值经过逐点相乘之后我们知道了 ()C(x) 在几个点上的值。但是最终需要的是 ()C(x) 的系数也就是要恢复()01223344C(x)c0​c1​xc2​x2c3​x3c4​x4这个过程就是插值。四次多项式可以由 5 个不同点上的取值唯一确定所以我们可以从C(0), C(1), C(-1), C(2), C(infinity)恢复出c0, c1, c2, c3, c4恢复出这些系数之后再把 xβ 代回去01223344ABc0​c1​βc2​β2c3​β3c4​β4这样就得到了原来两个大整数的乘积。这里的乘以 βi 主要对应把系数块放到正确的位置上并做进位归一化也不是重新做几次普通意义上的大整数乘法。伪代码如果把上面的流程写成伪代码大致是这样toom3_mul(A, B):split A into a0, a1, a2split B into b0, b1, b2evaluate:p0 a0p1 a0 a1 a2pm1 a0 - a1 a2p2 a0 2a1 4a2pinf a2q0 b0q1 b0 b1 b2qm1 b0 - b1 b2q2 b0 2b1 4b2qinf b2pointwise multiply:r0 p0 * q0r1 p1 * q1rm1 pm1 * qm1r2 p2 * q2rinf pinf * qinfinterpolate:recover c0, c1, c2, c3, c4 from r0, r1, rm1, r2, rinfrecombine:return c0 c1*beta c2*beta^2 c3*beta^3 c4*beta^4上面伪代码里真正昂贵的递归大整数乘法主要是这 5 行r0 p0 * q0r1 p1 * q1rm1 pm1 * qm1r2 p2 * q2rinf pinf * qinfToom-k 的一般形式更一般地如果把整数拆成 k 块那么对应的多项式次数是 −1k−1。两个 −1k−1 次多项式相乘乘积多项式的最高次数是2−22k−2因此它一共有2−12k−1个系数。为了恢复这个乘积多项式我们需要 2−12k−1 个取值点也就需要做 2−12k−1 次规模更小的子乘法。在 k 固定并且取值、插值、搬移、进位等额外操作都近似为线性级别的前提下可以粗略写出递推式()(2−1)(/)()T(n)(2k−1)T(n/k)O(n)根据主定理它的复杂度约为(log⁡(2−1))O(nlogk​(2k−1))几个常见算法可以粗略对比如下普通乘法O(n^2)Karatsuba拆成 2 块做 3 次子乘法指数约 1.585Toom-3拆成 3 块做 5 次子乘法指数约 1.465Toom-4拆成 4 块做 7 次子乘法指数约 1.404从理论上看k 越大复杂度指数越低。但是k 变大之后取值点更多插值更复杂中间结果也更多常数开销会明显增加。因此真实的大整数库不会盲目使用很大的 k而是会根据整数规模在普通乘法、Karatsuba、Toom-Cook、FFT 或 NTT 类算法之间切换。总结Toom-Cook 的整体思路可以概括为大整数太长直接乘很慢|v把大整数按 beta 拆成几块|v把这些块看成多项式系数|v在几个简单点上求值|v在相同点上逐点相乘|v通过插值恢复乘积多项式|v把 x beta 代回去得到最终整数结果其中最关键的洞察是大整数乘法可以转化成多项式乘法而多项式乘法可以通过“取值 逐点相乘 插值”来完成。Toom-Cook 并不是让乘法消失而是把许多直接的分块乘法变成更少次数的子问题乘法。