バイナリ法（繰り返し二乗法）による高速べき乗剰余計算のPython実装

繰り返し二乗法とは：まず答えとPythonの使い方

繰り返し二乗法は、指数を2進数に分解し、底の2乗を繰り返して必要な項だけを掛ける方法です。 指数 \(k > 0\) に対して乗算回数を \(O(\log k)\) にできます。剰余が不要な通常のべき乗にも使えますが、本記事では巨大な整数を扱いやすい \(a^k \bmod m\) を中心に説明します。

Pythonで結果だけを求めるなら、 組み込み関数 pow の3引数形式を使えます。^ はべき乗ではなくビット単位の排他的論理和なので、べき乗には ** または pow を使います。

print(pow(5, 21, 99))  # 71: 5**21 を99で割った余り
print(pow(5, 0, 1))    # 0: 指数0でも 1 % 1 = 0

以下では、計算例から実装、計算量の証明へ順に進みます。暗号への応用とタイミングの注意点は後半で扱います。

具体例

例えば、\(5^{21} \pmod{p}\) を計算する場合を考えます。

まず、指数 \(21\) を2進数で表現します。 \(21 = 16 + 4 + 1 = 1 \cdot 2^4 + 0 \cdot 2^3 + 1 \cdot 2^2 + 0 \cdot 2^1 + 1 \cdot 2^0\) したがって、\(21\) の2進数表記は \((10101)_2\) です。

これを利用して、\(5^{21}\) を以下のように分解します。

\(5^1, 5^2, 5^4, 5^8, 5^{16}, \dots\) といった、\(a^{2^i}\) の形のべき乗を順次計算（前の結果を2乗していくだけで）しておけば、あとは2進数表記でビットが立っている（1である）項だけを掛け合わせることで、最終的な結果を得ることができます。これにより、掛け算の回数を大幅に削減できます。

法を99に固定すると、途中の剰余も小さい整数だけで追えます。

したがって \((49 \times 31 \times 5) \bmod 99 = 71\) です。次の素直な実装では最後の不要な2乗も実行するため、5回の2乗と3回の掛け込みを行います。最後の2乗を省略すれば、上の表を作る2乗は4回で済みます。

アルゴリズムの実装

以下の自作3実装は、整数の底・非負整数の指数（exp >= 0）・正整数の法（mod > 0）を対象とします。負の指数と0以下の法には ValueError を送出します。

繰り返し二乗法には、指数の2進数表現を下位ビットから読む Right-to-Left 方式と、上位ビットから読む Left-to-Right 方式の2通りの実装があります。どちらも計算量は同じ \(O(\log k)\) ですが、その理由を天下り的に述べるのではなく、指数の2進展開から数学的に導出します。

なぜ乗算回数が \(O(\log k)\) で済むのか

指数 \(k\) を2進数で展開すると、ビット長を \(L = \lfloor \log_2 k \rfloor + 1\) として

と書けます。指数法則 \(a^{x+y} = a^x a^y\) より、

が成り立ちます。右辺に登場する \(a^{2^0}, a^{2^1}, \dots, a^{2^{L-1}}\) は、\(a^{2^0} = a\) から出発して直前の値を2乗するだけで

と、\(L-1\) 回の2乗算だけで全て求まります。あとはこのうちビット \(b_i = 1\) に対応する項だけを掛け合わせれば \(a^k\) が得られるので、必要な乗算はたかだか「\(L-1\) 回の2乗算＋\(L\) 回以下の掛け算」、合計 \(O(L) = O(\log k)\) 回です。この見積もりはビットをどちらの向きに走査しても変わりません。「\(2^i\) 乗を順に用意し、ビットが立っている位置だけ掛け合わせる」という構造自体が走査方向に依存しないからです。以下では、この構造を実際に「下位から」「上位から」の2通りのアルゴリズムとして具体化し、それぞれのループ不変量を示して正当性を証明します。

Right-to-Left方式（下位ビットから）

Right-to-Left方式は、指数の最下位ビットから順に読みながら、ビットが立っていればresにbase（そのときの \(a^{2^i}\) ）を掛け込み、その後baseを2乗しexpを1ビット右シフトします。

