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

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

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

エラトステネスの篩

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

コード

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桁の素数を判定できてしまう（驚）