第6章: モジュラー算術

6.1 基本 — なぜ除算を避けるか

モジュラー剰余 $a \bmod m$ は、定義上は除算 $\lfloor a / m \rfloor$ で決まる:

前章で見たように、多倍長除算は乗算の 2〜5 倍遅い。 ところが暗号や NTT では、同じ法 $m$で 数億回のモジュラー乗算を行うことが多い。 法 $m$ を固定して事前計算を行えば、実行時は乗算 1〜2 回で剰余が求まる。

6.2 Barrett 還元

Barrett (1986) のアイデアは、$\lfloor 2^k / m \rfloor$ を事前計算した「偽の逆数」 $\mu$ を使って、除算を 2 回の乗算に置き換えることである。

アルゴリズム

法 $m$ が $n$ ワードなら、$k = 2n \cdot 64$ として:

$\hat{q}$ は真の商 $\lfloor a/m \rfloor$ の近似で、誤差 $\le 2$。 よって $r$ は最大で $2m$ ほど余分に大きく、高々 2 回の減算で正規化できる。

なぜ速いか

$a \cdot \mu$ は $a$ が $2n$ ワード、$\mu$ が $n$ ワードなので $3n$ ワード乗算 → 実質 $2 \times$ 乗算コスト。 ただし上位 $2n$ ワードしか必要ないので半分のコストで済む (Karatsuba 式短縮乗算)。 除算 $O(n^2)$ より圧倒的に軽い。

6.3 Montgomery 乗算

Montgomery (1985) の発想は、$a \bmod m$ の代わりに $\bar{a} = a R \bmod m$ ($R = 2^{64n}$) という Montgomery 形式で データを保持し、この形式上で乗算 + 還元を「除算なし」で行うことである。

Montgomery 形式と REDC

REDC のアルゴリズム

ここが Montgomery の魔法である。$T$ (最大 $2n$ ワード) が与えられたとき:

1. $u = T \cdot (-m^{-1}) \bmod R$ を計算 (下位 $n$ ワードのみ)。
2. $T' = T + u \cdot m$ を計算。すると $T' \equiv T \pmod m$ かつ $T' \equiv 0 \pmod R$。
3. $T'' = T' / R$ は単なる右シフト (下位 $n$ ワード切り捨て)。
4. 必要なら $T''$ から $m$ を 1 回減算して正規化。

計算は乗算 2 回 ($u$ と $u \cdot m$) と加算、右シフト、高々 1 回の減算のみ。 除算は事前計算の $m^{-1} \bmod R$ に隠されている。 $m$ が固定なら $m^{-1}$ は 1 回だけ計算すればよい。

入出力のコスト

通常の整数 $a$ を Montgomery 形式に変換するには $\bar{a} = a R \bmod m$ を計算 (Montgomery 乗算 1 回)、 逆変換には $\bar{a} \cdot 1$ を REDC すればよい。 $k$ 回の乗算を Montgomery 空間で行うなら、変換コストは $2$ 回、内部乗算が $k$ 回。 $k \ge 3$ ほどから通常剰余より有利になる。

6.4 CIOS アルゴリズムと calx の実装

Montgomery の REDC を実装する際、$T = a \cdot b$ を一度に計算してから還元するSOS (Separated Operand Scanning) と、乗算と還元を交互に進めるCIOS (Coarsely Integrated Operand Scanning) がある。CIOS は中間結果の保存を節約でき、 暗号パラメータ (n=4〜64) で最も速い。

CIOS の構造

内側ループで $b$ の各ワード $b_i$ について:

1. 乗算段: $t \leftarrow t + a \cdot b_i$ (これで 1 ワード分の $T$ が確定)。
2. 還元段: $q \leftarrow t_0 \cdot m_{\text{inv}} \bmod \beta$ を求め、 $t \leftarrow (t + q \cdot m) / \beta$ (下位 1 ワードが 0 に)。

$n$ 回のループ終了時に $t$ が REDC 済みの結果になる。

calx の C++ ベースライン

// IntModular.cpp: mont_redc (汎用 CIOS)
// 前提: mpn::addmul_1 は「n ワード分の addmul 結果の最上位から溢れた
//       単一ワードのキャリー」を返す。tn > i+n のケースでは、このキャリー
//       を上位ワードへ伝搬させる必要がある (内側の while ループ)。
inline void mont_redc(uint64_t* r, uint64_t* t, size_t tn,
                      const uint64_t* m, size_t n, uint64_t m_inv) {
    for (size_t i = 0; i < n; ++i) {
        uint64_t q = t[i] * m_inv;                     // q = T[i] · (-m^{-1}) mod β
        uint64_t carry = mpn::addmul_1(t + i, m, n, q); // T += q·m·β^i の局所結果
        // addmul_1 から返るキャリーを上位ワードへ順次伝搬 (非ゼロの間のみ)
        for (size_t j = i + n; carry && j < tn; ++j) {
            uint64_t sum = t[j] + carry;
            carry = (sum < t[j]) ? 1 : 0;
            t[j] = sum;
        }
    }
    std::memcpy(r, t + n, n * sizeof(uint64_t));
    if (mpn::cmp(r, n, m, n) >= 0) {
        mpn::sub(r, r, n, m, n);                        // 正規化 (T < m に)
    }
}

