Day 103-Q4 — モンゴメリ乗算（Montgomery Multiplication）

2026-07-26 赤色 Master / Phase 8+ ★★★★★★★★☆ REDCアルゴリズムによる高速mod演算

問題

奇数の法 $M$ と $Q$ 個のクエリが与えられる。各クエリは整数 $A,E$（$0\le Aモンゴメリ乗算（Montgomery Multiplication）の考え方に基づくREDCアルゴリズムで実装する。

法$M$を「2のべき乗$R$」を使った特殊な形（モンゴメリ形式）$\tilde a=a\cdot R\bmod M$に変換しておくことで、以降の乗算をすべて「乗算とビットシフト」だけで行える。$R$を法とする剰余計算がシフトで済むため、除算命令を一切使わずに$a\times b\bmod M$を計算できる。低レベル言語での高速な多倍長mod演算や暗号ライブラリの内部実装で広く使われる実用的技法。

入力形式

M Q
A_1 E_1
:
A_Q E_Q

制約

$3 \le M \le 10^{18}$（$M$は奇数）
$1 \le Q \le 10^5$
$0 \le A_i < M$
$0 \le E_i \le 10^{18}$

入出力例

入力例1

1000000007 2
3 5
2 100

出力例1

243
976371285

$3^5=243$。$2^{100}\bmod 1000000007=976371285$。

概念図

通常表現 ⇄ モンゴメリ形式の変換とREDC a （通常表現） ×R² → redc ã = a·R mod M mont_mul (シフトのみ) 繰り返し二乗法 redc（逆変換） 結果 （通常表現） REDC(t) = ((t mod R)·M' mod R)·M を足して右シフト（除算命令を使わない） M' = -M⁻¹ mod R を事前計算しておくのが鍵

ヒント（段階的開示）

ヒント1: 方向性
$a\times b\bmod m$はPythonでは`(a*b)%m`で一瞬だが、内部的には大きな数の除算を行っている。法$M$を「2のべき乗」に関連づけて考えると、シフト演算だけで済む場面が作れることに気づけるかがカギ。
ヒント2: アプローチ
$M$より大きい2のべき乗$R=2^k$を選ぶ（$M$が奇数なので$\gcd(R,M)=1$は自動的に成立）。整数$a$を「モンゴメリ形式」$\tilde a=a\cdot R\bmod M$に変換しておくと、2つのモンゴメリ形式の積$\tilde a\cdot\tilde b$を「$R$で割って$\bmod M$する」処理（REDC）だけで正しく$\widetilde{a\cdot b}$が得られる。この「$R$で割る」操作がビットシフトのみで実現できる点が核心。$M'=-M^{-1}\bmod R$を前計算しておく。
ヒント3: 誘導（コード骨格） 
r_bits = M.bit_length() + 1
R = 1 << r_bits
mask = R - 1
M_inv = pow(M, -1, R)
M_prime = (-M_inv) & mask
R2 = (R * R) % M

def redc(t):
    s = ((t & mask) * M_prime) & mask
    u = (t + s * M) >> r_bits
    if u >= M:
        u -= M
    return u

def mont_mul(aR, bR):
    return redc(aR * bR)

繰り返し二乗法と組み合わせれば`pow(a,e,m)`と同じ結果を除算なしで計算できる。

模範解答 (Python)

import sys


def solve():
    data = sys.stdin.buffer.read().split()
    idx = 0
    m = int(data[idx]); idx += 1
    q = int(data[idx]); idx += 1

    r_bits = m.bit_length() + 1
    r = 1 << r_bits
    mask = r - 1
    shift = r_bits
    m_inv = pow(m, -1, r)          # m * m_inv ≡ 1 (mod R)
    m_prime = (-m_inv) & mask      # m * m_prime ≡ -1 (mod R)
    r2 = (r * r) % m               # R^2 mod M（モンゴメリ形式への変換用）

    def redc(t):
        # 0 <= t < R*M を仮定し、t*R^{-1} mod M を除算なしで計算する
        s = ((t & mask) * m_prime) & mask
        u = (t + s * m) >> shift
        if u >= m:
            u -= m
        return u

    def mont_mul(ar, br):
        return redc(ar * br)

    out = []
    for _ in range(q):
        a = int(data[idx]); idx += 1
        e = int(data[idx]); idx += 1

        ar = redc((a % m) * r2)         # aをモンゴメリ形式に変換
        result_r = redc(r2)             # 1のモンゴメリ形式 = R mod M = redc(R^2)
        while e > 0:
            if e & 1:
                result_r = mont_mul(result_r, ar)
            ar = mont_mul(ar, ar)
            e >>= 1
        out.append(str(redc(result_r)))  # モンゴメリ形式から通常表現に戻す

    print("\n".join(out))


solve()
計算量: $O(Q\log E)$（各クエリで$O(\log E)$回のモンゴメリ乗算、各乗算は$O(1)$の多倍長演算）。

Step-by-Step 解説

1パラメータの前計算
$R=2^k$（$M$より大きい）を選び、$M'=-M^{-1}\bmod R$を`pow(m,-1,r)`で求める（$M$に依存し使い回せる）。
2REDC関数
`redc(t)`は$0\le t
3変換と逆変換
$a$のモンゴメリ形式化は`redc(a*R2)`、逆変換は`redc(ar)`を1回呼ぶだけでよい。
4モンゴメリ形式での繰り返し二乗法
`mont_mul`を使って通常の繰り返し二乗法を回し、最後に結果を通常表現へ戻す。

よくあるミス

ミス原因正しい書き方
$R$を$M$より小さい2のべき乗にしてしまう$R>M$かつ$\gcd(R,M)=1$がREDCの前提条件$R$のビット数は`M.bit_length()+1`程度、確実に$M$を超えるようにする
$M$が偶数のケースを想定してしまう$R$（2のべき乗）と$M$が互いに素である必要がある本問では$M$は奇数と制約されているが、一般化時は奇数性を確認する
モンゴメリ形式への変換を忘れて生の値のまま演算する非モンゴメリ形式同士の積は正しい結果にならない冪乗の前に必ず$a$と$1$をモンゴメリ形式に変換する
最終結果をモンゴメリ形式のまま出力してしまう`redc`による逆変換をし忘れるループ終了後、必ず`redc(result_r)`で通常表現に戻す

次のステップ

  • 発展: 固定長整数レジスタ上（C/Rust等）でのモンゴメリ乗算実装によるさらなる高速化
  • 発展: バレット還元（Barrett Reduction）との比較実装・使い分け（$M$が毎回変わる場合はバレット還元が有利なことが多い）
  • 次回予告: 最大空き長方形（Largest Empty Rectangle・座標圧縮 + スイープライン）

自己評価

自分の回答

気づき・メモ