Miller–Rabin（ミラー・ラビン）テストの性能比較

1. はじめに

RSA暗号はランダムに選んだ2つの大きな素数（$=p, q$）をもとに、公開鍵と秘密鍵を生成する暗号方式である。

ランダムに生成した数が素数であるのかどうかを効率的に判定するためにMiller–Rabin（ミラー・ラビン）素数判定法と呼ばれる確率的アルゴリズムが広く使われている。

本記事では、Miller-Rabin テストについて解説し、ほかの素数判定アルゴリズムと性能を比較する。

2. Miller-Rabin テスト

Miller-Rabin テストは、次の性質を利用して素数であるかどうかを確率的に判定するアルゴリズムである。

はじめに、判定したい奇数 $n$ に対して次のように分解する。ここで $t$ は奇数。

$$
n - 1 = 2 ^ s \times t
$$

次に、$1 \le a \le n - 1$ の中から適当な基数 $a$ を選び、

$$
a^t, a^{2t}, a^{4t}, ..., a^{(s-1)t}
$$

を順に $n$ で割ったあまりを調べる。

もしこの中で、

• 最初に $n - 1$ が出て
• そのあとに $1$ が現れる

という並びになっていれば、$n$ は素数の可能性が高い。それ以外は、確実に合成数である。

Miller-Rabin テストは、合成数 $n$ に対してランダムに選んだ基数 $a$ が素数と誤って判定される確率（=誤判定）を起こす確率は高々 $\dfrac{1}{4}$ なので、テストを複数回繰り返すことで、誤判定の確率を指数関数的に小さくすることができる。実用上は50回程度のテスト回数で、ほぼ確実に素数とみなして問題はない。

3. Miller-Rabin の実装

ここでは、Python を用いて Miller-Rabin テストを実装する。

コードを以下に示す。

import random

def miller_rabin(n: int, k: int = 50) -> bool:
    """確率的素数判定（Miller–Rabin法）
    
    Args:
        n (int): 判定対象の整数
        k (int): テスト回数（大きいほど精度が高い）

    Returns:
        bool: n が素数と判定された場合 True、合成数の場合 False
    """
    if n in (2, 3):
        return True
    if n < 2 or n % 2 == 0:
        return False

    # n−1 = 2^s * t  に分解
    t, s = n - 1, 0
    while t % 2 == 0:
        t //= 2
        s += 1

    # k回テストを繰り返す
    for _ in range(k):
        a = random.randrange(2, n - 1)
        x = pow(a, t, n)  # a^t mod n
        flag = False

        if x == 1 or x == n - 1:
            continue  # 素数っぽい

        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                flag = True
                break

        if not flag:
            # 最後まで -1 (= n - 1) が出なかった
            return False
            
    return True

コードの流れ

1. $n$ が 2 または 3 のときは素数なので True を返す
2. $n$ が 2 未満、または偶数のときは素数ではないため False を返す
3. $n - 1 = 2^s \times t$ の形に分解する（$t$ が奇数になるまで 2 で割る）
4. $k$ 回のテストを繰り返す
5. 基数 $a$ をランダムに選び、最初の値 $a^t \bmod n$ を計算する
6. もし $a^t \equiv 1 \pmod{n}$ または $a^t \equiv n - 1 \pmod{n}$ なら、 そのラウンドでは「素数の可能性が高い」として次のテストに進む
7. そうでない場合は、$s - 1$ 回だけ平方を繰り返し、$x = a^{2^i t} \bmod n$ を順に計算する
8. 途中で $n - 1$ が現れた場合は、そのラウンドを合格とみなし、次のテストに進む
9. どの段階でも $n - 1$ が出なければ、その時点で $n$ は確実に合成数と判断し、False を返す
10. すべてのテスト（= $k$ 回のラウンド）を通過した場合、$n$ は高い確率で素数であるとみなし、True を返す

4. その他の素数判定法

ここまで Miller-Rabin テストを解説してきた。

素数を判定する手法はこれだけではなく、より単純な試し割り法や、フェルマーテストなどが存在する。

ここではそれらのアルゴリズムをいくつか取り上げ、Pythonでの実装を行う。

4.1. 試し割り法

試し割り法は、2 から $\sqrt{n}$ までの整数で $n$ を割り、どれでも割り切れなければ $n$ を素数とみなす。

理論的には単純であるが、計算量は $O(\sqrt{n})$ であり、大きな数になると非常に時間がかかる。

なぜ計算量が $O(\sqrt{n})$ なのか??

$n$ を合成数と仮定すると、ある整数 $a, b$（$a \ge b$）が存在して $n = a \times b$ と書ける。

ここで、$b > \sqrt{n}$ であると仮定すると、

