まさかの 5 年ぶりに更新です。今の所もうネタが無い。
- エラトステネスの篩の高速化
- エラトステネスの篩の高速化 (2)
- エラトステネスの篩の高速化 (3)
- エラトステネスの篩の高速化 (4)
- エラトステネスの篩の高速化 (5)
- エラトステネスの篩の高速化 (6)
- エラトステネスの篩の高速化 (7) ←今ココ
計算対象
今回は大きな数での篩を高速化する話になります。また、実装はすっ飛ばしてしまったのでやる気が出るまでは理論だけ。
パフォーマンス計測
パフォーマンス計測、といっても専用のツールを使ってどうこうする話ではありません。篩に使っている素数のうち、前回導入したループアンローリング部分をきちんと使っている(while 内を経由している)ものがどれだけあるでしょうか。ざっと見積もってみましょう。
素数 $p$ で篩っているとき、この展開したループを 1 回通過すると p バイト分 j が進むことになります。つまり区間篩のサイズ kSegmentSize が p より小さい場合(正確には k=0 への調整部分もあるのでもう少し多いケースで)は折角アンロールしたループを一切通過しないことが分かります。
サンプルコード内で設定している kSegmentSize = 2MB = $2^{21}$ バイトのことを考えると $p \leq \sqrt{x}$ で篩うので
$$
y > (2^{21})^2 = 2^{42} \sim 4\times 10^{12}
$$
では多くの素数でアンロール部分を通過しなくなります。 (割合いとしては $x$ が大きい方が素数は少ないが、絶対数としては多くなるので最大素数の大きさが計算時間において支配的になる)
これだけなら前回の準備が無駄になるだけなので構わないですが、もっと大きな、例えば $y>10^{13}$ での計算では区間篩の区間内で 1 ビットすら落としていかない素数が増えていくことになります。
対策
その様な素数のことを考えると、素数をイテレートしていく作業すら無駄になります。なので「今から篩う区間内で1ビットでも落とす可能性のある素数だけリスト」を作ります。
具体的には、素数の篩う位置を示す情報を持っていたリスト (コード内では indecies) について、それぞれの素数で次に篩い落とすビットの位置が
$$
30j+m_k
$$
となるわけですが、これは次の $\lfloor \frac{j}{S}\rfloor$ 番目の区間で篩う (0番目がこれからすぐ篩う区間) 、逆に言えば $\lfloor \frac{j}{S}\rfloor$ 個の区間では計算処理をしなくて良い、ということを意味しています。
このことを利用して、篩う位置リストを分割することを考えましょう。次に篩うビットのバイト位置が区間 $[Si, S(i+1))$ にある素数の位置情報を位置情報リストのリスト、プログラム的には 2 次元配列 index_list の i 番目 index_list[i] にストックしておきます。ただメモリも無限に使えるわけではないので、2 次元配列の 1 次元目をリングバッファの様に $N$ 個用意して index_list[i % N] に置いておくのが現実的かもしれません。
注意点
今回の記事では実装が含まれていないので注意点として書いておきます。
今回の提案ではこれまで定数リストのように使ってきた indecies 配列(配列の中身は変わっていくが、配列のサイズや各要素と対応する素数は変化しなかった)をキューの配列として使うことになります。そのため、一旦区間を篩い終わった素数が次に篩う区間に基づいてどのキューに追加していくかの計算を間違えないようにしましょう。特にキューの数を $N$ 個に固定した場合は「本当に今この区間を篩う素数なのか、あと$n$周した後で使う素数なのか」をきちんと区別できるようにしましょう。