★ 引子
最近在折腾 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
