Karatsuba 长乘法算法

• 下载源代码 - 73.9 KB

引言

在之前的系列文章中，我写了关于“小学算术”乘法算法的内容，并强调了它的重要性。

随着数字位数增长，小学算术乘法会被扩展性更好的算法超越。

1960 年代初，苏联数学家阿纳托利·卡拉楚巴发现的第一个比小学算术算法（以下简称 ES）“更快”的算法是卡拉楚巴算法。

ES 算法递归

卡拉楚巴算法可以从标准 ES 算法的一个递归变体开始学习。

示例

首先，将 A 和 B 分成两部分，长度大约是 A（最长操作数）的一半。

将分割的部分存储到变量中：a=932；b=8,225；c=39；d=9,103

A = 932e+8,225

B = 39e+9,103

其中 e = 10,000（我只需要符号 e 来提醒我们稍后需要移位）。

现在，因为 A = ae + b 且 B = ce + d ，所以 AxB = (ae + b)(ce + d) = ace² + ade + bce + bd = ace² + (ad + bc)e + bd

因此，ES 算法递归需要 4 次长乘法：ac、ad、bc、bd（将通过递归计算）。

ac = 932x39=36,348

ad = 932x9,103= 8,483,996

bc = 8,225x39=320,775

bd = 8,225x9,103= 74,872,175

计算完 4 次长乘法后，重新排列并累加部分结果。

卡拉楚巴算法

卡拉楚巴算法的精妙之处在于，它注意到，与其进行两次长乘法 ad 和 bc，我们可以通过一些额外的加减法来节省一次乘法。

我们有表达式

AxB = ace² + (ad + bc)e + bd

我们可以用以下方式替换灰色部分：

ace² + ((a + b)(c + d) - ac - bd)e + bd

“等一下，这有 5 次乘法而不是 4 次！”：是的；但我们已经有了 ac 和 bd，最终只需要 3 次长乘法。

算法步骤是：

1. 按照 ES 递归方式分割数字 A 和 B 并存储到临时变量。
   a=932； b=8,225； c=39； d=9,103
2. 递归计算并存储到其他临时变量，直到乘法可以由硬件完成。
   ac = 932x39=36,348
   bd = 8,225x9,103= 74,872,175
3. 现在，计算 a+b 和 c+d，并将结果存储在临时变量中。
   a_plus_b = a + b = 932 + 8225 = 9157
   c_plus_d = c + d = 39 + 9103 = 9142
4. 通过递归计算 a_plus_b x c_plus_d，并将结果存储在临时变量中。
   temp = a_plus_b x c_plus_d = 9157 x 9142 = 83,713,294
5. 现在从 temp 中减去 ac。
   temp = temp – ac = 83,713,294 - 36,348 = 83,676,946
6. 从 temp 中减去 bd。
   temp = temp – bd = 83,676,946 - 74,872,175 = 8,804,771
7. 现在进行移位和求和。

我们对 4 位数字进行了 3 次长乘法，4 次长加法，2 次减法，需要临时空间（分配）。

用 C 语言看看（已移除注释、断言和内存检查：请在附带的源代码中查找原始版本）

/*
KaratsubaRecursive: beware:
    simpleSum and simpleSub are required to allow operation in place
    (i.e the result of A - B stored in A)
*/

