読者です 読者をやめる 読者になる 読者になる

技術メモ

役に立てる技術的な何か、時々自分用の覚書。統計学、機械学習、神経科学、暗号理論...と幅広く興味があります。

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

Python 暗号理論 数学

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

ちなみに、暗号の分野で素数がどう使われるかはこちら
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秒程度で終わってしまう。

補足

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

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

なるほど、言われてみれば上のグラフは3次曲線っぽい。


エラトステネスの篩は遅い代わりに確実に素数を判定ができる。一方確率的な判定法では早さはずば抜けているがあくまで確率的な方法である。Wikipediaによると誤判定の確率は4^{-k}らしいが。

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

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



参考サイト:
http://d.hatena.ne.jp/pashango_p/20090704/1246692091
ミラー–ラビン素数判定法 - Wikipedia
素数判定 - Wikipedia

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

Python 機械学習

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

word2vecとは

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

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 のニューラルネットワーク学習過程を理解する · けんごのお屋敷