hnwの日記

PHPで素数を数えて落ち着いてみた

PHP

2,3,5,7,11,13,...と素数を順に列挙することで落ち着く人が世の中にはいるようです（参考：「素数を数えて落ち着くんだ…」）。とはいえ人力では素数を100個列挙するのさえつらいので、プログラミング言語に頼った方が落ち着けるはずです。PHPには、そんな状況で使えそうな関数が存在します。

gmp_nextprime ― 次の素数を見つける


PHP: gmp_nextprime - Manual


よし、この関数さえあれば落ち着けるぞ、と思いきや、マニュアルにはこんな記述もあります。

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


PHP: gmp_nextprime - Manual


えっ？「まずありません」ってことは少しくらいはあるんでしょうか。逆に不安で落ち着かなくなってしまいそうです。


本稿ではこの関数について詳しく調べてみます。

GNU MPと多倍長整数演算について

まず、gmp_nextprime関数の周辺知識を解説します。


この関数はGNU MPという多倍長整数演算ライブラリのPHPバインディングが提供する関数群の一つです。PHPではGMP 関数として分類されています。


多倍長整数演算とは、整数型の配列を一つの巨大な整数とみなして演算を行うことで、CPUの持つ整数レジスタよりも大きな整数を扱う技術です。ご存じの通り、PHPの整数型はCPUの整数レジスタと同じサイズであり、整数演算の結果がその上限を超えると勝手に浮動小数点数になってしまうという特徴があります。一方、これらのGMP関数では多倍長整数を扱うためのリソース型（PHPマニュアル上ではGMP数と書かれています）を提供しており、加減乗除やその他の演算について上限なしに多倍長整数演算を行うことができます。


GMP関数は多くのディストリビューションに含まれるPHPで有効になっているはずですが、環境によってはphp5-gmpのようなパッケージを追加でインストールしたり、自分でPHPをビルドしなおす必要があるかもしれません。


PHPにおけるGMP関数の知名度は現状でそれほど高くないような印象があります。とはいえ、PHP 5.6からはGMP数に対して加減乗除などの演算子が使えるようになるというアナウンスも出ていますので、今後はより身近な存在になってくるのではないでしょうか。（参考：「New features - GMP supports operator overloading」）

PHPのgmp_nextprime関数は何をしているのか

素数を求める関数の話題に戻りましょう。PHPのgmp_nextprime関数はGNU MPのmpz_nextprime関数に対応するものです。


このmpz_nextprime関数の説明にも、PHPのドキュメントと同様に「確率的アルゴリズムを使っている」「合成数を返す可能性は極めて小さい」などと書いてあります。


そこでGNU MPの実装に当たってみました。GNU MP 6.0.0aのソースコード中のmpz/nextprime.cを確認したところ、mpz_nextprime関数は簡易的なエラトステネスのふるいと25パスのMiller-Rabin素数判定法を組み合わせたようなものだとわかりました。

Miller-Rabin素数判定法とは

上の説明に出てきたMiller-Rabin素数判定法というのは、確率的素数判定アルゴリズムの一つです。ざっくり説明すると次のようになります。

  • 素数判定したい奇数をnとする
  • n-1を2^k*d（kは自然数、dは奇数の自然数）の形にする
  • 2以上n未満の自然数aを1個ランダムに選ぶ
    • a^d mod n が1ならnは「素数っぽい」
    • a^d mod n, a^(2*d) mod n, ..., a^(2^(k-1)*d) mod n のいずれかがn-1ならnは「素数っぽい」
    • いずれでもなければnは合成数
  • 素数っぽいと言われ続ける場合、気が済むまで異なるaについて上のテストを行う


ここで異なるaを何回まで選ぶかが、上で説明した「Miller-Rabin素数判定法のパス数」です。Miller-Rabin素数判定法を1パス行うと1/4以下の確率で合成数を「素数っぽい」と間違えて判定してしまいます。そこで、異なるaについて繰り返し試すことで間違う確率を減らせるというわけです。


素数かどうかは既に決まっていることなのに間違うかもしれない方法を使うのは不思議な気もしますが、このMiller-Rabin素数判定法は計算量のオーダーが他の方法より低い点がメリットです。僕の感覚値ですが、自身の平方根までの全素数で割るような試し割りアルゴリズムと比べると、10進8桁付近で互角、10進16桁付近では10000倍程度の速度差が出てしまいます。一定以上の大きさの数の素数判定をする場合、Miller-Rabin素数判定法は有力な選択肢の一つです。


詳しくはWikipediaの説明「ミラー-ラビン素数判定法 - Wikipedia」なども参照してください。

GNU MPのMiller-Rabin素数判定法が素数判定に失敗するとき

ところで、Miller-Rabin素数判定法での25パスというのは多いのでしょうか、少ないのでしょうか。合成数を間違って素数だと言ってしまう可能性はどれくらいあるのでしょうか。


まず確率の話だけをすると、間違える確率は(1/4)^25以下ですから、1000兆分の1以下の確率です。正直なところ良いのか悪いのか全くわからない数字です。


別の視点からも見てみましょう。世の中にはGNU MPのMiller-Rabin素数判定法がどれくらい間違うか調べるのが趣味の人がいるようです。それによれば、GNU MPの15パスのMiller-Rabin素数判定法で10進26桁の合成数10627945015121687221861591を素数だと判定してしまう例が見つかっているようです（参考：「GNU GMP mpz_probab_prime_p Pseudoprimes」、「MPZ_SPSP's under GMP 5.0.1」）。


また、GNU MPの挙動そのものではないのですが、Miller-Rabin素数判定法のaとして小さい方から素数n個を選んだ場合に、合成数なのに素数だと言ってしまう最小の整数を調べた資料も見つけました（参考：「The Prime Glossary: pseudoprime」）。これによれば、aとして2から71までの20個の素数を利用した場合、10^36までの全ての素数について正しく素数判定できるようです。


これらの結果から考えると、GNU MPの25パスのMiller-Rabin素数判定法であれば、かなり堅めに考えても10進35桁くらいまでは素数判定を間違わなさそうです。


なお、本稿ではGNU MP 6.0.0aを元に記事を書いています。GNU MP 5.0.2以前ではmpz_nextprimeのMiller-Rabinは10パスでしたので、10進19桁あたりから間違う例が見つかり始めます。

PHPで素数を数える

長々説明してきましたが、10進35桁くらいまでならgmp_nextprime関数で素数が数えられそうだ、という結論を得ました。素数を2から表示し続けるサンプルプログラムを貼り付けておきますので、皆様お役立てください。

<?php
$n = "1";
while (true) {
    $n = gmp_nextprime($n);
    echo gmp_strval($n), "\n";
}


$nの初期値を大きくすれば、かなり大きい素数を数えて落ち着くこともできます。

まとめ