大整數算法系列（十）：Comba 乘法高效實作三種模式詳解

編輯：關於C語言

Comba 乘法核心思想

逐列累加乘積，使用**三精度暫存器** (c2||c1||c0) 保存中間結果，僅在列結束時寫回記憶體，極大減少記憶體存取與進位傳播開銷。

演算法流程

c = a * b
for i = 0 .. a.used + b.used - 2:
    c2, c1, c0 = 0
    for j = max(0, i-b.used+1) .. min(i, a.used-1):
        k = i - j
        (c2||c1||c0) += a[j] * b[k]
    c[i] = c0
    // c1, c2 隱含進位到下一輪
壓縮前導零

三種實作模式

模式 1：雙精度型別可用 (`bn_udbl` 存在)

最常見、最高效。單精度乘法結果放入雙精度變數，再拆高低半部累加。

#define COMBA_INIT  { bn_udbl r; }
#define COMBA_MULADDC \
    r = (bn_udbl)(*px++) * (*py--) + c0; \
    c0 = (bn_digit)r; \
    r = (bn_udbl)c1 + (r >> biL); \
    c1 = (bn_digit)r; \
    c2 += (bn_digit)(r >> biL);
#define COMBA_STOP  }

模式 2：僅有單精度型別 (無 `bn_udbl`)

將每個數拆高/低半字 (half-word)，執行 4 次單精度乘法模擬雙精度乘法。

#define COMBA_INIT \
    { bn_digit a0,a1,b0,b1, t0,t1,r0,r1; }
#define COMBA_MULADDC \
    a0 = (*px & LO_MASK); a1 = *px++ >> HALF_BITS; \
    b0 = (*py & LO_MASK); b1 = *py-- >> HALF_BITS; \
    r0 = a0 * b0; r1 = a1 * b1; \
    t0 = a1 * b0; t1 = a0 * b1; \
    r1 += (t0 >> HALF_BITS) + (t1 >> HALF_BITS); \
    t0 = (t0 << HALF_BITS); t1 = (t1 << HALF_BITS); \
    r0 += t0; r1 += (r0 < t0); \
    r0 += t1; r1 += (r0 < t1); \
    c0 += r0; c1 += (c0 < r0); \
    c1 += r1; c2 += (c1 < r1);
#define COMBA_STOP  }

HALF_BITS = 位元寬度/2 (如 16 for 32-bit digit)。

模式 3：內聯組語 (x86 VC / GCC)

直接使用 MUL 指令產生雙精度乘積，ADC 累加進位。

/* VC 風格 */
#define COMBA_INIT \
    { __asm mov esi, px; __asm mov edi, py; }
#define COMBA_MULADDC \
    __asm lodsd; \
    __asm mov ebx, [edi]; \
    __asm mul ebx; \
    __asm add c0, eax; \
    __asm adc c1, edx; \
    __asm adc c2, 0; \
    __asm sub edi, 4;
#define COMBA_STOP \
    __asm mov px, esi; __asm mov py, edi; }

循環展開優化

內層迴圈按 32、16、8、4、1 步長展開，減少迴圈控制開銷：

for (j = k; j >= 32; j -= 32) { 32次 COMBA_MULADDC }
for (; j >= 16; j -= 16) { 16次 COMBA_MULADDC }
for (; j >= 8; j -= 8) { 8次 COMBA_MULADDC }
for (; j >= 4; j -= 4) { 4次 COMBA_MULADDC }
for (; j > 0; j--) { 1次 COMBA_MULADDC }

完整 `bn_mul_comba` 骨架

static void bn_mul_comba(bignum *z, const bignum *x, const bignum *y) {
    bn_digit c0=0, c1=0, c2=0;
    bn_digit *px, *py, *pz = z->dp;
    size_t nc = z->used;

    for (size_t i = 0; i < nc; ++i) {
        size_t ty = BN_MIN(i, y->used - 1);
        size_t tx = i - ty;
        size_t k = BN_MIN(x->used - tx, ty + 1);
        px = x->dp + tx; py = y->dp + ty;

        c0 = c1; c1 = c2; c2 = 0; // 右移一位 (進位傳遞)

        // 內層展開迴圈 (使用上述 COMBA_* 巨集)
        // ...

        *pz++ = c0;
    }
    bn_clamp(z); // 修剪前導零
}