技術メモ

役に立てる技術的な何か、時々自分用の覚書。幅広く色々なことに興味があります。

確率的素数判定法 vs. エラトステネスの篩

素数判定をするのは暗号の分野ではすごく重要なんだけど、そこでよく使われるのが確率的素数判定法というもの。
競技プログラミングの問題に出てくる程度の素数判定では、エラトステネスの篩と呼ばれる有名な素数判定アルゴリズムで充分なのだけど、暗号に使われるような数百桁の素数では遅すぎる。
じゃあどれくらい違うのというのを試してみたいと思う。

ちなみに、暗号の分野で素数がどう使われるかはこちら
RSA暗号を実装してみる(理論編) - 技術メモ

エラトステネスの篩

素数判定というよりは実際に素数を求める方法。
2から順に素数を求めていって、それまで求めてきた素数で割り切れないものが新しい素数だよねとする。
小さい数から反復的に素数を見つけていかなきゃいけない。

コード

def eratosthenes(n):
    def mark(primes,x):
        for i in xrange(2*x, n+1, x):
            primes[i-2] = False
    if n%2 == 0 or n%3 == 0 or n%5 == 0:
        return False
    else:
        primes = [(i%2!=0 and i%3!=0 and i%5!=0) for i in range(2,n+1)]
        primes[0] = primes[1] = primes[3] = True
        for x in xrange(7, int(np.sqrt(n+1))+1):
           if primes[x-2]: mark(primes, x)
        return primes[n-2]

確率的素数判定法

ある条件(フェルマーの小定理の対偶)を満たすなら素数じゃないねって判定を繰り返す方法。
いくつかあるが、今回使うのはミラーラビン素数判定法と呼ばれる方法。
どうして確率的と言うかというと、あくまで「素数じゃない」ことしかわからないから。
フェルマーの小定理は「素数ならAを満たす」ということしか言っていないので、その対偶の「Aを満たさないなら素数ではない」としか言い切れない。じゃあ「素数じゃなくてAを満たす」ものがあるかというと、実は存在する。
カーマイケル数と呼ばれ、求め方は未解決らしいけど無数に存在することだけは知られている。
解説を書きたかったけど長くなるので別の解説を読んでほしい。Wikipediaでも十分わかりやすかった。

コード

import random
def millerrabin(n):
    k = 50 # number of iteration times
    if n%2 == 0 or n%3 == 0 or n%5 == 0:
        return False
    else:
        s, d = 0, n-1
        while d%2==0: s,d = s+1, d/2
        while k > 0:
            k = k-1
            a = random.randint(1,n-1)
            t, y = d, pow(a,d,n)
            while t != n-1 and y != 1 and y != n-1:
                y = pow(y,2,n)
                t <<= 1
            if y != n-1 and t%2 == 0:
                return False
        return True

実験

※どちらも、明らかな合成数(2と3と5の倍数)は予め省いているので従来より少しだけ早くて正確なはず。
まず確認

print eratosthenes(11) # True
print eratosthenes(6499) # False
print eratosthenes(58427) # True
print millerrabin(11) #True
print millerrabin(6499) #False
print millerrabin(58427) #True

まあ合ってそう。(他にも色々試してみました)

3桁、4桁、…200桁の数をそれぞれ20個用意し、判定時間を見てみる。(CPU: corei5, RAM: 8GB)

from time import time
def testcase():
    for order in range(2,200):
        test_set = [random.randint(pow(10,order), pow(10,order+1)) for i in range(20)]
        start = time()
        for ts in test_set:
            millerrabin(ts) # or eratosthenes(ts)
        end = time()
        print "(%d)time: %f"%(order+1,end-start)

エラトステネスの篩

(3)time: 0.000000
(4)time: 0.004000
(5)time: 0.069000
(6)time: 1.431000
(7)time: 6.620000

8桁目を計算する時終わりそうな気配がなかったのでここで断念。
エラトステネスの篩で素数判定できるのは頑張って8桁って感じかな。

確率的素数判定法

(3)time: 0.000000
(4)time: 0.001000
(5)time: 0.000000
(6)time: 0.001000
…
(197)time: 0.008000
(198)time: 0.011000
(199)time: 0.009000
(200)time: 0.014000

200桁でもこの速さ。


ちなみに、純粋な方法(明らかな合成数を省く処理の最初の2行を除く)にして300桁までの判定時間をグラフにしてみると
f:id:swdrsker:20170226121932p:plain
こんな感じ。
300桁でも20個の素数を判定するのに3秒程度で終わってしまう。1個当たり0.15秒(!!)
エラトステネスの篩が6桁の素数を判定している間に、ミラーラビン素数判定では300桁の素数を判定できてしまう(驚)

補足

ミラーラビン素数判定についてWikipediaによると

平方の繰り返しをべき剰余を使って実現すると、このアルゴリズムの実行時間は O(k * log^3 n) となる。ここで、k は異なる a で判定を行う回数である。以上から、このアルゴリズム多項式時間の効率のよいアルゴリズムであることがわかる。

なるほど、言われてみれば上のグラフは3次曲線っぽい。
※上のグラフは横軸が「桁数」なのでlogスケールになっていることに注意

まとめ

6桁以下ならエラトステネスの篩で確実に判定するのが良さそう。7桁以上なら確率的な判定方法でやってみよう。
Wikipediaによると誤判定の確率は4^{-k}らしい。

雑感

結局、ハイブリッドにするのが一番いいんじゃないか。
「3桁くらいまでの素数は予め決定的な方法(エラトステネスの篩)でリスト化しいておいて、それ以上の数の判定の際はリストにある素数で割り切れるかをチェックした上で確率的な方法で判定する。」
これが時間的にも精度的にも"最強"の方法だと思う。

[追記]
提案方法を実装してみた→最強の素数チェッカーを作ってみた - 技術メモ



参考サイト:
Pythonを使って高速素数判定をしてみる - Pashango’s Blog
ミラー–ラビン素数判定法 - Wikipedia
素数判定 - Wikipedia