ループ不変量: whileループの第 \(i\) 回目の反復が始まる直前（\(i = 0, 1, 2, \dots\) ）、次が成り立ちます。

証明（数学的帰納法）: \(i=0\) では res\(=1=a^0\) 、base\(=a=a^{2^0}\) 、exp\(=k\) なので自明に成立します。\(i\) で成立すると仮定すると、反復では exp の最下位ビット \(b_i = \texttt{exp} \bmod 2\) を見て、\(b_i=1\) のときだけ res *= base とするので、更新後の res は

となり（2進展開の定義から \(k \bmod 2^{i+1} = (k \bmod 2^i) + b_i 2^i\) ）、続いて base を2乗し（\(a^{2^i} \to a^{2^{i+1}}\) ）、exp を右シフトする（\(\lfloor k/2^i \rfloor \to \lfloor k/2^{i+1} \rfloor\) ）ことで、\(i+1\) でも不変量が成立します。ループは exp \(=0\) 、すなわち \(i=L\) で終了し、\(k < 2^L\) より \(k \bmod 2^L = k\) なので res\(=a^k\) が得られます。∎

Pythonプログラム

def power(base, exp, mod):
    """
    繰り返し二乗法を用いて、(base^exp) % mod を効率的に計算する。

    :param base: 底
    :param exp: 指数
    :param mod: 法
    :return: 計算結果
    """
    if exp < 0 or mod <= 0:
        raise ValueError("exp must be nonnegative and mod must be positive")
    res = 1 % mod
    base %= mod
    while exp > 0:
        # 指数の最下位ビットが1の場合、結果に乗算する
        if exp % 2 == 1:
            res = (res * base) % mod

        # 底を2乗し、指数を右に1ビットシフト（2で割る）
        base = (base * base) % mod
        exp //= 2

    return res

# --- 使用例 ---
# 5^21 mod 99 を計算
k = 21
g = 5
p = 99
result = power(g, k, p)
print(f"{g}^{k} mod {p} = {result}") # -> 5^21 mod 99 = 71

# 巨大な数の計算例
k = 12345678901234567890
g = 987654321987654321
p = 1000000007
result = power(g, k, p)
print(f"巨大な数の計算結果: {result}")

プログラムの解説

1. res = 1 % mod：乗算の単位元1を法で正規化して初期化します。法が1なら0になります。
2. while exp > 0:：指数 exp が0になるまでループします。
3. if exp % 2 == 1:：指数 exp の最下位ビット（2で割った余り）が1であるかをチェックします。
4. res = (res * base) % mod：最下位ビットが1であれば、現在の base の値を結果 res に掛け合わせます。
5. base = (base * base) % mod：base を2乗します。これにより、\(a, a^2, a^4, a^8, \dots\) といった形で計算が進んでいきます。
6. exp //= 2：指数 exp を2で割り（整数除算）、1ビット右にシフトします。これにより、次のループでは1つ上のビットを評価できるようになります。
7. ループが終了すると、res に最終的な計算結果が格納されています。

各ステップで剰余を取った後の値は \([0,m)\) に収まります。ただし乗算直後・剰余を取る前の積は \(m\) 以上になりえます。Pythonの整数は任意精度ですが、固定長整数の言語では、この積のオーバーフロー対策が別途必要です。

Left-to-Right方式（上位ビットから）

Left-to-Right方式は、指数の最上位ビットから順に読みながら、毎回resを2乗し、ビットが立っていればbaseを掛け込みます。これは2進数を10進整数に変換するときの**ホーナー法（Horner’s method）**と全く同じ構造です。

上位から \(j\) 個のビットを読み終えた時点でのプレフィックス値を \(p_j = \lfloor k / 2^{L-j} \rfloor\) （\(j=0,\dots,L\) 、\(p_0=0\) 、\(p_L=k\) ）とすると、2進表記の性質から漸化式

が成り立ちます。

ループ不変量: 第 \(j\) 回目の反復が始まる直前（\(j=0,1,\dots,L-1\) ）、res\(\equiv a^{p_j} \pmod m\) が成り立ちます。

