前回までのあらすじ
詳しくは下記の記事を参照されたい。
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件 - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2) - Qiita
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(3) - Qiita
エラトステネスの篩（ふるい）のチューニング
前回の記事のコメント1より，重要な助言を頂いた。それは素数判定テーブル table[] への書き込み完了を待たずに素数リスト list[] の作成を開始できること。この結果，素数判定テーブル table[] を書き込んだ後，再度先頭から読み直す必要がなくなる。素数判定テーブルの大きさはキャッシュメモリから溢れているはずなので，この効果は大きい。
また記事2の実験結果より，効率的な配列アクセス方法に関する知見（TypedArray の使用）を得た。これらのアイディアを用いたチューニング後のコードを以下に示す。コードはエラトステネスの篩（ふるい）の奇数サーチ版をベースとしている。
//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
// 素数判定テーブル：table[i]，i 番目の奇数 2i + 1，素数は 0，合成数は 1 とする
// なお，初期化時には 0 で埋められることが保証されている
//------------------------------------------------------------------------------
function create_list( n ) {
var table = new Uint8Array( 1 + Math.floor( Math.sqrt( n ) / 2 ) );
var size = table.length;
var list = [];
list.push( 2 );
for( var i = 1; i < size; i++ ) {
if( table[i] != 0 ) continue;
var m = 2 * i + 1;
list.push( m );
for( var j = 2 * i * ( i + 1 ); j < size; j += m )
table[j] = 1;
}
return( list );
}
チューニング前後の比較を以下に示す。なお，今回より JavaScript エンジンは Chakra に統一することにした。また，与える素数の大きさを JavaScript で安全に扱える整数最大値 $2^{53} - 1$ に最も近い素数まで広げるようにした。これは，あっという間に計算が終わってしまうからである。
この結果，およそ 2.5～3 倍の高速化を達成することができた！
測定時間の詳細はコチラ
| 素数 | チューニング前 | チューニング後 |
|---|---|---|
| 1,250,000,000,111 | 0.065 | 0.050 |
| 2,500,000,000,009 | 0.069 | 0.051 |
| 5,000,000,000,053 | 0.082 | 0.056 |
| 10,000,000,000,037 | 0.103 | 0.062 |
| 20,000,000,000,021 | 0.126 | 0.071 |
| 40,000,000,000,001 | 0.174 | 0.085 |
| 80,000,000,000,027 | 0.245 | 0.100 |
| 160,000,000,000,069 | 0.317 | 0.126 |
| 320,000,000,000,029 | 0.458 | 0.158 |
| 640,000,000,000,033 | 0.697 | 0.224 |
| 1,250,000,000,000,111 | 0.884 | 0.315 |
| 2,500,000,000,000,043 | 1.345 | 0.479 |
| 5,000,000,000,000,023 | 1.922 | 0.707 |
| 9,007,199,254,740,881 | 2.349 | 0.960 |
まだ試し割り法は本気を出していない。
正しくは試し割り法の周期 8 の Wheel Factorization である。指導教官の指導3の下，オリジナルのコードに対して，以下の最適化を施した。
- 剰余演算
p % q == 0を(r = p / q) == Math.round(r)と分解した。 - 手動で loop unrolling を施した。
JavaScript の数値演算は浮動小数点で行われるので p % q = p - Math.floor(p / q) * q と展開されるはずで，処理内容は大きく変わらないはずであるが，なぜか大幅に高速化される。なお，２つめの loop unrolling 単体での効果は小さく，１つめの剰余演算の分解との合わせ技で効果を発揮する。
最適化後のコードを以下に示す。ちなみにこのコードでは for ループの条件判定 q * q <= p すら８回に１回と間引いている。この結果，for ループ実行中に q * q <= p の範囲を少し飛び越えてしまうが，それでも問題ないというのが恐ろしい。
素数判定の場合，与えられた自然数 $N$ に対して $\sqrt{N}$ までの範囲で試し割りのチェックを行えば良いが，別に $\sqrt{N}$ を超えてしまっても問題ない。もちろん超過した分の計算時間は余分にかかってしまうが，ループの条件判定を間引いた分のリターンのほうがはるかに大きいという目算である。
//------------------------------------------------------------------------------
// 素因数分解を行う：周期 8 の Wheel Factorization を適用したもの
//------------------------------------------------------------------------------
function prime_factorize( p ) {
var a = [];
var r;
while((r = p / 2) == Math.round(r)) { p = r; a.push(2); }
while((r = p / 3) == Math.round(r)) { p = r; a.push(3); }
while((r = p / 5) == Math.round(r)) { p = r; a.push(5); }
for( var q = 7; q * q <= p; ) {
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 4;
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 2;
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 4;
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 2;
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 4;
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 6;
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 2;
while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 6;
}
if( p > 1 ) a.push( p );
return( a );
}
最適化前後の比較を以下に示す。
うわぁぁぁぁ！コチラは一気に３倍以上高速化してしまった！
測定時間の詳細はコチラ
| 素数 | 最適前 | 最適後 |
|---|---|---|
| 1,250,000,000,111 | 0.056 | 0.043 |
| 2,500,000,000,009 | 0.062 | 0.042 |
| 5,000,000,000,053 | 0.069 | 0.041 |
| 10,000,000,000,037 | 0.084 | 0.043 |
| 20,000,000,000,021 | 0.103 | 0.042 |
| 40,000,000,000,001 | 0.128 | 0.045 |
| 80,000,000,000,027 | 0.165 | 0.046 |
| 160,000,000,000,069 | 0.215 | 0.070 |
| 320,000,000,000,029 | 0.305 | 0.109 |
| 640,000,000,000,033 | 0.372 | 0.137 |
| 1,250,000,000,000,111 | 0.544 | 0.175 |
| 2,500,000,000,000,043 | 0.777 | 0.230 |
| 5,000,000,000,000,023 | 1.032 | 0.305 |
| 9,007,199,254,740,881 | 1.455 | 0.397 |
高速化の理由ははっきり言ってよく分からない。厳密には乗算が減って，丸め演算の処理が増える。インテルのマニュアル4を見る限り，浮動小数点乗算はかなり高速で，むしろ丸め演算のほうが少し時間がかかりそうである。
除算の結果・商 r を再利用できるのは嬉しい。Intel の x86 プロッセッサの場合，整数除算命令では商と剰余を同時に得ることができるが，浮動小数点除算命令では商しか得られないからである。ただ，上記のテストではすべて素数を与えているため，商 r を一度も再利用していないはずで，高速化の理由は謎である。
ちなみに本来，剰余を求めるためには Math.floor() を使わなくてはならないが，本例のように割り切れるか否かの判定だけ行えばよい場合は Math.round() でも構わない。ただし Math.floor() よりも Math.round() のほうが有意に高速だったことから Math.round() を採用した。なお，インテルのマニュアルを参照する限り，丸め演算のモード round/floor/ceil による処理時間の違いは無いようである。
なお，最適化後において与える素数の大きさが $10^{14}$ を超えると急激に計算時間が増加し，その後，一律な増加レートに落ち着いている。おそらく $10^{14}$ を超える前後でインタプリタ型から JIT コンパイラ型へ切り替わっているものと推察される。
結論
悲報，差が広がってしまいました・・・
備考
- 測定マシン仕様 Core i5 10400 2.9GHz, 8GB, Windows10 22H2
- 計測時間は 10 回の平均値です。
次回
次回に続きます。
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(5) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(6) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(7) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(8) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(9) 最終回 - Qiita