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