More than 1 year has passed since last update.

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

Last updated at 2023-11-07Posted at 2023-09-09

前回までのあらすじ

詳しくは下記の記事を参照されたい。
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件 - Qiita)
Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2) - Qiita

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

２，３の倍数を除く場合

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

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

表１　２，３の倍数を除く自然数
$i$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$	$13$
$j$	$1$		$2$		$3$		$4$		$5$		$6$
$k$	$0$	$1$	$0$	$1$	$0$	$1$	$0$	$1$	$0$	$1$	$0$	$1$
$w_k$	$-1$	$+1$	$-1$	$+1$	$-1$	$+1$	$-1$	$+1$	$-1$	$+1$	$-1$	$+1$
$m_{j,k}$	$5$	$7$	$11$	$13$	$17$	$19$	$23$	$25$	$29$	$31$	$35$	$37$

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

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

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$ のときの位置の推移を示す。

\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$ のときの位置の推移を示す。

\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}

以上をまとめると以下の表になる。

表２　消去する合成数の位置 $(J_n,K_n)$ の推移
$k = 0$ のときの位置の推移				$k = 1$ のときの位置の推移
$n$	$J_n$	$K_n$	$\Delta J_n$	$n$	$J_n$	$K_n$	$\Delta J_n$
$0$	$6j^2 - 2j$	$1$	$2j$	$0$	$6j^2 + 2j$	$1$	$4j + 1$
$1$	$6j^2$	$0$	$4j-1$	$1$	$6j^2 + 6j + 1$	$0$	$2j$
$2$	$6j^2 + 4j - 1$	$1$	$2j$	$2$	$6j^2 + 8j + 1$	$1$	$4j + 1$
$3$	$6j^2 + 6j - 1$	$0$	$4j -1$	$3$	$6j^2 + 12j + 2$	$0$	$2j$

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

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

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

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

javascript エンジンを JScript9 に差し替えた場合を以下に示す。奇数版よりは速くなっているが，優位性は僅かである。

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

測定時間の詳細はコチラ

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

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

表５　Chakra エンジンに差し替えた場合
素数	試し割り法	エラトステネスの篩（ふるい）
素数	Wheel Factorization	全数サーチ	奇数サーチ	２，３の倍数除去
1,250,000,000,111	0.054	0.072	0.061	0.056
2,500,000,000,009	0.060	0.093	0.067	0.059
5,000,000,000,053	0.069	0.110	0.077	0.068
10,000,000,000,037	0.084	0.159	0.102	0.082
20,000,000,000,021	0.102	0.223	0.122	0.103
40,000,000,000,001	0.133	0.281	0.171	0.126
80,000,000,000,027	0.172	0.413	0.242	0.177
160,000,000,000,069	0.237	0.639	0.304	0.250
320,000,000,000,029	0.316	0.806	0.448	0.316
640,000,000,000,033	0.416	1.283	0.691	0.470

結論

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

次回

次回に続きます。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up