Wednesday, November 12, 2008

素数判定

ある数が素数かどうかを判定するときに、数をpとして、2からpの平方根までの数で割って確かめます。これだとpがおおきくなったときに時間がかかってしまいます。

ものの本を見ると、フェルマーの定理(小定理の方)で、pが素数の場合、
a^(p-1) = 1 (mod p)
であることを利用する。つまりいくつかのaを選び、この式に当てはめてみて、合格ならばpを素数の候補と考える、ということが書いてある。擬素数というのだそうです。なるほど、これはうまい手です。
pythonならば、

from random import randint

def is_pseudo_prime(p):
for i in xrange(0,10): #10個とって試してみる
a = randint(1, p-1)
if pow(a, p-1, p) != 1:
return False
return True

こんなもんです。大変簡単でよろしい。一方、ものの本の続きを見ると、Rabin Miller法というのが出ていて、これも似た考えのもの。GMPなどのライブラリにも入っています。

ところで、上の擬素数のコードをCで書いたら、非常に遅くてびっくりしました。ようするに

e = p-1;
a = randint(1, p-1);
r = 1;
while (e>0) {
r *=a;
r %= p;
e--;
}
return r;

といったことをしていたわけなのです。どうやら%=の処理が遅すぎる。というか多すぎる。pythonの組み込みのpowは便利で3項目のパラメタでmodの計算をしてくれているのです。じゃpythonはどうやってるのかとソースを見てみましたが、難しそうなことをやってます。

途方にくれていたのですが、wikiのページを見て、なるほどと思いました。上記のpow相当のwhileのところを、こう書くのです。

while (e > 0) {
if ((e & 1) == 1) r = (r * a) % p;
e >>= 1;
a = (a * a) % p;
}
return r;

詳しい説明はwikiの方を見てください。ようするにa同士を自乗の回数でかけるだけだから、もとのe回に比べてlog(e)に近いというわけです。これですばらしく速くなりました。

ただ、Cだと、変数を64bitにしても、かけ算が入るので、pはせいぜい32bitの範囲までしかできないのですよね。ちょっと残念です。
やっぱりGMPを使おうと思った次第です。

6 comments:

cobacco said...

おもしろい!

村松 said...

素数判定は P だという話をきいたことがあるよ。

Tsuchiya Yoshihiro said...

Pってpolynomialのこと?
割れるかどうかは
2,3,5,...,n^0.5を試せばいいので、nの大きさにたいしてO(n^0.5)では。

ぐれこ said...

計算量を言うときには、入力長に対してのオーダーなので、多項式時間というのは、log(n)の多項式オーダーです。発見者の頭文字をとったAKSという確定的素数判定法が何年か前に発見されて素数判定問題はPに属するということがわかって話題になりました。ルートnまで順に割っていくのは、最も素朴で最も効率が悪い方法で、計算量はlog(n)の指数時間オーダーです。

Tsuchiya Yoshihiro said...

なるほど。これのことですね。
http://en.wikipedia.org/wiki/AKS_primality_test
http://ja.wikipedia.org/wiki/AKS素数判定法
ありがとうございます。勉強してみます。

ぐれこ said...

ついでに言うと、RSA暗号のアルゴリズムは、べき乗剰余の計算そのものですが、1024ビットの数で割った余りを取るので、普通に書くと動かないため、高速化の研究が何十年も行われて、今では何気なくSSLで使われていたりしますので、こういう歴史を追ってみるのも面白いかもしれません。