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

# 前回までのあらすじ

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

今回のテーマは打倒！試し割り法（Wheel Factorization）である。試し割り法の Wheel Factorization と同様のアイディアを用いてエラトステネスの篩（ふるい）をより省メモリ化，効率化，高速化する。

# ２，３の倍数を除く場合

奇数と異なり，数値が等間隔に並んでいないので取り扱いが格段に難しくなる。

２，３の倍数を除いた $(j，k)$ 番目の自然数を $m_{j,k} = 6j + w_k$ とおく。
ただし $j \ge 1$，$w_k = \lbrace -1, +1 \rbrace$ である。

<table>
<caption>表１　２，３の倍数を除く自然数</caption>
<thead>
<tr><th style="text-align:center;">$i$</th>
    <th style="text-align:right;">$2$</th>
    <th style="text-align:right;">$3$</th>
    <th style="text-align:right;">$4$</th>
    <th style="text-align:right;">$5$</th>
    <th style="text-align:right;">$6$</th>
    <th style="text-align:right;">$7$</th>
    <th style="text-align:right;">$8$</th>
    <th style="text-align:right;">$9$</th>
    <th style="text-align:right;">$10$</th>
    <th style="text-align:right;">$11$</th>
    <th style="text-align:right;">$12$</th>
    <th style="text-align:right;">$13$</th></tr>
</thead>
<tbody>
<tr><td style="text-align:center;">$j$</td>
    <td colspan=2 style="text-align:center;">$1$</td>
    <td colspan=2 style="text-align:center;">$2$</td>
    <td colspan=2 style="text-align:center;">$3$</td>
    <td colspan=2 style="text-align:center;">$4$</td>
    <td colspan=2 style="text-align:center;">$5$</td>
    <td colspan=2 style="text-align:center;">$6$</td></tr>
<tr><td style="text-align:center;">$k$</td>
    <td style="text-align:right;">$0$</td>
    <td style="text-align:right;">$1$</td>
    <td style="text-align:right;">$0$</td>
    <td style="text-align:right;">$1$</td>
    <td style="text-align:right;">$0$</td>
    <td style="text-align:right;">$1$</td>
    <td style="text-align:right;">$0$</td>
    <td style="text-align:right;">$1$</td>
    <td style="text-align:right;">$0$</td>
    <td style="text-align:right;">$1$</td>
    <td style="text-align:right;">$0$</td>
    <td style="text-align:right;">$1$</td></tr>
<tr><td style="text-align:center;">$w_k$</td>
    <td style="text-align:right;">$-1$</td>
    <td style="text-align:right;">$+1$</td>
    <td style="text-align:right;">$-1$</td>
    <td style="text-align:right;">$+1$</td>
    <td style="text-align:right;">$-1$</td>
    <td style="text-align:right;">$+1$</td>
    <td style="text-align:right;">$-1$</td>
    <td style="text-align:right;">$+1$</td>
    <td style="text-align:right;">$-1$</td>
    <td style="text-align:right;">$+1$</td>
    <td style="text-align:right;">$-1$</td>
    <td style="text-align:right;">$+1$</td></tr>
<tr><td style="text-align:center;">$m_{j,k}$</td>
    <td style="text-align:right;">$5$</td>
    <td style="text-align:right;">$7$</td>
    <td style="text-align:right;">$11$</td>
    <td style="text-align:right;">$13$</td>
    <td style="text-align:right;">$17$</td>
    <td style="text-align:right;">$19$</td>
    <td style="text-align:right;">$23$</td>
    <td style="text-align:right;">$25$</td>
    <td style="text-align:right;">$29$</td>
    <td style="text-align:right;">$31$</td>
    <td style="text-align:right;">$35$</td>
    <td style="text-align:right;">$37$</td></tr>
</tbody>
</table>

$m_{j,k}^2$ 以上の $m_{j,k}$ の倍数を合成数として素数判定テーブルから除去していく。

初回位置 $m_{j,k}^2 = m_{J_0,K_0}$ とおく。

```math
m_{j,k}^2 = (6j + w_k)^2 = 6(6j^2 + 2j\cdot w_k) + w_k^2
```
より $J_0 = 6j^2 + 2j\cdot w_k$，$K_0 = w_k^2 = 1$ となる。すなわち $m_{j,k}^2$ は２，３の倍数を除いた自然数において $(J_0，K_0)$ 番目の数である。

$k = 0$ のときの位置の推移を示す。

