0. 前回までのあらすじ
自作の素因数分解プログラムを Qiita で公開した。プログラムに添削が入り，30倍以上速くなった。ちなみに教えて頂いた方法は「試し割り法」と呼ばれている。教えられた方法をそのまま使うのもアレなので，もっと速い方法を編み出せないか考えてみた。
現代のプロセッサでも割り算は比較的重いとされるので，与えられた自然数に対してエラトステネスの篩（ふるい）と呼ばれる手法を用いて 以下の素数を列挙し，それらの素数だけで試し割りをすれば除算回数は最小化され，格段の高速化が期待できるはずだ！しかし，エラトステネスの篩（ふるい）は大量にメモリを使用するせいか，意図したほど速くならなかった。
詳しくは下記の記事を参照されたい。
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件 - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(3) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(4) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(5) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(6) - Qiita
なので，今回は区分篩（ふるい）に挑戦する。
1. エラトステネスの篩（ふるい）に必要なメモリ容量
JavaScript で扱える安全な整数の最大値 $2^{53} - 1 = 9,007,199,254,740,991$ まで素数判定を行うためには，$2$ から最大で $94,906,249$（$2^{53} - 1$ の平方根以下の最大素数）までの全ての素数で割り切れないことを示さなくてはならない。
このため $2$ から $94,906,249$ までの素数リストをエラトステネスの篩（ふるい）を用いて作成するとしよう。エラトステネスの篩（ふるい）は素数判定のためのテーブルを用意する必要があり，大量のメモリを必要とする。$94,906,249$ までの自然数を格納するとして，以下のメモリ（単位はバイト）を必要とする。素数の判定値 true/false を１バイトで表した場合と，１ビットで表した場合を示す。
| 格納対象 | バイト単位 | ビット単位 |
|---|---|---|
| 自然数を格納する場合 | $94,906,250$ | $11,863,282$ |
| 奇数を格納する場合 | $47,453,125$ | $5,931,641$ |
| ２・３の倍数を除く場合 | $31,635,418$ | $3,954,428$ |
| ２・３・５の倍数を除く場合 | $25,308,334$ | $3,163,542$ |
近年のインテルプロセッサはキャッシュ容量を増やしているが，それでも L2 キャッシュに収まり切らないことが分かる。
2. Segmented Sieve 区分篩（ふるい）
この問題自体は古くから知られており，対策として Segmented Sieve，日本語だと区分篩（ふるい）という手法が知られている。残念ながら，入門者向けの平易な日本語サイトは見当たらなかった。強いて挙げるなら脚注1のサイトであろうか？
2.1. 素数リスト
$2$ から $94,906,249$ までの素数判定テーブルを区分篩（ふるい）を用いて作成するとしよう。まず，区分篩（ふるい）を行うためには $94,906,249$ の平方根以下の最大素数 $9,739$ 以下の素数リストを必要とする。$9,739$ 以下の素数は $1,201$ 個あり，これらは符号なし 16bit 整数で収まるので $2,402$ バイト必要となる。
2.2. 素数判定テーブル
$9,739$ 以下の領域は篩（ふるい）にかける必要がないので，$9,740$ から $9,740$ 個ずつ区分篩（ふるい）にかけるものとする。true/false を１バイトで表すとして $9,740$ バイト必要となる。
2.3. 篩（ふるい）をかける初期位置
通常のエラトステネスの篩（ふるい）の場合，素数 $p$ に対して $p^2$ 以上の $p$ の倍数を合成数として消していく。すなわち初期位置は $p^2$ である。一方，区分篩（ふるい）では $p^2 \ge 9740$ の場合と $p^2 < 9740$ の場合に分けて考える必要がある。
- $p^2 \ge 9740$ の場合，初期位置は $p^2$ で良い。
- $p^2 < 9740$ の場合，初期位置は $\lceil 9740/p \rceil \cdot p$ となる。
初期位置は符号なし 32bit 整数で収まるので $4 \times 1,201 = 4,804$ バイト必要となる。
初期位置を $p^2$ の大小に関わらず常に $9,740$ 以上の $p$ の倍数とする方法もある。この場合，初期位置は $\lceil 9740/p \rceil \cdot p$ となるため，符号なし 16bit 整数で収まる。こうすると必要なメモリ容量は半減する一方，篩（ふるい）をかける回数が僅かに（約2.4%）増加することになる。実際に両者を試作して比較してみたところ，$p^2 \ge 9740$ のときは $p^2$ から始めたほうが僅かに速いという結果が得られた。
2.4. 必要なメモリ容量
これで合計 $16,946$ バイト必要となるが，大半のプロセッサの L1 データキャッシュに収まるはずだ。
| 内容 | 型 | 個数 | 容量 |
|---|---|---|---|
| 素数リスト | Uint16Array | $1,201$ | $2,402$ |
| 初回位置 | Uint32Array | $1,201$ | $4,804$ |
| 素数判定テーブル | Uint8Array | $9,740$ | $9,740$ |
| 合計 | $16,946$ |
3. 区分篩（ふるい）の実装コード
区分篩（ふるい）の実装コードを以下に示す。
配列の要素数に直値 immediate number を使用するという極めて行儀の悪いスタイルだが，パファーマンスを上げるために敢えてこうしたと理解されたい。
配列 start[i] には，list1[i] の示す素数がマーキングを行う初期位置が格納されている。この初期位置は素数判定テーブル table[] における相対位置，インデクスを示すようにした。このため初期位置 $x$ とおくと
- $p^2 \ge 9740$ のとき $x = p^2 - 9740$
- $p^2 < 9740$ のとき
- $9740 \bmod p = 0$ のとき $x = 0$
- $9740 \bmod p \ne 0$ のとき $x = \lceil 9740 / p \rceil \times p - 9740$
となる。ここで $9740 \bmod p \ne 0$ のとき
\begin{aligned}
x &= \lceil 9740 / p \rceil \times p - 9740 \\
&= \left( \lfloor 9740 / p \rfloor + 1 \right) \times p - 9740 \\
&= p - \left( 9740 - \lfloor 9740 / p \rfloor \times p \right) \\
&= p - (9740 \bmod p)
\end{aligned}
と変形できる。通常 JavaScript では数値を倍精度浮動小数点として扱うが，素数 $p$ も TypedArray に格納された整数値であることから，整数剰余命令を使うことで JIT コンパイラが整数演算に変換してくれることを期待しての式変形である。
なお配列 start[] および table[] は，JavaScriptの TypedArray が初期化時に必ずゼロクリアすることを利用している。table[] は区分を変えるたびにゼロクリアする必要があるが，まとめてゼロクリアするよりも，読んだ直後に同じアドレスをクリアしたほうが速い。
//------------------------------------------------------------------------------
// 9,740 以下の素数リスト ※要素数 1,201個
//------------------------------------------------------------------------------
var list1 = new Uint16Array([
2,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,
103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,
/* 中略 */
9547,9551,9587,9601,9613,9619,9623,9629,9631,9643,9649,9661,9677,9679,9689,
9697,9719,9721,9733,9739]);
//------------------------------------------------------------------------------
// 与えられた自然数 n の平方根 √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list(n) {
var nsqrt = Math.floor(Math.sqrt(n));
var list2 = []; // 得られた素数を格納するリスト
var start = new Uint32Array(1201); // 初期位置
var table = new Uint8Array(9740); // 区分篩（ふるい）を行う素数判定テーブル
//--------------------------------------------------------------------------
// 初期位置の計算
//--------------------------------------------------------------------------
for(var i = 0; i < 1201; i++) {
var p = list1[i];
list2.push(p);
var r, s;
if((r = p * p - 9740) >= 0)
start[i] = r;
else if((s = 9740 % p) != 0)
start[i] = p - s;
}
//--------------------------------------------------------------------------
// 区分篩（ふるい）の実行
//--------------------------------------------------------------------------
for(var offset = 9740; offset < nsqrt; offset += 9740) {
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; // 次回位置を更新
}
for(var j = 0; j < 9740; j++) {
if(table[j] == 0) list2.push(offset + j);
table[j] = 0; // 次回に備えてゼロクリア
}
}
return list2;
}
4. パフォーマンス比較
JavaScript において安全に使用できる整数最大値 Number.MAX_SAFE_INTEGER 以下の最大素数 $9,007,199,254,740,881$ を入力し，素数判定に要した時間を計測した。なお JavaScript エンジンは Chakra エンジンを使用した。TypedArray を使用してるので，Internet Explorer の JavaScript エンジンでは動作しない。
区分篩（ふるい）によるパフォーマンス改善効果は著しく，これまでのエラトステネスの篩（ふるい）に対して，あらゆるプロセッサでほぼ倍速を達成している。試し割り法にはまだ及ばないものの，これまでと同様に奇数探索，２・３の倍数除去などを行っていけば試し割り法を超える見込みがある。
5. 次回
次回に続きます。
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(8) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(9) 最終回 - Qiita