前提とする知識

合同式に関する知識が必要ですが、正直割った余り程度に考えてくれれば構いません。

趣旨

暗号を含むを多くの応用では巨大でランダムな素数を発見する必要があります。
RSA公開鍵暗号系などが有名な例です。

ランダムに整数$n$を選択し、それが素数であるかどうかを判定するのに必要な試行回数の期待値はほぼ$\ln n$となります。(メンドウなので理由は省略します。素数定理とかでググれば多分証明あると思う。)

このことから$n$と同じ長さの素数を見つけるには、$n$の周辺からランダムに選択したおおよそ$\ln n$個の整数を調べれば良いと期待できます。

たとえば、2048ビットの素数を見つけるためには、ランダムに選択したおおよそ$\ln 2^{2048}\approx$1420個の2048ビット数を調べればよいと期待できることになります。

以下では大きな整数$n$が素数であるか判定する問題を考えていきます。

素数判定法

試し割り

素数性を判定する簡単な方法として各整数$2, 3, ..., \lfloor \sqrt n \rfloor$による$n$の除算を試行する試し割り法がある。
$n$が素数であることの必要十分条件はどの試行でも$n$が割り切れないことである。
それぞれの試行が定数時間で実行できると仮定すると、最悪実行時間は$\Theta(\sqrt n)$であり、これは$n$の長さに関して指数関数的に伸びる。

愚直に実装すると以下のようになる。

def is_prime(n):
    for i in range(2, int(n**0.5) + 1):
        if n % i == 0:
            return False
    return n != 1

もちろん愚直な実装なので改善点はたくさんあるのだが、実際のところこれで十分だろう。

正直なところ$n=2^{32}$程度であれば計算できてしまっているようだが、$n=2^{1024}$などとなると流石に無理と思われる。

ではどうすればよいのだろうか

疑似素数判定法

ところで、フェルマーの小定理をご存知でしょうか。

$n$を素数、$a$を$n$と互いに素な整数とすると、
$$a^{n-1} \equiv 1 \pmod{n}$$
が成り立つ。
というやつです。このことから
$$a^{n-1} \not\equiv 1 \pmod{n}$$
となるような$a$を発見できたら、$n$は合成数であることがわかります。
そして、ほとんどの場合、この逆も成り立ちます。($a^{n-1} \equiv 1 \pmod{n}$を満たしてしまうような合成数$n$は$a$を底とする疑似素数と呼ばれます)

というわけで以下に$a=2$とした場合に$n$が素数であるかどうかを擬似的に判定するコードを記します。

# 入力nは2より大きい奇数であることを仮定します。
def pseudo_prime(n):
    if pow(2, n - 1, n) != 1:
        return "間違いなく合成数"
    else:
        return "素数...と思われる"

ちなみにpow(a, b, n)は冪剰余$a^{b}\mod{n}$を計算します。
冪剰余についてはWikiにも簡単に書いてありますので気になったら参照ください。

$10,000$未満の$n$の値でpseudo_prime(n)が間違えるのは$22$個であり、
$341, 561, 645,1105$が最初の4つの値となります。

実際に上のプログラムを少し手を加えたものに$5, 9, 19, 341$を与えると以下のように出力されました。

5は素数...と思われる
9は間違いなく合成数
19は素数...と思われる
341は素数...と思われる

最後341でばっちし間違ってくれています。

ところで、$gcd(a, n)=1$を満たすような任意の整数$a$に対して
$$a^{n-1} \equiv 1 \pmod{n}$$
を満たす合成数$n$が存在することが知られています。（カーマイケル数とか言われてる）

なので、別の底の候補、たとえば、$a=3$を選んでテストするだけではすべての間違いを排除することはできません。¹

ではカーマイケル数に負けないようにするにはどうすればよいのだろうか

Miller-Rabinの乱択素数判定法

合同式$x^2 \equiv 1 \pmod n$の解$x$が、法$n$の下でこの合同式の2つの”自明”な平方根$1$と$-1$
のどちらとも合同ではないとき、法$n$の下での$1$の非自明な平方根といいます。

それで証明は省略するのですが、

法$n$の下で1の非自明な平方根が存在すれば、$n$は合成数です。

気になる方はWikiにある程度詳細があるようなのでご参照ください。

Miller-Rabinの乱択素数判定法は前に紹介した疑似素数判定法に以下２つの修正を盛り込んだ内容となっています。

1. 底$a$を１つの値だけではなく、ランダムに選択した複数の値を試す。
2. 各冪剰余を計算している間、法$n$の下で非自明な$1$の平方根を見つけたら合成数とする。

2についてですが$n-1=2^{t}u$とすると、$a^{n-1} \equiv (a^u)^{2^t} \pmod n$なので、
$a^{n-1}\mod n$を計算するには、$a^u \mod n$を計算して、その結果を2乗する操作をt回続ければいいです。

というわけで本を参考にPythonでそれっぽく実装してみました。

# 疑似素数判定pseudo_primeの拡張関数
def witness(a, n, t, u):
    x = pow(a, u, n)  # x0 = a^u mod nを計算

    for _ in range(0, t):
        y = (x * x) % n

        # 法nの下での1の非自明な平方根が存在するのはnが合成数であるときに限る
        if y == 1 and x != 1 and x != (n - 1):
            return True
        x = y
    return y != 1  # xt ≡ a^(n-1) (mod n)が1でなければTrueそうでなければFalseを返す

# 確率的素数判定を行う.
# 入力はn>2を満たし、witnessesは証拠となるaの配列となる。
def is_probably_prime(n, witnesses):
    t = 1
    u = n >> 1
    while u & 1 == 0:
        t = t + 1
        u >>= 1
    assert(2**t * u == n - 1)

    for a in witnesses:
        if a < n and witness(a, n, t, u):
            return False
    return True

ちなみに証人となる$a$をどれだけ用意すればよいかというのはある程度の範囲でわかっているようだ。
詳しい$a$の値についてはここに記載されている。

とりあえず、$n < 2^{32}$の範囲でならば$a=2, 7, 61$で試してみたらよいみたい。

# n < 4,759,123,141を想定
def is_definitely_prime(n):
    if ((not (n & 1)) and n != 2) or (n < 2) or ((n % 3 == 0) and (n != 3)):
        return False
    elif n <= 3:
        return True
    else:
        return is_probably_prime(n, [2, 7, 61])

AOJの素数判定問題で動作を確認できたので大丈夫そう。
カーマイケル数にも負けていませんでした。

実際、どのくらいの精度なのかという話はそもそも$a$の値をランダムに決定した場合とそうでない場合で証明できる範囲が違うっぽいです。
まあ$a$で試す回数が多いほど間違える確率は下がるので実用的にはたくさん試しているのでしょう。
実際は試す回数が少なくても間違えている可能性はかなり低いようですが。

参考

アルゴリズムイントロダクション

さいごに

最近、プログラミングに対するモチベーションが低く、リハビリとしてこの記事を書いた。
なのであまり初心者向けというよりか自分が思い出すために書いた部分も大きい。
正直なところ需要がわからんし、こんな記事書いてもカネになるわけでもないし。

少しでもこの記事に関心したらモチベ回復のためにいいねボタンを押してくれると助かる。

また、わかりづらかったらコメント欄でアドバイスなどあればうれしいです。
編集リクエストでも構いません。

それと、「へぇーこんな記事書くやつがいるんだ」みたいなこと少しでも思ったら私の書いた他の記事も一度見てみてほしいです。
もっとも、書いている記事のレベルがバラバラすぎるので参考になるとは限りませんが、、

( ｀・∀・´)ﾉﾖﾛｼｸ