繰り返し二乗法による高速べき乗剰余計算のPython実装

$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である）項だけを掛け合わせることで、最終的な結果を得ることができます。これにより、掛け算の回数を大幅に削減できます。

アルゴリズムの実装

繰り返し二乗法には、指数のビットを右から見るか（Right-to-Left）、左から見るか（Left-to-Right）で少し異なる実装がありますが、どちらも計算量は同じです。ここでは、よりシンプルに実装できるRight-to-Left方式を紹介します。

Pythonプログラム

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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}")

プログラムの解説

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

この方法では、計算の各ステップで剰余を取るため、扱う数値が法 mod より大きくなることがなく、巨大な数のべき乗計算でもオーバーフローを防ぎつつ高速に処理できます。