★ 引子
前面兩篇介紹了 Comba 乘法，最後提到當輸入的規模很大時，所需的計算時間會急劇增長，因為 Comba 乘法的時間復雜度仍然是 O(n^2)。想要打破乘法中 O(n^2) 的限制，需要從一個完全不同的角度來看待乘法。在下面的乘法算法中，需要使用 x 和 y 這兩個大整數的多項式基表達式 f(x) 和 g(x) 來表示。
令 f(x) = a * x + b，g(x) = c * x + d，h(x) = f(x) * g(x)。這裡的 x 相當於一個基，比如十進制下，123456 可以表示成 123 * 100 + 456，這時 x 就是 100 了。
★ Karatsuba 乘法原理
既然當輸入規模 n 增加時計算量會以 n^2 增加，那麼考慮使用分治的方式把一個規模較大的問題分解成若干個規模較小的問題，這樣解決這幾個規模較小的問題就會比較容易了。
對於 x 和 y 的乘積 z = x * y，可以把 x 和 y 都拆成兩段：x 的左半段和右半段分別是 a 和 b，y 的左半段和右半段分別是 c 和 d。
現在考慮乘積 h(x) = f(x) * g(x)：
h(x) = (a * x + b) * (c * x + d)
= a * c * (x^2) + (a * d + b * c) * x + b * d
可以看到，計算 h(x)，需要計算 4 個一半大小的乘法，3 次加法，乘以 x 或者乘以 x^2 可以通過移位實現，所有的加法和移位共用 O(n) 次計算。
設 T(n) 是兩個 n 位整數相乘所需的運算總數，則遞歸方程有：
T(n) = 4 * T(n / 2) + O(n) 當 n > 1
T(n) = O(1) 當 n = 1
解這個遞歸方程，得到 T(n) = O(n^2)，即時間復雜度和 Comba 乘法是一樣的，沒有什麼改進。要想降低計算的復雜度，必須減少乘法的計算次數。
注意到 a * d + b * c 可以用 a * c 和 b * d 表示： (a + b) * (c + d) - a * c - b * d。
原式 = (a + b) * (c + d) - a * c - b * d
= a * c + a * d + b * c + b * d - a *c - b *d
= a * d + b * c
上邊的換算，說明計算 a * d + b * c 可以通過兩次加法一次乘法搞定，這樣總的乘法次數就減少到 3 次，列出新的遞歸方程有：
T(n) = 3 * T(n / 2) + O(n) 當 n > 1
T(n) = O(1) 當 n = 1
解遞歸方程得到 T(n) = O(n^log3)，注意這裡的 log 是以 2 為底，近似計算，時間復雜度為：T(n) = O(n^1.585)，比 Comba 乘法的 O(n^2) 要小。
以上就是使用分治的方式計算乘法的原理。上面這個算法，由 Anatolii Alexeevitch Karatsuba 於1960年提出並於1962年發表，所以也被稱為 Karatsuba 乘法。
★ 實現思路
原理弄明白了，現在整理一下思路：
1. 拆分輸入：
計算分割的基：B = MIN(x->used, y->used) / 2
x0，y0 分別存儲 x 和 y 的低半部分，x1，y1 分別存儲 x 和 y 的高半部分。
2. 計算三個乘積：
x0y0 = x0 * y0 //遞歸調用乘法 bn_mul_bn
x1y1 = x1 * y1
t1 = x0 + x1
x0 = y0 + y1
t1 = t1 * x0
3. 計算中間項：
x0 = x0y0 + x1y1
t1 = t1 - x0
4. 計算最終乘積：
t1 = t1 * (2^(n * B)) //左移 B個數位
x1y1 = x1y1 * (2^(2 * B)) //左移 2 * B 個數位
t1 = x0y0 + t1
z = t1 + x1y1 //最終結果
★ 實現
根據上面的思路，Karatsuba乘法的實現代碼如下：
static int bn_mul_karatsuba(bignum *z, const bignum *x, const bignum *y)
{
int ret;
size_t i, B;
register bn_digit *pa, *pb, *px, *py;
bignum x0[1], x1[1], y0[1], y1[1], t1[1], x0y0[1], x1y1[1];
B = BN_MIN(x->used, y->used);
B >>= 1;
BN_CHECK(bn_init_size(x0, B));
BN_CHECK(bn_init_size(x1, x->used - B));
BN_CHECK(bn_init_size(y0, B));
BN_CHECK(bn_init_size(y1, y->used - B));
BN_CHECK(bn_init_size(t1, B << 1));
BN_CHECK(bn_init_size(x0y0, B << 1));
BN_CHECK(bn_init_size(x1y1, B << 1));
x0->used = y0->used = B;
x1->used = x->used - B;
y1->used = y->used - B;
px = x->dp;
py = y->dp;
pa = x0->dp;
pb = y0->dp;
for(i = 0; i < B; i++)
{
*pa++ = *px++;
*pb++ = *py++;
}
pa = x1->dp;
pb = y1->dp;
for(i = B; i < x->used; i++)
*pa++ = *px++;
for(i = B; i < y->used; i++)
*pb++ = *py++;
bn_clamp(x0);
bn_clamp(y0);
BN_CHECK(bn_mul_bn(x0y0, x0, y0));
BN_CHECK(bn_mul_bn(x1y1, x1, y1));
BN_CHECK(bn_add_abs(t1, x0, x1));
BN_CHECK(bn_add_abs(x0, y0, y1));
BN_CHECK(bn_mul_bn(t1, x0, t1));
BN_CHECK(bn_add_abs(x0, x0y0, x1y1));
BN_CHECK(bn_sub_abs(t1, t1, x0));
BN_CHECK(bn_lshd(t1, B));
BN_CHECK(bn_lshd(x1y1, B << 1));
BN_CHECK(bn_add_abs(t1, x0y0, t1));
BN_CHECK(bn_add_abs(z, t1, x1y1));
clean:
bn_free(x0);
bn_free(x1);
bn_free(y0);
bn_free(y1);
bn_free(t1);
bn_free(x0y0);
bn_free(x1y1);
return ret;
}
上面的代碼，需要很多臨時的 bignum 變量，但由於一開始就知道各個 bignum 的大小，所以使用 bn_init_size 函數初始化並且分配指定的數位，避免後面再進行內存的重新分配了，節約了時間。
算法一開始將輸入的 x 和 y 拆分成兩半，使用三個循環搞定。為了提高效率，這裡在指針的前邊加上了 register 關鍵字來暗示在執行過程中盡量把這幾個指針變量放到 CPU 的寄存器中，以此來加快變量的訪問速度。
乘法的操作是遞歸調用 bn_mul_bn 函數，這個是有符號數乘法的計算函數，後面會講，當遞歸調用到某一個臨界點後，乘法的計算會直接調用 Comba 方法進行計算，而不是一直使用 Karatsuba 遞歸下去。 bn_mul_bn 的函數原型是：int bn_mul_bn(bignum *z, const bignum *x, const bignum *y);
需要注意的是，每個計算操作（加法，減法，乘法和移位）在執行過程中都有可能出錯，所以必須加上 BN_CHECK 宏進行錯誤檢查，一旦函數調用出錯，調到 clean 後面執行內存清理操作。
★ 分割點
雖然 Karatsuba 乘法執行時所需的單精度乘法比 Comba 方法少，但是也多了一項 O(n) 級別的開銷來解一個方程組，用於計算中間項以及合並最後的結果，這就使得 Karatsuba 乘法在應對輸入比較小的數字時所需的計算時間會更多。因此在實際操作中，遞歸計算到一定大小後，就應該改用 Comba 方法計算了。在 bn_mul_bn 函數中，分割點大小是 80（bn_digit 字長為 32 bit 時）或 64（bn_digit 字長為 64 bit 時），當輸入兩個數的規模有一個小於分割點時，就應該改用 Comba 方法計算乘法，只有當兩個數的規模大於或等於分割點才使用 Karatsuba 方法遞歸計算。
★ 總結
Karatsuba 算法是比較簡單的遞歸乘法，把輸入拆分成 2 部分，不過對於更大的數，可以把輸入拆分成 3 部分甚至 4 部分。拆分為 3 部分時，可以使用 Toom-Cook 3-way 乘法，復雜度降低到 O(n^1.465)。拆分為 4 部分時，使用 Toom-Cook 4-way 乘法，復雜度進一步下降到 O(n^1.404)。對於更大的數字，可以拆成 100 段，使用快速傅裡葉變換，復雜度接近線性，大約是 O(n^1.149)。可以看出，分割越大，時間復雜度就越低，但是所要計算的中間項以及合並最終結果的過程就會越復雜，開銷會增加，因此分割點上升，對於公鑰加密，暫時用不到太大的整數，所以使用 Karatsuba 就合適了，不用再去弄更復雜的遞歸乘法。
下一篇將使用 Comba 方法和 Karatsuba 方法構建有符號數的乘法。
【回到本系列目錄】
版權聲明
原創博文，轉載必須包含本聲明，保持本文完整，並以超鏈接形式注明作者Starrybird和本文原始地址：http://www.cnblogs.com/starrybird/p/4445566.html