証明: \(j=0\) では res\(=1=a^{p_0}=a^0\) で成立。\(j\) で成立すると仮定すると、反復ではまず res を2乗し（\(a^{p_j} \to a^{2p_j}\) ）、次にビット \(b_{L-1-j}\) が1であれば base（\(=a\) ）を掛けるので、更新後の res は

となり、\(j+1\) でも不変量が成立します。\(j=L\) でループが終わると res\(=a^{p_L}=a^k\) が得られます。∎

def power_left_to_right(base, exp, mod):
    """
    Left-to-Right方式の繰り返し二乗法。指数の最上位ビットから走査する。
    """
    if exp < 0 or mod <= 0:
        raise ValueError("exp must be nonnegative and mod must be positive")
    if exp == 0:
        return 1 % mod
    base %= mod
    res = 1
    for bit in bin(exp)[2:]:  # "0b..." を除いた、MSBから始まる2進文字列
        res = (res * res) % mod
        if bit == "1":
            res = (res * base) % mod
    return res

掲載したRight-to-Left方式・Left-to-Right方式は、\(k > 0\) のビット長を \(L\) 、ハミング重みを \(w(k) = \sum_i b_i\) とすると、どちらも2乗算を \(L\) 回、掛け込みを \(w(k)\) 回行います。合計は \(L+w(k)\) 回で、\(L+1 \le L+w(k) \le 2L\) です。理論上必要な \(L-1\) 回の2乗に対し、Right-to-Leftには最後の不要な2乗、Left-to-Rightには最初の \(1^2\) が含まれます。これらを省く最適化でも定数項が変わるだけで、いずれも \(\Theta(\log k)\) です。これは乗算の回数であり、任意精度整数の乗算1回の計算時間が一定という意味ではありません。

サイドチャネル攻撃と定数時間性：モンゴメリラダー

ここまでの2方式には共通する弱点があります。if bit == 1: 乗算する というデータ依存の分岐です。指数が公開情報（RSA暗号化の \(e\) など）であれば無害ですが、指数が秘密情報の場合——RSA復号の秘密指数 \(d\) （ RSA暗号の記事 参照）や、楕円曲線スカラー倍算の秘密スカラー（ 楕円曲線暗号の記事 の double-and-add 参照）——では、この分岐が処理時間・消費電力・電磁波などの物理的に観測可能な副作用（サイドチャネル）を通じて秘密情報を漏洩させる危険があります。

タイミング攻撃の原理

固定ビット長 \(L\) の指数に対し、Right-to-Left方式は常に \(L\) 回（最後の不要な2乗を含む）の2乗算を行いますが、掛け算の回数はハミング重み \(w(k)\) に一致します。同じ秘密指数を使い回す処理（例えばRSA秘密鍵 \(d\) による繰り返しの復号リクエスト）に対して攻撃者が実行時間を多数回測定できれば、統計的に「掛け算が多かった＝ハミング重みが大きい」ことを推定でき、さらに測定を工夫すればビットを1つずつ復元することも可能です（Paul Kocher が1996年にRSA/DH/DSSの素朴な実装に対して実証した攻撃で、以降のRSA/ECC実装は指数計算を定数時間化することが標準になりました）。本記事末尾の実行検証で、この「掛け算回数の違いが実行時間の違いとして現れる」ことを実測します。

モンゴメリラダーの導出

対策の要点は「分岐をなくす」ことではなく、「分岐の先で行う演算の個数・種類が秘密ビットの値に依存しないようにする」ことです。モンゴメリラダーは、常に2つのレジスタ \(R_0, R_1\) を保持し、どのビットでも必ず1回の乗算と1回の2乗算を行い、ビットの値によって変わるのは「どちらのレジスタにどちらの結果を書き込むか」だけ、という構成でこれを実現します。

Left-to-Right方式と同じプレフィックス \(p_j = \lfloor k/2^{L-j} \rfloor\) を使い、次のループ不変量を課します。

すなわち「\(R_1\) は常に \(R_0\) よりちょうど1乗だけ先にいる」という不変量です。次のビット \(b_{L-1-j}\) を処理して \(p_{j+1} = 2p_j + b_{L-1-j}\) に進めるとき、

