Go to Qiita Advent Calendar 2024 Top
LoginSignup
54
37

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@TETSURO1999

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

Last updated at Posted at 2023-09-07

前回のあらすじ

  1. 自作の素因数分解プログラムを Qiita で公開した。
  2. プログラムに添削が入り，30倍以上速くなった。
  3. なぜ速くなったのか，逆になぜ自作プログラムは遅かったのか考えた。
  4. そのまま使うのもアレなので添削プログラムを改良して少しだけ速くした。

詳しくは下記の記事を参照されたい。

エラトステネスの篩（ふるい）

素数判定では片っ端から自然数で割って確認してみるのだが（試し割りという），除算は現代の最新プロセッサでも比較的コストのかかる計算であり，試し割りの回数はできるだけ抑制したい。

このため当初の自作プログラムでは，試し割りをする回数を抑えるために試し割りをする自然数の素数判定を行っていたのだが，素数判定ルーチンの中で試し割りをしていたので意味がなかった。というか逆に計算量が大幅に増えてしまっていた。親会社は経営合理化したつもりで，単に子会社に面倒を押し付けただけっていうアレよ。

だがしかし，世の中には「除算」を行わなくても素数の選別ができる方法がある。それがエラトステネスの篩（ふるい）だ。詳しくは脚注1あるいは脚注2などを参照されたい。

素数判定を行う自然数 $N$ とおくと $\sqrt{N}$ までの素数を求めればよい。こうして厳選された素数のみによって試し割りを行えば除算回数は最小化される。また，以下に示すように素数判定テーブルおよび素数リストの計算には高コストな除算は使っていない。

//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list( n ) {
	//--------------------------------------------------------------------------
	// エラトステネスのふるいにより，素数判定テーブルを作成する
	// 素数判定テーブル，true：素数，false：合成数
	//--------------------------------------------------------------------------
	var	table = Array( 1 + Math.floor( Math.sqrt( n ) ) );
	var	size = table.length;
	table[0] = false;	// 0 は素数ではない
	table[1] = false;	// 1 は素数ではない
	for( var i = 2; i < size; i++ ) table[i] = true;
	for( var i = 2; i < size; i++ ) {
		if( !table[i] ) continue;
		for( var j = i * i; j < size; j += i ) table[j] = false;
	}
	//--------------------------------------------------------------------------
	// 素数リストを作成する
	//--------------------------------------------------------------------------
	var	list = [];
	for( var i = 2; i < size; i++ ) if( table[i] ) list.push( i );
	return( list );
}
//------------------------------------------------------------------------------
// 素因数分解を行う
//------------------------------------------------------------------------------
function prime_factorize( n, list ) {
	var	a = [];
	for( var i = 0, p = n; i < list.length && ( q = list[i],  q * q <= p ); i++ ) {
		while( p % q == 0 ) {
			p /= q; a.push( q );
		}
	}
	if( p > 1 ) a.push( p );
	return( a );
}
長いので全コードはコチラ
CHECKPRIME.JS
//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list( n ) {
	//--------------------------------------------------------------------------
	// エラトステネスのふるいにより，素数判定テーブルを作成する
	// 素数判定テーブル，true：素数，false：合成数
	//--------------------------------------------------------------------------
	var	table = Array( 1 + Math.floor( Math.sqrt( n ) ) );
	var	size = table.length;
	table[0] = false;	// 0 は素数ではない
	table[1] = false;	// 1 は素数ではない
	for( var i = 2; i < size; i++ ) table[i] = true;
	for( var i = 2; i < size; i++ ) {
		if( !table[i] ) continue;
		for( var j = i * i; j < size; j += i ) table[j] = false;
	}
	//--------------------------------------------------------------------------
	// 素数リストを作成する
	//--------------------------------------------------------------------------
	var	list = [];
	for( var i = 2; i < size; i++ ) if( table[i] ) list.push( i );
	return( list );
}
//------------------------------------------------------------------------------
// 素因数分解を行う
//------------------------------------------------------------------------------
function prime_factorize( n, list ) {
	var	a = [];
	for( var i = 0, p = n; i < list.length && ( q = list[i],  q * q <= p ); i++ ) {
		while( p % q == 0 ) {
			p /= q; a.push( q );
		}
	}
	if( p > 1 ) a.push( p );
	return( a );
}
//------------------------------------------------------------------------------
// メイン関数
//------------------------------------------------------------------------------
function main( args ) {
	//--------------------------------------------------------------------------
	// オプション解析
	//--------------------------------------------------------------------------
	if( args.Count == 0 ) {
		WScript.Echo( "素数のチェック，合成数の場合は素因数分解を行います。" );
		WScript.Echo( "" );
		WScript.Echo( "CHECKPRIME(.JS) [整数]" );
		WScript.Echo( "" );
		WScript.Echo( "素数の場合は 0，合成数の場合は 1 を返します。" );
		return( -1 );
	}
	var	n = parseInt( args(0) );
	if( n < 2 ) {
		WScript.Echo( "2 以上の整数を入力して下さい!!" );
		return( -1 );
	}
	//--------------------------------------------------------------------------
	// 素因数分解を行う
	//--------------------------------------------------------------------------
	var	list = create_list( n );
	var	a = prime_factorize( n, list );
	if( a.length == 1 ) {
		WScript.Echo( "整数 " + n + " は素数です。" );
		return( 0 );
	} else {
		WScript.Echo( "整数 " + n + " は合成数です。" );
		WScript.Echo( a.join(",") + " を素因数に持ちます。" );
		return( 1 );
	}
}
var	ret = main( WScript.Arguments.Unnamed );
try {
	WScript.Quit( ret );
} catch(e) {
	/* 何もしない */
}

