第7章: モジュラー算術

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

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

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

7.2 Barrett 還元

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

アルゴリズム

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

$\hat{q}$ は真の商 $\lfloor a/m \rfloor$ の近似で、誤差は高々 $1$。 よって $r < 2m$ であり、高々 1 回の減算で正規化できる。

なぜ速いか

$m$ が $n$ ワードなら $2^k = \beta^{2n}$ (ここで $\beta = 2^{64}$) であり、 $\mu = \lfloor \beta^{2n} / m \rfloor \ge \beta^{n}$ となるので、$\mu$ は一般に $n+1$ ワードである ($n$ ワードではない)。したがって $a \cdot \mu$ は $2n$ ワード $\times$ $(n{+}1)$ ワードの乗算になる。

ただし商の近似に必要なのは $\lfloor a\mu / \beta^{2n} \rfloor$、すなわち積のうち $\beta^{2n}$ より上の部分だけである。 全積を作らずに高位部分積 (high product) だけを求める短縮乗算を使えばその分安くなるが、 実際にどこまで削れるかは採用する短縮乗算のアルゴリズムと実装に依存する。 いずれにせよ $O(n^2)$ の除算より軽い。

7.3 Montgomery 乗算

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

Montgomery 形式と REDC

REDC のアルゴリズム

ここが Montgomery の魔法である。$0 \le T < Rm$ を満たす $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 回だけ計算すればよい。

入出力のコスト

あらかじめ $R^2 \bmod m$ を用意しておけば、通常の整数 $a$ の Montgomery 形式は $\text{REDC}(a \cdot (R^2 \bmod m)) = aR \bmod m$ として 1 回の Montgomery 乗算で得られる。 逆変換には $\bar{a} \cdot 1$ を REDC すればよい。 $k$ 回の乗算を Montgomery 空間で行うなら、変換コストは $2$ 回、内部乗算が $k$ 回。 通常剰余より有利になるのは $k \ge 3$ ほどからというのが目安だが、実際の分岐点は実装と法のサイズに依存する。

7.4 CIOS アルゴリズムと sangi の実装

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

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 済みの結果になる。

sangi の汎用 REDC 実装

次に挙げるのは sangi の汎用の還元ルーチンである。名前のとおり REDC そのもので、 上で述べた CIOS の「乗算段」は含んでいない — 積 $T = a \cdot b$ は呼び出し側が先に作り、 この関数はその $T$ を還元するだけである (構成としては SOS 側にあたる)。 CIOS の形にするには、$a b_i$ の累積と $q_i m$ の累積を同じ外側ループの中で交互に行う必要がある。

// IntModular.cpp: mont_redc (汎用 REDC。積 T は呼び出し側が用意する)
// 前提: T は 2n+1 ワード分の領域を持ち、0 ≤ T < Rm を満たす。
// 前提: mpn::addmul_1 は「n ワード分の addmul 結果の最上位から溢れた
//       単一ワードのキャリー」を返す。tn > i+n のケースでは、このキャリー
//       を上位ワードへ伝搬させる必要がある (内側の for ループ)。
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;
        }
    }
    // 結果は t[n..2n-1] だけでなく、その上の桁上げワード t[2n] との組である。
    // t[2n] != 0 なら値は β^n > m なので減算は必須。
    const uint64_t top = (tn > 2 * n) ? t[2 * n] : 0;
    std::memcpy(r, t + n, n * sizeof(uint64_t));
    if (top != 0 || mpn::cmp(r, n, m, n) >= 0) {
        mpn::sub(r, r, n, m, n);                        // 正規化 (r < m に)
    }
}

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

7.5 Montgomery n=8 特化 — sangi のアセンブリ

512 bit (8 ワード) は多倍長モジュラー演算のベンチマークとして代表的なサイズである。 （実際の暗号方式の法はこれとは異なる: RSA-4096 の法 $n$ は 4096 bit = 64 ワード、CRT 内部の $p, q$ でも 2048 bit = 32 ワード。ECDSA P-521 は 521 bit なので 64-bit ワードでは $\lceil 521/64 \rceil = 9$ ワード必要である。） この 1 サイズに特化したアセンブリカーネルを書くと、 汎用の C++ 実装より 30-50% 速くなる (sangi の実測値。効き幅は CPU と比較対象の実装に依存する)。sangi の 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;
}