```math
\begin{align}
m_{J_0,K_0} &= m_{j,0}^2 = (6j - 1)(6j - 1) = 6(6j^2 - 2j) + 1 \\
m_{J_1,K_1} &= m_{j,0}\cdot m_{j,1} = (6j - 1)(6j + 1) = 6(6j^2) - 1  \\
m_{J_2,K_2} &= m_{j,0}\cdot m_{j+1,0} = (6j - 1)(6j + 5) = 6(6j^2 + 4j -1) + 1  \\
m_{J_3,K_3} &= m_{j,0}\cdot m_{j+1,1} = (6j - 1)(6j + 7) = 6(6j^2 +6j - 1) - 1
\end{align}
```

$k = 1$ のときの位置の推移を示す。

```math
\begin{align}
m_{J_0,K_0} &= m_{j,1}^2 = (6j + 1)(6j + 1) = 6(6j^2 + 2j) + 1 \\
m_{J_1,K_1} &= m_{j,1}\cdot m_{j+1,0} = (6j + 1)(6j + 5) = 6(6j^2 + 6j + 1) - 1 \\
m_{J_2,K_2} &= m_{j,1}\cdot m_{j+1,1} = (6j + 1)(6j + 7) = 6(6j^2 + 8j + 1) + 1 \\
m_{J_3,K_3} &= m_{j,1}\cdot m_{j+2,0} = (6j + 1)(6j + 11) = 6(6j^2 +12j + 2) - 1
\end{align}
```
以上をまとめると以下の表になる。

<table>
<caption>表２　消去する合成数の位置 $(J_n,K_n)$ の推移</caption>
<tr><th colspan=4 style="text-align:center;">$k = 0$ のときの位置の推移</th>
    <th rowspan=6></th>
    <th colspan=4 style="text-align:center;">$k = 1$ のときの位置の推移</th></tr>
<tr><th style="text-align:center;">$n$</th>
    <th style="text-align:center;">$J_n$</th>
    <th style="text-align:center;">$K_n$</th>
    <th style="text-align:center;">$\Delta J_n$</th>
    <th style="text-align:center;">$n$</th>
    <th style="text-align:center;">$J_n$</th>
    <th style="text-align:center;">$K_n$</th>
    <th style="text-align:center;">$\Delta J_n$</th></tr>
<tr><td style="text-align:center;">$0$</td>
    <td style="text-align:center;">$6j^2 - 2j$</td>
    <td style="text-align:center;">$1$</td>
    <td style="text-align:center;">$2j$</td>
    <td style="text-align:center;">$0$</td>
    <td style="text-align:center;">$6j^2 + 2j$</td>
    <td style="text-align:center;">$1$</td>
    <td style="text-align:center;">$4j + 1$</td></tr>
<tr><td style="text-align:center;">$1$</td>
    <td style="text-align:center;">$6j^2$</td>
    <td style="text-align:center;">$0$</td>
    <td style="text-align:center;">$4j-1$</td>
    <td style="text-align:center;">$1$</td>
    <td style="text-align:center;">$6j^2 + 6j + 1$</td>
    <td style="text-align:center;">$0$</td>
    <td style="text-align:center;">$2j$</td></tr>
<tr><td style="text-align:center;">$2$</td>
    <td style="text-align:center;">$6j^2 + 4j - 1$</td>
    <td style="text-align:center;">$1$</td>
    <td style="text-align:center;">$2j$</td>
    <td style="text-align:center;">$2$</td>
    <td style="text-align:center;">$6j^2 + 8j + 1$</td>
    <td style="text-align:center;">$1$</td>
    <td style="text-align:center;">$4j + 1$</td></tr>
<tr><td style="text-align:center;">$3$</td>
    <td style="text-align:center;">$6j^2 + 6j - 1$</td>
    <td style="text-align:center;">$0$</td>
    <td style="text-align:center;">$4j -1$</td>
    <td style="text-align:center;">$3$</td>
    <td style="text-align:center;">$6j^2 + 12j + 2$</td>
    <td style="text-align:center;">$0$</td>
    <td style="text-align:center;">$2j$</td></tr>
</tbody>
</table>

以上より，$K_n$ は $1 \rightarrow 0 \rightarrow 1 \rightarrow 0$ のように必ず $1$ からスタートして，$0$ と $1$ の値を交互にとる。一方，$\Delta J_n$ は以下のようにまとめることができる。

```math
\Delta J_n = \cases{
2j & \(k = K_n\) のとき \\
4j + w_k & \(k \ne K_n\) のとき
}
```

以上をまとめると下記のコードになる。奇数版に比べるとずいぶん込み入ったコードになってしまった。