計算時間の実測結果を以下に示す。

あれれ？？おかしいぞ？？めっちゃ遅いぞ！？

ちなみに素数の大きさが $10^{14}$ を超えると急激に遅くなるのは素数の大きさの平方根に比例して使用メモリ容量が増大し，キャッシュから溢れたためだと思われる。※根拠は後で述べる。

ちなみに javascript エンジンを JScript9 に入れ替えると下記の通り。全般的に速くなるが，相対的な関係は変わらない。またメモリ管理方法が改善されたためか，急激な低速化は見られない。上記のコードで用いている素数判定テーブルは bool 型で保持すれば十分キャッシュ内に収まるはずだが，おそらくデフォルトの javascript エンジンでは double 型で保持しているものと考える。

さらに Chakra だと下記の通り。差が縮まったように見えるのは何故か Chakra だと Wheel Factorization が JScript9 よりも少し遅くなるからである。

測定時間の詳細はコチラ
表１　デフォルトの JScript エンジンの場合
素数 計算時間（秒）
オリジナル Wheel
Factorization		 エラトステネス
の篩（ふるい）
1,250,000,000,111 3.940 0.098 0.560
2,500,000,000,009 6.050 0.117 0.730
5,000,000,000,053 9.460 0.143 1.010
10,000,000,000,037 14.930 0.179 1.460
20,000,000,000,021 22.790 0.234 2.040
40,000,000,000,001 36.020 0.307 2.940
80,000,000,000,027 56.620 0.411 6.680
160,000,000,000,069 91.540 0.562 32.850
320,000,000,000,029 146.120 0.775 92.640
640,000,000,000,033 234.660 1.078 221.080
表２　JScript9 エンジンに差し替えた場合
素数 計算時間（秒）
オリジナル Wheel
Factorization		 エラトステネス
の篩（ふるい）
1,250,000,000,111 0.230 0.044 0.078
2,500,000,000,009 0.280 0.047 0.085
5,000,000,000,053 0.420 0.053 0.112
10,000,000,000,037 0.620 0.062 0.175
20,000,000,000,021 0.930 0.073 0.216
40,000,000,000,001 1.480 0.089 0.318
80,000,000,000,027 2.380 0.112 0.463
160,000,000,000,069 3.720 0.143 0.668
320,000,000,000,029 5.660 0.189 0.862
640,000,000,000,033 8.720 0.253 1.251
表３　Chakra エンジンに差し替えた場合
素数 計算時間（秒）
オリジナル Wheel
Factorization		 エラトステネス
の篩（ふるい）
1,250,000,000,111 0.180 0.054 0.075
2,500,000,000,009 0.210 0.060 0.091
5,000,000,000,053 0.310 0.069 0.110
10,000,000,000,037 0.480 0.084 0.156
20,000,000,000,021 0.720 0.102 0.223
40,000,000,000,001 1.140 0.133 0.276
80,000,000,000,027 1.770 0.172 0.418
160,000,000,000,069 2.760 0.237 0.672
320,000,000,000,029 4.210 0.316 0.820
640,000,000,000,033 6.890 0.416 1.273

現代のプロセッサは除算が高速化された一方，メモリアクセスは相対的に遅くなっており，エラトステネスの篩（ふるい）のような大量のメモリを使用するアルゴリズムは苦手になっているのかもしれない。

であれば使用するメモリ容量を減らす工夫をしてみる。