• \(b_{L-1-j}=0\) なら \(p_{j+1}=2p_j\) 。必要なのは \(R_0' = a^{2p_j} = R_0^2\) と \(R_1' = a^{2p_j+1} = a^{p_j}\cdot a^{p_j+1} = R_0 \cdot R_1\) 。
• \(b_{L-1-j}=1\) なら \(p_{j+1}=2p_j+1\) 。必要なのは \(R_0' = a^{2p_j+1} = R_0 \cdot R_1\) と \(R_1' = a^{2p_j+2} = (a^{p_j+1})^2 = R_1^2\) 。

どちらの場合も「\(R_0 \cdot R_1\) という積を1回」と「どちらかの2乗を1回」——つまり演算の個数・種類は完全に同一で、違いは計算結果 \(R_0 \cdot R_1\) と (2乗結果) をどちらの変数に代入するかだけです。初期状態は \(j=0\) （\(p_0=0\) ）で \(R_0=a^0=1\) 、\(R_1=a^1=a\) とし、\(L\) ビット処理し終えた \(j=L\) （\(p_L=k\) ）で \(R_0 = a^k\) を返します。

def power_montgomery_ladder(base, exp, mod):
    """
    モンゴメリラダーによる、ビットパターンに依存しない演算回数のべき乗剰余計算。
    どのビットでも「乗算1回＋2乗算1回」を必ず実行し、代入先だけが変わる。
    """
    if exp < 0 or mod <= 0:
        raise ValueError("exp must be nonnegative and mod must be positive")
    if mod == 1:
        return 0
    base %= mod
    if exp == 0:
        return 1 % mod

    bitlen = exp.bit_length()
    R0, R1 = 1, base
    for i in range(bitlen - 1, -1, -1):
        bit = (exp >> i) & 1
        if bit == 0:
            R1 = (R0 * R1) % mod
            R0 = (R0 * R0) % mod
        else:
            R0 = (R0 * R1) % mod
            R1 = (R1 * R1) % mod
    return R0

重要な注意点として、このif bit == 0: ... else: ...というPythonレベルの分岐そのものは、CPUの分岐予測やキャッシュタイミングを通じてなお情報を漏らしうるため、真に定数時間な実装（X25519 の Curve25519 上のモンゴメリラダーなど）では、分岐の代わりにビット演算によるconditional swap（cswap）——mask = -bit のようなビットマスクを使い、R0, R1 = R0 ^ (mask & (R0 ^ R1)), R1 ^ (mask & (R0 ^ R1)) のように無分岐で入れ替える——を用います。本記事のモンゴメリラダーが解決しているのは、あくまで**「演算の回数・種類」が秘密ビットの値に依存する問題**（＝素朴な二進法べき乗算がハミング重みに応じて掛け算回数が変わる問題）であり、これは次の実行検証で実測するタイミング差の主要因です。実務でRSA/ECCの秘密鍵演算を実装する際は、cryptographyライブラリやOpenSSLなど、cswapレベルまで定数時間化された実装を使うべきです。

エッジケースと実務上の注意

• 指数が0の場合: 戻り値は 1 % mod です。Right-to-Left方式は初期値を 1 % mod にすることで、ループが0回でも正しい結果を返します。Left-to-Right方式とモンゴメリラダーは早期returnで処理しています。pow(0, 0, m) も正整数の法に対して 1 % m を返します。

• 負の指数: \(a^{-n} \bmod m\) を計算するには、\(a\) の \(m\) を法とする乗法逆元 \(a^{-1}\) （\(a \cdot a^{-1} \equiv 1 \pmod m\) 、\(\gcd(a,m)=1\) のときのみ存在）を拡張ユークリッド互除法で求め、\(a^{-n} \equiv (a^{-1})^n \pmod m\) として通常の繰り返し二乗法に帰着させます。これは RSA暗号の記事のモジュラ逆元の節 で説明している考え方そのものです。Pythonでは3引数powが負の指数を拡張しており、pow(3, -1, 11)は拡張ユークリッド互除法により4（\(3 \times 4 = 12 \equiv 1 \pmod{11}\) ）を返すことを確認しました。一方 \(\gcd(a,m) \ne 1\) の場合は逆元が存在せず、例えばpow(2, -1, 4)はValueError: base is not invertible for the given modulusを送出することも確認済みです。

