★ 引子
上一篇文章讲了 Comba 乘法的原理,这次来讲讲如//代码效果参考:http://www.lyjsj.net.cn/wx/art_23737.html
何实现。为了方便移植和充分发挥不同平台下的性能,暂时用了三种不同的实现方式:1、单双精度变量都有的情况。
2、只有单精度变量的情况。
3、可以使用内联汇编的情况。
前面已经介绍了 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。
上面所说的三种不同实现方式,主要体现在 3.5 中的循环。 Comba 乘法的实现代码如下:(暂//代码效果参考:http://www.lyjsj.net.cn/wx/art_23735.html
时不考虑 x = x y 这种输入和输出是同一个变量的情况)?1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495static void bn_mul_comba(bignum z, const bignum x, const bignum y){ bn_digit c0, c1, c2; bn_digit px, py, pz; size_t nc, i, j, k, tx, ty; pz = z->dp; nc = z->used; c0 = c1 = c2 = 0; for(i = 0; i < nc; i++) { ty = BN_MIN(i, y->used - 1); tx = i - ty; k = BN_MIN(x->used - tx, ty + 1); px = x->dp + tx; py = y->dp + ty; c0 = c1; c1 = c2; c2 = 0U; //Comba 32 for(j = k; j >= 32; j -= 32) { COMBA_INIT COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_STOP } //Comba 16 for(; j >= 16; j -= 16) { COMBA_INIT COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_STOP } //Comba 8 for(; j >= 8; j -= 8) { COMBA_INIT COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC //代码效果参考:http://www.lyjsj.net.cn/wz/art_23733.html
COMBA_STOP } //Comba 4 for(; j >= 4; j -= 4) { COMBA_INIT COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_MULADDC COMBA_STOP } //Comba 1 for(; j > 0; j--) { COMBA_INIT COMBA_MULADDC COMBA_STOP } pz++ = c0; } //for i bn_clamp(z);}在 3.5 步中的内循环,乘法的计算按照 1,4,8,16,32 的步进展开,展开的过程需要比较大的空间,但是减少了循环控制的开销,所以可以节省大量时间。执行乘法的关键代码都定义在 COMBA_INIT,COMBA_MULADDC 和 COMBA_STOP 这三个宏内。COMBA_INIT 用于变量定义或者初始化,COMBA_MULADDC 用于计算单精度乘法和累加,COMBA_STOP 用于存储当前计算结果。另外,这里假设 z 的精度已经足够。
★ 单双精度变量都有的情况
这种情况下,因为 bn_digit 和 bn_udbl 同时有定义,所以是最容易实现的一种。三个宏的定义如下:
?1234567891011121314#define COMBA_INIT { &