[PHP] 1,000,000番目の素数を求める

素数を見つけろ!

まずは100万番目を計算する気にもならなかった愚直コードたちを紹介する!

全ての数を全ての数で割るコード(13.119s)

<?php
define("N", 10000);
$time=microtime(true);
for($i=3, $p=1, $flg=true; $p<N; $i++, $flg=true){
    for($k=2; $k<=$i; $k++) {
        if($i%$k==0) {
            $flg=false;
            break;
        }
    }
    if($flg==true) $p++;
    if($p==N) {
        print(round(microtime(true) - $time, 3)." ");
        print($i." is ".($p+1)." th prime number\n");
    }
}
?>

2は最初の素数で唯一偶数の素数なので取り除く。

あとは2から順番に増やしていってその数までの何かの数で割れたら素数ではないから出力するというコード。

10,000番目の素数を見つけるのにも13.119秒もかかったので10万番目すらチャレンジする気が起きなかった。

奇数だけを奇数で割るコード(7.881s)

<?php
define("N", 10000);
$time=microtime(true);
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=3; $k<=$i; $k+=2) {
        if($i%$k==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true) $p++;
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

よく考えれば2以上の偶数はそもそも素数ではないので奇数だけ素数判定をするように変えたコード。

チェックするべき数が半分に減ったので計算量が半分になり、結果ほぼ計算時間を半分にすることができた。

半分までを奇数で割るコード(6.543s)

<?php
define("N", 10000);
$time=microtime(true);
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=3; $k<=$1/2; $k+=2) {
        if($i%$k==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true) $p++;
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

1以上で割り切れるということは少なくとも2で割り切れるということである。

つまり、9999が5001で割り切れるはずがないのだから最初から考えなくていいということである。

1/3までを奇数で割るコード(4.331s)

<?php
define("N", 10000);
$time=microtime(true);
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=3; $k<=$i/3; $k+=2) {
        if($i%$k==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true) $p++;
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

よく考えると2で割り切れる数は最初から省いているので、半分ではなく1/3まで考えれば良いことがわかる。

除算は遅いぞ!

除算命令は加算・減算命令に比べてとてつもなく遅い。

乗算命令も遅いが除算命令に比べれば相当速いので、ループ中に何度も除算をする場合は乗算命令に変えたほうが良い。

// Before
for($k=3; $k<=$i/3; $k+=2)
// After
for($k=3; $k*3<=$i; $k+=2) 

たったこれだけ変えるだけで相当速くなるぞ!

除算命令を乗算に変更(2.852s)

<?php
define("N", 10000);
$time=microtime(true);
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=3; $k*3<=$i; $k+=2) {
        if($i%$k==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true) $p++;
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

数学の知識を使っていけ!

さて、全部の数を全部で割るというところから奇数を1/3までの数で割るというコードに書き換えることに成功した。

しかし、数学に詳しい人ならたぶん最初からわかっていたと思う。

「それ、平方根までで良くね?」

詳しい解説はおいておくとして、要するにある数Pが素数ではなかったとする。

するとなにかの数で割れるわけなのだが、小さい順にどんどん試していって割り切れた数が仮にMだったとしよう。

つまり、P÷M=Nといった風に書けるわけである。これは当然、P=M*Nと書き直すことができる。

そして、このMは高々Pの平方根くらいまでしか考えなくていいのである。

何故なら例えば91が何で割り切れるか考えたいとして、小さい数から順に割っていったときに13で割り切れることに気付いたとしよう。

実はこういうことは起こりえないのだが、わかるだろうか?

そう、13で割り切れることが分かる前により小さい7で割り切れることに気付いていなければおかしいからである。

これは極めて計算量を減らすことができる重要なテクニックであるのでしっかり覚えておくように。

平方根までで除算する(0.075s)

<?php
define("N", 10000);
$time=microtime(true);
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=3; $k<=sqrt($i); $k+=2) {
        if($i%$k==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true) $p++;
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

なんといきなり計算速度が40倍になってしまった。

これ以上高速化するととても10,000番目では差がつかないので、以後100,000番目の素数を求めるコードに変更する。

ちなみにこのコードだと100,000番目を計算するのにかかった時間は2.518秒である。

平方根を求めるのはPHPだとsqrt()関数である。

しかし、懸命な諸君なら気付いたであろう。

そう、このsqrt()関数は除算を伴うため非常に遅いのである。

gmp_sqrt()っていう数学専門っぽいライブラリ関数使ったらなんと4倍くらい遅くなりました。

なんじゃそりゃ!!

平方までで除算する(1.679s)

<?php
define("N", 100000);
$time=microtime(true);
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=3; $k*$k<=$i; $k+=2) {
        if($i%$k==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true) $p++;
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

ちなみにPHPには指数を計算するpow()という関数もあるが、これも使えば使うだけ遅くなるのでオススメしない。

pow()を使ったところ、計算時間が2.616sになってしまった。これでは除算をしているのと変わらないではないか!

三乗を計算したいなら$k*$k*$kと書くこと。愚直だけどこっちの方が速いのである。

奇数で割る必要ないぞ

素数かどうか調べたい数の平方根までを調べることでおよそ40倍の高速化を果たしたが、これではまだ遅い。

例えば、ある数が素数かどうか調べるのに9で割る必要があるだろうか?ないよね??

何故なら、9で割り切れるならそれよりも前に3で割り切れていることに気付いていなければおかしいからである。

つまり、調べたい数の平方までの素数で割り切れるかどうか調べれば良いことになる。

素数リストで除算(0.988s)

<?php
define("N", 100000);
$time=microtime(true);
$prime = array();
array_push($prime, 2);
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=0; $prime[$k]*$prime[$k]<=$i; $k++) {
        if($i%$prime[$k]==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true){
        array_push($prime, $i);
        $p++;
    }
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

新しい素数が見つかったらどんどん素数リストに追加してそれで割れるかどうかを試していくわけである。

これによってとうとう100,000番目の素数を1秒以内に見つけることに成功した!

ちなみに1,000,000番目の素数も27.031秒で見つけることができた!

array_push()は遅いぞ!

ライブラリ関数を使うと遅いというのはpow()で述べたが、array_push()も憎まれるくらい遅い。

配列の末尾に数を追加するのはもっと高速な方法があるのでそれで実装しよう。

// Before
array_push($prime, $i);
// After
$prime[] = $i;

これだけで平気で10%くらいパフォーマンスに差が出るぞ!

演算子で末尾に代入(24.686s)

<?php
define("N", 1000000);
$time=microtime(true);
$prime = array();
$prime[] = 2;
for($i=3, $p=1, $flg=true; $p<N; $i+=2, $flg=true){
    for($k=0; $prime[$k]*$prime[$k]<=$i; $k++) {
        if($i%$prime[$k]==0) {
            $flg=false;
            break;
        }
    }
    if($flg == true){
        $prime[] = $i;
        $p++;
    }
    if($p==N) {
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

そもそも除算がナンセンス

除算を乗算にしたり工夫はしているものの、結局その数が素数かどうかは除算して計算している。

「えっ、でも結局除算しなきゃ素数かどうか確かめられなくないですか?」と思ったそこのキミはまだまだ甘い。

見つけた素数の倍数を消してしまえば除算の代わりに乗算を使えるのだ。

ただし、此処から先は素数リストとは別に数リストなるものが必要となる。また、予め探したい(N番目の)素数がどのくらいの数までにあるかわかっていなければいけない。

エラトステネスの篩

<?php
define("N", 16000000);
$time=microtime(true);
$n=floor(sqrt(N));
$number=array_fill(2, N-1, true);
for($i=2; $i<=$n; $i++){
	if(isset($number[$i])){
		for($j=$i*2; $j<=N; $j+=$i){
			unset($number[$j]);
		}
	}
}
$number=array_keys($number);
print(round(microtime(true)-$time, 3) . " ");
print($number[999999]." is "."10000000"." th prime number\n");
?>

2からNまでの配列をあらかじめつくっておき、unset()を使って配列を削除していくわけである。

なお、unset()は除算を伴わないうため非常に高速である。

平方根を求めるsqrt()も使用するが、これは一回しか呼び出さないので全く気にしなくて良い。

今回の問題でいうとこの方法は予めN番目の素数が1600万以下だとわかっているから使える方法である。

これがエラトステネスの篩と呼ばれる素数発見法である。

エラトステネスの篩の実装にあたってビット演算がよくわからなかったので下記サイトを参考にさせていただきました。

https://tech.basicinc.jp/articles/52

注意点としてはエラトステネスの篩は最初に配列を確保するのでメモリをバカ食いします。

ライブラリ使えばええんや!

ミラーラビン素数判定法(7.402s)

<?php
define("N", 1000000);
$time=microtime(true);
for($i=3, $p=1; $p<1000000; $i+=2){
    if(gmp_prob_prime($i)==(1||2)) $p++;
    if($p==N){
        print(round(microtime(true)-$time, 3) . " ");
        print($i." is ".$p." th prime number\n");
    }
}
?>

ミラーラビン素数判定法のライブラリを使って計算してみました。

普通に素数リストで割るよりは4倍くらい速いですね。

next_prime()(15.12s)

<?php
define("N", 1000000);
$time=microtime(true);

for($i=1, $p=1; $p<=N; $p++){
    $i=gmp_nextprime($i);
}
print(round(microtime(true)-$time, 3) . " ");
print($i." is ".($p-1)." th prime number\n");
?>

next_prime()は次の素数を返す関数です。

はいはい、どうせ内部ではミラーラビン判定法が動いてるだけやろって思っていたら倍近く遅くなってびっくりしました、なんじゃこりゃ?

決定的ミラーラビン判定法(—)

<?php
define("N", 10000);
$time = microtime(true);

function is_prime($n)
{
    if ($n == 2) return true;
    if ($n == 3) return true;
    if ($n == 5) return true;
    if ($n < 2 || $n % 2 == 0) return false;

    $d = $n - 1;
    $s = 0;

    while (($d &amp; 1) == 0) {
        $d = $d >> 1;
        $s++;
    }

    foreach ($a as $k) {
        $x = bcpowmod($k, $d, $n);
        if ($x == 1 || $x == $n - 1) continue;
        for ($r = 1; $r < $s; $r++) {
            $x = bcmod(bcmul($x, $x), $n);
            if ($x == 1) return false;
            if ($x == $n - 1) continue 2;
        }
        return false;
    }
    return true;
}

for ($i = 1, $p = 1; $p <= N; $i++) {
    if (is_prime($i)) {
        $p++;
    }
    if ($p == N + 1) {
        print(round(microtime(true) - $time, 3) . " ");
        print($i . " is " . ($p - 1) . " th prime number\n");
    }
}
?>

PHPのライブラリのミラーラビンアルゴリズムは確率的なので決定的アルゴリズムに直してみました。

テストするaとして2, 3, 5を選べば25,326,001までの数については素数かどうか決定的に判定することができます。

結果発表

10,000100,0001,000,000
コード113.119s
コード27.881s
コード36.543s
コード44.331s
コード52.852s
sqrt()0.075s2.518s80.231s
gmp_sqrt()0.293s9.933s
$k*$k0.045s1.679s55.65s
素数リスト0.04s0.988s27.031s
素数リスト改0.037s0.943s24.686s
エラトステネスの篩0.007s0.16s2.383s
ミラーラビン判定法0.024s0.559s7.402s
next_prime()0.119s 1.364s15.12s
決定ミラーラビン法1.371s23.535s

まあ上の方は放っておくとして、エラトステネスの篩のぶっちぎりの速さが際立ちました。

1000万くらいまでの素数を調べるのであればこれが一番手っ取り早い気がします。

1000万ビットくらいメモリがあれば足りるので、そこらのパソコンでは十分でしょう。ちなみにPHPの仕様なのかあんまり大きい配列を予約しようとするとダメでした…

あと、多分ちゃんと書けば自作の決定ミラーラビン法はこれより良くなるはずなので誰かコード書いてくださいお願いします。

PHPはコードを走らせるマシンのスペックにも依りますが、概ねこういう順番で速いコードということになりそうです。

次は別の言語で勝負してみたいですね。

シェアする