汎用版は $n$ によらない。特定サイズ (n=4, 8, 16, 32) では個別のアセンブリカーネルに ディスパッチする。

6.5 Montgomery n=8 特化 — calx のアセンブリ

RSA-4096 や ECDSA P-521 といった用途では、法 $m$ が 512 bit (8 ワード) である。 この 1 サイズに特化したアセンブリカーネルを書くと、 汎用 CIOS より 30-50% 速くなる。calx の mpn_mont_mul_8 は 3 段構成で実装されている。

3 フェーズ構成

累算器のレジスタ常駐

8 ワードの累算器 $r_8, r_9, \ldots, r_{15}$ を GPR (r8-r15) に常駐させる。 これでメモリアクセスを大幅に減らせる:

; mpn_x64_mont.asm: Phase 1 の骨格
xor r8d, r8d          ; 8 個の累算器を 0 クリア
xor r9d, r9d
...                   ; r15 までゼロ
mov rdx, [rcx]        ; b[0] を rdx へ (MULX の暗黙オペランド)
ADDMUL_FIRST_8        ; 種乗算: r8 = a[0]·b[0]_lo, r9 に桁上げ
mov QWORD PTR [rdi], r8  ; r8 は確定したので product buffer へ格納
ADDMUL_REST_8         ; r9-r15 に残りの乗算を累積
; b[1]..b[7] で繰り返し

ADDMUL マクロ

; ADDMUL_FIRST_8: b[j] との乗算の最初の 1 ワード (rbp:rbx に a[0]·rdx)
ADDMUL_FIRST_8 MACRO
    xor eax, eax
    mulx rbx, rbp, [rsi]    ; rbp:rbx = a[0] · rdx (b[j])
    adox r8, rbp            ; r8 += rbp (OF を残す)
ENDM

; ADDMUL_REST_8: 残り a[1]..a[7] の積を r9..r15 に累積 (2 本のキャリーチェーン)
ADDMUL_REST_8 MACRO
    mulx rbp, rax, [rsi+8]   ; a[1]·rdx の低半
    adcx r8, rbx             ; 前回の rbx を CF 鎖で r8 に
    adox r9, rax             ; 低半を OF 鎖で r9 に
    mulx rbx, rax, [rsi+16]  ; a[2]·rdx
    adcx r9, rbp
    adox r10, rax
    ; ... 同じパターンで a[3]..a[7] を r11..r15 まで展開
    ; 最後に adcx r15, rbx; adox r15, 0 で両キャリーを合流
ENDM

; REDC_ITER: 還元 1 ラウンド (iter_idx = 0..7)
REDC_ITER MACRO iter_idx
    mov rax, r8
    imul rax, QWORD PTR [rsp+168]   ; rax = r8 · m_inv
    mov rdx, rax                     ; q = rax (MULX の暗黙オペランド)
    ADDMUL_FIRST_8                   ; T += q · m (先頭ワード)
    ADDMUL_REST_8                    ; T += q · m (残り 7 ワード)
    add r15, QWORD PTR [rdi + (8 + iter_idx)*8]  ; 上位 product buffer の合流
    jnc ri_nc
ENDM

MULX / ADCX / ADOX の 2 並列キャリーチェーン

BMI2 の MULX は乗算で CF/OF を変えない。 ADX 拡張の ADCX / ADOX は CF と OF を別々のキャリーフラグとして扱う。これにより 2 本のキャリーチェーンを同時に走らせ、 依存性を半分にできる。Intel/AMD の現代 CPU で Montgomery 乗算の典型的な実装技法である。

ディスパッチ

// IntModular.cpp: n==8 の分岐
if (n == 8 && mpn::detail::has_bmi2_adx()) {
    mpn_mont_mul_8(r, a, b, m, m_inv);
    return;
}

6.6 Chinese Remainder Theorem

互いに素な法 $m_1, m_2, \ldots, m_k$ が与えられたとき、中国剰余定理 (CRT) により

が 1 対 1 対応する。この写像は 群準同型 (加減乗で可換) なので、 大きな法 $M$ での計算を $k$ 個の小さな法での計算に分割できる。 calx では NTT の素数合成にも使うが、 整数 CRT 単独でも RSA 鍵生成や多項式補間で利用する。