• 法が1の場合: 指数0を含め、結果は常に0です。Right-to-Left方式は res = 1 % mod、Left-to-Right方式は指数0の早期returnと各ステップの剰余、モンゴメリラダーは mod == 1 の早期returnで処理します。

• 整数の大きさと計算量: 正規化後の値は \([0,m)\) 、乗算直後の積は \(m^2\) 未満なので、法が \(B\) ビットなら積は最大約 \(2B\) ビットです。\(a^k\) 全体を先に作る方法より必要メモリを抑えられます。ただしPythonの整数演算コストは桁数に依存するため、\(O(\log k)\) は乗算・剰余の回数についての評価です。

実行検証：正当性とタイミング測定

正当性の検証

Right-to-Left方式・Left-to-Right方式・モンゴメリラダーの3実装が、Python組み込みのpow(base, exp, mod)と常に同じ結果を返すことを検証しました。ビット長 8〜512 の乱数ケース500件に加え、\(k=0\) 、\(0^0\) 、\(a=0\) 、\(k=1\) 、\(m=1\) 、20桁の巨大な指数を含むケースなどのエッジケース14件を合わせた514件のテストケース全てで、4つの実装（pow・Right-to-Left・Left-to-Right・モンゴメリラダー）の出力が完全に一致し、不一致は0件でした（負の指数・モジュラ逆元の挙動は前節で個別に確認済みです）。

import random

random.seed(42)
cases = []
for _ in range(500):
    bits = random.choice([8, 16, 64, 256, 512])
    mod = random.randrange(1, 2**bits)
    base = random.randrange(0, 2**bits)
    exp = random.randrange(0, 2**bits)
    cases.append((base, exp, mod))
cases += [(5, 0, 99), (0, 0, 99), (0, 5, 99), (5, 1, 99), (5, 21, 1),
          (5, 0, 1), (0, 0, 1), (-5, 21, 99), (-5, 0, 99),
          (5, 21, 99), (5, 65537, 99),
          (12345678901234567890, 1, 1000000007), (2, 3, 5),
          (987654321987654321, 12345678901234567890, 1000000007)]

mismatches = [
    (b, e, m) for b, e, m in cases
    if not (pow(b, e, m) == power(b, e, m)
            == power_left_to_right(b, e, m)
            == power_montgomery_ladder(b, e, m))
]
print(len(cases), len(mismatches))  # -> 514 0

タイミング測定：ビットパターンへの依存性

Right-to-Left方式（通常の二進法べき乗算）とモンゴメリラダーで、指数のビットパターンによって実行時間がどう変わるかを実測しました。RSA-2048相当の2048ビットの法のもとで、いずれもビット長は2048で揃えつつハミング重みだけが異なる3種類の指数を用意しています。

• all_ones：全ビットが1（ハミング重み2048、\(2^{2048}-1\) ）
• sparse：最上位ビットと最下位ビットのみ1（ハミング重み2）
• alternating：1と0を交互（ハミング重み1024）

各条件につき150回試行し、さらに熱ドリフトやOSジッタが特定の条件だけに偏らないよう、3条件をラウンドロビンでインターリーブして計測しました（全条件をまとめて先に測ると、後半になるほどCPUスロットリング等で系統的に遅くなるバイアスが生じるため）。

条件間の最大・最小の差を「(最大中央値－最小中央値)／最小中央値」で見ると、Right-to-Left方式は121.6%（sparseの27.22msに対しall_onesは60.31msと、ハミング重みが1024多いだけで実行時間が2倍以上）もの差が生じたのに対し、モンゴメリラダーはわずか0.1%（3条件とも59.7〜59.8ms付近）とほぼ完全にビットパターンから独立していました。これは前節の理論（Right-to-Left方式の乗算回数がハミング重みに比例して増える一方、モンゴメリラダーは常に固定回数の演算を行う）と定量的に整合する結果であり、「ハミング重みが大きい秘密指数ほど処理時間が長くなる」というタイミング攻撃の前提が、実測でも明確に確認できたことになります。

import time, random, statistics

