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

繰り返し二乗法（バイナリ法）による高速べき乗剰余計算の理論と Python 実装。O(log n) の計算量解析、Right-to-Left 方式の power(base, exp, mod) 自作実装、組み込み関数 pow(base, exp, mod) との使い分け、RSA 暗号のべき乗剰余計算や楕円曲線スカラー倍算（double-and-add）への応用まで解説します。

17 min read

\(a^k \pmod{p}\) のような、大きな数のべき乗を法 \(p\) で計算する場合、単純に \(a\) を \(k\) 回掛けてから剰余を求めると、計算途中の数値が非常に大きくなり、計算時間やメモリの点で非効率です。

繰り返し二乗法 (Exponentiation by Squaring)、またはバイナリ法 (Binary Method) は、この計算を効率的に行うためのアルゴリズムです。RSA暗号などで巨大な数のべき乗剰余を計算する際に不可欠な技術です。

このアルゴリズムの基本は、指数 \(k\) を2進数表現で捉え、計算量を \(O(\log k)\) に削減することにあります。

具体例

例えば、\(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^{21} = 5^{16+4+1} = 5^{16} \cdot 5^4 \cdot 5^1 \]

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

アルゴリズムの実装

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

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

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

\[ k = \sum_{i=0}^{L-1} b_i 2^i, \qquad b_i \in \{0, 1\} \]

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

\[ a^k = a^{\sum_{i=0}^{L-1} b_i 2^i} = \prod_{i=0}^{L-1} \left(a^{2^i}\right)^{b_i} = \prod_{i \,:\, b_i = 1} a^{2^i} \]

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

\[ a^{2^{i+1}} = \left(a^{2^i}\right)^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\) ）、次が成り立ちます。

\[ \texttt{res} \equiv a^{\,k \bmod 2^i} \pmod{m}, \qquad \texttt{base} \equiv a^{2^i} \pmod{m}, \qquad \texttt{exp} = \left\lfloor \frac{k}{2^i} \right\rfloor \]

証明（数学的帰納法）: \(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 は

\[ a^{k \bmod 2^i} \cdot a^{2^i \cdot b_i} = a^{(k \bmod 2^i) + b_i 2^i} = a^{k \bmod 2^{i+1}} \]

となり（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: 計算結果
    """
    res = 1
    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 = 20

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

プログラムの解説

  1. res = 1：結果を格納する変数を1で初期化します。
  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 に最終的な計算結果が格納されています。

この方法では、計算の各ステップで剰余を取るため、扱う数値が法 mod より大きくなることがなく、巨大な数のべき乗計算でもオーバーフローを防ぎつつ高速に処理できます（詳細は後述の「エッジケースと実務上の注意」を参照）。

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

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

\[ k = b_{L-1} 2^{L-1} + b_{L-2} 2^{L-2} + \cdots + b_0 = \Big( \big( (b_{L-1}) \cdot 2 + b_{L-2} \big) \cdot 2 + b_{L-3} \Big) \cdot 2 + \cdots + b_0 \]

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

\[ p_{j+1} = 2 p_j + b_{L-1-j} \]

が成り立ちます。

ループ不変量: 第 \(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 は

\[ a^{2p_j} \cdot a^{b_{L-1-j}} = a^{2p_j + b_{L-1-j}} = a^{p_{j+1}} \]

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

def power_left_to_right(base, exp, mod):
    """
    Left-to-Right方式の繰り返し二乗法。指数の最上位ビットから走査する。
    """
    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方式のいずれも、2乗算は必ず \(L-1\) 回、掛け算はビットが立っている個数（ハミング重み \(w(k) = \sum_i b_i\) ）だけ発生します。つまり乗算回数は \(L-1 \le (\text{乗算回数}) \le 2L-1\) の範囲に収まり、いずれの走査方向でも \(\Theta(\log k)\) であることが確認できます。実務上は、Right-to-Left方式はbaseの2乗とexpの右シフトが独立に進むためループ本体が単純になり、Left-to-Right方式は指数を文字列やビット列として扱う実装（後述のモンゴメリラダー含む）と相性が良い、という使い分けがあります。

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

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

タイミング攻撃の原理

固定ビット長 \(L\) の指数に対し、Right-to-Left方式は常に \(L-1\) 回の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_0 \equiv a^{p_j} \pmod m, \qquad R_1 \equiv a^{p_j + 1} \pmod m \]

すなわち「\(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 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の場合: \(a^0 = 1\) です（\(a=0\) でも慣習的に \(0^0=1\) と定義することが多く、実際 Python の pow(0, 0, m) も 1 % m を返します）。Right-to-Left方式は exp=0 だとループが1度も実行されず res=1 を返すので自然に正しく、Left-to-Right方式・モンゴメリラダーの実装ではこのケースを明示的な early return として扱っています（後者は bin(0) や bit_length() がビット列を持たないため必須です）。

  • 負の指数: \(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の場合: 任意の \(a^k \bmod 1\) は常に \(0\) です。モンゴメリラダーの実装では mod == 1 を明示的に特別扱いしていますが、これがなくても base %= mod で base が常に0になるため結果的に0が返ります。境界条件として明示しておくのが安全です。

  • オーバーフロー回避: 素朴に \(a^k\) を先に計算してから剰余を取る方式では、\(a^k\) 自体の桁数がおよそ \(k \log_{10} a\) 桁にもなり、RSA-2048のように \(a, k\) が数百〜2048ビットの場合は非現実的な時間・メモリを要します。本記事のコードのように各ステップで必ず法 \(m\) を取ることで、扱う数値は常に \([0, m)\) に収まり、法のビット長が固定されたまま \(O(\log k)\) 回の「固定長同士の掛け算＋剰余」で計算が完了します。RSA・DHなどの実装で剰余を計算し忘れると、正しさは変わらなくても実行時間・メモリが桁違いに悪化する点に注意してください。

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

正当性の検証

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桁の巨大な指数を含むケースなどのエッジケース8件を合わせた508件のテストケース全てで、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),
          (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))  # -> 508 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-Leftall_ones（2048）60.31 ms7.21 ms
Right-to-Leftalternating（1024）43.73 ms5.95 ms
Right-to-Leftsparse（2）27.22 ms3.31 ms
モンゴメリラダーall_ones（2048）59.74 ms4.27 ms
モンゴメリラダーalternating（1024）59.79 ms4.92 ms
モンゴメリラダーsparse（2）59.76 ms3.00 ms

条件間の最大・最小の差を「(最大中央値－最小中央値)／最小中央値」で見ると、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暗号では、暗号化・復号の両方で大きな整数のべき乗剰余演算が必要になります。

\[c = m^e \bmod n \quad \text{（暗号化）}\]

\[m = c^d \bmod n \quad \text{（復号）}\]

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

計算量の比較

方法計算量\(e = 2^{16} + 1\) での乗算回数秘密指数に対する安全性
素朴なべき乗\(O(e)\)65,537回（非現実的な速度のため対象外）
繰り返し二乗法（分岐あり）\(O(\log e)\)17回ハミング重みが実行時間に漏れる
モンゴメリラダー\(O(\log e)\)固定 \(2\lceil \log_2 e\rceil\) 回演算回数が指数のビットパターンに依存しない

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

関連記事

関連ツール

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション DevToolBox Regex Tester 正規表現をリアルタイムテスト

関連記事