$$
a > b > \sqrt{n}
$$

であるから、

$$
a \times b > \sqrt{n} \times \sqrt{n} = n
$$

となって矛盾する。

したがって、少なくとも一方の因数は $\sqrt{n}$ 以下である。

これは、$\sqrt{n}$ までの範囲で割り算を試せば、
もし $n$ が合成数であれば必ずどこかで割り切れることを意味する。

よってループを最大で $\sqrt{n}$ 回まで回せば十分であり、これが試し割り法の時間計算量が $O(\sqrt{n})$ となる理由である。

コードを以下に示す。

import math

def is_prime(n):
    if n < 2:
        return False
        
    if n == 2:
        return True

    for i in range(3, int(math.sqrt(n)) + 1):
        if n % i == 0:
            return False

    return True

コードの流れ（試し割り法）

1. $n$ が 2 以下であるとき、素数ではないので False を返す
2. $n$ が 2 の時は True を返す
3. ループを 3 から $\sqrt{n}$ まで回す
4. $n$ を $i$ で割ったときの剰余が 0 のとき、これは合成数であるので False を返す
5. いずれの数（$= i$）でも割り切れなかったとき、素数であるので True を返す

4.2. 試し割り法の改良

4.1. では試し割り法について説明をした。ここでは、その計算量を少しでも減らすための改良版を紹介する。

1. 偶数時の早期リターン
   $n$ が 2 で割り切れるとき（$=$ 偶数の時）、2 以外は必ず合成数になる。したがって、偶数を先に判定してしまえば、以降は奇数を調べるだけでよくなる。
2. 素数は $6k \pm 1$ の形で表される
   3 以上の素数はすべて、$6k \pm 1$ の形で表されるので、この形だけを調べれば十分である。

なぜ素数は $6k \pm 1$ の形であらわされるのか??

ある整数 $n$ を 6 で割ったときのあまりで分類すると、$k$ を整数として次のようにあらわされる。

\begin{align}
n &= 6k, \\
n &= 6k + 1, \\
n &= 6k + 2, \\
n &= 6k + 3, \\
n &= 6k + 4, \\
n &= 6k + 5
\end{align}

このうち、$6k, \ 6k+2, \ 6k+3, \ 6K+4$ は合成数である。それぞれ、6、2、3、2 で割り切れる。

\begin{align}
n &= 6k = 6 \times k, \\
n &= 6k + 2 = 2(3k + 1), \\
n &= 6k + 3 = 3(2k + 1), \\
n &= 6k + 4 = 2(3k + 2), \\
\end{align}

残るのは $6k + 1$ と $6k + 5$ ($= 6 \pm 1$) のみである。

つまり、3以上の全ての素数は $6k \pm 1$ の形で表される。

コードを以下に示す。

def is_prime_6k(n):
    if n < 2:
        return False
        
    if n in (2, 3):
        return True

    if n % 2 == 0 or n % 3 == 0:
        return False

    i = 5
    while i * i <= n:
        if n % i == 0 or n % (i + 2) == 0:
            return False
        i += 6
    return True

コードの流れ（試し割り法・改良版）

1. $n$ が 2 未満のとき、素数ではないので False を返す
2. $n$ が 2 または 3 のとき、素数なので True を返す
3. $n$ が 2 または 3 で割り切れるとき、合成数なので False を返す（このような数は $6k \pm 1$ の形で表せない）
4. この時点で素数候補の最小値は 5 なので、5 からループを開始する（試し割り法と同様に、$\sqrt{n}$ まで調べる）
5. $n$ が $i$ または $i + 2$ で割り切れる場合、これは合成数なので False を返す（ここで $i = 6k - 1$, $i + 2 = 6k + 1$ の形を同時にチェックしている）
6. どの $i$ に対しても割り切れなかった場合、$n$ は素数と判定され、True を返す

4.3. フェルマーテスト

ここでは、フェルマーテストについて説明をする。この手法は、フェルマーの小定理に基づく、確率的な素数判定アルゴリズムの一種である。

フェルマーの小テストとは、素数 $p$ に対して、$p$ と互いに異なる任意の整数 $a$ について次が成り立つ。

$$
a ^ {p-1} \equiv 1 \pmod{p}
$$

もし $n$ が素数であれば、上の式が $n$ に対しても成り立つはずである。したがって、ランダムに $a$ を選び、次の条件を満たすかどうかを調べる。

$$
a^{n-1} \bmod n = 1
$$

これが成立すれば、$n$ は素数の可能性が高く、成立しなければ $n$ は確実に合成数である。

フェルマーテストは、単純で高速に判定することができる一方、カーマイケル数のような特殊な合成数を誤って素数と判定してしまうことがある。