def make_exponents(bitlen):
    all_ones = (1 << bitlen) - 1
    sparse = (1 << (bitlen - 1)) | 1
    alternating = int("10" * (bitlen // 2), 2) | (1 << (bitlen - 1))
    return {"all_ones": all_ones, "sparse": sparse, "alternating": alternating}

def time_interleaved(func, base, exps, mod, repeats):
    samples = {label: [] for label in exps}
    for _ in range(repeats):
        for label, exp in exps.items():          # ラウンドロビン
            t0 = time.perf_counter()
            func(base, exp, mod)
            samples[label].append(time.perf_counter() - t0)
    return {l: (statistics.median(v), statistics.stdev(v)) for l, v in samples.items()}

random.seed(7)
bitlen = 2048
mod = random.getrandbits(bitlen) | 1 | (1 << (bitlen - 1))
base = random.getrandbits(bitlen) % mod
exps = make_exponents(bitlen)

r2l = time_interleaved(power, base, exps, mod, repeats=150)
mont = time_interleaved(power_montgomery_ladder, base, exps, mod, repeats=150)

RSA暗号における応用

繰り返し二乗法の最も重要な応用の一つが RSA暗号 です。RSA暗号では、暗号化・復号の両方で大きな整数のべき乗剰余演算が必要になります。

秘密指数 \(d\) は通常大きな整数であり、高速なべき乗剰余計算が必要です。公開指数 \(e\) には、比較的小さい65,537がよく使われます。特に復号で使う秘密指数 \(d\) は文字通り秘密情報であり、前述のサイドチャネル攻撃の標的になります。実運用のRSA実装（OpenSSL、cryptographyライブラリなど）は、単純な二進法べき乗算ではなく、モンゴメリラダーやモンゴメリ乗算（Montgomery multiplication、剰余計算自体を除算なしで高速化する別の技法）と定数時間の条件選択を組み合わせて、\(d\) のビットパターンが実行時間に漏れないようにしています。

計算量の比較

Pythonの組み込み関数pow(base, exp, mod)は内部で繰り返し二乗法を使用しており、RSA暗号の実装で直接利用できます。ただしCPython のpow実装自体は必ずしも全ケースで定数時間性を保証するものではなく、秘密鍵を扱う本番システムではcryptographyのようにOpenSSLへ委譲するライブラリを使うのが安全です。

関連記事

• RSA暗号の原理とPythonによる実装 - べき乗剰余計算がRSA暗号の暗号化・復号の核心的な処理として使われていることを解説しています。
• 楕円曲線ElGamal暗号の原理とPythonによる簡易実装 - べき乗剰余計算を楕円曲線上のスカラ倍算に応用した公開鍵暗号を解説しています。
• Diffie-Hellman鍵交換プロトコル：理論とPython実装 - べき乗剰余計算を基盤とするDiffie-Hellman鍵交換の原理と実装を解説しています。
• セキュリティ資格比較：CISSP vs 情報処理安全確保支援士 - RSA暗号を含む暗号技術が試験範囲に含まれるセキュリティ資格の比較を解説しています。
• 暗号ロードマップ：対称鍵・公開鍵・ハッシュ・鍵交換の体系ハブ — 繰り返し二乗法を中核計算とする RSA・DH・ElGamal を一枚地図で俯瞰する暗号メタハブ。
• RSA 暗号：理論と Python 実装 — 繰り返し二乗法の代表的応用先。\(c = m^e \bmod n\) の高速計算が本アルゴリズム抜きには不可能であることを示します。
• Diffie-Hellman 鍵交換の理論と Python 実装 — \(g^a \bmod p\) の高速計算が DH の実用化を支える点を解説。X25519 の定数時間モンゴメリラダーにも言及しています。
• 楕円曲線暗号（ECC）の数学と Python 実装 — 繰り返し二乗法の「2倍して足す」構造が、楕円曲線上のスカラー倍算（double-and-add）にそのまま対応することを詳解しています。

関連ツール

• 進数変換ツール（CalcBox） - 2進数・8進数・10進数・16進数の相互変換
• 進数変換ツール（DevToolBox） - 開発者向け進数変換ツール

参考資料

• Python公式ドキュメント：pow — 3引数形式と負の指数・逆元の条件。