2019-06-10に更新

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

PHP

image.png

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

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


100万番目の素数かぁ…


こうかな?


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

<?php
$prime = 1;

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

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

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

100万番目の素数の答え


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

15485863

15485863!!

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

実コードで実装よろぴく


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

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


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

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

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

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


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

素数判定
image.png

ミラー–ラビン素数判定法
image.png
image.png

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

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


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

こうだっ!


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

<?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万回やってみようとした


(´・ω・`)諦めた

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


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

image

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


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

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

直した


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

<?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)さん!

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は組み込み系言語の常識が通用しないと聞いたことあるけど、それなんだろか…。

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


ウラル

Splatoonの二次創作サイト「スプランプ」の管理人です。サーモンラン研究所やオクトチャット、フェス速報などを作りました。

Crieitは個人で開発中です。 興味がある方は是非記事の投稿をお願いします! どんな軽い内容でも嬉しいです。
なぜCrieitを作ろうと思ったか

また、「こんな記事が読みたいけど見つからない!」という方は是非記事投稿リクエストボードへ!

こじんまりと作業ログやメモ、進捗を書き残しておきたい方はボード機能をご利用ください!

ボードとは?

関連記事

コメント