確率的素数判定法 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桁までの判定時間をグラフにしてみると

こんな感じ。
300桁でも20個の素数を判定するのに3秒程度で終わってしまう。1個当たり0.15秒(!!)
エラトステネスの篩が6桁の素数を判定している間に、ミラーラビン素数判定では300桁の素数を判定できてしまう（驚）

補足

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

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

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

まとめ

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

雑感

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

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

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

word2vecで「単語の足し算引き算」をしてみる

word2vecを試してみたいけど使ったことがないという人が対象
とりあえず動かしてみるのが目的。
※python2.7.12で動かしているけど、3系でも大丈夫なはず。

word2vecとは

word2vecとは、ざっくり言えば
「単語をベクトル表現にでき、意味の足し算や引き算ができるニューラルネット」
これまでの言語処理ではone hotという表現で単語の数だけ次元を持つベクトルで処理していたけれど、次元圧縮したベクトルで単語を意味（的なもの）で表せるようになったってのが凄いよねって話。
よく使われる有名な例がこれ

king - man + woman = queen

とりあえずこれを動かしてみたい！というのが今回の記事
詳しい解説は下の参考サイトがめっちゃ分かりやすかったのでそちらを参照して欲しい。

必要なライブラリのインストール

pipで簡単にインストールできる

pip install numpy
pip install gensim

学習用のデータセット

日本語でもできるが、簡単に手に入るデータセットがないので英語で試してみる。
英語では簡単に手に入るデータセットが用意されているのでそれを使う。

wget http://mattmahoney.net/dc/text8.zip
unzip text8.zip

学習用のスクリプト

learning.pyという名前で保存

from gensim.models import word2vec
import logging

logging.basicConfig(format='%(asctime)s : %(levelname)s', level=logging.INFO)

sentences = word2vec.Text8Corpus('text8')
model = word2vec.Word2Vec(sentences, size=200, min_count=20, window=15)
model.save("sample.model")
print("Finish!")

モデルを試してみるスクリプト

result.pyという名前で保存

from gensim.models import word2vec
import logging
import sys

model = word2vec.Word2Vec.load("sample.model")

def neighbor_word(posi, nega=[], n=10):
    count = 1
    result = model.most_similar(positive = posi, negative = nega, topn = n)
    for r in result:
        print(str(count)+" "+str(r[0])+" "+str(r[1]))
        count += 1


def calc(equation):
    if "+" not in equation or "-" not in equation:
        neighbor_word([equation])
    else:
        posi,nega = [],[]
        positives = equation.split("+")
        for positive in positives:
            negatives = positive.split("-")
            posi.append(negatives[0])
            nega = nega + negatives[1:]
        neighbor_word(posi = posi, nega = nega)

if __name__=="__main__":
    equation = sys.argv[1]
    calc(equation)

使い方

python leraning.py
python result.py king
python result.py king-man+woman

learning.pyはモデルを作るためのものなので1回だけで充分。結構(3分くらい)時間がかかった。
足し算引き算はスペースを入れないで書くこと。
確かに第一候補が

1 queen 0.660548448563

になった！嬉しい！

参考サイト：
Word2VecをPythonでやってみる | Foolean
Word2vec Tutorial | RaRe Technologies

word2vecについてわかりやすい記事：
Word2Vec のニューラルネットワーク学習過程を理解する · けんごのお屋敷

pythonの仮想環境に関する話

pythonのvirtualenvで仮想環境の作成について。windows版

どうして仮想環境を使うのか

チームで開発する時など、環境を統一する必要がある。
ネットで公開する時は、ゼロから環境構築を再現して試してみたい時がある。
仮想環境ではライブラリは何もないのでゼロから環境構築を再現することができる。
加えてpythonは特に、2系と3系で異なるところが多いので使い分けられたら嬉しいことが多い。

使い方

準備

virtualenv自体はpipでinstallできる

pip install virtualenv

環境を作る

新しく環境を作るには、作業ディレクトリでこのコマンドを実行

virtualenv test1

これでtest1という名前の環境用のフォルダができる。

2系と3系を使い分けたい場合pythonのフォルダがあるパスを指定する

virtualenv --python=PATH test1

環境への出入り

環境に入る場合は
環境用のフォルダがある場所で（virtualenvを実行した場所で）

source activate test1

他の場所からこの環境を使う場合はそのフォルダへのパスにしてやる

環境から出る場合は

source deactivate

環境を共有する

今の環境を他のマシンや新しい環境でも再現したい時はどうするか。
virtualenvではpipでパッケージを管理すると思うけれど、

pip freeze > requirements.txt

として、パッケージを書き出すことができる。
ここに書き出したパッケージを一括で落としてくることができる。

pip install -r requirements.txt

パソコンを買い替える時とかにも使えるね。

補足

anacondaでも簡単に仮想環境を作ることができる。
特に異なるバージョンを動かしたい時は、いちいち別バージョンのpythonを予め持っている必要はないので手軽。
例えば、python3.6の環境を作るコマンドは

conda create -n test1 python=3.6

環境用のフォルダが
…/Anaconda2/envs/test1
のようにできるので(上の場合はデフォルトで2系を落としている時)
Anacondaのパスを通しているとどこからでも

source activate test1

として環境に入ることができる。

関連：