Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(4)

# 前回までのあらすじ
詳しくは下記の記事を参照されたい。
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件 - Qiita](https://qiita.com/TETSURO1999/items/ffc1cae68c16a5c26fde)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2) - Qiita](https://qiita.com/TETSURO1999/items/28ca00f958628423d64f?1)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(3) - Qiita](https://qiita.com/TETSURO1999/items/38ddcc0f4f3411924fe9)

# エラトステネスの篩（ふるい）のチューニング

前回の記事のコメント[^1]より，重要な助言を頂いた。それは素数判定テーブル `table[]` への書き込み完了を待たずに素数リスト `list[]` の作成を開始できること。この結果，素数判定テーブル `table[]` を書き込んだ後，再度先頭から読み直す必要がなくなる。素数判定テーブルの大きさはキャッシュメモリから溢れているはずなので，この効果は大きい。

また記事[^2]の実験結果より，効率的な配列アクセス方法に関する知見（`TypedArray` の使用）を得た。これらのアイディアを用いたチューニング後のコードを以下に示す。コードはエラトステネスの篩（ふるい）の奇数サーチ版をベースとしている。

[^1]:[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2)のコメント](https://qiita.com/TETSURO1999/items/28ca00f958628423d64f#comment-fead1338bc929c25ce72)

[^2]:[JavaScriptによる大規模配列初期化のパフォーマンス検証](https://qiita.com/TETSURO1999/items/2621a4d6c00adba85ff6)

```javascript
//------------------------------------------------------------------------------
// 与えられた自然数 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$ に最も近い素数まで広げるようにした。これは，あっという間に計算が終わってしまうからである。

```plantuml
@startuml
<style>
root [
  Padding 0
  Margin 0
}
header {
  fontColor black
  fontSize 14
}
footer {
  fontColor black
  fontSize 16
}
</style>
left header <font color=red>▲</font>チューニング前　<font color=cyan>■</font>チューニング後
footer
図１　エラトステネスの篩（ふるい）奇数サーチ版
チューニング前後の比較
end footer
title <img https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/886d02c2-dc66-bc70-9473-e984329aa726.png>
@enduml
```

この結果，およそ 2.5～3 倍の高速化を達成することができた！

<details><summary>測定時間の詳細はコチラ</summary>

<table>
<caption>表１　エラトステネスの篩（ふるい）奇数サーチ版<br>チューニング前後の比較（単位は秒）</caption>
<thead>
<tr><th style="text-align:center;">素数</th>
    <th style="text-align:center;">チューニング前</th>
    <th style="text-align:center;">チューニング後</th></tr>
</thead>
<tbody>
<tr><td style="text-align:right;">1,250,000,000,111</td>
    <td style="text-align:right;">0.065</td>
    <td style="text-align:right;">0.050</td></tr>
<tr><td style="text-align:right;">2,500,000,000,009</td>
    <td style="text-align:right;">0.069</td>
    <td style="text-align:right;">0.051</td></tr>
<tr><td style="text-align:right;">5,000,000,000,053</td>
    <td style="text-align:right;">0.082</td>
    <td style="text-align:right;">0.056</td></tr>
<tr><td style="text-align:right;">10,000,000,000,037</td>
    <td style="text-align:right;">0.103</td>
    <td style="text-align:right;">0.062</td></tr>
<tr><td style="text-align:right;">20,000,000,000,021</td>
    <td style="text-align:right;">0.126</td>
    <td style="text-align:right;">0.071</td></tr>
<tr><td style="text-align:right;">40,000,000,000,001</td>
    <td style="text-align:right;">0.174</td>
    <td style="text-align:right;">0.085</td></tr>
<tr><td style="text-align:right;">80,000,000,000,027</td>
    <td style="text-align:right;">0.245</td>
    <td style="text-align:right;">0.100</td></tr>
<tr><td style="text-align:right;">160,000,000,000,069</td>
    <td style="text-align:right;">0.317</td>
    <td style="text-align:right;">0.126</td></tr>
<tr><td style="text-align:right;">320,000,000,000,029</td>
    <td style="text-align:right;">0.458</td>
    <td style="text-align:right;">0.158</td></tr>
<tr><td style="text-align:right;">640,000,000,000,033</td>
    <td style="text-align:right;">0.697</td>
    <td style="text-align:right;">0.224</td></tr>
<tr><td style="text-align:right;">1,250,000,000,000,111</td>
    <td style="text-align:right;">0.884</td>
    <td style="text-align:right;">0.315</td></tr>
<tr><td style="text-align:right;">2,500,000,000,000,043</td>
    <td style="text-align:right;">1.345</td>
    <td style="text-align:right;">0.479</td></tr>
<tr><td style="text-align:right;">5,000,000,000,000,023</td>
    <td style="text-align:right;">1.922</td>
    <td style="text-align:right;">0.707</td></tr>
<tr><td style="text-align:right;">9,007,199,254,740,881</td>
    <td style="text-align:right;">2.349</td>
    <td style="text-align:right;">0.960</td></tr>
</tbody>
</table>

</details>

# まだ試し割り法は本気を出していない。

正しくは試し割り法の周期 8 の Wheel Factorization である。指導教官の指導[^3]の下，オリジナルのコードに対して，以下の最適化を施した。

[^3]:[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2)のコメント](https://qiita.com/TETSURO1999/items/28ca00f958628423d64f#comment-4bfdde2d40becd097f64)

- 剰余演算 `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` の範囲を少し飛び越えてしまうが，それでも問題ないというのが恐ろしい。

:::note
素数判定の場合，与えられた自然数 $N$ に対して $\sqrt{N}$ までの範囲で試し割りのチェックを行えば良いが，別に $\sqrt{N}$ を超えてしまっても問題ない。もちろん超過した分の計算時間は余分にかかってしまうが，ループの条件判定を間引いた分のリターンのほうがはるかに大きいという目算である。
:::

```javascript
//------------------------------------------------------------------------------
// 素因数分解を行う：周期 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 );
}
```

最適化前後の比較を以下に示す。

```plantuml
@startuml
<style>
root [
  Padding 0
  Margin 0
}
header {
  fontColor black
  fontSize 14
}
footer {
  fontColor black
  fontSize 16
}
</style>
left header <font color=magenta>●</font>最適化前　<font color=black>×</font>最適化後
footer 図２　試し割り法（Wheel Factorization）最適化前後の比較
title <img https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/bbc96ec0-ea02-2ef5-c4f6-c362a8151003.png>
@enduml
```

うわぁぁぁぁ！コチラは一気に３倍以上高速化してしまった！

<details><summary>測定時間の詳細はコチラ</summary>

<table>
<caption>表２　試し割り法（Wheel Factorization）<br>最適化前後の比較（単位は秒）</caption>
<thead>
<tr><th style="text-align:center;">素数</th>
    <th style="text-align:center;">最適前</th>
    <th style="text-align:center;">最適後</th></tr>
</thead>
<tbody>
<tr><td style="text-align:right;">1,250,000,000,111</td>
    <td style="text-align:right;">0.056</td>
    <td style="text-align:right;">0.043</td></tr>
<tr><td style="text-align:right;">2,500,000,000,009</td>
    <td style="text-align:right;">0.062</td>
    <td style="text-align:right;">0.042</td></tr>
<tr><td style="text-align:right;">5,000,000,000,053</td>
    <td style="text-align:right;">0.069</td>
    <td style="text-align:right;">0.041</td></tr>
<tr><td style="text-align:right;">10,000,000,000,037</td>
    <td style="text-align:right;">0.084</td>
    <td style="text-align:right;">0.043</td></tr>
<tr><td style="text-align:right;">20,000,000,000,021</td>
    <td style="text-align:right;">0.103</td>
    <td style="text-align:right;">0.042</td></tr>
<tr><td style="text-align:right;">40,000,000,000,001</td>
    <td style="text-align:right;">0.128</td>
    <td style="text-align:right;">0.045</td></tr>
<tr><td style="text-align:right;">80,000,000,000,027</td>
    <td style="text-align:right;">0.165</td>
    <td style="text-align:right;">0.046</td></tr>
<tr><td style="text-align:right;">160,000,000,000,069</td>
    <td style="text-align:right;">0.215</td>
    <td style="text-align:right;">0.070</td></tr>
<tr><td style="text-align:right;">320,000,000,000,029</td>
    <td style="text-align:right;">0.305</td>
    <td style="text-align:right;">0.109</td></tr>
<tr><td style="text-align:right;">640,000,000,000,033</td>
    <td style="text-align:right;">0.372</td>
    <td style="text-align:right;">0.137</td></tr>
<tr><td style="text-align:right;">1,250,000,000,000,111</td>
    <td style="text-align:right;">0.544</td>
    <td style="text-align:right;">0.175</td></tr>
<tr><td style="text-align:right;">2,500,000,000,000,043</td>
    <td style="text-align:right;">0.777</td>
    <td style="text-align:right;">0.230</td></tr>
<tr><td style="text-align:right;">5,000,000,000,000,023</td>
    <td style="text-align:right;">1.032</td>
    <td style="text-align:right;">0.305</td></tr>
<tr><td style="text-align:right;">9,007,199,254,740,881</td>
    <td style="text-align:right;">1.455</td>
    <td style="text-align:right;">0.397</td></tr>
</tbody>
</table>
</details>

高速化の理由ははっきり言ってよく分からない。厳密には乗算が減って，丸め演算の処理が増える。インテルのマニュアル[^4]を見る限り，浮動小数点乗算はかなり高速で，むしろ丸め演算のほうが少し時間がかかりそうである。

除算の結果・商 `r` を再利用できるのは嬉しい。Intel の x86 プロッセッサの場合，整数除算命令では商と剰余を同時に得ることができるが，浮動小数点除算命令では商しか得られないからである。ただ，上記のテストではすべて素数を与えているため，商 `r` を一度も再利用していないはずで，高速化の理由は謎である。

ちなみに本来，剰余を求めるためには `Math.floor()` を使わなくてはならないが，本例のように割り切れるか否かの判定だけ行えばよい場合は `Math.round()` でも構わない。ただし `Math.floor()` よりも `Math.round()` のほうが有意に高速だったことから `Math.round()` を採用した。なお，インテルのマニュアルを参照する限り，丸め演算のモード round/floor/ceil による処理時間の違いは無いようである。

なお，最適化後において与える素数の大きさが $10^{14}$ を超えると急激に計算時間が増加し，その後，一律な増加レートに落ち着いている。おそらく $10^{14}$ を超える前後でインタプリタ型から JIT コンパイラ型へ切り替わっているものと推察される。

[^4]:[インテル64およびIA-32アーキテクチャー最適化リファレンス・マニュアル参考訳](https://www.intel.co.jp/content/dam/www/central-libraries/jp/ja/documents/2022-04/ia-optimization-reference-manual-2022-02-jp.pdf)

# 結論

<font color=magenta><b>悲報，差が広がってしまいました・・・</b></font>

# 備考

- 測定マシン仕様 Core i5 10400 2.9GHz, 8GB, Windows10 22H2
- 計測時間は 10 回の平均値です。

# 次回
次回に続きます。

- [Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(5) - Qiita](https://qiita.com/TETSURO1999/items/4957acd1614aa319f787)
- [Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(6) - Qiita](https://qiita.com/TETSURO1999/items/a51747cfac9f89b51a11)
- [Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(7) - Qiita](https://qiita.com/TETSURO1999/items/392fd3a8f999a018edfd)
- [Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(8) - Qiita](https://qiita.com/TETSURO1999/items/599d4b7ccb2b4e407a34)
- [Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(9) 最終回 - Qiita](https://qiita.com/TETSURO1999/items/5d69e183af27487de7bc)