```javascript
//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list( n ) {
	//--------------------------------------------------------------------------
	// エラトステネスのふるいにより，素数フラグのテーブルを作成する
	//--------------------------------------------------------------------------
	var	table = Array( 2 + Math.floor( Math.sqrt( n ) / 3 ) );
	var	size = table.length;
	table[0] = false;	// -1 は素数ではない
	table[1] = false;	//  1 は素数ではない
	for( var i = 2; i < size; i++ ) table[i] = true;
	var	w = [-1, 1];
	for( var i = 2, j = 1, k = 0; i < size; i++, j += k, k ^= 1 ) {
		if( !table[i] ) continue;
		var	ww = w[k];
		for( var jj = 6 * j * j + 2 * j * ww, kk = 1; ( ii = 2 * jj + kk ) < size; kk ^= 1 ) {
			table[ii] = false;
			jj += ( kk == k ) ? 4 * j + ww : 2 * j;
		}
	}
	//--------------------------------------------------------------------------
	// 素数リストを作成する
	//--------------------------------------------------------------------------
	var	list = [];
	list.push( 2 );
	list.push( 3 );
	for( var i = 2, j = 1, k = 0; i < size; i++, j += k, k ^= 1 )
		if( table[i] ) list.push( 6 * j + w[k] );
	return( list );
}
```
計算時間の実測結果を以下に示す。

デフォルトの javascript エンジンの場合を以下に示す。使用するメモリ容量は，エラトステネスの篩（ふるい）の１／３，奇数版に対しても２／３になっているが，複雑になってしまったコードと相殺して優位性が見えない。ただ，使用するメモリ容量が減ったことにより，キャッシュメモリ不足による計算時間の急増が抑制されている。

```plantuml
@startuml
<style>
header {
  fontColor black
  fontSize 14
  horizontalAlignment left
}
footer {
  fontColor black
  fontSize 16
}
title {
  fontStyle normal
}
</style>
header
<font color=black>×</font>Wheel Factorization <font color=magenta>●</font>全数サーチ　<font color=cyan>◆</font>奇数サーチ　<font color=lime>■</font>２，３の倍数除去
end header
footer 図１　デフォルトの javascript エンジンの場合
title <img https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/782736f4-6b37-55aa-04f0-5d841b8e54e8.png>
@enduml
```
javascript エンジンを JScript9 に差し替えた場合を以下に示す。奇数版よりは速くなっているが，優位性は僅かである。

```plantuml
@startuml
<style>
header {
  fontColor black
  fontSize 14
  horizontalAlignment left
}
footer {
  fontColor black
  fontSize 16
}
title {
  fontStyle normal
}
</style>
header
<font color=black>×</font>Wheel Factorization <font color=magenta>●</font>全数サーチ　<font color=cyan>◆</font>奇数サーチ　<font color=lime>■</font>２，３の倍数除去
end header
footer 図２　JScript9 エンジンに差し替えた場合
title <img https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/dadfe693-b101-5147-3fe8-869dad4fd696.png>
@enduml
```

javascript エンジンを Chakra に差し替えた場合を以下に示す。JIT コンパイラが優秀なせいか，複雑なコードをうまく最適化して，エラトステネスの篩（ふるい）勢の計算時間は使用するメモリ容量に比例するようになっている。また試し割り法（Wheel Factorization）にも肉薄している。

```plantuml
@startuml
<style>
header {
  fontColor black
  fontSize 14
  horizontalAlignment left
}
footer {
  fontColor black
  fontSize 16
}
title {
  fontStyle normal
}
</style>
header
<font color=black>×</font>Wheel Factorization <font color=magenta>●</font>全数サーチ　<font color=cyan>◆</font>奇数サーチ　<font color=lime>■</font>２，３の倍数除去
end header
footer 図３　Chakra エンジンに差し替えた場合
title <img https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/aa2735d9-797f-1aea-30dc-c9e91cac36cc.png>
@enduml
```

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

<table>
<caption>表３　デフォルトの Javascript エンジンの場合</caption>
<thead>
<th rowspan=2 style="text-align:center;">素数</th>
    <th style="text-align:center;">試し割り法</th>
    <th colspan=3 style="text-align:center;">エラトステネスの篩（ふるい）</th></tr>
<tr><th style="text-align:center;">Wheel<br>Factorization</th>
    <th style="text-align:center;">全数サーチ</th>
    <th style="text-align:center;">奇数サーチ</th>
    <th style="text-align:center;">２，３の<br>倍数除去</th></tr>
</thead>
<tbody>
<tr><td style="text-align:right;">1,250,000,000,111</td>
    <td style="text-align:right;">0.098</td>
    <td style="text-align:right;">0.560</td>
    <td style="text-align:right;">0.310</td>
    <td style="text-align:right;">0.310</td></tr>