例：561, 1105, 1729

コードは以下の通り。

import random

def fermat_test(n, k=10):
    if n in (2, 3):
        return True
        
    for _ in range(k):
        a = random.randrange(2, n - 1)
        if pow(a, n - 1, n) != 1:
            return False
            
    return True

コードの流れ（フェルマーテスト）

1. $n$ が 2 あるいは 3 のとき素数であるので True を返す
2. それ以外の整数に対して、$k$ 回テストを繰り返す
3. 各テストで、基数 $a$ をランダムに、$[2, n-2]$ の範囲から選ぶ
4. フェルマーの小定理に基づき、$a^{n-1} \bmod n$ を計算して、結果が1でなければ、$n$ は合成数なので False を返す
5. もしすべてのテストで 1 となった場合、$n$ は高い確率で素数とみなして True を返す

5. 素数判定法の比較

ここでは、4. その他の素数判定法で説明した手法と Miller-Rabin テスト都の性能の違いについて比較・検討する。

比較・検討するために、以下のテストを行う。

1. 指定されたサイズの素数を生成する
2. 結果を返すまでの時間を計測する

なお、テストの実行時間が長くなりすぎる場合は、5秒を上限として計測を打ち切る。

また、Miller-Rabin テストとフェルマーテストはそれぞれテスト回数は50回に設定してある。

注意

試し割り法が確定的素数判定法であるのに対して、フェルマーテストおよび Miller-Rabin テストは確率的素数判定法に分類される点に注意する。そのため、今回のテスト結果は厳密なものではない。

また今回のテストでは、カーマイケル数のような擬素数に対する特別な判定も行っていないため、一部のケースでは誤判定が起こる可能性がある。

コードを次に示す。

import time, random, math, signal
import numpy as np
import matplotlib.pyplot as plt
from sympy import nextprime
from functools import wraps

LIMIT = 5

class TimeoutException(Exception):
    """時間超過時に送出される例外"""
    pass


def measure_time(limit=None):
    """
    実行時間を計測し、limit秒を超えたら中断するデコレータ
    例: @measure_time(limit=2)
    """
    def decorator(func):
        def _handle_timeout(signum, frame):
            raise TimeoutException(f"{func.__name__}: timeout (> {limit}s)")

        @wraps(func)
        def wrapper(*args, **kwargs):
            signal.signal(signal.SIGALRM, _handle_timeout)
            if limit:
                signal.alarm(limit)
            start = time.time()
            try:
                result = func(*args, **kwargs)
                elapsed = time.time() - start
                print(f"{func.__name__}: {elapsed:.3e} 秒")
                return time.time() - start

            except TimeoutException as e:
                print(e)
                return None

            except Exception as e:
                print(e)
                return None

            finally:
                signal.alarm(0)  # アラーム解除
        return wrapper
    return decorator


@measure_time(LIMIT)
def MillerRabin(n: int, k: int = 50) -> bool:
    """確率的素数判定（Miller–Rabin法）
    
    Args:
        n (int): 判定対象の整数
        k (int): テスト回数（大きいほど精度が高い）

    Returns:
        bool: n が素数と判定された場合 True、合成数の場合 False
    """
    if n in (2, 3):
        return True

    if n < 2 or n % 2 == 0:
        return False

    # n−1 = 2^s * t  に分解
    t, s = n - 1, 0
    while t % 2 == 0:
        t //= 2
        s += 1

    # k回テストを繰り返す
    for _ in range(k):
        a = random.randrange(2, n - 1)
        x = pow(a, t, n)  # a^t mod n
        flag = False

        if x == 1 or x == n - 1:
            continue  # 素数っぽい

        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                flag = True
                break

        if not flag:
            # 最後まで -1 (= n - 1) が出なかった
            return False

    return True

@measure_time(LIMIT)
def is_prime(n):
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    for i in range(3, int(n**0.5) + 1, 2):
        if n % i == 0:
            return False
    return True

@measure_time(LIMIT)
def is_prime_6k(n):
    if n < 2:
        return False

    if n in (2, 3):
        return True

    if n % 2 == 0 or n % 3 == 0:
        return False

    i = 5
    while i * i <= n:
        if n % i == 0 or n % (i + 2) == 0:
            return False
        i += 6
    return True

@measure_time(LIMIT)
def fermat_test(n, k=50):
    if n <= 3:
        return n > 1
    for _ in range(k):
        a = random.randrange(2, n - 1)
        if pow(a, n - 1, n) != 1:
            return False
    return True


sizes = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]
funcs = [is_prime, is_prime_6k, MillerRabin, fermat_test]
results = {f.__name__: [] for f in funcs}

