★ 引子
最近在折腾 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
"
"movl %4, %%edi
"
"movl %5, %%ecx
"
#define MULADDC_WORD_CORE
"lodsl
"
"mull %6
"
"addl %%ecx, %%eax
"
"adcl $0, %%edx
"
"movl %%edx, %%ecx
"
"stosl
"
#define MULADDC_WORD_STOP
"movl %%esi, %0
"
"movl %%edi, %1
"
"movl %%ecx, %2
"
:"=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
"
"movq %4, %%rdi
"
"movq %5, %%rcx
"
#define MULADDC_WORD_CORE
"lodsq
"
"mulq %6
"
"addq %%rcx, %%rax
"
"adcq $0, %%rdx
"
"movq %%rdx, %%rcx
"
"stosq
"
#define MULADDC_WORD_STOP
"movq %%rsi, %0
"
"movq %%rdi, %1
"
"movq %%rcx, %2
"
:"=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