挑戦！100万番目の素数を求めろ

きっかけは@tkglingさんのツイートでした。

100万番目の素数を返すプログラム、どのくらいで計算できるだろう？

100万番目の素数を返すプログラム、どのくらいで計算できるだろう？

— nme (@tkgling) 2019年6月6日

100万番目の素数を返すコードをC++, Python, Javascriptの三つで対決させよう。

— nme (@tkgling) 2019年6月6日

100万番目の素数かぁ…

100万番目の素数かぁ…

— ウラル (@barley_ural) 2019年6月6日

こうかな？

GMP(GNU MP)関数を使ってみました。これは、PHPに標準で含まれるライブラリです。
PHPはお呼びでないって？うるせーーー！知らねーーー！（PHPの苦手分野だけどな！）

$prime = 1;for($i=0;$i<1000000;$i++){$prime = gmp_nextprime($prime);}echo $prime; // 100万番目の素数

— ウラル (@barley_ural) 2019年6月6日

<?php
$prime = 1;

for($i=0;$i<1000000;$i++){
$prime = gmp_nextprime($prime);
}

echo $prime;  //  100万番目の素数

gmp_nextprime()は、引数の次の素数を返す関数です。これなら、100万回繰り返すだけで求まりますね。

100万番目の素数の答え

上記のコードを実行します。

15485863

15485863！！

64秒かかった…

— ウラル (@barley_ural) 2019年6月6日

set_time_limit(0);を書かないとタイムアウトするくらいにはかかりますね。
PCのスペックと実行中の負荷によってもかなり変わります。僕のPCは
記事執筆中、3回実行して1番早かったのは72.1秒でした。

実コードで実装よろぴく

実コードで実装よろぴく。

— nme (@tkgling) 2019年6月6日

ですよねーーー。GMPライブラリは、標準で切られている環境が多いです。それにGMP関数使うの、ちょっとずるい。

gmp_nextprime()はどうやって素数を求めてるの？

ところで、さっきのgmp_nextprime()はどうやって素数を求めているんでしょうか。

公式リファレンスの一番下に、こう書かれています。

この関数は素数を識別するのに確率的アルゴリズムを使用します。 誤って合成数を取得してしまうことは、まずありません。

「確率的アルゴリズム」…？

さっきの「確率的アルゴリズム」とは何なのでしょうか。Wikipediaで調べてみました。

素数判定

ミラー–ラビン素数判定法

説明を読んでもサッパリですね。
ざっくり言うと、確率的アルゴリズムとは、「おそらく素数」を高速で求める方法らしいです。なぜ「おそらく」なのかと言うと、この方法では（多分）フェルマーの小定理を用いた合同式の対偶を利用して計算していて、その中で乱数を使用するからです。ごく稀に、合成数なのに奇数として出してしまうことがあるということですね。実用上は問題ないらしいです。

決定的アルゴリズムで考えよう

確率的アルゴリズムは難しいからやめましょう。地道に確実な方法で求めようよ、ね？
このような方法を確率的アルゴリズムと比較して「決定的アルゴリズム」と言います。

こうだっ！

ファーストトライとりあえず動いた… pic.twitter.com/h0P6m6m4BL

— ウラル (@barley_ural) 2019年6月7日

少し考えて作ってみました。

<?php

set_time_limit(0);
$debug_start = microtime(true);

$a = 2;//~無限
$primes = [2];

for($cnt=0;$cnt<10000;$cnt++){
    //素数が見つかるまでの処理

    $flg = false;

    while($flg == false){
        $a++;
        $flg = true;

        for($i=0;$i<=$cnt;$i++){
            if($a % $primes[$i] == 0){
                $flg = false;
                break;
            }
        }
    }
    $primes[$cnt + 1] = $a;
}

echo $primes[$cnt - 1] . '<br>' . (microtime(true) - $debug_start) . 's';

とりあえず1万まで。

最初に、$debug_start = microtime(true);で開始時間を記録します。
$a = 2;は、素数かどうか判定される自然数です。これを増やしていきます。
$primes = [2];は素数リストです。素数を見つけたらここに追加していきます。
for($cnt=0;$cnt<10000;$cnt++)で、素数リストが1万個になるまで素数判定を繰り返します。
ここからは素数判定の方法です。
$flg = false;で最初に$aが素数かどうかのフラグを偽にします。
while($flg == false)で、「$aが素数でなければ繰り返す」ようにします。
$a++;します。これは単に順番の問題でここにあります。
これから合成数かどうか判定するので、$flg = true;で最初にこれが素数であると仮定します。
for($i=0;$i<=$cnt;$i++)は、素数リストの手前から順番に取り出すための繰り返しです。
if($a % $primes[$i] == 0)で、「もし$aを素数リストから取り出した数で割った余りが0」になったら、$aは合成数だとわかります。
合成数なので$flg = false;で素数フラグは偽です。
素数リストを辿る繰り返しを抜け、$a+1して次のループに入ります。
$aがそれまで発見したどの素数でも割れなかったら、それは素数です。
$flg = true;だと仮定されたまま、while($flg == false)に入らずに進みます。
$primes[$cnt + 1] = $a;で素数リストに新たにその数が追加されます。
1万回繰り返したら、出力します。

