約分、素因数分解、素数判定

約分、素因数分解、素数判定

この３つは実質同じことをする。

ある数について、約数があるかどうかを確かめる。

たとえば、「3007/1067 という分数を約分しなさい」という問題。

パっとみ、割れる数が思いつかない。2, 3, 5あたりはダメ。

すぐにできる約数を見つける方法は、

末尾が偶数... 2

各桁の数を足すと３の倍数... 3

末尾が5か0... 5

そして、２桁だったら素数は列挙してもそんなに時間はかからない。

2, 3, 5, 7, 11, 13, 17, 19, 23, 

くらいまではすぐ出る。

29, 31, 37, 41, 43, 47, 53, 

この辺から怪しくなってくるか

でも２桁だったら、「素数を列挙してその中にあるか」ではなく、

いくつかの数で割ってみればわかる。１分もかからないだろう。

また、桁数は１０桁とか２０桁とか、「どんな桁数であっても使える素数判定法」を見つけようというのでもない。

まあ、４、５桁くらい。

先ほど例にあげた、3007, 1067 といったような数字で使える、そんなに厳密でなくてもいい方法。

この例では、２つの数字がともに末尾が7となっている。

もしこれらの数字に約数があるならば、２つの数値をかけて結果の末尾（１の位）が7になる数字、ということである。

それはどういう数かと考えて、まず思ったのは、「末尾が7である数と、末尾が1である数の積」である。そして、これだけだと思った。

しかし、それだけではなく、3x9=27 でもある。

ただ、これだけわかればだいぶ楽だ。

3007について。

2, 3, 5 では割れないことはわかっている。

7は、微妙だが割ってみればすぐわかる。割れない。

8, 9, 10 はそれぞれその約数で割れないのですでに対象外。

２桁になっていく。

４桁の数で、「約分しなさい」という問題だから、おそらく２桁x２桁になるのだろう。

候補は下１桁が1, 3, 7, 9

これくらいなら順番に割ってもよさそうだが、

目的はこの問題を解くことではなく簡易的でもいいから判定法を見つけることなので、

もう少し考える。

3007という数の大きさ。

3000は、10x300, 15x200, 30x100 などである。（ほかにもあるが、たとえば。）

ここに登場する、10, 15, 30, 100, 200, 300 などをいじって、

先ほど候補として絞った「下一桁が1と7、あるいは3と9」にしたらどうだろうか？

10x300→ 11x299, 13x299, 17x291

15x200→13x209, 17x191, 19x193

30x100→31x97, 33x99, 29x93

適当にやっただけなので、ざっと見ると

13x209, 29x93は3000に満たないので除外

11x299, 13x299, 17x291

17x191, 19x193

31x97, 33x99

13x299は明らかに大きい

19x193も大きい

33x99も大きい

11x299, 17x291

17x191, 31x97,

おっと、末尾が7にならないのが混じっていた。

17x291, 17x191, 31x97,

17x291=17x(300-9)=5100-153

17x191=17x(200-9)=3400-153

31x97=31x(100-3)=3100-93=3007 ... これ！

1067は...

同様にやってもいいのだが

すでに3007=31x97が判明しているので

どっちかで割れるであろう。

31x37かな

(30+1)x37 ... 違う

31x27=(30+1)x27=810+27... 違う

97x11は...

(100-3)x11=1100-33=1067... これ！

ちなみに1067の方を先にやった場合の候補は

1067だと、1000との差がやや大きい。

1000=10x100, 20x50, 500x2, 330x3+α, 250x4

で、末尾が7になるから

11x97, 21x47, 

500x2, 330x3, 250x4は除外

もう答えを知ってるから 11x97がすぐ出てきちゃったけど...

bigint、素数判定

というものがあることを知った。cpanからインストールして

use bigint;

を書くと1023乗より大きい数を計算できた。

2のベキ乗を計算させると

1.28421286658896e+207

という表示ではなく、全桁表示する。

たとえば２の1247乗

2423285551989543969259886147306320615721694717012975552426444448158985017722789267546553034738712987127346362442309271495645764807314487385596126924659433020959638410571315406303196994043324038030803068668509323897700215957022383545283810557899591580902307500436706661372076856211818662627186819885605667216486349283517459646495626188985295134800963771933733229691796170211328

ちなみにこんなのもある。

なんでこれで判定できるのかさっぱりわからない。

perl -le 'print "PRIME" if (1 x shift) !~ /^(11+)\1+$/' 19

http://www.drdobbs.com/web-development/tpj-one-liners/184416234




「フェルマー小定理」方式C言語版

判定できる最大の素数は 16381である。

#include <stdio.h>
#include <math.h>

long double gojo(long double a, long double b);

int main(int argc, char *argv[]){
long double d;
long double  a, m, p, x;
int count, max;

count = 0;
a = 2;

if(argc > 1) {
    max = atoi(argv[1]);
}else{
    printf("please specify max number.\n");
    return 1;
}

printf("start\n");

for(p=3;p<max;p++){
    x = gojo(p, 2L);
    if(x == 1){
        d = powl(a, p-1);
        m = fmodl(d, p);
        if(m == 1){
            printf ("%.1Lf is PRIME\n", p);
            count++;
        }
    }
}

printf("end\n");
printf("count:%d\n", count);
}

long double gojo(long double a, long double b){
    long double c;

    while(b>0) {
        c = fmodl(a, b);
        a = b;
        b = c;
    }
    return a;
}

素数判定

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乗までしかできない。