Day 122-Q1 | モンゴメリ乗算（Montgomery Multiplication）

問題

奇数 $\mathrm{MOD}$ が与えられる。$Q$ 個のクエリが与えられるので、各クエリ $(a, e)$ に対して $a^e \bmod \mathrm{MOD}$ を出力せよ。ただし、繰り返し二乗法の内部で行う「掛け算のたびに mod を取る」処理を、通常の除算命令（%）を使わずに、モンゴメリ乗算（REDCアルゴリズム）を用いて実装すること。

入力形式

MOD Q
a_1 e_1
...
a_Q e_Q

制約

$\mathrm{MOD}$ は奇数

$3 \le \mathrm{MOD} \le 10^9$

$1 \le Q \le 10^5$

$0 \le a_i < \mathrm{MOD}$

$0 \le e_i \le 10^{18}$

入出力例

入力例1

出力例1

54
1
30

入力例2

999999937 2
123456789 999999936
999999936 5

出力例2

1
999999936

概念図: REDCアルゴリズムの流れ

ヒント（段階的開示）

ヒント1: 方向性

Pythonの pow(a, e, mod) は内部でCの最適化された剰余演算を使うので実務上はそれで十分速い。しかし「除算命令を使わずに掛け算のたびに mod を取る」という制約を守るとどうなるか——という視点に立つと、除算の代わりにシフト演算とビットマスクだけでmod演算を実現する仕組みが必要になる。ハードウェアレベルで除算が乗算より何倍も遅い環境（組み込み・暗号処理系）では、この置き換えが実際に高速化に直結する。

ヒント2: アプローチ

鍵となる発想は「MODで割る代わりに、$2^k$（ビット演算で扱える数）で割る」ことである。$R=2^k$（$R>\mathrm{MOD}$、$\mathrm{MOD}$が奇数なので $\gcd(R,\mathrm{MOD})=1$ は自動成立）を選び、すべての数を「モンゴメリ表現」$\tilde{x}=xR\bmod \mathrm{MOD}$ に変換してから演算する。この表現同士の掛け算を $R$ で割って $\mathrm{MOD}$ で mod する処理（REDC）が、$R$ が2の冪であることを利用してシフトとANDだけで書ける。

ヒント3: 誘導（コード骨格）

R_BITS = MOD.bit_length() + 2
R = 1 << R_BITS
R_MASK = R - 1
mod_inv_neg = (-pow(MOD, -1, R)) % R  # MOD * mod_inv_neg ≡ -1 (mod R)

def redc(t):
    m = ((t & R_MASK) * mod_inv_neg) & R_MASK
    t2 = (t + m * MOD) >> R_BITS
    return t2 - MOD if t2 >= MOD else t2

模範解答 (Python)

import sys

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

    R_BITS = MOD.bit_length() + 2
    R = 1 << R_BITS
    R_MASK = R - 1
    mod_inv_neg = (-pow(MOD, -1, R)) % R  # MOD * mod_inv_neg ≡ -1 (mod R)

    def redc(t):
        m = ((t & R_MASK) * mod_inv_neg) & R_MASK
        t2 = (t + m * MOD) >> R_BITS
        return t2 - MOD if t2 >= MOD else t2

    def to_mont(x):
        return (x * R) % MOD

    def mont_mul(a, b):
        return redc(a * b)

    out = []
    for _ in range(Q):
        a = int(data[idx]); idx += 1
        e = int(data[idx]); idx += 1
        a_m = to_mont(a)
        res_m = to_mont(1)
        while e > 0:
            if e & 1:
                res_m = mont_mul(res_m, a_m)
            a_m = mont_mul(a_m, a_m)
            e >>= 1
        out.append(str(redc(res_m)))
    print('\n'.join(out))

solve()

計算量: モンゴメリ乗算1回あたりO(1)（多倍長乗算のコストのみ）。繰り返し二乗法のループはO(log e)回の乗算なので1クエリあたり O(log e)、全体で O(Q log(max e))。奇数MOD 3≤MOD<10^9 の範囲で3000通りの (MOD,a,e) をランダム生成し、Pythonの組み込み pow(a,e,MOD) と本実装の出力が完全一致することを確認済み。

Step-by-Step 解説

1なぜ「mod を取らない mod 演算」が要るのか
a×bのあとmodを取る処理は除算命令（DIV）を必要とし乗算命令（MUL）より数倍〜数十倍遅いことが多い。繰り返し二乗法では掛け算のたびにmodを取るためこの差が積み重なる。

2モンゴメリ表現への変換
R=2^R_BITS（MODより大きい2の冪）を選び、数xを $\tilde{x}=xR\bmod\mathrm{MOD}$ という表現に変換する。変換自体には通常の%を使ってよい（最初と最後の1回だけ）。

3REDC（モンゴメリ簡約）の仕組み
$T=\tilde a\tilde b$ に対し $m=(T\bmod R)\cdot\mathrm{MOD}^{-1}\bmod R$ を計算し、$t=(T+m\cdot\mathrm{MOD})/R$ を求める（必ず割り切れるので右シフトでよい）。Rが2の冪なので除算命令は一切使わない。

4繰り返し二乗法との組み合わせ
to_montでaをモンゴメリ表現に変換し、ループ内の掛け算をすべてmont_mulに置き換える。ループを抜けたら最後にredcを1回呼んで通常表現に戻す。

5正しさの検証
3000通りの(MOD,a,e)をランダム生成しPython組み込みpowと突き合わせて全一致を確認済み。

よくあるミス

ミス	原因	正しい書き方
モンゴメリ表現のまま出力してしまう	ループの最後にredcを呼び忘れる	res_mはredc(res_m)を通してから出力する
RがMOD以下になってしまう	bit_length()だけでR_BITSを決め余裕を持たせない	R>MODを確実にするため+2などの余裕を持たせる
MODが偶数の入力にも使ってしまう	gcd(R,MOD)=1が前提であることを忘れる	MODが奇数であることを確認する（本問は制約で保証）
mod_inv_negの符号を逆にしてしまう	「MOD·MOD^-1≡1」で計算しredcの加算をそのまま使う	REDCの定義に合わせ「≡-1」になるよう符号反転した値を使う

次のステップ

発展: RSA等の暗号処理ではモンゴメリ乗算はハードウェア実装の標準技法。バレット簡約（Barrett Reduction）との使い分けも調べてみるとよい。
次回予告: 橋・関節点列挙（Low-Link法 / Tarjanのアルゴリズム）

自己評価

完全理解おおむね理解要復習

Day 122-Q1 — モンゴメリ乗算（Montgomery Multiplication）

問題