⼤整数算法[09]Comba乘法（原理）
★引⼦
原本打算⼀篇⽂章讲完，后来发现篇幅会很⼤，所以拆成两部分，先讲原理，再讲实现。实现的话相对复杂，要⽤到内联汇编，要考虑不同平台等等。
在⼤整数计算中，乘法是⾮常重要的，因为在公钥密码学中模幂运算要频繁使⽤乘法，所以乘法的性能会直接影响到模幂运算的效率。下⾯将会介绍两种乘法：基线乘法和 Comba 乘法，尽管他们的原理和计算看起来⼗分类似，⽽且算法的时间复杂度都是 O(n^2)，但是他们的效率差别是很⼤的。
★基线乘法（Baseline Multiplication）
基线乘法是 Bignum Math ⼀书中的⼀个术语，这⾥直接拿来⽤。实际上基线乘法就是通常所说的笔算乘法。先来看看计算 123 * 456 ⽤笔算乘法应怎么做：
1          2          3
x      4          5          6
-----------------------------------------------
7          3          8
6          1          5
4          9            2
------------------------------------------------
5            6            0          8          8
在笔算算法中，第⼆个数的每⼀位分别与第⼀个数的每⼀个位相乘，并且将进位传递到下⼀位，⽐如 3 * 6 = 18 保留本位 8，进位 1 传递到下⼀位。按照笔算的原理，很容易得到基线乘法的算法思路：
设⼤整数 x 和 y 分别有 x->used 和 y->used 个数位，进位 c，相乘结果放在 z 中，算法的步骤是：
1. 对于下标 i，从 0 到 x->used + y->used - 1 循环，将 z 中每⼀个数位清 0。其实就是把 z 设为 0。
2. 对于下标 i，从 0 到 y->used - 1 循环，执⾏如下⼏个操作：
2.1  c = 0
2.2  对于下标 j，从 0 到 x->used - 1循环计算：（注：i，j 都是从 0 开始的，r 是⼀个双精度变量）
r  = z(i + j) + x(j) * y(i) + c。其中 z(i + j)，x(j)，y(i) 分别表⽰ z 的第 i + j + 1个数位，x 的第 j + 1 个数位，y 的第 i +1 个数位。
z(i + j) = r & BN_MASK          //取 r 的低半部分作为 z(i + j) 的本位。
c = r >> biL                                  //取 r 的⾼半部分作为进位 c。
2.3  z(i + x->used) = c
3. 返回结果 z。
根据上⾯的算法思路，便可以写出基线乘法，C 代码如下：
int bn_baseline_mul(bignum * z, const bignum *x, const bignum *y)
{
int ret;
bn_udbl r;
size_t i, j, c;
size_t *px, *py, *pz;
bignum ta[1], tb[1];
bn_init(ta);
bn_init(tb);
if(x == z)
{
BN_CHECK(bn_copy(ta, x));
x = ta;
}
if(y == z)
{
BN_CHECK(bn_copy(tb, y));
y = tb;
}
BN_CHECK(bn_grow(z, x->used + y->used));
BN_CHECK(bn_set_word(z, 0));
z->used = x->used + y->used;
for(i = 0; i < y->used; i++)
{
c = 0;
px = x->dp;
py = y->dp + i;
pz = z->dp + i;
for(j = 0; j < x->used; j++)
{
r = (*pz) + (bn_udbl)(*px++) * (*py) + c;
*pz++ = (bn_digit)r;
c = (bn_digit)(r >> biL);
}
*pz = c;
}
z->sign = x->sign * y->sign;
bn_clamp(z);
clean:
bn_free(ta);
bn_free(tb);
booth算法乘法例题讲解return ret;
}
上⾯的代码中 ta，tb 临时变量是为了处理类似 x = x * y，y = x * y，x = x * x 之类的输⼊和输出是同⼀个变量的情况，因为在计算中，z 会先被置 0，如果不⽤临时变量存储 x 或 y，则 x 或 y 的值就会被改变。要注意的是这个乘法并没有考虑到没有双精度变量的情况，如果在 64 位的环境下，没有 128 bit 的双精度类型，则处理起来要复杂很多，这⾥暂时不讲，下⼀篇会慢慢讨论。
★ Comba 乘法（Comba  Multiplication）原理
Comba 乘法以（在密码学⽅⾯）不太出名的得名。上⾯的笔算乘法，虽然⽐较简单，但是有个很⼤的问题：在 O(n^2) 的复杂度上进⾏计算和向上传递进位，看看前⾯的那个竖式，每计算⼀次单精度乘法都要计算和传递进位，这样的话就使得嵌套循环的顺序性很强，难以并⾏展开和实现。Comba 乘法则
⽆需进⾏嵌套进位来计算乘法，所以虽然其时间复杂度和基线乘法⼀样，但是速度会快很多。还是以计算 123 * 456 为例：
1            2            3
x        4            5            6
-----------------------------------------------
6          12          18
5          10        15
4          8          12
------------------------------------------------
4          13          28        27        18
4          13          28        28          8
4          13          30          8
4          16          0
5          6
0        5
------------------------------------------------------
5          6            0            8            8
和普通的笔算乘法很类似，只是每⼀次单精度乘法只是单纯计算乘法，不计算进位，进位留到每⼀列累加后进⾏。所以原来需要 n * n 次进位，现在最多只需要 2n 次即可。
以上就是 Comba 乘法的原理，不过这⾥有个⽐较严重的问题：如何保证累加后结果不溢出。上⾯的例⼦，假设单精度数 1  位数，双精度是两位数，那万⼀累加后的结果超过两位数则么办？那没办法，只能⽤三精度变量了。在⼤整数算法中，单精度能表⽰的最⼤整数是 2^n - 1（n 是单精度变量的⽐特数），⽤三个单精度变量 c2，c1，c0 连在⼀起作为⼀个三精度变量（⾼位在左，低位在右），则 c2 || c1 || c0 能表⽰的最⼤整数是
2^(3n) - 1，最多能存放 (2^(3n) - 1) / ((2^n - 1)^2) 个单精度乘积结果。当 n = 32 时，能够存放 42949
67298 个单精度乘积结果；当 n = 64 时，能够存放约 1.845 * 10^19 个单精度乘积结果，⽽我⼀开始规定 bignum 不能超过 25600 个数位，这样使⽤三精度变量就可以保证累加结果不会溢出了。
有了上⾯的铺垫，下⾯就把 Comba 乘法的思路列出来：
计算 c = a * b，c0，c1，c2 为单精度变量。
1. 增加 c 到所需要的精度，并且令 c = 0，c->used = a->used + b->used。
2. 令 c2 = c1 = c0 = 0。
3. 对于 i 从 0 到 c->used - 1 循环，执⾏如下操作：
3.1  ty = BN_MIN(i, b->used - 1)
3.2  tx = i - ty
3.3  k = BN_MIN(a->used - tx, ty + 1)
3.4  三精度变量右移⼀个数位：(c2 || c1 || c0)  =  (0 || c2 || c1)
3.5  对于 j 从 0 到 k - 1 之间执⾏循环，计算：
(c2 || c1 || c0) = (c2 || c1 || c0) + a(tx + j) * b(ty - j)
3.6  c(i) = c0
4. 压缩多余位，返回 c。
BN_MIN 是宏定义：#define BN_MIN(x, y)                  (((x) < (y)) ? (x) : (y))
以上就是 Comba 乘法的思路，先计算每⼀列，然后求和累加，将累加结果求余数得到本位，进位传递到下⼀列。
第三步中，分别计算 tx， ty 和 k 的值，⽤于控制列的计算，光这么说⽐较抽象，所以还是举个例⼦⽐较直观。假设要计算 12345 * 678。
Index i:                    6          5          4            3            2            1            0
-----------------------------------------------------------------------
1          2            3            4            5
x                                  6            7            8
-----------------------------------------------------------------------
8        16        24        32        40
7        14        21        28        35
6        12        18        24        30
2          4          6          8            7            4          0
------------------------------------------------------------------------
8          3          6          9            9            1          0
过程：a->used = 5， b->used = 3
i = 0，ty = 0，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 1，k = 1，循环计算第⼀列：5 * 8 = 40
i = 1，ty = 1，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 2，k = 2，循环计算第⼆列：5 * 7 = 35，4 * 8 = 32
i = 2，ty = 2，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 3，k = 3，循环计算第三列：5 * 6 = 30，4 * 7 = 28，3 * 8 = 24
i = 3，ty = 2，tx = i - ty = 1，a->used - tx = 4，ty + 1 = 3，k = 3，循环计算第四列：4 * 6 = 24，3 * 7 = 21，2 * 8 = 16
i = 4，ty = 2，tx = i - ty = 2，a->used - tx = 3，ty + 1 = 3，k = 3，循环计算第五列：3 * 6 = 18，2 * 7 = 14，1 * 8 = 8
i = 5，ty = 2，tx = i - ty = 3，a->used - tx = 2，ty + 1 = 3，k = 2，循环计算第六列：2 * 6 = 12，1 * 7 = 7
i = 6，ty = 2，tx = i - ty = 4，a->used - tx = 1，ty + 1 = 3，k = 1，循环计算第七列：1 * 6 = 6
可以看出，Comba 乘法每⼀次都是计算⼀列，所以如果想进⾏并⾏计算的话，那会⽅便很多：同时计算所有列并累加，最后进⾏⼀次进位传递。
前⾯提到基线乘法和 Comba 乘法的时间复杂度是⼀样的，这是对于单精度乘法⽽⾔，Comba 乘法之所以要⽐基线乘法快，主要是因为减少了进位传递的次数，所以减少了加法的计算量。
★总结
篇幅好像有点长了，所以暂时写到这⾥。说了这么多，其实关键就是：都是使⽤笔算乘法，基线乘法计算每⼀⾏的结果，在累加的同时计算进位，Comba 乘法则按列计算，先累加，然后传递进位，减少了计算量，所以速度快。这次先把 Comba 乘法的原理讲清楚，讲讲讲如何实现 Comba 乘法，包括在没有双精度的环境下和在 x86 环境下使⽤内联汇编优化速度。
【】

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。