この方法なら、奇数だけに絞る意味はあまりありません。なぜなら、素数リストの1番目が2なので一瞬で合成数だと分かるからです。

これで100万回やってみようとした

これワンチャン30分くらいかかるな　やべえ

— ウラル (@barley_ural) 2019年6月7日

(´・ω・｀)諦めた

平方根まででいいのでは…？

ちょっと待ってよ。素数か判定するために素数リストを全部当てはめて確かめたけど、合成数か確かめるには「平方根×平方根」までで良くない？
「11」が素数か判定したいとき、5以上の素数で割れないのは明らかだよね。
ということで

平方根までにするだけでめちゃくちゃ早くなった pic.twitter.com/ZqWbRMtKDz

— ウラル (@barley_ural) 2019年6月7日

これで100万番目の素数は求まるだろ！

それでも100万番目の素数を出すのにうちのPCだと206秒かかる

— ウラル (@barley_ural) 2019年6月7日

求まった。でも206秒かかる…！

for文の中に関数を入れると毎回呼び出すので（しかもsqrt()は遅い）、外で一回変数定義しておくか$prime[$k]**みたいに平方根でなくて平方を使うのがオススメです。

— nme (@tkgling) 2019年6月7日

forループの中にいるsqrt($a)を変数に入れるとちょっと早くなるんじゃないかな

— もぎゃ🔌(daisuke furukawa) (@mogya) 2019年6月7日

for文は毎回判定されるので、確かに関数入れたら毎回実行されちゃう…！

直した

平方根を出す処理をfor文の外に出しました。
ついでにインクリメントを前置きにすると早いと聞いたのでそこも直しました。（実測すると意外と逆だったりするのでおまじない程度です）

アドバイスを受けて pic.twitter.com/1ySLeLUcOI

— ウラル (@barley_ural) 2019年6月8日

<?php

set_time_limit(0);
$debug_start = microtime(true);

$a = 2;//~無限
$primes = [2];

for($cnt=0;$cnt<1000000;++$cnt){
    //素数が見つかるまでの処理

    $flg = false;

    while($flg == false){
        ++$a;
        $flg = true;
        $sqrted = sqrt($a);
        for($i=0;$primes[$i]<=$sqrted;++$i){
            if($a % $primes[$i] == 0){
                $flg = false;
                break;
            }
        }
    }
    $primes[$cnt + 1] = $a;
}

echo $primes[$cnt - 1] . "\n" . (microtime(true) - $debug_start) . 's';

85.7秒！めちゃくちゃ早くなった！！

これにてウラルはターンエンド！！

（補足：記事執筆中、3回実行して1番早かったのは90.6秒でした）

他の人の実装紹介

まずはネタ元@tkglingさん！
記事にされていたので、ここで書くよりずっと分かりやすいので読んで下さい。

[PHP] 1,000,000番目の素数を求める
https://tkgstrator.work/?p=16234

細かい部分は違いますが、同様に素数リストを使った方法でした。
記事執筆中、素数リストの最終形を3回実行して1番早かったのは115.4秒でした。

エラトステネスの篩やライブラリを用いたミラーラビン素数判定法のほか、決定的ミラーラビン判定法についても書かれてる…！！

次に、組み込み系のくわい(@21692fa)さん！

100万番目の素数計算ページ100万番目の素数と計算時間が出るよ（※読み込みに数秒かかります）https://t.co/dOdg754j0a

— くわい👎@赤眼鏡 (@21692fa) 2019年6月8日

https://kunow.ml/PrimeNum.html

jsはっやい。うちのPCだと14.5秒でした。iPhoneだと数秒です。

嫌々PHPでも書いてくれました。申し訳ない…

<?php
set_time_limit(0);

$PrimeNumbers  = [2];
$PrimeNumCount = 1;
$CurrentNumber = 3;
$RoundUpSqRoot = 2;
$RoundUpSqNum  = $RoundUpSqRoot * $RoundUpSqRoot;
$IsPrimeNumber = true;

$ExecuteTime = microtime(true);

for (; $PrimeNumCount < 1000000; $CurrentNumber += 2) {
  if ($CurrentNumber > $RoundUpSqNum) {
    $RoundUpSqRoot++;
    $RoundUpSqNum = $RoundUpSqRoot * $RoundUpSqRoot;
  }

  for ($i = 0; $i < $PrimeNumCount; ++$i) {
    if ($CurrentNumber % $PrimeNumbers[$i] === 0) {
      $IsPrimeNumber = false;
      break;
    }
    if ($PrimeNumbers[$i] > $RoundUpSqRoot) {
      $IsPrimeNumber = true;
      break;
    }
  }

  if ($IsPrimeNumber) {
    $PrimeNumCount++;
    $PrimeNumbers[] = $CurrentNumber;
  }
}

$ExecuteTime = microtime(true) - $ExecuteTime;

echo $PrimeNumbers[$PrimeNumCount - 1] . "\n" . $ExecuteTime . "s";

記事執筆中、3回実行して1番早かったのは99.1秒でした。
関数も何も使ってないのに思ったほど早くない…という不思議な結果に。PHPは組み込み系言語の常識が通用しないと聞いたことあるけど、それなんだろか…。

ということで、楽しかったです。定期的にこういうのやりたい。