0. 前回までのあらすじ
詳しくは下記の記事を参照されたい。
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件 - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2) - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(3) - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(4) - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(5) - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(6) - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(7) - Qiita
1. 区分篩（ふるい）のさらなる高速化
前回の記事より，区分篩（ふるい）自然探索版のコア部分を抽出すると以下のようになる。素数リスト list1[] および 素数判定テーブル table[] は L1 データキャッシュ内に十分収まるはずである。
for(var i = 0; i < 1201; i++) {
var p = list1[i];
for(var j = start[i]; j < 9740; j += p)
table[j] = 1;
start[i] = j - 9740;
}
一方，奇数探索版は以下のようになる。素数リスト list1[] は $2$ が抜けるのでサイズが一つ減る。また素数判定テーブル table[] のサイズが半分になることから，倍速化が期待できる。
for(var i = 0; i < 1200; i++) {
var p = list1[i];
for(var j = start[i]; j < 4870; j += p)
table[j] = 1;
start[i] = j - 4870;
}
2. 奇数探索版の設計
2.1 素数リストの定義
$9739$ 以下の奇数の素数は $1200$ 個ある。素数 $p$ の素数判定テーブルおける位置を $q$ とする。
| $i$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $\cdots$ | $1199$ |
|---|---|---|---|---|---|---|---|---|---|---|
| $q$ | $1$ | $2$ | $3$ | $5$ | $6$ | $8$ | $9$ | $11$ | $\cdots$ | $4869$ |
| $p$ | $3$ | $5$ | $7$ | $11$ | $13$ | $17$ | $19$ | $23$ | $\cdots$ | $9739$ |
素数の値 $p$ と位置 $q$ は相互変換可能である。
$$
q = \frac{p - 1}{2} \iff p = 2 q + 1
$$
2.2 素数判定テーブル（奇数のみ）の定義
素数判定テーブル（奇数のみ）を定義する。位置 $j$ における値 $n$ とおく。
| $j$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $\cdots$ | $3246$ | $4869$ |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $n$ | $1$ | $3$ | $5$ | $7$ | $9$ | $11$ | $13$ | $15$ | $\cdots$ | $9737$ | $9739$ |
素数判定テーブルの値 $n$ と位置 $j$ も相互変換可能である。
j = \frac{n - 1}{2} \iff n = 2 j + 1
2.3 初回位置
通常のエラトステネスの篩（ふるい）では $p^2$ 以上の $p$ の倍数を合成数として消していく。すなわち初回位置は $p^2$ である。一方，今回の区分篩（ふるい）の場合，$9739$ 以下の素数リストは既に得られており，初回位置は $9741$ 以上の $p$ の倍数となる。すなわち $p^2$ と $9741$ の大小関係によって場合分けが必要である。
区分篩（ふるい）を行う最初の値 $n_0$，素数判定テーブルにおける位置 $j_0$ とおく。
-
$p^2 \ge 9741$ のとき（$\sqrt{9741} = 98.7 \cdots$ なので等号が成立することはない）
\begin{aligned} n_0 &= p^2 = (2q + 1)^2 = 4q^2 + 4q + 1 = 2 \cdot \left\lbrace 2q \cdot (q + 1) \right\rbrace + 1 \\ &= 2j_0 + 1 \\ j_0 &= 2q \cdot (q + 1) \end{aligned} -
$p^2 < 9741$ のとき，$(2m + 1) \cdot p \ge 9741$ となる最小の自然数 $m$ とおく。
\begin{aligned} (2m + 1) \cdot p &= 2mp + p = 2mp + 2q + 1 \ge 9741 \\ mp + q &\ge 4870 \end{aligned}より
m = \left\lceil \frac{4870 - q}{p} \right\rceilとなる。なお $\lceil x \rceil $ は $x$ 以上の最小の整数という意味である。※いわゆる天井関数
\begin{aligned} n_0 &= (2m + 1) \cdot p \\ &= (2m + 1) \cdot (2q + 1) \\ &= 2\cdot(2mq + m + q) + 1 \\ &= 2j_0 + 1 \\ j_0 &= 2mq + m + q \\ \end{aligned}- $(4870 - q) \bmod p = 0$ のとき，$(2m + 1) \cdot p = 9741$ となるので
\begin{aligned} n_0 &= (2m + 1) \cdot p \\ &= 9741 \\ &= 2j_0 + 1 \\ j_0 &= 4870 \end{aligned} - $(4870 - q) \bmod p \ne 0$ のとき
なお，$\lfloor x \rfloor$ は $x$ 以下の最大の整数という意味である。※いわゆる床関数
\begin{aligned} n_0 &= (2m + 1 ) \cdot p \\ &= 2mp + 2q + 1 \\ &= 2 \left\lbrace \left\lceil \frac{4870 - q}{p} \right\rceil \cdot p + q \right\rbrace + 1 \\ &= 2j_0 + 1 \\ j_0 &= \left\lceil \frac{4870 - q}{p} \right\rceil \cdot p + q \\ &= \left\lbrace \left\lfloor \frac{4870 - q}{p} \right\rfloor + 1 \right\rbrace \cdot p + q \\ &= (4870 - q) - \lbrace (4870 - q) \bmod p \rbrace + p + q \\ &= p + 4870 + \lbrace (4870 - q) \bmod p \rbrace \\ \end{aligned}
- $(4870 - q) \bmod p = 0$ のとき，$(2m + 1) \cdot p = 9741$ となるので
2.4 移動幅
区分篩（ふるい）を行う $k$ 番目の値 $n_k$ および位置 $j_k$ とする。$n_k$ は素数 $p$ の奇数倍であるから，適当な自然数 $x$ を用いて
\begin{aligned}
n_k &= (2x + 1) \cdot p \\
&= (2x + 1) \cdot (2q + 1) \\
&= 2\cdot(2xq + x + q) + 1 \\
&= 2j_k + 1 \\
j_k &= 2xq + x + q
\end{aligned}
と表せる。一方，$k + 1$ 番目の値 $n_{k+1}$ および位置 $j_{k+1}$ は
\begin{aligned}
n_{k+1} &= \left\lbrace 2(x + 1) + 1\right\rbrace \cdot p \\
&= (2x + 3) \cdot (2q + 1) \\
&= 2\cdot(2xq + x + 3q + 1) + 1 \\
&= 2j_{k+1} + 1 \\
j_{k+1} &= 2xq + x + 3q + 1 \\
\end{aligned}
となる。これより素数判定テーブルにおける位置の移動幅 $\Delta j$ は $k$ に依存せず
\Delta j = j_{k+1} - j_k = 2q + 1 = p
一定となる。
3. 実装コード
実装コードを以下に示す。
//------------------------------------------------------------------------------
// 3 以上 9,740 以下の素数リスト ※要素数 1,200個
//------------------------------------------------------------------------------
var list1 = new Uint16Array([
3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,
/* 中略 */
9697,9719,9721,9733,9739]);
//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list(n) {
var nsqrt = Math.floor(Math.sqrt(n));
var list2 = []; // 得られた素数を格納するリスト
var start = new Uint32Array(1200); // 初期位置
var table = new Uint8Array(4870); // 区分篩（ふるい）を行う素数判定テーブル
//--------------------------------------------------------------------------
// 初期位置の計算
//--------------------------------------------------------------------------
list2.push(2);
for(var i = 0; i < 1200; i++) {
var p = list1[i];
list2.push(p);
var q = (p - 1) >>> 1;
var r, s;
if((r = 2 * q * (q + 1) - 4870) >= 0)
start[i] = r;
else if((s = (4870 - q) % p) != 0)
start[i] = p - s;
}
//--------------------------------------------------------------------------
// 区分篩（ふるい）の実行
//--------------------------------------------------------------------------
for(var offset = 9740; offset < nsqrt; offset += 9740) {
for(var i = 0; i < 1200; i++) {
var p = list1[i];
for(var j = start[i]; j < 4870; j += p)
table[j] = 1;
start[i] = j - 4870;
}
for(var j = 0; j < 4870; j++) {
if(table[j] == 0) list2.push(offset + 2 * j + 1);
table[j] = 0;
}
}
return list2;
}
4. 結果
JavaScript において安全に使用できる整数最大値 Number.MAX_SAFE_INTEGER 以下の最大素数を入力し，素数判定に要した時間を計測した。なお JavaScript エンジンは Chakra エンジンを使用した。TypedArray を使用してるので，Internet Explorer の JavaScript エンジンでは動作しない。
残念ながら自然数探索版の２倍速とはならなかったが，どの世代のプロセッサでも著しい改善効果が見られた。試し割り法に迫る勢いである。あとちょいでアレである。
5. 次回
次回に続きます・・・