Pythonを使って高速素数判定をしてみる

みなさん、素数を数えてますか？
『素数』は１と自分の数でしか割ることのできない孤独な数字。
暗号化できたり、乱数を作れたり、心を落ち着いたりして、私達に勇気を与えてくれます。
素数といえば「エラトステネスのふるい」ですが、あれは大きい桁の素数を生成しようとすると、とんでもなく時間が掛ります。
今回は、どんな大きな桁の素数でも高速で素数判定するプログラムを作ってみます。

基本は「フェルマーの小定理」

素数判定の基本は「フェルマーの小定理」です、数式は１行だけのごく簡単なものです。

a^(p-1) mod p の答えが1以外ならpは合成数である
ただし、aとpが素の関係（最大公約数が1)であること


2つの数を「べき剰余算」して答えが1以外なら合成数(not 素数)という事です。
aに2を代入してqが素数なら答えが1になる、たったこれだけです簡単でしょ？

def is_prime(q):
    q = abs(q)
    if q == 2: return True
    if q < 2 or q&1 == 0: return False
    return pow(2, q-1, q) == 1


引数に素数を代入するとTrueを返します。
Pythonのべき乗関数pow()ですが、第3引数を代入するとべき剰余算してくれます！なんてめんこい言語なのかしら。
（きっと中では高速なモンゴメリ乗算をしてくれているはず！）
たったこれだけです、みなさん思いついた素数を入れて遊んでみてください。

>>> is_prime(5)
True
>>> is_prime(9)
False
>>> is_prime(19)
True
>>> is_prime(45877878467)
True


こんな短いコードなのに、ちゃんと出来てますね！！
ちなみにフェルマーの小定理を使った素数判定を「フェルマーテスト」と言います。

改良版フェルマーテスト

お手軽＆高速なフェルマーテストですが、重大な欠陥があります。


じつは、たまに素数じゃない数もフェルマーテスト通っちゃいます、テヘ☆


フェルマーテストが通っちゃう合成数を「疑素数」といいます。
a=2のとき2〜10,000まで間で素数を判定すると1,250個の素数候補を検出しますが、うち22個が疑素数ですので約1.76%の確率で誤判定をします。
以下a=2のときの10,000の疑素数です。


341,561,645,1105,1387,1729,1905,2047,2465,2701,2821,3277,
4033,4369,4371,4681,5461,6601,7957,8321,8481,8911


さて、a=2の疑素数「341」ですが、a=3の時はどうでしょうか？

>>> pow(3, 341-1, 341)
56

ちゃんと合成数と判定していますね、aの値が変わると判定結果も変わるんですね。
ならば、aの値を変えながら100回くらいテストしちゃえば、正確な素数が出るんじゃないでしょうか。
以下が改良版のフェルマーテストです。

def is_prime2(q,k=100):
    q = abs(q)
    if q == 2: return True
    if q < 2 or q&1 == 0: return False
    for i in xrange(3,k):
        x,y = q,i
        while y:
            x, y =  y, x % y
        if x != 1: continue
        if pow(i, q-1, q) != 1:
            return False
    return True


aの値を第2引数kの分だけ回しますのですが、aとqは最大公約数が1の関係でなければいけないので「ユークリッドの互除法」を使ってチェックをしています。
さっそく、前回のフェルマーテストで漏れた疑素数をテストしてみましょう。

pseudoprime = [341,561,645,1105,1387,1729,1905,2047,2465,2701,
               2821,3277,4033,4369,4371,4681,5461,6601,7957,8321,
               8481,8911]

for x in pseudoprime:
    if is_prime2(x):
        print(x)

#---出力結果
561
1105
1729
2465
2821
6601
8911


・・・ゲェー！！これもすり抜ける疑素数がある！！

ミラー・ラビン(Miller-Rabin)テスト

・・・ゲェー！！これもすり抜ける疑素数がある！！
(先週からのつづき）


うそうそ、実はこうなるの知ってた。
561とかの数字は「カーマイケル数」って言って、aにどんな値が入ってもフェルマーテストが合格しちゃうやつらです。
カーマイケル数は無限にあることも知られていて、簡単な判別方法は存在しません。


せっかく高速なフェルマーテストですが、カーマイケル数のせいで使えなくなるのは惜しいですよね。
そこで登場するのがフェルマーテストの改良版でもある「ミラー・ラビンテスト」です。

ミラー・ラビンテストを簡単に説明すると、判定したい正の奇数をqとした場合、

q - 1 = 2^s * d
(但しsは正整数、dは奇数とする）

の式を解き、sとdを求めます。


そして、1〜q-1までの間の数を適当に選びaとし、「pow(a,d,q) != 1」がTrueだった場合、0〜s-1までのカウンタをrとし、すべての「pow(a, 2**(r*d), d)!=-1」がTrueになると、qは合成数と判断できます。


時々、誤った判断を下すaが現れますが、何回かaをランダムで選ぶことにより誤る確率を下げます。
以下、ミラー・ラビンテストの関数です。

#coding:utf8
import random

def is_prime3(q,k=50):
    q = abs(q)
    #計算するまでもなく判定できるものははじく
    if q == 2: return True
    if q < 2 or q&1 == 0: return False

    #n-1=2^s*dとし（但しaは整数、dは奇数)、dを求める
    d = (q-1)>>1
    while d&1 == 0:
        d >>= 1
    
    #判定をk回繰り返す
    for i in xrange(k):
        a = random.randint(1,q-1)
        t = d
        y = pow(a,t,q)
        #[0,s-1]の範囲すべてをチェック
        while t != q-1 and y != 1 and y != q-1: 
            y = pow(y,2,q)
            t <<= 1
        if y != q-1 and t&1 == 0:
            return False
    return True


それでは試してみましょう。

>>> is_prime3(13)
True
>>> is_prime3(561)
False
>>> is_prime3(8911)
False
>>> is_prime3(45877878467)
True


できました！
疑素数とカーマイケル数も、ちゃんと見破ってますね。


でも、ミラー・ラビンテストも結局は確率的素数判定でありますので、運が悪いと判定を間違えてしまいます。
しかし、is_prime3()の第2引数kに大きい数を設定する事で、間違う確率を下げる事ができます。


ちなみに間違える確率は「4^-k」となっています、kに50を設定したとすると1/1267650600228229401496703205376となります。
逆に間違った場面に遭遇した方がラッキーだと思います。

まとめ

確率的素数判定のミラー・ラビンテストですが、そのぶんかなり高速です。
真面目に素数判定したら一生終わらない桁数まで、一瞬で判断することができるのです。


例えば、以下のように300桁の素数も簡単に判別可能です。


1049261498049476504856419608432108674896402341867409604
1684060461674968406126850469654006353403205485749261498
0494765048564196084321086748964023418674096041684060461
6749684061268504696540063534032054857492614980494765048
5641960843210867489640234186740960416840604616749684061
2685046965400635340329729


それでは次回は素数を使って、いろいろ楽しい事をしてみたいとおもいます。

RSA暗号で「ふっかつのじゅもん」を作る(1) - Pashango’s Blogへ続きます。