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を使おうと思った次第です。