Garner のアルゴリズム

残差 $r_i = a \bmod m_i$ から $a$ を復元する効率的な方法:

最後に $a = v_1 + v_2 m_1 + v_3 m_1 m_2 + \cdots$ を計算する。 Garner の利点は $M = \prod m_i$ を明示的に構築せず、各ステップで $m_1 m_2 \cdots m_{i-1}$ を累積しながら逐次合成する点にある (小さな $m_i$ で モジュラ逆元を計算すれば済む)。

// Garner のアルゴリズム (逐次構築版, 擬似コード)
template <std::integral T>
T garner_crt(const std::vector<T>& remainders,
             const std::vector<T>& moduli) {
    std::vector<T> v(moduli.size());
    v[0] = remainders[0];
    for (size_t i = 1; i < moduli.size(); ++i) {
        T t = remainders[i];
        T prod = 1;  // m_0 * m_1 * ... * m_{i-1} mod moduli[i]
        for (size_t j = 0; j < i; ++j) {
            t = ((t - v[j]) * mod_inverse(moduli[j], moduli[i])) % moduli[i];
            if (t < 0) t += moduli[i];
        }
        v[i] = t;
    }
    // 最後に a = v[0] + v[1]·m_0 + v[2]·m_0·m_1 + ... を多倍長で合成
    T a = v.back();
    for (size_t i = moduli.size() - 1; i-- > 0;) a = a * moduli[i] + v[i];
    return a;
}

なお、$M$ を先に計算する単純な Lagrange 型 CRT ($a = \sum r_i \cdot M_i \cdot (M_i^{-1} \bmod m_i) \bmod M$, ここで $M_i = M / m_i$) も 実装上は有効で、小さな $k$ では差が小さい。両方の長所短所を理解して用途に応じて選ぶ。

6.7 高速冪剰余 — スライディングウィンドウ

$m^e \bmod n$ を計算する最も簡単な方法は、 $e$ を 2 進展開して $e$-ビットを左から右 (あるいは右から左) に処理することである (バイナリ法)。これだけで $O(\log e)$ 回の乗算で済む。

$k$ ビットウィンドウ法

$e$ の $k$ 連続ビットをまとめて処理すると、加算連鎖がさらに短くなる。 前処理で $m^{1}, m^{3}, m^{5}, \ldots, m^{2^k - 1}$ (奇数冪) を計算し、 実行時はウィンドウ $w$ について:

calx のウィンドウ幅選択は次の通り:

// IntModular.cpp: choose_window_width
int w = 1;                          // e ≤ 24 bit
if (expBits > 24)  w = 3;           // ≤ 256 bit
if (expBits > 256) w = 4;           // ≤ 1024 bit
if (expBits > 1024) w = 5;          // ≤ 4096 bit
if (expBits > 4096) w = 7;          // 超大サイズ

スライディングウィンドウ

通常の固定ウィンドウでは、ウィンドウが $0$ で始まる場合も乗算を行う必要がある。 スライディングウィンドウでは、ウィンドウを常に $1$ で始まるようにずらし、 連続する $0$ はスキップして平方のみ行う。平均的に $20\text{-}30\%$ の乗算削減。

// IntModular.cpp: powerMod (スライディングウィンドウ, 概略)
int oddTableSize = 1 << (w - 1);   // 奇数冪のみプレコンピュート
uint64_t g[oddTableSize];
g[0] = baseR;                       // m^1 (Montgomery 形式)
if (oddTableSize > 1) {
    uint64_t base2R = mont_mul(baseR, baseR);    // m^2
    for (int i = 1; i < oddTableSize; ++i)
        g[i] = mont_mul(g[i-1], base2R);         // m^3, m^5, m^7, ...
}
// 本体: 指数のビット列を左から走査、ウィンドウ検出で乗算

Montgomery 空間で全計算を行うので、各 mont_mul は除算なし。 最後に結果を通常形式に戻す。

6.8 定時間実装 — powerModSec

暗号用途では、計算時間が秘密鍵の値に依存すると タイミング攻撃で秘密鍵を推定されうる。 1996 年の Kocher の攻撃以来、このリスクは定期的に実用攻撃として報告されている (CVE-2018-0737 など)。

何がタイミングに影響するか

• 分岐: if 文の分岐予測ミスで数サイクル変わる。
• キャッシュアクセスパターン: ルックアップテーブルの index が秘密なら L1 ヒット/ミスで 10 倍変わる。
• データ依存の命令: 除算や乗算が一部 CPU ではオペランドでサイクル数が変わる。

定時間冪剰余 (calx の powerModSec)

calx の powerModSec は次の 3 点を通常の powerMod と変えている:

1. 固定ウィンドウ (スライディングなし): 常に同じ回数の乗算を行う。 ウィンドウが 0 でも $m^0 = 1$ との乗算を実行。
2. マスク選択のテーブル参照: すべてのテーブルエントリを読み、 マスク演算で所望のものを選ぶ。キャッシュアクセスパターンを秘密依存にしない。
3. 条件分岐の除去: if の代わりに算術マスクを使う。

// IntModular.cpp: powerModSec の constant-time table select
int tableSize = 1 << w;  // 奇数のみでなく 0..2^w-1 すべて
std::vector<uint64_t> g_buf(tableSize * n);
// 注: スライディングなし、全エントリ常時準備

auto ct_select = [&](uint64_t* dst, int idx) {
    std::memset(dst, 0, n * sizeof(uint64_t));
    for (int i = 0; i < tableSize; ++i) {
        // i==idx なら 1, それ以外なら 0 → int64_t キャスト後の単項マイナスで
        // 0xFFFF...FFFF (all-ones) または 0x0000...0000 に (2 の補数表現の利用)
        uint64_t mask = -(static_cast<int64_t>(i == idx));
        for (size_t j = 0; j < n; ++j)
            dst[j] |= g_buf[i * n + j] & mask;              // マスク OR で分岐なし選択
    }
};

さらなる対策

完全な定時間実装には、以下も考慮する必要がある:

• ブラインディング: 乱数 $r$ で $m \leftarrow m r^e \bmod n$ としてから計算。出力は $\cdot r^{-1}$ で戻す。
• メモリアクセスのキャッシュライン整列: テーブルエントリをキャッシュライン境界に配置し、1 ラインで収まるよう設計。
• 電力解析 (DPA) 対策: 物理攻撃まで視野に入れるなら、マスク化計算 (secret sharing) が必要。

calx は暗号ライブラリではないため、ブラインディングや DPA 対策は実装していない。 用途に応じて OpenSSL や libsodium 等の確立された暗号ライブラリを使うのが原則である。

6.9 まとめ

• Barrett: 事前計算 $\mu = \lfloor 2^k / m \rfloor$ で除算を 2 回の乗算に置換。単発の剰余計算に最適。
• Montgomery: Montgomery 形式 $\bar{a} = a R \bmod m$ で乗算連鎖を除算なしに。REDC が核。暗号・NTT のホットループで支配的。
• CIOS: 乗算と還元を交互に進める実装パターン。calx のベースライン。n=4, 8, 16, 32 では特化アセンブリにディスパッチ。
• mpn_mont_mul_8: 512 bit 特化。8 累算器を r8-r15 レジスタ常駐。MULX/ADCX/ADOX の 2 キャリーチェーンで並列化。
• CRT: 互いに素な法での並列計算を統合。Garner で $O(n^2)$ 合成。RSA-CRT は標準技法。
• スライディングウィンドウ冪剰余: 奇数冪のプレコンピュート + ウィンドウ検出で乗算を 20-30% 削減。
• 定時間実装: powerModSec は固定ウィンドウ + マスク選択。タイミング攻撃対策は暗号実装の必須要件。

次章 (第7章 GCD) では、除算を繰り返して最大公約数を求める Euclid 法とその高速化 (binary GCD, half-GCD) を扱う。

参考文献

• Barrett, P. "Implementing the Rivest Shamir and Adleman Public Key Encryption Algorithm on a Standard Digital Signal Processor", CRYPTO '86, LNCS 263, pp. 311-323, 1987.
• Montgomery, P.L. "Modular Multiplication without Trial Division", Mathematics of Computation, 44(170), pp. 519-521, 1985.
• Menezes, A.J., van Oorschot, P.C., Vanstone, S.A. Handbook of Applied Cryptography, CRC Press, 1996. 第14章。
• Brent, R.P., Zimmermann, P. Modern Computer Arithmetic, Cambridge University Press, 2010. 第2章。
• Koç, Ç.K., Acar, T., Kaliski, B.S. "Analyzing and Comparing Montgomery Multiplication Algorithms", IEEE Micro, 16(3), pp. 26-33, 1996. CIOS / SOS / FIOS の比較。
• Kocher, P.C. "Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems", CRYPTO '96, LNCS 1109, pp. 104-113.
• Bernstein, D.J. "Curve25519: New Diffie-Hellman Speed Records", PKC '06, LNCS 3958, pp. 207-228. 定時間実装の実例。

関連章

• 第5章 高速除算アルゴリズム — モジュラー剰余の素朴な実装
• 第7章 整数 GCD アルゴリズム — 拡張 Euclid と $m^{-1}$ の計算
• 第9章 NTT 実装比較 — NTT での Montgomery 乗算
• 第10章 素数判定 — Miller-Rabin の冪剰余
• 第11章 整数因数分解 — Pollard rho, 楕円曲線法のモジュラー演算