[转]大整数算法[11] Karatsuba乘法

★引子

前面两篇介绍了 Comba 乘法，最后提到当输入的规模很大时，所需的计算时间会急剧增长，因为 Comba 乘法的时间复杂度仍然是 O(n^2)。想要打破乘法中 O(n^2) 的限制，需要从一个完全不同的角度来看待乘法。在下面的乘法算法中，需要使用 x 和 y 这两个大整数的多项式基表达式 f(x) 和 g(x) 来表示。

令 f(x) = a * x + b，g(x) = c * x + d，h(x) = f(x) * g(x)。这里的 x 相当于一个基，比如十进制下，123456 可以表示成 123 * 100 + 456，这时 x 就是 100 了。

★ Karatsuba 乘法原理

既然当输入规模 n 增加时计算量会以 n^2 增加，那么考虑使用分治的方式把一个规模较大的问题分解成若干个规模较小的问题，这样解决这几个规模较小的问题就会比较容易了。

对于 x 和 y 的乘积 z = x * y，可以把 x 和 y 都拆成两段：x 的左半段和右半段分别是 a 和 b，y 的左半段和右半段分别是 c 和 d。

现在考虑乘积 h(x) = f(x) * g(x)：

h(x) = (a * x + b) * (c * x + d)

= a * c * (x^2) + (a * d + b * c) * x + b * d

可以看到，计算 h(x)，需要计算 4 个一半大小的乘法，3 次加法，乘以 x 或者乘以 x^2 可以通过移位实现，所有的加法和移位共用 O(n) 次计算。

设 T(n) 是两个 n 位整数相乘所需的运算总数，则递归方程有：

T(n) = 4 * T(n / 2) + O(n) 当 n > 1

T(n) = O(1) 当 n = 1

解这个递归方程，得到 T(n) = O(n^2)，即时间复杂度和 Comba 乘法是一样的，没有什么改进。要想降低计算的复杂度，必须减少乘法的计算次数。

注意到 a * d + b * c 可以用 a * c 和 b * d 表示：(a + b) * (c + d) - a * c - b * d。

原式 = (a + b) * (c + d) - a * c - b * d

= a * c + a * d + b * c + b * d - a *c - b *d

= a * d + b * c

上边的换算，说明计算 a * d + b * c 可以通过两次加法一次乘法搞定，这样总的乘法次数就减少到 3 次，列出新的递归方程有：

T(n) = 3 * T(n / 2) + O(n) 当 n > 1

T(n) = O(1) 当 n = 1

解递归方程得到 T(n) = O(n^log3)，注意这里的 log 是以 2 为底，近似计算，时间复杂度为：T(n) = O(n^1.585)，比 Comba 乘法的 O(n^2) 要小。

以上就是使用分治的方式计算乘法的原理。上面这个算法，由Anatolii Alexeevitch Karatsuba于1960年提出并于1962年发表，所以也被称为 Karatsuba 乘法。

★实现思路

原理弄明白了，现在整理一下思路：

1. 拆分输入：

计算分割的基：B = MIN(x->used, y->used) / 2

x0，y0 分别存储 x 和 y 的低半部分，x1，y1 分别存储 x 和 y 的高半部分。

2. 计算三个乘积：

x0y0 = x0 * y0 //递归调用乘法 bn_mul_bn

x1y1 = x1 * y1

t1 = x0 + x1

x0 = y0 + y1

t1 = t1 * x0

3. 计算中间项：

x0 = x0y0 + x1y1

t1 = t1 - x0

4. 计算最终乘积：

t1 = t1 * (2^(n * B)) //左移 B个数位

x1y1 = x1y1 * (2^(2 * B)) //左移 2 * B 个数位

t1 = x0y0 + t1

z = t1 + x1y1 //最终结果

★实现

根据上面的思路，Karatsuba乘法的实现代码如下：