numsize_t KaratsubaRecursive(
    reg_t* A /* first operand */,
    numsize_t m /* size of first operand */,
    reg_t* B /* second operand */,
    numsize_t n /* size of second operand */,
    reg_t* R /* pointer to pre-allocated area where result will be stored */,
    operation simpleMul, /* pointer to primitive multiplication function (ES) */
    operation simpleSum, /* pointer to a long sum function */
    operation simpleSub, /* pointer to a long subtraction function */
    numsize_t simpleMulThreshold /* if operand short than this, use primitive long 
                                    multiplication and stop recursion 
                                    this is useful to create hybrid algorithm as 
                                    discussed later in article
                                  */)
{
        
    numsize_t min = m;
    numsize_t max = n;
    if (m > n) {
        min = n;
        max = m;
    }
   
    if (min == 0)
        return 0;

    if (min <= simpleMulThreshold)
        return simpleMul(A, m, B, n, R);

    numsize_t split_point = (max + 1) >> 1;

    if (split_point > min)
        split_point = min;

    reg_t* a = A + split_point;
    reg_t* b = A;
    reg_t* c = B + split_point;
    reg_t* d = B;

    numsize_t len_a = m - split_point,        
        len_c = n - split_point;

    numsize_t allocsize = ((max - split_point + 2)
        + (max - split_point + 2) +
        ((max - split_point + 2) << 1)
        )
        * sizeof(reg_t);

    reg_t* a_plus_b;
    
    a_plus_b = (reg_t*)ALLOC(allocsize);

    reg_t* c_plus_d = a_plus_b + (max - split_point + 2) ;
    reg_t* z1 = c_plus_d + (max - split_point + 2);  
 
    numsize_t len_z0 = KaratsubaRecursive(a, len_a, c, len_c, R + (split_point << 1), 
                       simpleMul, simpleSum, simpleSub, simpleMulThreshold);   
    numsize_t len_z2 = KaratsubaRecursive(b, split_point, d, split_point, R, 
                       simpleMul, simpleSum, simpleSub, simpleMulThreshold);

    numsize_t len_a_plus_b = simpleSum(a, len_a, b, split_point, a_plus_b);
    numsize_t len_c_plus_d = simpleSum(c, len_c, d, split_point, c_plus_d);

    numsize_t z1_len = KaratsubaRecursive(a_plus_b, len_a_plus_b, c_plus_d, 
           len_c_plus_d, z1, simpleMul, simpleSum, simpleSub, simpleMulThreshold);

    z1_len = simpleSub(z1, z1_len, R, len_z2, z1);
    z1_len = simpleSub(z1, z1_len, R + (split_point << 1), len_z0, z1);
    z1_len = simpleSum(R + split_point, len_z0 + split_point, 
                       z1, z1_len, R+split_point);

    DEALLOC(a_plus_b);

    return z1_len+split_point;
}

卡拉楚巴成本

读者须知：文章到此结束，以下是附加内容。

tl;dr：为什么 ES 在小数字上优于卡拉楚巴？在本章中，我们将计算复杂度并得出混合卡拉楚巴/ES 算法优于两者的结论。

ES 的渐近复杂度为 O(N²)，卡拉楚巴约为 O(N^1.58)。

关于大 O 记法的哲学

在研究卡拉楚巴算法时，我重新发现了自己知道但低估的东西；大 O 测量的是 N 趋向无穷大的算法性能，而不是算法的实际速度。

所有 N 都相等，但有些 N 比其他 N 更相等。

在许多实际应用场景中，常数和较小的因子（当 N 趋向无穷大时变得可忽略不计）并非那么可忽略。

卡拉楚巴成本为 N^1.58 对比 N²，意味着当 N 趋向无穷大时，卡拉楚巴优于 ES。

为什么我的卡拉楚巴实现比 ES 算法慢？我曾假设对于 500 字大小的操作数，卡拉楚巴应该更快，但测量时间表明我低估了较低指数的成本（并且使用 malloc 分配非常小的内存块，这极大地影响了性能，但这并不是重点）。

卡拉楚巴完成任务的步骤更少，但每一步的成本都比 ES 高。

这意味着卡拉楚巴在 N 超过某个我预料之外的阈值后才开始变快。

计算卡拉楚巴算法的复杂度

我们开始计算单个递归步骤的成本，暂时忽略嵌套调用。

• 6 次大小为 N 的长加法，总成本为 6Ns，其中 s 是 CPU 单次加法的成本。
• 3 次大小为 N/2 的递归调用。
• 一些常数成本 K（堆栈准备、分配、样板代码）。

因此，操作数长度为 N 的单个递归步骤的成本是：

Cost(N) = 6Ns + 3Cost(N/2) + K

尝试得到一个非递归公式。

上述公式没有帮助，因为它具有递归性，我们必须对每次递归调用累加成本。

举个例子，假设常数 k 的成本是 100，单次加法成本是 1，单次乘法成本是 4，N 是 8。

