どうも、吟遊詩人です。今回も素数探索、やっていきますよ。
今回で吟遊詩人の素数探索第一部完となります。
前回gmpyを用いて、リュカ・レーマーテストを高速化し、35番目のメルセンヌ素数p=1,398,269を素数判定しました。
(前回：吟遊詩人の素数探索その4〜gmpy使ったらメッチャ速くなった話〜)
今回はさらなる高速化を目指して剰余算をビット演算化します。（なんか改良みたいに書いていますが、本来ビット演算ができるからこそのリュカ・レーマーテストなんですけどね。）
リュカ–レーマー・テストは二進計算機用のアルゴリズムに向いており、コンピュータによるメルセンヌ素数の発見には、この判定法が用いられてきた。例えば、2p ≡ 1 (mod Mp) より、A·2p + B ≡ A + B (mod Mp) が成り立つので、Mp で割る割り算の代わりに、二進法で p 桁のシフト演算と足し算だけで計算できる。
アルゴリズムは以下の擬似コードで表される。
入力: p:奇素数であるテスト対象の整数
出力: PRIME:素数の場合, COMPOSIT:合成数の場合
Lucas_Lehmer_Test_FAST(p):
var s = 4
var M = 2p − 1
for n in range(2, p):
var sqrt = s × s
s = (sqrt & m) + (sqrt >> p)
if s >= m then
s = s − m
s = s − 2
if s == 0 then
return PRIME
else
return COMPOSIT
出典: フリー百科事典『ウィキペディア（Wikipedia）』：メルセンヌ数
という感じです。
上位ビットと下位ビットの足し算として表すことができるということですが、理解のため例をあげてみましょう。
Mn＝１２７のときに出てくる4489(=67^2) mod 127を例に考えます。
4489と127とついでに128を２進数で表すと以下のとおりです。
ここで128 mod 127は1です。127は7ビットなので、4489を下位７ビットと上位ビットに分けます。
ここで、上位ビット（赤）と下位ビット（青）での余りをそれぞれ足していくのですが、下位ビットの方は127(1111111)より少ないので、これがそのまま余りになります。
プログラムでは下位７ビットを抽出するために、1111111とのANDを取ります。
次に上位ビットの方を見ていきます。上位ビットは1000110000000ですが、まずこれを128で割ると、８ビット以上の部分の数になります。
当たり前といえば当たり前ですが、次に127で割ると、先程128で割った答えと同じ余りが出てきます。
128 mod 127 = 1なので、128で割った数に対して１ずつ余りが生じるからなんですね。（わかりにくい言い方かな？）
ですので、上位ビット（赤の部分）の剰余は、結局もとの数を7ビット右にシフトした数になります。下位ビットと上位ビットを足すと、
となり、4489%127=44が正しく導けていることがわかります。
(今回は足した余りが127以下でしたが、１２７以上になった場合、１２７をひいた数が剰余となります。)
結局、Mn=2^p-1の形なので、p桁より上のビットの剰余がビットシフトで持ってこれるのが大きいんですね。
では、wikiのコードをそのままコピーして実装しませう。
実行結果はコチラ↓
赤いドットが前回で、今回ビット演算を使ったのが緑のドットです。
だいぶ早くなりましたね。
最終的に８時間ほどで36番目のメルセンヌ素数（p=2,976,221, 1997年GIMPS / Gordon Spenceによって発見）まで素数判定できるようになりました。
この曲線によれば、現在最大の素数（p = 82,589,933, 2018年12月にGIMPS / Patrick Larocheによって発見）も半年ほどで素数判定できることになります。
ただし、吟遊詩人予想の最小のp=2147483647だと、約５００年かかります。
おわりに
さて、これまでリュカ・レーマーテストの高速化を行ってきましたが、とりあえず今の所できるのはこの程度といったところです。
PCのスペックを上げていけばもう少し早くなるかもしれませんが、数１００倍の高速化は難しいところです。
マルチスレッド化ができれば良いのですが、いろいろ考えてみましたが難しそうです。ループのところが前の計算結果を用いるため、並列処理化の良いアイデアが浮かびません。計算結果があとから反映されるような仕組みでもあればいいんですが（それって量子コンピューター）
あとはアルゴリズム自体を変えるしか無いですが、リュカ・レーマーテストより高速なアルゴリズムを作るのは至難の業だと思います。
そのためには理論をもっと勉強しないと。（するとは言っていない。）
気長にやっていきましょう。（でも早く早期リタイアしたい。）
ひとまず吟遊詩人の素数探索第一部完となります。
Cを使ったリュカレーマーとかは気が向いたらやります。
ここまでご視聴ありがとうございました。
またの機会に会いましょう。
では。
第１部 完
そーめんブログ 