<tr><td style="text-align:right;">2,500,000,000,009</td>
    <td style="text-align:right;">0.117</td>
    <td style="text-align:right;">0.730</td>
    <td style="text-align:right;">0.430</td>
    <td style="text-align:right;">0.400</td></tr>
<tr><td style="text-align:right;">5,000,000,000,053</td>
    <td style="text-align:right;">0.143</td>
    <td style="text-align:right;">1.010</td>
    <td style="text-align:right;">0.550</td>
    <td style="text-align:right;">0.550</td></tr>
<tr><td style="text-align:right;">10,000,000,000,037</td>
    <td style="text-align:right;">0.179</td>
    <td style="text-align:right;">1.460</td>
    <td style="text-align:right;">0.810</td>
    <td style="text-align:right;">0.770</td></tr>
<tr><td style="text-align:right;">20,000,000,000,021</td>
    <td style="text-align:right;">0.234</td>
    <td style="text-align:right;">2.040</td>
    <td style="text-align:right;">1.100</td>
    <td style="text-align:right;">1.080</td></tr>
<tr><td style="text-align:right;">40,000,000,000,001</td>
    <td style="text-align:right;">0.307</td>
    <td style="text-align:right;">2.940</td>
    <td style="text-align:right;">1.530</td>
    <td style="text-align:right;">1.530</td></tr>
<tr><td style="text-align:right;">80,000,000,000,027</td>
    <td style="text-align:right;">0.411</td>
    <td style="text-align:right;">6.680</td>
    <td style="text-align:right;">2.160</td>
    <td style="text-align:right;">2.160</td></tr>
<tr><td style="text-align:right;">160,000,000,000,069</td>
    <td style="text-align:right;">0.562</td>
    <td style="text-align:right;">32.850</td>
    <td style="text-align:right;">3.110</td>
    <td style="text-align:right;">3.070</td></tr>
<tr><td style="text-align:right;">320,000,000,000,029</td>
    <td style="text-align:right;">0.775</td>
    <td style="text-align:right;">92.640</td>
    <td style="text-align:right;">6.840</td>
    <td style="text-align:right;">4.410</td></tr>
<tr><td style="text-align:right;">640,000,000,000,033</td>
    <td style="text-align:right;">1.078</td>
    <td style="text-align:right;">221.080</td>
    <td style="text-align:right;">32.920</td>
    <td style="text-align:right;">6.560</td></tr>
</tbody>
</table>

<table>
<caption>表４　JScript9 エンジンに差し替えた場合</caption>
<thead>
<th rowspan=2 style="text-align:center;">素数</th>
    <th style="text-align:center;">試し割り法</th>
    <th colspan=3 style="text-align:center;">エラトステネスの篩（ふるい）</th></tr>
<tr><th style="text-align:center;">Wheel<br>Factorization</th>
    <th style="text-align:center;">全数サーチ</th>
    <th style="text-align:center;">奇数サーチ</th>
    <th style="text-align:center;">２，３の<br>倍数除去</th></tr>
</thead>
<tbody>
<tr><td style="text-align:right;">1,250,000,000,111</td>
    <td style="text-align:right;">0.044</td>
    <td style="text-align:right;">0.072</td>
    <td style="text-align:right;">0.055</td>
    <td style="text-align:right;">0.057</td></tr>
<tr><td style="text-align:right;">2,500,000,000,009</td>
    <td style="text-align:right;">0.047</td>
    <td style="text-align:right;">0.084</td>
    <td style="text-align:right;">0.061</td>
    <td style="text-align:right;">0.063</td></tr>
<tr><td style="text-align:right;">5,000,000,000,053</td>
    <td style="text-align:right;">0.053</td>
    <td style="text-align:right;">0.112</td>
    <td style="text-align:right;">0.077</td>
    <td style="text-align:right;">0.074</td></tr>
<tr><td style="text-align:right;">10,000,000,000,037</td>
    <td style="text-align:right;">0.062</td>
    <td style="text-align:right;">0.174</td>
    <td style="text-align:right;">0.088</td>
    <td style="text-align:right;">0.093</td></tr>
<tr><td style="text-align:right;">20,000,000,000,021</td>
    <td style="text-align:right;">0.073</td>
    <td style="text-align:right;">0.217</td>
    <td style="text-align:right;">0.119</td>
    <td style="text-align:right;">0.109</td></tr>
<tr><td style="text-align:right;">40,000,000,000,001</td>
    <td style="text-align:right;">0.089</td>
    <td style="text-align:right;">0.316</td>
    <td style="text-align:right;">0.185</td>
    <td style="text-align:right;">0.146</td></tr>