メモリの書き込み回数について
本記事で計測している素数の大きさ $N = 10^{12} \thicksim 10^{15}$ なので，エラトステネスの篩（ふるい）を行っている範囲は，その平方根である $\sqrt{N} = 10^6 \thicksim 10^7$ になる。このあたりの素数密度は素数定理を用いて計算すると $0.06 \thicksim 0.07$ 程度しかないので，素数判定テーブルへの書き込みは初期化時を除いて $\sqrt{N} \times 0.93 \thicksim 0.94$ 回行っているものと考えられた。しかし，実際に書き込みは重複して行われており，書き込み回数をカウントしてみると $\sqrt{N} \times 2.2$ 回と書き込み密度が $2$ を超えてしまった。これは考えてみれば当然の話で $\sqrt{N}$ 以下の範囲に $2^2$ 以上の $2$ の倍数は約 $\sqrt{N}/2$ 個，$3^2$ 以上の $3$ の倍数も同様に約 $\sqrt{N}/3$ 個あると考えられる。こうして素数判定テーブルへの書き込み回数を書き込み密度係数 $\mu$ を用いて表すと以下のようになる。

\frac{\sqrt{N}}{2} + \frac{\sqrt{N}}{3} + \frac{\sqrt{N}}{5} + \frac{\sqrt{N}}{7} + \frac{\sqrt{N}}{11} + \cdots + \frac{\sqrt{N}}{p_{max}} \approx \mu \times \sqrt{N}

※$p_{max}$ は $\sqrt{N}$ 以下の最大の素数とする。

$\sqrt{N}$ が $1000$ を超えた時点で $\mu > 2$ となるが，そこからの上昇は極めて緩やかである。なお，$\sqrt{N} \rightarrow \infty$ とすると $\mu \rightarrow \infty$ となることが知られている。

素数判定テーブルを奇数だけにする

$j$ 番目の奇数 $m_j = 2j + 1$ とおく。

$m_j^2$ 以上の $m_j$ の倍数を合成数として素数判定テーブルのフラグを消していく。

初回位置 $m_j^2 = m_{J_0}$ とおく。

\begin{align}
m_j^2 &= 4j^2 + 4j + 1 \\
&= 2\lbrace 2j\cdot(j + 1) \rbrace + 1 \\
&= 2J_0 + 1
\end{align}

より $J_0 = 2j\cdot(j + 1)$ 番目の奇数である。

２番目の位置 $m_j \cdot m_{j+1} = m_{J_1}$ とおく。

\begin{align}
m_j \cdot m_{j+1} &= (2j + 1)\lbrace 2j \cdot(j + 1) + 1\rbrace \\
&= (2j + 1)(2j + 3) \\
&= 4j^2 + 8j + 3 \\
&= 2(2j^2 + 4j + 1) + 1 \\
&= 2J_1 + 1
\end{align}

３番目の位置 $m_j \cdot m_{j+2} = m_{J_2}$ とおく。

\begin{align}
m_j \cdot m_{j+1} &= (2j + 1)\lbrace 2j \cdot(j + 2) + 1\rbrace \\
&= (2j + 1)(2j + 5) \\
&= 4j^2 + 12j + 5 \\
&= 2(2j^2 + 6j + 2) + 1 \\
&= 2J_2 + 1
\end{align}

以上より合成数の位置の移動幅 $\Delta j$ は

\Delta j = J_2 - J_1 = J_1 - J_0 = 2j + 1 = m_j

一定となる。

具体的に考えると分かり易い。
奇数版のエラトステネスの篩（ふるい）では，$j = 1$ 番目の素数 $m_j = 3$ に対し，$3^2$ 以上の $3$ の倍数，すなわち $3 \cdot 3，3 \cdot 5，3 \cdot 7 \cdots$ を合成数として消していくが，奇数テーブル内において $3 \cdot 3 $ の位置は $2j\cdot(j+1)=4$ 番目にある。続く $3 \cdot 5$ の位置は $m_j = 3$ だけ離れた $7$ 番目にある。以降，$3$ の倍数は $m_j = 3$ だけ離れた位置（等間隔）にある。

表４　奇数テーブル
$\!j$ $\!0$ $\!1$ $\!2$ $\!3$ $\!4$ $\!5$ $\!6$ $\!7$ $\!8$ $\!9$ $\!10$ $\!11$ $\!12$ $\!13$
$\!m_j$ $\!1$ $\!3$ $\!5$ $\!7$ $\!9$ $\!11$ $\!13$ $\!15$ $\!17$ $\!19$ $\!21$ $\!23$ $\!25$ $\!27$