7.6 Chinese Remainder Theorem

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

が 1 対 1 対応する。この対応は 環同型 $\mathbb{Z}/M\mathbb{Z} \cong \prod_i \mathbb{Z}/m_i\mathbb{Z}$ であり、 加算・減算・乗算のいずれも法ごとに独立に行えるので、 大きな法 $M$ での計算を $k$ 個の小さな法での計算に分割できる。 sangi では 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 のアルゴリズム (逐次構築版, 擬似コード)
// std::integral は C++20 concept。C++17 以前では template <typename T> に置換。
// 注: 混合基数係数 v[i] は m_i 未満なので T に収まるが、最後の再構成
//     a = v[0] + v[1]·m_0 + ... は M = Π m_i の大きさになるため、組込み整数型では
//     容易にオーバーフローする。実装では戻り値と再構成部を多倍長整数型にすること。
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$ では差が小さい。両方の長所短所を理解して用途に応じて選ぶ。

7.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$ について:

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

// IntModular.cpp: choose_window_width (sangi 2026-05 時点)
inline int choose_window_width(size_t expBits) {
    if (expBits <= 24)   return 1;
    if (expBits <= 64)   return 3;
    if (expBits <= 256)  return 4;
    if (expBits <= 1024) return 5;
    if (expBits <= 4096) return 6;
    return 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 は除算なし。 最後に結果を通常形式に戻す。

7.8 定時間化を意識した実装 — powerModSec

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

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

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

定時間冪剰余 (sangi の powerModSec)

sangi の 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 で分岐なし選択
    }
};

さらなる対策

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

• ブラインディング: RSA 復号 $c^d \bmod n$ を例にとると、乱数 $r$ で $c^{\prime} = c\,r^e \bmod n$ と目隠ししてから $m^{\prime} = (c^{\prime})^d \bmod n$ を計算し、$m = m^{\prime} r^{-1} \bmod n$ として戻す。秘密指数 $d$ が毎回異なる入力に作用するので、実行時間から $d$ を推定しにくくなる（本章では $m$ は一般に法を表すが、この項に限り $m$ は平文、法は $n$ である）。
• メモリアクセスパターンの固定: テーブルを規則的に配置し、秘密値に依存せず同じキャッシュライン群を同じ順序で走査する。512 bit なら 1 エントリがちょうど 64 バイト = 1 キャッシュラインに収まるが、それより大きい法では 1 ラインには収まらない。本質は「1 ラインに収める」ことではなく「アクセス列を秘密に依存させない」ことである。
• 電力解析 (DPA) 対策: 物理攻撃まで視野に入れるなら、マスク化計算 (secret sharing) などの対策が必要になる場合がある。

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

7.9 まとめ

• Barrett: 事前計算 $\mu = \lfloor 2^k / m \rfloor$ で除算を 2 回の乗算に置換。単発の剰余計算に最適。
• Montgomery: Montgomery 形式 $\bar{a} = a R \bmod m$ で乗算連鎖を除算なしに。REDC が核。暗号・NTT のホットループで支配的。
• CIOS: 乗算と還元を交互に進める実装パターン。sangi の汎用経路は積を先に作ってから還元する形 (REDC 単独) で、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 は固定ウィンドウ + マスク選択で骨格を定時間化する。ただし下位の Montgomery 乗算にはデータ依存の処理が残っており、全体の定時間性は機械語レベルの検証を経て初めて主張できる。

次章 (第8章 数論関数の高速計算) では、Euler の $\phi$、 Möbius の $\mu$、約数和 $\sigma_k$ といった乗法的関数と、素数計数関数 $\pi(x)$ の高速計算を扱う。

参考文献

• 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. 定時間実装の実例。

関連章

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