★ 引子
原本打算一篇文章講完，後來發現篇幅會很大，所以拆成兩部分，先講原理，再講實現。實現的話相對復雜，要用到內聯匯編，要考慮不同平台等等。
在大整數計算中，乘法是非常重要的，因為在公鑰密碼學中模冪運算要頻繁使用乘法，所以乘法的性能會直接影響到模冪運算的效率。下面將會介紹兩種乘法：基線乘法和 Comba 乘法，盡管他們的原理和計算看起來十分類似，而且算法的時間復雜度都是 O(n^2)，但是他們的效率差別是很大的。
★ 基線乘法 （Baseline Multiplication）
基線乘法是 Bignum Math 一書中的一個術語，這裡直接拿來用。實際上基線乘法就是通常所說的筆算乘法。先來看看計算 123 * 456 用筆算乘法應怎麼做：
1 2 3
x 4 5 6
-----------------------------------------------
7 3 8
6 1 5
4 9 2
------------------------------------------------
5 6 0 8 8
在筆算算法中，第二個數的每一位分別與第一個數的每一個位相乘，並且將進位傳遞到下一位，比如 3 * 6 = 18 保留本位 8，進位 1 傳遞到下一位。按照筆算的原理，很容易得到基線乘法的算法思路：
設大整數 x 和 y 分別有 x->used 和 y->used 個數位，進位 c，相乘結果放在 z 中，算法的步驟是：
1. 對於下標 i，從 0 到 x->used + y->used - 1 循環，將 z 中每一個數位清 0。其實就是把 z 設為 0。
2. 對於下標 i，從 0 到 y->used - 1 循環，執行如下幾個操作：
2.1 c = 0
2.2 對於下標 j，從 0 到 x->used - 1循環計算：（注：i，j 都是從 0 開始的，r 是一個雙精度變量）
r = z(i + j) + x(j) * y(i) + c。其中 z(i + j)，x(j)，y(i) 分別表示 z 的第 i + j + 1個數位，x 的第 j + 1 個數位，y 的第 i +1 個數位。
z(i + j) = r & BN_MASK //取 r 的低半部分作為 z(i + j) 的本位。
c = r >> biL //取 r 的高半部分作為進位 c。
2.3 z(i + x->used) = c
3. 返回結果 z。
根據上面的算法思路，便可以寫出基線乘法，C 代碼如下：
int bn_baseline_mul(bignum * z, const bignum *x, const bignum *y)
{
int ret;
bn_udbl r;
size_t i, j, c;
size_t *px, *py, *pz;
bignum ta[1], tb[1];
bn_init(ta);
bn_init(tb);
if(x == z)
{
BN_CHECK(bn_copy(ta, x));
x = ta;
}
if(y == z)
{
BN_CHECK(bn_copy(tb, y));
y = tb;
}
BN_CHECK(bn_grow(z, x->used + y->used));
BN_CHECK(bn_set_word(z, 0));
z->used = x->used + y->used;
for(i = 0; i < y->used; i++)
{
c = 0;
px = x->dp;
py = y->dp + i;
pz = z->dp + i;
for(j = 0; j < x->used; j++)
{
r = (*pz) + (bn_udbl)(*px++) * (*py) + c;
*pz++ = (bn_digit)r;
c = (bn_digit)(r >> biL);
}
*pz = c;
}
z->sign = x->sign * y->sign;
bn_clamp(z);
clean:
bn_free(ta);
bn_free(tb);
return ret;
}
上面的代碼中 ta，tb 臨時變量是為了處理類似 x = x * y，y = x * y，x = x * x 之類的輸入和輸出是同一個變量的情況，因為在計算中，z 會先被置 0，如果不用臨時變量存儲 x 或 y，則 x 或 y 的值就會被改變。要注意的是這個乘法並沒有考慮到沒有雙精度變量的情況，如果在 64 位的環境下，沒有 128 bit 的雙精度類型，則處理起來要復雜很多，這裡暫時不講，下一篇會慢慢討論。
★ Comba 乘法 （Comba Multiplication）原理
Comba 乘法以（在密碼學方面）不太出名的 Paul G. Comba 得名。上面的筆算乘法，雖然比較簡單， 但是有個很大的問題：在 O(n^2) 的復雜度上進行計算和向上傳遞進位，看看前面的那個豎式，每計算一次單精度乘法都要計算和傳遞進位，這樣的話就使得嵌套循環的順序性很強，難以並行展開和實現。Comba 乘法則無需進行嵌套進位來計算乘法，所以雖然其時間復雜度和基線乘法一樣，但是速度會快很多。還是以計算 123 * 456 為例：
1 2 3
x 4 5 6
-----------------------------------------------
6 12 18
5 10 15
4 8 12
------------------------------------------------
4 13 28 27 18
4 13 28 28 8
4 13 30 8
4 16 0
5 6
0 5
------------------------------------------------------
5 6 0 8 8
和普通的筆算乘法很類似，只是每一次單精度乘法只是單純計算乘法，不計算進位，進位留到每一列累加後進行。所以原來需要 n * n 次進位，現在最多只需要 2n 次即可。
以上就是 Comba 乘法的原理，不過這裡有個比較嚴重的問題：如何保證累加後結果不溢出。上面的例子，假設單精度數 1 位數，雙精度是兩位數，那萬一累加後的結果超過兩位數則麼辦？那沒辦法，只能用三精度變量了。在大整數算法中，單精度能表示的最大整數是 2^n - 1（n 是單精度變量的比特數），用三個單精度變量 c2，c1，c0 連在一起作為一個三精度變量（高位在左，低位在右），則 c2 || c1 || c0 能表示的最大整數是 2^(3n) - 1，最多能存放 (2^(3n) - 1) / ((2^n - 1)^2) 個單精度乘積結果。當 n = 32 時，能夠存放 4294967298 個單精度乘積結果；當 n = 64 時，能夠存放約 1.845 * 10^19 個單精度乘積結果，而我一開始規定 bignum 不能超過 25600 個數位，這樣使用三精度變量就可以保證累加結果不會溢出了。
有了上面的鋪墊，下面就把 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。
BN_MIN 是宏定義：#define BN_MIN(x, y) (((x) < (y)) ? (x) : (y))
以上就是 Comba 乘法的思路，先計算每一列，然後求和累加，將累加結果求余數得到本位，進位傳遞到下一列。
第三步中，分別計算 tx， ty 和 k 的值，用於控制列的計算，光這麼說比較抽象，所以還是舉個例子比較直觀。假設要計算 12345 * 678。
Index i: 6 5 4 3 2 1 0
-----------------------------------------------------------------------
1 2 3 4 5
x 6 7 8
-----------------------------------------------------------------------
8 16 24 32 40
7 14 21 28 35
6 12 18 24 30
2 4 6 8 7 4 0
------------------------------------------------------------------------
8 3 6 9 9 1 0
過程：a->used = 5， b->used = 3
i = 0，ty = 0，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 1，k = 1，循環計算第一列：5 * 8 = 40
i = 1，ty = 1，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 2，k = 2，循環計算第二列：5 * 7 = 35，4 * 8 = 32
i = 2，ty = 2，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 3，k = 3，循環計算第三列：5 * 6 = 30，4 * 7 = 28，3 * 8 = 24
i = 3，ty = 2，tx = i - ty = 1，a->used - tx = 4，ty + 1 = 3，k = 3，循環計算第四列：4 * 6 = 24，3 * 7 = 21，2 * 8 = 16
i = 4，ty = 2，tx = i - ty = 2，a->used - tx = 3，ty + 1 = 3，k = 3，循環計算第五列：3 * 6 = 18，2 * 7 = 14，1 * 8 = 8
i = 5，ty = 2，tx = i - ty = 3，a->used - tx = 2，ty + 1 = 3，k = 2，循環計算第六列：2 * 6 = 12，1 * 7 = 7
i = 6，ty = 2，tx = i - ty = 4，a->used - tx = 1，ty + 1 = 3，k = 1，循環計算第七列：1 * 6 = 6
可以看出，Comba 乘法每一次都是計算一列，所以如果想進行並行計算的話，那會方便很多：同時計算所有列並累加，最後進行一次進位傳遞。
前面提到基線乘法和 Comba 乘法的時間復雜度是一樣的，這是對於單精度乘法而言，Comba 乘法之所以要比基線乘法快，主要是因為減少了進位傳遞的次數，所以減少了加法的計算量。
★ 總結
篇幅好像有點長了，所以暫時寫到這裡。說了這麼多，其實關鍵就是：都是使用筆算乘法，基線乘法計算每一行的結果，在累加的同時計算進位，Comba 乘法則按列計算，先累加，然後傳遞進位，減少了計算量，所以速度快。這次先把 Comba 乘法的原理講清楚，下一篇文章講講講如何實現 Comba 乘法，包括在沒有雙精度的環境下和在 x86 環境下使用內聯匯編優化速度。
【回到本系列目錄】
版權聲明
原創博文，轉載必須包含本聲明，保持本文完整，並以超鏈接形式注明作者Starrybird和本文原始地址：http://www.cnblogs.com/starrybird/p/4419444.html
