1. 前回までのあらすじ
自作の素因数分解プログラムを Qiita で公開した。プログラムに添削が入り，30倍以上速くなった。ちなみに教えて頂いた方法は「試し割り法」と呼ばれている。教えられた方法をそのまま使うのもアレなので，もっと速い方法を編み出せないか考えてみた。
現代のプロセッサでも割り算は比較的重いとされるので，与えられた自然数に対してエラトステネスの篩（ふるい）と呼ばれる手法を用いて必要最小限の素数を列挙し，それらの素数だけで試し割りをすれば除算回数は最小化され，格段の高速化が期待できるはずだ！しかし，エラトステネスの篩（ふるい）は大量にメモリを使用するせいか，意図したほど速くならなかった。
この問題自体は古くから知られており，対策として Segmented Sieve，日本語だと区分篩（ふるい）という手法が知られている。区分篩（ふるい）の効果は著しく，前回の記事までに試し割り法に肉薄するところまできた。
詳しくは下記の記事を参照されたい。
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
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(8) - Qiita
今回，区分篩（ふるい）を行う素数判定テーブルを前回記事の奇数のみからさらに２・３の倍数を除く。これにより最大で $1.5$ 倍の高速化が期待できるはずである。
2. 素数リストの定義
$5$ 以上 $9743$ 以下の素数は $1200$ 個ある。素数 $p$ の素数判定テーブルおける位置を $q$ とする。
| $i$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $\cdots$ | $1199$ |
|---|---|---|---|---|---|---|---|---|---|---|
| $q$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $9$ | $\cdots$ | $3247$ |
| $p$ | $5$ | $7$ | $11$ | $13$ | $17$ | $19$ | $23$ | $29$ | $\cdots$ | $9743$ |
素数の値 $p$ と位置 $q$ は相互変換可能である。
q = \frac{p - (p \bmod 3)}{3} = \left\lfloor\frac{p}{3}\right\rfloor \iff p = 3 q + 1 + (q \bmod 2)
3. 素数判定テーブルの定義
素数判定テーブル（２・３の倍数を除く）を定義する。位置 $j$ における値 $n$ とおく。
| $j$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $\cdots$ | $3246$ | $3247$ |
|---|---|---|---|---|---|---|---|---|---|---|---|
| $n$ | $1$ | $5$ | $7$ | $11$ | $13$ | $17$ | $19$ | $23$ | $\cdots$ | $9739$ | $9743$ |
素数判定テーブルの値 $n$ と位置 $j$ も相互変換可能である。
j = \frac{n - (n \bmod 3)}{3} = \left\lfloor\frac{n}{3}\right\rfloor \iff n = 3 j + 1 + (j \bmod 2)
4. 基本方針
比較のため前回の区分篩（ふるい）奇数探索版のコア部分を以下に示す。
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;
}
しかしながら自然数探索版，奇数探索版と異なり，２・３の倍数を除くと数列が等間隔では無くなってしまうため，非常に取り扱いが難しくなる。例えば素数 $p$ の倍数は以下のように非等間隔に並ぶ。
$5p \rightarrow 7p \rightarrow 11p \rightarrow 13p \rightarrow 17p \rightarrow 19p \rightarrow \cdots$
しかし，過去の記事で述べたように数列をひとつ置きに見ると綺麗に等間隔に並ぶことが分かる。
$5p \rightarrow 11p \rightarrow 17p \rightarrow 23p \rightarrow \cdots$
$7p \rightarrow 13p \rightarrow 19p \rightarrow 25p \rightarrow \cdots$
このため，区分篩（ふるい）のコア部分を次のように２つのループに分ける。ループ回数の総和はちょうど $2/3$ になるが，ループを２つに分けたことによる効率低下を考えると，どこまでの高速化を果たせるであろうか？
for(var i = 0; i < 1200; i++) {
var p = list1[i];
for(var j = start1[i]; j < 3248; j += 2 * p)
table[j] = 1;
start1[i] = j - 3248;
for(var j = start2[i]; j < 3248; j += 2 * p)
table[j] = 1;
start2[i] = j - 3248;
}
5. 代案（没案）
素数判定テーブル上での移動幅は $\Delta j = \lbrace 2p - (2q + 1), 2q + 1\rbrace$ の二値を交互に取ること，二つを足すとちょうど $2p$ になることから，以下のようにループを一つにまとめることもできる。しかし，実際にコーディングして比較してみるとパフォーマンスに劣ることが分かり，不採用となった。ループの中身の処理が僅かに増えるだけだが，それが致命的のようだ。初回位置の配列 start[i] は上記の start1[i] と start2[i] の最小値を与える。使用するメモリ容量は少なくなるが，その効果は薄かったようである。
for(var i = 0; i < 1200; i++) {
var p = list1[i];
var q = Math.floor(p / 3);
for(var j = start[i], dj = ((j & 1) == (q & 1)) ? 2 * p - (2 * q + 1) : 2 * q + 1;
j < 3248; j += dj, dj = 2 * p - dj)
table[j] = 1;
start[i] = j - 3248;
}
6. 初回位置
先に述べたように２つの初回位置 $j_1$，$j_2$ を持ち，それぞれの値 $n_1$，$n_2$ とする。
- $p^2 \ge 9745$ のとき
これより初回位置 $j_1$，$j_2$ は以下のように求められる。
\begin{aligned} n_1 &= p^2 \\ &= \left\lbrace 3q + 1 + (q \bmod 2)\right\rbrace^2 \\ &= 9q^2 + 6q + 6q \cdot (q \bmod 2) + 3 (q \bmod 2) + 1 \\ &= 3 \left\lbrace 3q^2 + 2q + (2q + 1) \cdot (q \bmod 2) \right\rbrace + 1 \\ &= 3 j_1 + 1 \\ j_1 &= 3q^2 + 2q + (2q + 1) \cdot (q \bmod 2) \\ n_2 &= \left\lbrace 3q + 1 + (q \bmod 2)\right\rbrace \cdot \left\lbrace 3(q + 1) + 2 - (q \bmod 2)\right\rbrace \\ &= 9q^2 + 18q + 3 (q \bmod 2) + 5 \\ &= 3 \left\lbrace 3q^2 + 6q + (q \bmod 2) + 1 \right\rbrace + 2 \\ &= 3 j_2 + 2 \\ j_2 &= 3q^2 + 6q + (q \bmod 2) + 1 \end{aligned}\begin{aligned} j_1 &= \begin{cases} 3q^2 + 2q & (q \bmod 2 = 0) \\ 3q^2 + 4q + 1 & (q \bmod 2 \ne 0) \end{cases} \\ j_2 &= \begin{cases} 3q^2 + 6q + 1 & (q \bmod 2 = 0) \\ 3q^2 + 6q + 2 & (q \bmod 2 \ne 0) \end{cases} \end{aligned} - $p^2 < 9745$ のとき
ここで $m_1$，$m_2$ は，
\begin{aligned} n_1 &= p + 6 m_1 \cdot p \\ &= 3q + 1 + (q \bmod 2) + 6 m_1 \cdot p \\ &= 3 (q + 2 m_1 \cdot p) + 1 + (q \bmod 2) \\ &= 3 j_1 + 1 + (q \bmod 2) \\ j_1 &= q + 2 m_1 \cdot p \\ n_2 &= 5p + 6 m_2 \cdot p \\ &= 5\left\lbrace 3q + 1 + (q \bmod 2) \right\rbrace + 6 m_2 \cdot p \\ &= 15q + 5 + 5(q \bmod 2) + 6 m_2 \cdot p \\ &= 3 \left\lbrace 5q + 1 + 2(q \bmod 2) + 2 m_2 \cdot p \right\rbrace + 2 - (q \bmod 2) \\ &= 3 j_2 + 2 - (q \bmod 2) \\ j_2 &= 5q + 1 + 2(q \bmod 2) + 2 m_2 \cdot p \\ &= 2p - q - 1 + 2 m_2 \cdot p \end{aligned}を満たす最小の自然数であり，天井関数 $\lceil x \rceil$ を用いて次のように表せる。\begin{aligned} j_1 &= q + 2 m_1 \cdot p \ge 3248 \\ j_2 &= 2p - q - 1 + 2 m_2 \cdot p \ge 3248 \end{aligned}これより初回位置 $j_1$，$j_2$ は以下のように求められる。\begin{aligned} m_1 &= \left\lceil \frac{3248 - q}{2p} \right\rceil \\ m_2 &= \left\lceil \frac{3248 - 2p + q + 1}{2p} \right\rceil \end{aligned}\begin{aligned} j_1 &= \begin{cases} 3248 & (3248 - q \bmod 2p = 0) \\ 3248 + 2p - (3248 - q) \bmod 2p & (3248 - q \bmod 2p \ne 0) \end{cases} \\ j_2 &= \begin{cases} 3248 & (3248 - 2p + q + 1 \bmod 2p = 0) \\ 3248 + 2p - (3248 - 2p + q + 1) \bmod 2p & (3248 - 2p + q + 1 \bmod 2p \ne 0) \end{cases} \end{aligned}
7. 実装コード
実装コードを以下に示す。なお，素数判定テーブル table[] のサイズは一次データキャッシュに収まるぎりぎりの範囲まで広げている。
| 内容 | 型 | 個数 | 容量 |
|---|---|---|---|
| 素数リスト | Uint16Array | $1,200$ | $2,400$ |
| 初回位置 | Uint32Array | $2 \times 1,200$ | $9,600$ |
| 素数判定テーブル | Uint8Array | $6 \times 3,248$ | $19,488$ |
| 合計 | $31,488$ |
//------------------------------------------------------------------------------
// 5 以上 9,744 以下の素数リスト ※要素数 1,200個
//------------------------------------------------------------------------------
var list1 = new Uint16Array([
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,
/* 中略 */
9697,9719,9721,9733,9739,9743]);
//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list(n) {
var nsqrt = Math.floor(Math.sqrt(n));
var list2 = []; // 得られた素数を格納するリスト
var start1 = new Uint32Array(1200); // 初回位置１
var start2 = new Uint32Array(1200); // 初回位置２
var table = new Uint8Array(3248*6); // 区分篩（ふるい）を行う素数判定テーブル
//--------------------------------------------------------------------------
// 初回位置の計算
//--------------------------------------------------------------------------
list2.push(2);
list2.push(3);
for(var i = 0; i < 1200; i++) {
var p = list1[i];
list2.push(p);
var q = Math.floor(p / 3);
if(p * p >= 9745) {
if(q & 1) {
start1[i] = (3 * q + 4) * q + 1 - 3248;
start2[i] = (3 * q + 6) * q + 2 - 3248;
} else {
start1[i] = (3 * q + 2) * q - 3248;
start2[i] = (3 * q + 6) * q + 1 - 3248;
}
} else {
var p2 = 2 * p;
var r = (3248 - q) % p2;
var s = (3248 - p2 + q + 1) % p2;
if(r != 0) start1[i] = p2 - r;
if(s != 0) start2[i] = p2 - s;
}
}
//--------------------------------------------------------------------------
// 区分篩（ふるい）の実行
//--------------------------------------------------------------------------
for(var offset = 9744; offset < nsqrt; offset += 9744*6) {
for(var i = 0; i < 1200; i++) {
var p = list1[i];
var p2 = 2 * p;
for(var j = start1[i]; j < 3248*6; j += p2)
table[j] = 1;
start1[i] = j - 3248*6;
for(var j = start2[i]; j < 3248*6; j += p2)
table[j] = 1;
start2[i] = j - 3248*6;
}
for(var j = 0; j < 3248*6; j++) {
if(table[j] == 0) list2.push(offset + 3 * j + 1 + (j & 1));
table[j] = 0;
}
}
return list2;
}
8. 結果
JavaScript において安全に使用できる整数最大値 Number.MAX_SAFE_INTEGER 以下の最大素数を入力し，素数判定に要した時間を計測した。なお JavaScript エンジンは Chakra エンジンを使用した。TypedArray を使用してるので，Internet Explorer の JavaScript エンジンでは動作しない。
公平のため前回までの記事で述べてきた区分篩（ふるい）の自然数探索版および奇数探索版についても，素数判定テーブルの大きさを一次データキャッシュに収まるぎりぎりの範囲まで広げており，前回の記事よりも少しだけパフォーマンスが高まっている。
| 内容 | 自然数探索版 | 奇数探索版 | ２・３の倍数を除く |
|---|---|---|---|
| 素数リスト | $2 \times 1,201$ | $2 \times 1,200$ | $2 \times 1,200$ |
| 初回位置 | $4 \times 1,201$ | $4 \times 1,200$ | $2 \times 4 \times 1,200$ |
| 素数判定テーブル | $2 \times 9,740$ | $4 \times 4,870$ | $6 \times 3,248$ |
| 合計 | $26,686$ | $26,680$ | $31,488$ |
9. 結論
ついに試し割り法（２・３・５の倍数を除く Wheel Factorization）に追いつくことが出来た！厳密には第8～11世代では上回ることは出来たのだが，最新の第13世代のみ上回ることができなかった。
今回の結果より，おそらく２・３・５の倍数を除く場合を実装してもこれ以上パフォーマンスを向上させることは期待できないと考えられる。つまり，今回が最終回である。
なお，この区分篩（ふるい）による素因数分解法（あるいは素数判定法）は与えられた自然数 $N$ が素数の場合に限り試し割り法よりも僅かに速いか同等レベルだが，合成数のときは試し割り法のほうが圧倒的に速いことを断っておく。これはいったん $\sqrt{N}$ 以下の素数リストを全て揃えてから試し割りを開始するためである・・・あれれ？ちょっと待て・・・その必要あるかな？
10. 延長戦
よくよく考えると，本題はあくまで与えられた自然数 $N$ の素数判定もしくは素因数分解を行えば良いのであって，もちろん $\sqrt{N}$ 以下の素数は全て必要だが，それらを素数リストとして保存しておく必要はないはずだ。ということで新たな素数を見出したらすぐに試し割りを行うよう書き換えた。
//------------------------------------------------------------------------------
// 5 以上 9,744 以下の素数リスト ※要素数 1,200個
//------------------------------------------------------------------------------
var list = new Uint16Array([
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,
/* 中略 */
9697,9719,9721,9733,9739,9743]);
//------------------------------------------------------------------------------
// 素因数分解を行う
//------------------------------------------------------------------------------
function prime_factorize(n) {
var nsqrt = Math.floor(Math.sqrt(n));
var start1 = new Uint32Array(1200); // 初回位置１
var start2 = new Uint32Array(1200); // 初回位置２
var table = new Uint8Array(3248*6); // 区分篩（ふるい）を行う素数判定テーブル
var a = []; // 素因数リスト
var m = n;
while(m % 2 == 0) {m /= 2; a.push(2);}
while(m % 3 == 0) {m /= 3; a.push(3);}
//--------------------------------------------------------------------------
// 初回位置の計算
//--------------------------------------------------------------------------
for(var i = 0; i < 1200; i++) {
var p = list[i];
while(m % p == 0) {m /= p; a.push(p);}
var q = Math.floor(p / 3);
if(p * p >= 9745) {
if(q & 1) {
start1[i] = (3 * q + 4) * q + 1 - 3248;
start2[i] = (3 * q + 6) * q + 2 - 3248;
} else {
start1[i] = (3 * q + 2) * q - 3248;
start2[i] = (3 * q + 6) * q + 1 - 3248;
}
} else {
var p2 = 2 * p;
var r = (3248 - q) % p2;
var s = (3248 - p2 + q + 1) % p2;
if(r != 0) start1[i] = p2 - r;
if(s != 0) start2[i] = p2 - s;
}
}
//--------------------------------------------------------------------------
// 区分篩（ふるい）の実行
//--------------------------------------------------------------------------
for(var offset = 9744; offset < nsqrt; offset += 9744*6) {
for(var i = 0; i < 1200; i++) {
var p = list[i];
var p2 = 2 * p;
for(var j = start1[i]; j < 3248*6; j += p2)
table[j] = 1;
start1[i] = j - 3248*6;
for(var j = start2[i]; j < 3248*6; j += p2)
table[j] = 1;
start2[i] = j - 3248*6;
}
for(var j = 0; j < 3248*6; j++) {
if(table[j] == 0) {
var p = offset + 3 * j + 1 + (j & 1);
while(m % p == 0) {m /= p; a.push(p);}
}
table[j] = 0;
}
}
if(m > 1) a.push(m);
return a;
}
上記のプログラムは突貫で作成したため，非常に作りが甘い。与えられた自然数が合成数の場合，もっと速くループを抜けることができる。
比較のため前回までの記事で述べてきた区分篩（ふるい）の自然数探索版および奇数探索版についても同様に書き換えた。ご覧のように全世代で試し割り法を上回ることができた！
しかし，自然数探索版は速くならず，奇数探索版および２・３の倍数を除く場合が高速化される理由は何だろう？区分篩（ふるい）では各テーブルが L1 データキャッシュ内に収まるよう工夫していたつもりだが，素数リスト list2[] は L1 データキャッシュに収まり切らなかったはず。つまり，素数リスト list2[] を保存しないことによって L1 データキャッシュの効率が上がる恩恵は自然数探索版でも等しく受けるはずなのだ。