素数判定いろいろ - フェルマーテスト・ミラーラビン素数判定法

数学好きのエンジニアしてます @srtk86 です。本日の日付 1219 は素数です。

素数かどうかを判定するアルゴリズムについて色々見てみたので、その話をしようと思います。

前回は素数の定義と、そこから自然に生まれるアルゴリズムについて書きました。
今回は確率的な素数判定について、裏側の理論に触れながら追っていきたいと思います。

目次

シンプルな素数判定と、素数の分布
フェルマーテスト・ミラーラビン素数判定法 - 今ここ！
AKS素数判定法

フェルマーテスト

合同式の定義

本格的にフェルマーテストの話をする前に、あまりの計算を表す形式とその意味について、簡単に触れておきます。
「$n$を$d$でわると$a$あまる」といった性質を、以下のような式で表します（一般的に合同式と呼ばれるものです）

n \equiv a \pmod{d}

合同式の性質などに関しては以降の話でそれほど重要ではないため、割愛します。
「この記号が出てくるところでは、割り算のあまりを計算しているんだなぁ」と思っていただければ。

フェルマーの小定理

フェルマーテストは、次の「フェルマーの小定理」の対偶を利用する判定方法です。

フェルマーの小定理

$p$を素数とするとき、$p$と互いに素な整数$a$に対して以下が成り立つ。

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

---

たとえば $p = 7$ の場合は

$a = 2$ のとき $a^{6} = 64 \equiv 1 \pmod{7}$
$a = 3$ のとき $a^{6} = 729 \equiv 1 \pmod{7}$
$a = 4$ のとき $a^{6} = 4096 \equiv 1 \pmod{7}$
...

という具合に、あまりが必ず$1$になってくれます。

$p$が合成数の場合は、上記の等式が成り立たない場合があります。
たとえば $p = 4$ とすると
$a = 3$ のとき $a^{3} = 27 \equiv 3 \pmod{4}$

フェルマーテストは、この性質を利用して確実に合成数のものを弾こう、というアプローチを取るものです。

フェルマーテストのプログラム

先ほど紹介したフェルマーの小定理の対偶を書くと、こうなります。

フェルマーの小定理の対偶

$p$を整数とする。$p$と互いに素な整数$a$に対して

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

となることがあれば、$p$は合成数。

---

これを活かして、確率的に素数判定を行うのが、フェルマーテストです。

フェルマーテストの流れ

以下の試行を1回分として、何回か試行を繰り返します（判定したい数の大きさなどに応じて適当に選択します）

• $a$を、$2$から$n - 1$までの数からランダムに選ぶ
• $a$が$n$と互いに素でなければ、$n$は合成数 (この時点で$n$より小さい素因数が発見されたことになるので)
• $a^{p-1} \not\equiv 1 \pmod{p}$ なら、$n$は合成数

プログラムで表現すると、こんな感じになります。（ここでは試行回数を100回としています）

fermat_prime_test.py

import math
import random

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

    for k in range(100):
        a = random.randint(2, n - 1)

        if gcd(n, a) != 1:
            return False

        if pow(a, n - 1, n) != 1:
            return False

    return True

def gcd(a, b):
    while b > 0:
        a, b = b, a % b
    return a

if __name__ == '__main__':
    print(is_prime(2 ** 89 - 1))

gcd(a, b) は最大公約数を計算する関数を切り出していて、ユークリッドの互除法と呼ばれるアルゴリズムに基づいています。
また、$a^{n-1} \pmod{p}$ の計算は、Pythonに組み込まれている pow 関数を使っています。普通に a ** (n - 1) % n と書きたくなるところですが、これだと$n$の桁が多い場合の処理が劇的に遅くなってしまいます。
※ちなみに main 関数内でテストしている数($2^{89} - 1$)は10番目のメルセンヌ素数です。

フェルマーテストの精度

このテストは、運が悪いと合成数を誤って素数と判定します。
たとえば$341$という数は $341 = 11 \times 31$ と素因数分解される合成数ですが、
$2^{340} \equiv 1 \pmod{341}$
が成り立ってしまいます。($a = 2, p = 341$として、$p$が合成数なのに、合同式が成り立ってしまう状態)
このように、合成数でも、判定を通過してしまうような場合があるわけです。こういった数を擬素数と呼んだりします。