static int bn_mul_karatsuba(bignum *z, const bignum *x, const bignum *y)
{
    int ret;
    size_t i, B;
    register bn_digit *pa, *pb, *px, *py;
    bignum x0[1], x1[1], y0[1], y1[1], t1[1], x0y0[1], x1y1[1];
 
    B = BN_MIN(x->used, y->used);
    B >>= 1;
 
    BN_CHECK(bn_init_size(x0, B));
    BN_CHECK(bn_init_size(x1, x->used - B));
    BN_CHECK(bn_init_size(y0, B));
    BN_CHECK(bn_init_size(y1, y->used - B));
    BN_CHECK(bn_init_size(t1, B << 1));
    BN_CHECK(bn_init_size(x0y0, B << 1));
    BN_CHECK(bn_init_size(x1y1, B << 1));
 
    x0->used = y0->used = B;
    x1->used = x->used - B;
    y1->used = y->used - B;
 
    px = x->dp;
    py = y->dp;
    pa = x0->dp;
    pb = y0->dp;
 
    for(i = 0; i < B; i++)
    {
        *pa++ = *px++;
        *pb++ = *py++;
    }
 
    pa = x1->dp;
    pb = y1->dp;
 
    for(i = B; i < x->used; i++)
        *pa++ = *px++;
 
    for(i = B; i < y->used; i++)
        *pb++ = *py++;
 
    bn_clamp(x0);
    bn_clamp(y0);
 
    BN_CHECK(bn_mul_bn(x0y0, x0, y0));
    BN_CHECK(bn_mul_bn(x1y1, x1, y1));
 
    BN_CHECK(bn_add_abs(t1, x0, x1));
    BN_CHECK(bn_add_abs(x0, y0, y1));
    BN_CHECK(bn_mul_bn(t1, x0, t1));
 
    BN_CHECK(bn_add_abs(x0, x0y0, x1y1));
    BN_CHECK(bn_sub_abs(t1, t1, x0));
 
    BN_CHECK(bn_lshd(t1, B));
    BN_CHECK(bn_lshd(x1y1, B << 1));
 
    BN_CHECK(bn_add_abs(t1, x0y0, t1));
    BN_CHECK(bn_add_abs(z, t1, x1y1));
 
clean:
    bn_free(x0);
    bn_free(x1);
    bn_free(y0);
    bn_free(y1);
    bn_free(t1);
    bn_free(x0y0);
    bn_free(x1y1);
 
    return ret;
}

上面的代码，需要很多临时的 bignum 变量，但由于一开始就知道各个 bignum 的大小，所以使用 bn_init_size 函数初始化并且分配指定的数位，避免后面再进行内存的重新分配了，节约了时间。

算法一开始将输入的 x 和 y 拆分成两半，使用三个循环搞定。为了提高效率，这里在指针的前边加上了 register 关键字来暗示在执行过程中尽量把这几个指针变量放到 CPU 的寄存器中，以此来加快变量的访问速度。

乘法的操作是递归调用 bn_mul_bn 函数，这个是有符号数乘法的计算函数，后面会讲，当递归调用到某一个临界点后，乘法的计算会直接调用 Comba 方法进行计算，而不是一直使用 Karatsuba 递归下去。bn_mul_bn 的函数原型是：int bn_mul_bn(bignum *z, const bignum *x, const bignum *y);

需要注意的是，每个计算操作（加法，减法，乘法和移位）在执行过程中都有可能出错，所以必须加上 BN_CHECK 宏进行错误检查，一旦函数调用出错，调到 clean 后面执行内存清理操作。

★分割点

虽然 Karatsuba 乘法执行时所需的单精度乘法比 Comba 方法少，但是也多了一项 O(n) 级别的开销来解一个方程组，用于计算中间项以及合并最后的结果，这就使得 Karatsuba 乘法在应对输入比较小的数字时所需的计算时间会更多。因此在实际操作中，递归计算到一定大小后，就应该改用 Comba 方法计算了。在 bn_mul_bn 函数中，分割点大小是 80（bn_digit 字长为 32 bit 时）或 64（bn_digit 字长为 64 bit 时），当输入两个数的规模有一个小于分割点时，就应该改用 Comba 方法计算乘法，只有当两个数的规模大于或等于分割点才使用 Karatsuba 方法递归计算。

★总结

Karatsuba 算法是比较简单的递归乘法，把输入拆分成 2 部分，不过对于更大的数，可以把输入拆分成 3 部分甚至 4 部分。拆分为 3 部分时，可以使用 Toom-Cook 3-way 乘法，复杂度降低到 O(n^1.465)。拆分为 4 部分时，使用 Toom-Cook 4-way 乘法，复杂度进一步下降到 O(n^1.404)。对于更大的数字，可以拆成 100 段，使用快速傅里叶变换，复杂度接近线性，大约是 O(n^1.149)。可以看出，分割越大，时间复杂度就越低，但是所要计算的中间项以及合并最终结果的过程就会越复杂，开销会增加，因此分割点上升，对于公钥加密，暂时用不到太大的整数，所以使用 Karatsuba 就合适了，不用再去弄更复杂的递归乘法。