この結果，以下のようにコードをあまり複雑にしないで使用するメモリ容量を半減できる。

//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list( n ) {
	//--------------------------------------------------------------------------
	// エラトステネスのふるいにより，素数判定テーブル（奇数のみ）を作成する
	// 素数判定テーブル，true：素数，false：合成数
	//--------------------------------------------------------------------------
	var	table = Array( 1 + Math.floor( Math.sqrt( n ) / 2 ) );
	var	size = table.length;
	table[0] = false;	// 1 は素数ではない
	for( var i = 1; i < size; i++ ) table[i] = true;
	for( var i = 1; i < size; i++ ) {
		if( !table[i] ) continue;
		var	m = 2 * i + 1;
		for( var j = 2 * i* ( i + 1 ); j < size; j += m ) table[j] = false;
	}
	//--------------------------------------------------------------------------
	// 素数リストを作成する
	//--------------------------------------------------------------------------
	var	list = [];
	list.push( 2 );
	for( var i = 1; i < size; i++ ) if( table[i] ) list.push( 2 * i + 1 );
	return( list );
}

計算時間の実測結果を以下に示す。

デフォルトの javascript エンジンの場合は以下の通り。奇数版の効果を確認できた。また，使用するメモリ容量が半減したため，計算時間が急増する素数の大きさが４倍になっている。※使用するメモリ容量は素数の大きさの平方根に比例するため。

javascript エンジンを JScript9 に入れ替えた場合は以下の通り。全般的に速くなっているが，相対的な関係は変わらない。また，計算時間の急増が無くなっている。

javascript エンジンを Chakra に入れ替えた場合は以下の通り。何故か Chakra だと Wheel Factorization が JScript9 よりも少し遅くなるため，差が縮まったように見える。しかし，それでもまだ Wheel Factorization とは２倍近い差が開いている。

測定時間の詳細はコチラ
表５　デフォルトの JScript エンジンの場合
素数 計算時間（秒）
Wheel
Factorization		 エラトステネス
の篩		 エラトステネス
の篩（奇数版）
1,250,000,000,111 0.098 0.560 0.310
2,500,000,000,009 0.117 0.730 0.430
5,000,000,000,053 0.143 1.010 0.550
10,000,000,000,037 0.179 1.460 0.810
20,000,000,000,021 0.234 2.040 1.100
40,000,000,000,001 0.307 2.940 1.530
80,000,000,000,027 0.411 6.680 2.160
160,000,000,000,069 0.562 32.850 3.110
320,000,000,000,029 0.775 92.640 6.840
640,000,000,000,033 1.078 221.080 32.920
>表６　JScript9 エンジンに差し替えた場合
素数 計算時間（秒）
Wheel
Factorization		 エラトステネス
の篩		 エラトステネス
の篩（奇数版）
>1,250,000,000,111 0.044 0.072 0.055
2,500,000,000,009 0.047 0.084 0.061
5,000,000,000,053 0.053 0.112 0.077
10,000,000,000,037 0.062 0.174 0.088
20,000,000,000,021 0.073 0.217 0.119
40,000,000,000,001 0.089 0.316 0.185
80,000,000,000,027 0.112 0.468 0.233
160,000,000,000,069 0.143 0.672 0.336
320,000,000,000,029 0.189 0.854 0.492
640,000,000,000,033 0.253 1.254 0.702
>表７　Chakra エンジンに差し替えた場合
素数 計算時間（秒）
Wheel
Factorization		 エラトステネス
の篩		 エラトステネス
の篩（奇数版）
1,250,000,000,111 0.054 0.072 0.061
2,500,000,000,009 0.060 0.093 0.067
5,000,000,000,053 0.069 0.110 0.077
10,000,000,000,037 0.084 0.159 0.102
20,000,000,000,021 0.102 0.223 0.122
40,000,000,000,001 0.133 0.281 0.171
80,000,000,000,027 0.172 0.413 0.242
160,000,000,000,069 0.237 0.639 0.304
320,000,000,000,029 0.316 0.806 0.448
640,000,000,000,033 0.416 1.283 0.691

結論

なかなか試し割り法に勝てないな・・・
あと２倍速くなれば Wheel Factorization に追いつけるのに・・・

謝辞

@sio-funmatsu さん，ご指摘ありがとうございました。
誤：エラストテネス
正：エラトステネス

次回

次回に続きます。

  1. エラトステネスの篩の高速化 (Qiita) ↩

  2. エラトステネスの篩の活用法を総特集！ 〜 高速素因数分解・メビウスの反転公式 〜 (Qita) ↩

54
37
13

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
54
37

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?