このブログを検索

2013/04/04

素数判定

youtubeのおすすめ動画に、NHKスペシャルの「リーマン予想」が出てきたので見てみた。

「リーマン予想」で検索してみると、専門家から言わせるとこの番組にはいろいろボロがあるようであるが、リーマン予想というのは、素数の出現に関する規則性のようなものだということはわかった。素数の並びそのものではなく、ゼータ関数と呼ばれる数式のある値に関する規則性だそうである。

この規則性については、正しい値が多数存在することは証明されても、その例外がみつかっていないのでほぼ正しいだろうとは考えられているものの、正しいという証明がいまだにされておらず、100万ドルの懸賞金までかかっているそうである。

私からすると、素数というものになぜ数学者が魅了されるのかわからない。そんなものに規則性なんかあるわけがないというのが私の直感的な考えであるが、どうやら規則性があるようなのである。

というわけで、とりあえず素数判定プログラムを書いてみようと思った。 

use strict;
use warnings;

my $end = shift || 10000;
my $count = 0;

my $starttime = localtime();

for (2..$end) {
    if(isprime($_) == 0){
        $count++;
        print $_;
        print "\t";
    }
}

print "\ncount: $count\n";

my $endtime = localtime();

print $starttime." - ".$endtime;

sub isprime {
    my ($num) = @_;

    for (my $i=2;$i<$num;$i++){
        if(($num % $i) == 0){
            return 1;
        }
    }
    return 0;
}

これだと100000以下の素数を列挙するのに1分半くらいかかる。

2から順に割っているのだが、判定対象の数-1まで割るのはムダである。たとえば1600が800で割れるというのは、20で割ったときにわかる。では、いくつまで割ればいいのか？

これは昔数学の授業で習った記憶がある。最初、「1/2かな」と思って、やってみると素数の個数からして正しく判定できたようだ。

時間は20秒くらい短くなった。だが、なぜ1/2かと言われると説明できない。

高校の数学Bの教科書を開いてみた。

「自然数Nが合成数ならば、必ず√N以下の約数をもつ」

ということであった。だから、100なら10まで、10000なら100まで試し割をすればよいことになる。 

for (my $i=2;$i<=(sqrt $num);$i++){

このように修正すると、3秒くらいになった。

「エラトステネスのふるい」方式。

(2013/04/04修正） 

use strict;
use warnings;

my @array =();
my @new_array =();

my $max = shift || 1000;


my $start_time = localtime();

my $count=0;

for(0..$max-1){
    push @array, $count;
    $count++;
}

$array[1] = 0;

for(my $i=2;$i0){
        $count++;
        print;
        print "\t";
    }
}

print "\n $count \n";


my $end_time = localtime();

print "$start_time - $end_time\n";

「倍数をふるいにかける」ところを、mod (%)でやっているのはちょっとずるいかな。 （後記：ずるいというかムダ。$jを$iずつ増分した値を無条件に消していけばよい）

あと、「ふるいにかける」は、配列の要素を削除してやりたかったのだが、 spliceを使ったり別の配列にpushしたりしてみたがうまくできず、 値をゼロにして表示するときはゼロを飛ばすという方法にした。 100000まででやると、20秒かかる。試し割り方式より遅い。 （後記：修正後は約1秒になった）

さらに、素数判定では「フェルマーの小定理」を使うのが一番速いそうだ。 その定理とは「pを素数とし、aをpの倍数でない整数とするときにaのp-1乗をpで割った余りは1となる(wikipediaより）」 である。

実際に確認してみようか。

p=5とする。aは・・・3にしようか。

3^(5-1) / 5 の余りが 1になるということだね。

3^4 = 81

5でわると16...1だね。

p=23, a=17でやってみよう。

17^22 = 1174562876521148458974062689

1174562876521148458974062689 mod 23 は、電卓で計算したら1になった。

aとpの二つの数字が必要で、どちらかが判定対象だとして、もう1個の数字は互いに素になればなんでもいいのか?

p=100, a=3

3^99 を100でわった余りは67、だから素数でない。 これをプログラムする場合、判定対象をpとし、それと互いに素となるaを選ばねばならない。 a=2 としておいて、それが対象と互いに素でなければその時点で素数ではない。

プログラミングすると・・・ 

$b = 2;
$max = shift;

print "2\t";
$count = 1;

for($p=3;$p<$max;$p++){
    if (&gojo($p, $b) ==1){
        if($b**($p-1) % $p == 1){
            print "$p\t";
            $count++;
        }
    }
    $b = 2;
}
print "\ncount: $count\n";

sub gojo{
    my ($a,$b) =@_;

    while ($b>0) {
        $c = ($a % $b);
        $a = $b;
        $b = $c;
    }
    return $a;
}

100までは正しく判定できているようだ。 1000までにすると3個、余計に素数と判定したものが混じる。

誤判定した数字は、341, 561, 645。 電卓で計算すると、2^(p-1) mod p は皆1になる。 561は「カーマイケル数」だが、あとの二つは違う。 645は5で割れる。341は11で割れる。561も11で割れる。

これらはa=2とした場合に誤判定される数値のようだ。 1021までしか判定できないのだが、2のべき乗を計算すると1023乗までしかできない。

at
Labels: ,