请看下图，我们的递归树将有深度 log2(8) 即 3 层 + 根层（第 0 层）。

第 3 层：27 次调用，操作数大小 1：支付 27*100 + 27*(6 次大小为 1 的加法) + 27*4（乘法成本）= 27*(100+6+4)=2970
第 2 层：9 次调用，操作数大小 2：支付 9*100 + 9*(6 * 2 次加法) = 1008
第 1 层：3 次调用，操作数大小 4：支付 3*100 + 3*(6 * 4 次加法) = 372
第 0 层：1 次调用，操作数大小 8：支付 100 + (6 * 8 次加法) = 148

总成本 = 2970+1008+372+148=4498（度量单位是假想的时钟周期，并非真实值）。

泛化

在第 i 层

但我们必须对每个递归层累加上述总和（并且记住要加上叶子节点的乘法成本，设 CPU 乘法成本为 \(\mu\)）。
由于我们有 log2(N) 层，为简化表示，设 \( log_2(N) = M \)，因此

等于

由于 \( 3^i \) 和 \(\left(\frac{3}{2} \right)^i\) 是等比数列，我们知道如何计算它们前 M 项的和。

稍微简化一下。

现在我们可以将成本公式重写为：

稍微简化一下。

代入一些数字来测试公式是否确实有效。

与手动计算结果一致。

请记住，我们必须根据测试来确定常数 k 和 \(\mu\) 的值，但我们有足够的信息来模拟卡拉楚巴何时比 ES 更方便使用。

比较卡拉楚巴和 ES 的成本

如果我们知道 k 和 μ 的值，就能定义卡拉楚巴的成本函数，我将使用一些合理的数字（度量单位是时钟周期）。

ES 成本函数（这是根据实际实验得出的）

卡拉楚巴在大约 N > 1200 后开始更快。现在很清楚为什么我最初的测试显示卡拉楚巴比 ES 慢，因为我需要使用更长的输入大小（并且需要降低 K 的成本……我在叶子节点分配了大约 4 个 int 使用 malloc）。

每次调用的常数成本 K 使该算法在小数字上比 ES 慢。

再次可视化，您在哪里支付递归成本最多？运用直觉：您能看出在哪个深度支付的常数成本最高吗？

正如您所见，约 1/2 的固定成本在叶子节点，与树的大小无关，通过消除叶子节点，我们消除了总 K 成本的 1/2。

所以……我们可以修改算法，当递归树的某个点，如果操作数大小小于 X，则使用小学算术算法并停止递归。

实际算法之间的速度比较

卡拉楚巴 vs ES

真实数据表明卡拉楚巴的常数时间远低于我的估计，我当时试图估算 malloc 的平均成本，但在实现时，由于我使用了 C 语言，我通过使用“alloca”函数在堆栈上动态分配内存，只有在需要超过 4KB 内存时才使用 malloc（但我认为我可以在接近 1KB 时触发 malloc），如果你只需要 16 字节这样的小量内存，调用 malloc 就有点过分了，所以分配的成本比我预期的要低得多。结果是，卡拉楚巴开始优于 ES 的点在 512 到 1024 字之间，而不是像模拟成本那样在 1200 之后。

在那 1024 字之后，它变得越来越好，在 32K 字时快了 5 倍。

卡拉楚巴 vs 卡拉楚巴/ES 混合

现在我们将卡拉楚巴算法的阈值设置为 > 1，这样在递归过程中，如果操作数大小小于或等于阈值，就切换到 ES 算法。

正如您所见，当 N 大于 16（阈值 12）时，卡拉楚巴开始优于 ES，对于更大的 N，甚至无法将 ES 与带阈值的卡拉楚巴进行比较。

我们在 32K 字时快了 30 倍，另外请注意，在某个阈值之后与卡拉楚巴结合使用的 ES 使其在 32K 字时比纯卡拉楚巴快 4 倍。

在讨论了卡拉楚巴（+小学算术，N=12 的阈值）对于大于 16 字的数字很有用之后，我预告一下，这并不是终点，对于更大的数字还有比卡拉楚巴更好的算法……迟早会到来。

感谢您一直跟到这里。

历史

• 2023 年 1 月 12 日：初版