for bits in sizes:
    print(f"\n== {bits}-bit number ==")
    n = random.getrandbits(bits) | 1
    n = nextprime(n)
    for func in funcs:
        t = func(n)
        results[func.__name__].append(t)

実行結果を次に示す。

== 1-bit number ==
is_prime: 2.384e-06 秒
is_prime_6k: 1.192e-06 秒
MillerRabin: 7.153e-07 秒
fermat_test: 9.537e-07 秒

== 2-bit number ==
is_prime: 1.192e-06 秒
is_prime_6k: 7.153e-07 秒
MillerRabin: 7.153e-07 秒
fermat_test: 7.153e-07 秒

== 4-bit number ==
is_prime: 1.001e-05 秒
is_prime_6k: 1.192e-06 秒
MillerRabin: 8.154e-05 秒
fermat_test: 4.768e-05 秒

== 8-bit number ==
is_prime: 5.007e-06 秒
is_prime_6k: 1.907e-06 秒
MillerRabin: 8.941e-05 秒
fermat_test: 6.962e-05 秒

== 16-bit number ==
is_prime: 1.121e-05 秒
is_prime_6k: 9.060e-06 秒
MillerRabin: 1.640e-04 秒
fermat_test: 1.366e-04 秒

== 32-bit number ==
is_prime: 1.707e-03 秒
is_prime_6k: 1.593e-03 秒
MillerRabin: 3.250e-04 秒
fermat_test: 3.202e-04 秒

== 64-bit number ==
is_prime: timeout (> 5s)
is_prime_6k: timeout (> 5s)
MillerRabin: 9.813e-04 秒
fermat_test: 1.017e-03 秒

== 128-bit number ==
is_prime: timeout (> 5s)
is_prime_6k: timeout (> 5s)
MillerRabin: 6.738e-03 秒
fermat_test: 4.385e-03 秒

== 256-bit number ==
is_prime: timeout (> 5s)
is_prime_6k: timeout (> 5s)
MillerRabin: 9.523e-03 秒
fermat_test: 9.005e-03 秒

== 512-bit number ==
is_prime: timeout (> 5s)
is_prime_6k: timeout (> 5s)
MillerRabin: 4.312e-02 秒
fermat_test: 4.230e-02 秒

== 1024-bit number ==
is_prime: timeout (> 5s)
is_prime_6k: timeout (> 5s)
MillerRabin: 2.498e-01 秒
fermat_test: 2.564e-01 秒

== 2048-bit number ==
int too large to convert to float
is_prime_6k: timeout (> 5s)
MillerRabin: 1.566e+00 秒
fermat_test: 1.553e+00 秒

これをグラフにしたものを次に示す。1枚目は試し割り法とその改良、2枚目は Miller-Rabin テストとフェルマーテストの結果を図示したものである。

1枚目の試し割り法とその改良は、32 bit 以降計測上限の5秒に達してしまい、以降すべて timeout となっている。そのためグラフも横軸の上限は 32 bit までしかない。グラフから改良版はより高速であることがわかる。

2枚目の Miller-Rabin テストとフェルマーテストは、フェルマーテストのほうがわずかに高速であることがわかる。また、試し割りとはことなり、2048 bit でも制限時間（$=5$秒以内）に処理を終えることができている。

この結果から処理時間という観点からみると、フェルマーテストが最適な素数判定法であるように見える。

しかし、特定の合成数（$=$ カーマイケル数）を素数と誤判定してしまう可能性があり、また、一般の合成数に対しても誤判定確率はおおむね $\dfrac{1}{2}$ 以下とされている。これは、Miller–Rabin テストの誤判定確率（高々 $\dfrac{1}{4}$）よりも高い。

なお、この誤判定確率はテスト回数を重ねることで無視できるレベルにまで抑えることができる。

したがって、速度だけみればフェルマー・テストが有利だが、信頼性の面では Miller-Rabin テストが優れている。

6. まとめ

本記事では、代表的な素数判定法である 試し割り法・フェルマー・テスト・Miller–Rabin テスト の概要と性能を比較した。

手法	種類	誤判定の有無	計算量の目安	特徴
試し割り法	確定的	なし	$O(\sqrt{n})$	最も単純で確実だが、大きな数では実用的でない
フェルマー・テスト	確率的	あり（カーマイケル数で常に誤判定）	$O(k (\log n)^3)$	高速だが信頼性に欠ける
Miller–Rabin テスト	確率的	あり（$\dfrac{1}{4}$ 以下）	$O(k (\log n)^3)$	速度と信頼性のバランスが良く、実用上の標準

7. 参考資料

[1]. Pythonで学ぶ暗号理論

[2]. 中学生でも解ける｢素数はほぼ6の倍数±1｣の証明