标签：
★ 引子
最近在折腾 wxWidgets，同时拖延症又犯了，所以中断了好久。这次来讲讲单数位乘法，前面讲到 Comba 和 Karatsuba 乘法，这两个算法适合用来处理比较大的整数，但是对于一个大整数和一个单精度数相乘，其效果反而会不好，因为计算量过多。实际上单数位乘法只是基线乘法的一个特例，不存在嵌套循环进位，因此可以通过优化减少计算量。另外与完整的乘法不同的是，单数位乘法不需要什么临时变量存储和内存分配（目标精度增加除外）。
★ 算法思路
单数位乘法类似于计算 1234567890 * 8 这种计算，第二个数只有一位，在大整数中就是一个单精度变量，只需要执行 O(n) 次单精度乘法就可以完成主要的计算。每一次单精度乘法计算完后，执行进位传递。具体的实现思路如下：
计算 z = x * y，其中 z 和 x 是 bignum，y 是无符号的单精度数。
1. olduse = z->used 记录当前 z 使用了多少数位，用于辅助后面的高位清零。
2. 目标精度增加 1 ： z->used = x->used +1
3. z->sign = x->sign （因为 y 为无符号数，所以结果符号只和 x 有关）
4. 初始化进位： u = 0
5. 对于 i 从 0 到 x->used - 1 之间进行循环：
5.1 r = u + x(i) * y //r 是双精度变量
5.2 z(i) = r & BN_MASK0 //取低半部分作为本位 （32 bit 下，BN_MASK0 = 0xFFFFFFFF，其他情况同理，为 2^n - 1）
5.3 u = r >> biL //取高半部分作为进位
6. z(x->used) = u //传递最后一个进位
7. 剩余的高位清零
8. 压缩多余位
这里 y 是无符号数，如果想计算和 -y 的乘积，在计算完后将 z 的符号取反即可。
★ 实现
和 Comba 乘法一样，先给出总体实现，具体细节后面讲。考虑到计算效率和可移植性的问题，第五步的循环关键代码还是写在宏定义里面，然后按照 1，4，8，16 和 32 的步进展开乘法器，减少循环控制的开销。
int bn_mul_word(bignum *z, const bignum *x, const bn_digit y)
{
int ret;
size_t i, olduse;
bn_digit u, *px, *pz;
olduse = z->used;
z->sign = x->sign;
BN_CHECK(bn_grow(z, x->used + 1));
u = 0;
px = x->dp;
pz = z->dp;
for(i = x->used; i >= 32; i -= 32)
{
MULADDC_WORD_INIT
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_STOP
}
for(; i >= 16; i -= 16)
{
MULADDC_WORD_INIT
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_STOP
}
for(; i >= 8; i -= 8)
{
MULADDC_WORD_INIT
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_STOP
}
for(; i >= 4; i -= 4)
{
MULADDC_WORD_INIT
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_CORE MULADDC_WORD_CORE
MULADDC_WORD_STOP
}
for(; i > 0; i--)
{
MULADDC_WORD_INIT
MULADDC_WORD_CORE
MULADDC_WORD_STOP
}
*pz++ = u;
for(i = x->used + 1; i < olduse; i++)
*pz++ = 0;
z->used = x->used + 1;
bn_clamp(z);
clean:
return ret;
}
以上是单数位乘法的总体实现，关键的地方都在宏定义中，下面将讲讲不同环境下的实现方式。
★ 单双精度变量都有的情况
此情况下， bn_digit 和 bn_udbl 同时有定义，最容易实现的。三个宏的定义如下：
#define MULADDC_WORD_INIT { bn_udbl r;
#define MULADDC_WORD_CORE r = u + (bn_udbl)(*px++) * y; *pz++ = (bn_digit)r; u = (bn_digit)(r >> biL);
#define MULADDC_WORD_STOP }
这种情况完全是按照思路实现的，具体原理就不多说了。
★ 只有单精度变量的情况
如果遇到这种情况，则 bn_udbl 无定义，单精度乘法需要转换成 4 个 半精度的乘法来计算，相对比较复杂。具体的实现原理和 Comba 的乘法器类似，参考此处：http://www.cnblogs.com/starrybird/p/4441022.html 。 具体的宏实现如下：
#define MULADDC_WORD_INIT { bn_digit a0, a1, b0, b1; bn_digit t0, t1, r0, r1;
#define MULADDC_WORD_CORE a0 = (*px << biLH) >> biLH; b0 = ( y << biLH) >> biLH; a1 = *px++ >> biLH; b1 = y >> biLH; r0 = a0 * b0; r1 = a1 * b1; t0 = a1 * b0; t1 = a0 * b1; r1 += (t0 >> biLH); r1 += (t1 >> biLH); t0 <<= biLH; t1 <<= biLH; r0 += t0; r1 += (r0 < t0); r0 += t1; r1 += (r0 < t1); r0 += u; r1 += (r0 < u); *pz++ = r0; u = r1;
#define MULADDC_WORD_STOP }
★ 使用内联汇编的情况
C 的内联汇编细节就不多说了，如果你不会可以跳过。
VC x86:
#define MULADDC_WORD_INIT { __asm mov esi, px __asm mov edi, pz __asm mov ecx, u
#define MULADDC_WORD_CORE __asm lodsd __asm mul y __asm add eax, ecx __asm adc edx, 0 __asm mov ecx, edx __asm stosd
#define MULADDC_WORD_STOP __asm mov px, esi __asm mov pz, edi __asm mov u, ecx }
#endif
GCC x86:
#define MULADDC_WORD_INIT { asm ( "movl %3, %%esi \n\t" "movl %4, %%edi \n\t" "movl %5, %%ecx \n\t"
#define MULADDC_WORD_CORE "lodsl \n\t" "mull %6 \n\t" "addl %%ecx, %%eax \n\t" "adcl $0, %%edx \n\t" "movl %%edx, %%ecx \n\t" "stosl \n\t"
#define MULADDC_WORD_STOP "movl %%esi, %0 \n\t" "movl %%edi, %1 \n\t" "movl %%ecx, %2 \n\t" :"=m"(px),"=m"(pz),"=m"(u) :"m"(px),"m"(pz),"m"(u),"m"(y) :"%eax","%ecx","%edx","%esi","%edi" ); }
GCC x64:
#define MULADDC_WORD_INIT { asm ( "movq %3, %%rsi \n\t" "movq %4, %%rdi \n\t" "movq %5, %%rcx \n\t"
#define MULADDC_WORD_CORE "lodsq \n\t" "mulq %6 \n\t" "addq %%rcx, %%rax \n\t" "adcq $0, %%rdx \n\t" "movq %%rdx, %%rcx \n\t" "stosq \n\t"
#define MULADDC_WORD_STOP "movq %%rsi, %0 \n\t" "movq %%rdi, %1 \n\t" "movq %%rcx, %2 \n\t" :"=m"(px),"=m"(pz),"=m"(u) :"m"(px),"m"(pz),"m"(u),"m"(y) :"%rax","%rcx","%rdx","%rsi","%rdi" ); }
★ 总结
算法也很简单，按照 Baseline Multiplication 的方法做就行了，注意关键的地方优化一下。下一篇讲讲平方的计算。
【回到本系列目录】
版权声明
原创博文，转载必须包含本声明，保持本文完整，并以超链接形式注明作者Starrybird和本文原始地址：http://www.cnblogs.com/starrybird/p/4489859.html
标签：
原文地址：http://www.cnblogs.com/starrybird/p/4489859.html