<tr><td style="text-align:right;">80,000,000,000,027</td>
    <td style="text-align:right;">0.112</td>
    <td style="text-align:right;">0.468</td>
    <td style="text-align:right;">0.233</td>
    <td style="text-align:right;">0.217</td></tr>
<tr><td style="text-align:right;">160,000,000,000,069</td>
    <td style="text-align:right;">0.143</td>
    <td style="text-align:right;">0.672</td>
    <td style="text-align:right;">0.336</td>
    <td style="text-align:right;">0.281</td></tr>
<tr><td style="text-align:right;">320,000,000,000,029</td>
    <td style="text-align:right;">0.189</td>
    <td style="text-align:right;">0.854</td>
    <td style="text-align:right;">0.492</td>
    <td style="text-align:right;">0.391</td></tr>
<tr><td style="text-align:right;">640,000,000,000,033</td>
    <td style="text-align:right;">0.253</td>
    <td style="text-align:right;">1.254</td>
    <td style="text-align:right;">0.702</td>
    <td style="text-align:right;">0.567</td></tr>
</tbody>
</table>

<table>
<caption>表５　Chakra エンジンに差し替えた場合</caption>
<thead>
<th rowspan=2 style="text-align:center;">素数</th>
    <th style="text-align:center;">試し割り法</th>
    <th colspan=3 style="text-align:center;">エラトステネスの篩（ふるい）</th></tr>
<tr><th style="text-align:center;">Wheel<br>Factorization</th>
    <th style="text-align:center;">全数サーチ</th>
    <th style="text-align:center;">奇数サーチ</th>
    <th style="text-align:center;">２，３の<br>倍数除去</th></tr>
</thead>
<tbody>
<tr><td style="text-align:right;">1,250,000,000,111</td>
    <td style="text-align:right;">0.054</td>
    <td style="text-align:right;">0.072</td>
    <td style="text-align:right;">0.061</td>
    <td style="text-align:right;">0.056</td></tr>
<tr><td style="text-align:right;">2,500,000,000,009</td>
    <td style="text-align:right;">0.060</td>
    <td style="text-align:right;">0.093</td>
    <td style="text-align:right;">0.067</td>
    <td style="text-align:right;">0.059</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.110</td>
    <td style="text-align:right;">0.077</td>
    <td style="text-align:right;">0.068</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.159</td>
    <td style="text-align:right;">0.102</td>
    <td style="text-align:right;">0.082</td></tr>
<tr><td style="text-align:right;">20,000,000,000,021</td>
    <td style="text-align:right;">0.102</td>
    <td style="text-align:right;">0.223</td>
    <td style="text-align:right;">0.122</td>
    <td style="text-align:right;">0.103</td></tr>
<tr><td style="text-align:right;">40,000,000,000,001</td>
    <td style="text-align:right;">0.133</td>
    <td style="text-align:right;">0.281</td>
    <td style="text-align:right;">0.171</td>
    <td style="text-align:right;">0.126</td></tr>
<tr><td style="text-align:right;">80,000,000,000,027</td>
    <td style="text-align:right;">0.172</td>
    <td style="text-align:right;">0.413</td>
    <td style="text-align:right;">0.242</td>
    <td style="text-align:right;">0.177</td></tr>
<tr><td style="text-align:right;">160,000,000,000,069</td>
    <td style="text-align:right;">0.237</td>
    <td style="text-align:right;">0.639</td>
    <td style="text-align:right;">0.304</td>
    <td style="text-align:right;">0.250</td></tr>
<tr><td style="text-align:right;">320,000,000,000,029</td>
    <td style="text-align:right;">0.316</td>
    <td style="text-align:right;">0.806</td>
    <td style="text-align:right;">0.448</td>
    <td style="text-align:right;">0.316</td></tr>
<tr><td style="text-align:right;">640,000,000,000,033</td>
    <td style="text-align:right;">0.416</td>
    <td style="text-align:right;">1.283</td>
    <td style="text-align:right;">0.691</td>
    <td style="text-align:right;">0.470</td></tr>
</tbody>
</table>

</details>

# 結論
- すっごく苦労してようやく試し割り法（Wheel Factorization）に追いついた。
- ２，３，５の倍数除去までやれば，試し割り法（Wheel Factorization）を逆転できるかもしれないが，コードが複雑になり過ぎる。

# 次回

次回に続きます。

- [Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(4) - Qiita](https://qiita.com/TETSURO1999/items/a8321c80405a70755c41)
- [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)