1本目 からの続きになります。その他の高速化は
- エラトステネスの篩の高速化 (1)
- エラトステネスの篩の高速化 (2) ← 今ココ
- エラトステネスの篩の高速化 (3)
- エラトステネスの篩の高速化 (4)
- エラトステネスの篩の高速化 (5)
- エラトステネスの篩の高速化 (6)
辺りを参照してください。
3,5 の倍数の省略
前回 2 の倍数を省いて倍速化したように、3 の倍数を省けば 3/2 倍、5 の倍数を省けば 5/4 倍の高速化が期待されます。が、それ以上に 5 の倍数まで除くと候補として残る数が
30k+\{1,7,11,13,17,19,23,29\}
と全体の 8/30 になるので、コンピュータで扱うことがとても便利になります。具体的には $k$ バイト目の $i$ ビット目を $30k+m_i$ と対応させることができます。なお $\{m_i\}$ は上記の式の {} の中身に対応します。
数とビットとの対応が簡単になったので、今度は篩う位置の計算をどのようにするか検討しましょう。
篩の初期位置
まず素数 $p=30k+m$ で篩う初期ビットは、$p^2$ の数を求める式を展開することで
\begin{eqnarray}
(30k+m)^2 &=& 30^2k^2 + 2\cdot30km + m^2\\
&=& 30\times \left\{k(30k+2m)+\left\lfloor\frac{m^2}{30}\right\rfloor\right\} + (m^2 \bmod 30)
\end{eqnarray}
に対応するビットだということがわかります。つまり $k(30k+2m)+\lfloor m^2/30\rfloor$ バイト目の $m^2\bmod 30$ に対応するビットです。前者についてはそのまま計算することになりますが、後者については $m$ から求められる表を用意しましょう。ただ、素直に考えると 8 要素の配列で大丈夫なのですが、後の再利用を兼ねて 8x8 の表を作っておくと楽です。詳細は後述。
篩う位置の移動
ついで移動処理についてですが、これについては
- 今どの素数 ($p=30i_1+m_{j_1}$) で篩っているのか
- 今どの数 ($q=30i_2+m_{j_2}$) を掛けた合成数を篩ったのか
の情報が必要になります。まずはバイト位置を計算しましょう。ちょっと面倒ですが、
\begin{eqnarray}
&& \left\lfloor\frac{1}{30}(30i_1+m_{j_1})(30i_2+m_{j_2+1})\right\rfloor
- \left\lfloor\frac{1}{30}(30i_1+m_{j_1})(30i_2+m_{j_2})\right\rfloor \\
&=&i_1(m_{j_2+1}-m_{j_2}) + \left\lfloor\frac{m_{j_1}m_{j_2+1}}{30}\right\rfloor - \left\lfloor\frac{m_{j_1}m_{j_2}}{30}\right\rfloor
\end{eqnarray}
となります。$j_2$ が定まれば $m_{j_2}$ も $m_{j_2+1}$ も決まるので、要素 8 の配列と、8x8 の 2 次元配列を用意すれば効率的に計算できることがわかります。
ビット位置については初期位置で使いまわそうとした 8x8 のテーブルが便利です。このテーブルは要するに $m_im_j \bmod 30$ を表せばいいだけのテーブルですが、計算処理のことを考えると $m_im_j \bmod 30$ の数を直接持つよりも該当ビットを外すビットマスクにする方が効率的です。
以上を踏まえた今回のコードは github の Version 2 に該当します。
計算時間
目安として、計算時間を記録しておきます。単位は秒。
| $x$ | $\pi(x)$ | Ver. 0 | Ver. 1 | Ver. 2 |
|---|---|---|---|---|
| $10^7$ | 664579 | 0.047 | 0.017 | 0.007 |
| $10^8$ | 5761455 | 1.039 | 0.352 | 0.189 |
| $10^9$ | 50847534 | - | 6.180 | 3.377 |
補足：ビット処理
立っている最下位ビットを取得
n & (-n)
立っている最下位ビットを削除
n & (n - 1)
ビットカウント（実装しなくても専用命令が多くの CPU の場合あります）
uint64 PopCnt(uint64 n) {
n = (n & 0x5555555555555555ULL) + ((n >> 1) & 0x5555555555555555ULL);
n = (n & 0x3333333333333333ULL) + ((n >> 2) & 0x3333333333333333ULL);
n = (n & 0x0f0f0f0f0f0f0f0fULL) + ((n >> 4) & 0x0f0f0f0f0f0f0f0fULL);
return (n * 0x0101010101010101ULL) >> (64 - 8);
}