であれば、試行回数を増やして $a$ を何回か変化させれば確実だよね、と言いたいところなのですが、実はテストを必ずパスしてしまう合成数の存在が知られています。たとえば$561$はこの代表例で、

a^{560} \equiv 1 \pmod{561} \quad (2 \leq a \leq 560)

が成り立ちますが、$561 = 3 \times 11 \times 17$なので、$561$は合成数です。これは、何回フェルマーテストを試行しても、誤って素数と検出されてしまうような合成数が存在する、ということを意味します。やっかいですね。

この「フェルマーテスト絶対合格するマン」な合成数のことをカーマイケル数(あるいは絶対擬素数)と呼びます。ラマヌジャンのタクシー数のエピソードで(一部で)有名な$1729$も、カーマイケル数の1つです。

カーマイケル数は無限に存在することも知られていて、最初からこれらを弾くといった手法も使えないことが分かっています。本当にやっかいですね。
ですが比較的最近になって、カーマイケル数を導出する方法が中国で発見されたというニュースも飛び交ったりしています。この発見が今後どう影響してくるのか、カーマイケル数絶対不合格にしたいマンの一人として非常に楽しみです。

割合としては、1回の試行に合格すると、それが素数である確率は$\frac{1}{2}$以上であることが知られています。（群論を使って証明できるようです） 試行回数を繰り返すことで、結構な精度で判定できるため実用性はそれなりにありますが、これから紹介する判定法の方が精度として確実であるため、実際はそちらが使われることが多いです。

ミラーラビンの素数判定法

ミラーラビンの判定法も、素数が満たす合同式をもとに組まれるアルゴリズムです。まずはその性質を見てみることにします。

ミラーラビンの素数判定法を支える素数の性質

$p$を奇素数（奇数の素数 = 2以外の素数） とします。このとき、$p - 1$は偶数なので、$2$で何回か割り切ることができます。その回数を$s$とすると、以下のように表現することができます。

p - 1 = 2^{s} \times d \quad (s \in \mathbb{N}, d: 奇数)

このとき、$1 \leq a \leq n - 1$ を満たす$a$に対して

a^{d} \equiv 1 \pmod{p}

または、ある $0 \leq r \leq s - 1$ に対して

a^{2^{r} \cdot d} \equiv -1 \pmod{p}

が成り立ちます。

---

これは、$p$が素数の場合、$x^2 \equiv 1 \pmod{p}$ を満たす$x$が$1$または$-1$以外には存在しない、という性質をもとにしています。
先ほどのフェルマーの小定理は、$p$が素数のときに

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

が成り立つ、というものでしたが、上の性質とあわせて考えると、これが成立するためには下の赤枠で囲われた合同式のうちのどれかが成り立つ必要があります。

ならば実際にこれらの合同式を全て判定してしまって、いずれも成り立たなければ合成数と判定できるよね、というのがミラーラビン素数判定法の方針です。

ミラーラビンの素数判定法のプログラム

実際のプログラムは下記のようになります。(Wikipedia にもRuby版のコードが掲載されています)

miller_rabin_prime_test.py

import random

def is_prime(n):
    if n == 2: return True
    if n == 1 or n & 1 == 0: return False

    d = (n - 1) >> 1
    while d & 1 == 0:
        d >>= 1

    for k in range(100):
        a = random.randint(1, n - 1)
        t = d
        y = pow(a, t, n)

        while t != n - 1 and y != 1 and y != n - 1:
            y = (y * y) % n
            t <<= 1

        if y != n - 1 and t & 1 == 0:
            return False

    return True

if __name__ == '__main__':
    print(is_prime(2 ** 89 - 1))

試していただくとわかりますが、100桁以上の数についても、かなり高速で判定できます。マシンスペックにもよりますが、200桁程度であれば数秒以内に判定できちゃいます。おそろしいですね...！

ミラーラビン素数判定法の精度

上記のやり方で判定すると、1回の判定につき$3/4$以上の精度で合成数をふるい落すことができることが知られています。試行回数を増やすほど、指数関数的に判定精度は上がるそうです。
隕石に遭遇する確率の方が高いのではないか、と思える精度が100回ぐらいの試行で実現できるので、実用性はバツグンです。
このあたりの証明までは追えていないので、そのうち理解できるといいな。。。

まとめ

• フェルマーテストすごい
• ミラーラビンの判定法はもっとすごい
• フェルマーの小定理はえらい