Go to Qiita Advent Calendar 2024 Top
LoginSignup
20
21

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@TETSURO1999

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

Last updated at Posted at 2023-09-01

はじめに

以下の自分の記事で自作の素数判定（素因数分解）プログラムを公開したのだが，こうしたほうがもっと速いんじゃね？という添削指導が入った。

【京大2021】3ⁿ－2ⁿが素数ならnも素数になることを示せ - Qiita
【京大2021】pが素数ならばp⁴＋14が素数ではないことを示せ - Qiita

まずは自作のコード（一部抜粋）を示す。
このプログラムのキモは，与えられた自然数の約数を検索する search_divisor() 関数である。与えられた自然数に対し既知の素数リスト内に割り切れるものがあればそれが約数であり，ここで検索が終了するが，見つからなかった場合，素数リストの最大値（奇数）に２を加えた奇数を新たな約数の候補とし，この候補が素数であることを search_divisor() を呼び出してチェックする。もしも約数候補が素数でなければ，試し割りを行わない。素数であれば素数リストに末尾に追加して，与えられた自然数が割り切れるかどうか確認する。割り切れればそれが約数であり，検索を終了する。割り切れなければ，さらに２を加えた奇数を新たな約数の候補として繰り返す。

//------------------------------------------------------------------------------
// 100 以下の素数リスト
//------------------------------------------------------------------------------
var	list = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97];
//------------------------------------------------------------------------------
// 約数を検索する
// ・素数の場合は 1 を返す
// ・合成数の場合は 1 より大きい最小の約数を返す
//------------------------------------------------------------------------------
function search_divisor( n ) {
	var	m;
	//--------------------------------------------------------------------------
	// 素数リストをチェックする
	//--------------------------------------------------------------------------
	for( var i = 0; i < list.length; i++ ) {
		m = list[i];
		if( m * m > n ) return( 1 );
		if( n % m == 0 ) return( m );
	}
	//--------------------------------------------------------------------------
	// 約数の候補を探す
	//--------------------------------------------------------------------------
	for(;;) {
		m += 2;
		if( m * m > n ) return( 1 );
		if( search_divisor( m ) != 1 ) continue;
		list.push( m );
		if( n % m == 0 ) return( m );
	}
}
//------------------------------------------------------------------------------
// 素因数分解を行う
//------------------------------------------------------------------------------
function prime_factorize( p ) {
	var	a = [];
	for( var q; ( q = search_divisor( p ) ) > 1; p = p / q )
		a.push( q );
	a.push( p );
	return( a );
}
長いので全コードはコチラ WSH（Windows Script Host）で動かす javascript プログラムなので Windows ユーザなら誰でもコマンドプロンプト上で動かせる。
checkprime.js
//------------------------------------------------------------------------------
// 素数チェックおよび素因数分解を行います。
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// 100 以下の素数リスト
//------------------------------------------------------------------------------
var	list = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97];
//------------------------------------------------------------------------------
// 約数を検索する
// ・素数の場合は 1 を返す
// ・合成数の場合は最小の約数を返す
//------------------------------------------------------------------------------
function search_divisor( n ) {
	var	m;
	//--------------------------------------------------------------------------
	// 素数リストをチェックする
	//--------------------------------------------------------------------------
	for( var i = 0; i < list.length; i++ ) {
		m = list[i];
		if( m * m > n ) return( 1 );
		if( n % m == 0 ) return( m );
	}
	//--------------------------------------------------------------------------
	// 約数の候補を探す
	//--------------------------------------------------------------------------
	for(;;) {
		m += 2;
		if( m * m > n ) return( 1 );
		if( search_divisor( m ) != 1 ) continue;
		list.push( m );
		if( n % m == 0 ) return( m );
	}
}
//------------------------------------------------------------------------------
// 素因数分解を行う
//------------------------------------------------------------------------------
function prime_factorize( p ) {
	var	a = [];
	for( var q; ( q = search_divisor( p ) ) > 1; p = p / q )
		a.push( q );
	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	a = prime_factorize( n );
	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) {
	/* 何もしない */
}

どこが遅いのか？

このプログラムはもともと大学入試問題で取り扱う程度の素数を判定（もしくは素因数分解）するために作られたので $100$ 以下の素数リストを内蔵しており，$100^2 = 10000$ 以下の自然数までならば高速に素数判定できる。

しかし，これらを超えてしまうと一気に計算時間がかかるようになる。計算時間を実測した結果を以下に示す。素数の大きさが1000倍になると計算時間は100倍になり，綺麗に $O(N^{2/3})$ のラインに乗っている。普通1に作っても $O(N^{1/2})$ なのでかなり頭の悪いアルゴリズムといえよう。

測定時間の詳細はコチラ
表１　素数判定に要した時間
素数 時間（秒）
1,250,000,000,111 3.940
2,500,000,000,009 6.050
5,000,000,000,053 9.460
10,000,000,000,037 14.930
20,000,000,000,021 22.790
40,000,000,000,001 36.020
80,000,000,000,027 56.620
160,000,000,000,069 91.540
320,000,000,000,029 146.120
640,000,000,000,033 234.660

さて，どこが問題点なのかというと，試し割りの回数を節約するため，約数候補の素数チェックを行っていることである。素数チェックの中では約数候補 $m$ が素数であることを確認するため素数リストにある $\sqrt{m}$ 以下の素数で試し割りを行っており，すべて割り切れなかった場合に $m$ が素数と判断される。

素数定理2より，$\sqrt{m}$ 以下の素数の数はおよそ $\sqrt{m}/\log{\sqrt{m}}$ であるから，試し割りを１回節約するつもりが，逆に大幅に増えていたことになる。

具体的な添削内容

具体的な添削内容を以下に示す。元のコードに比べて随分とすっきりしている。

これは約数の候補 $q$ が素数であろうと合成数であろうと気にせず，どんどん試し割りしていこうという方式だ。さすがに $4$ とか $6$ とか偶数でも試し割りするのは如何なものかと思うが，コチラのほうが桁違いに速くなる。

なおオリジナルの添削コードではループ条件を q <= Math.sqrt(p) としていたが，時間のかかる平方根計算を避けて q * q <= p と修正した。※ちょっとだけ速くなった。

//------------------------------------------------------------------------------
// 素因数分解を行う（添削指導後）
//------------------------------------------------------------------------------
function prime_factorize( p ) {
	var	a = [];
	for( var q = 2; q * q <= p; q++ ) {
		while( p % q == 0 ) {
			p /= q; a.push( q );
		}
	}
	if( p > 1 ) a.push( p );
	return( a );
}

プログラムがすっきりすると改善点が分かり易くなる。先に $2$ で試し割りを終えて，その後は奇数のみ試し割りを行うようにするとさらに２倍速くなる。

//------------------------------------------------------------------------------
// 素因数分解を行う（添削指導後＋奇数サーチ）
//------------------------------------------------------------------------------
function prime_factorize( p ) {
	var	a = [];
	while( p >= 2 && p % 2 == 0 ) {
		p /= 2; a.push( 2 );
	}
	for( var q = 3; q * q <= p; q += 2 ) {
		while( p % q == 0 ) {
			p /= q; a.push( q );
		}
	}
	if( p > 1 ) a.push( p );
	return( a );
}

計算時間の比較を示す。添削指導後および奇数サーチは桁違いに速くなっているほか，計算量も理論通り $O(N^{1/2})$ となっていることが分かる。

計算時間の比較詳細はコチラ
表２　計算時間の比較
素数 時間（秒）
オリジナル 添削指導後 奇数サーチ
1,250,000,000,111 3.940 0.146 0.100
2,500,000,000,009 6.050 0.180 0.117
5,000,000,000,053 9.460 0.233 0.145
10,000,000,000,037 14.930 0.305 0.183
20,000,000,000,021 22.790 0.408 0.236
40,000,000,000,001 36.020 0.554 0.311
80,000,000,000,027 56.620 0.764 0.420
160,000,000,000,069 91.540 1.061 0.570
320,000,000,000,029 146.120 1.480 0.785
640,000,000,000,033 234.660 2.075 1.094

さらに改善するには？

Wheel Factorization 法3を使えば，奇数サーチのさらに２倍速くなる。ただ，WSH の javascript エンジンはアホの子なので Wheel Factorization 法を使っても速くはならなかった。奇数サーチより僅かに速い程度。ただし，JScript9，Chakra など JIT コンパイラ型のエンジンに切り替えると高速化の効果がはっきり見える。

//------------------------------------------------------------------------------
// 素因数分解を行う（Wheel Factorization）
//------------------------------------------------------------------------------
function prime_factorize( p ) {
	var	a = [];
	while( p >= 2 && p % 2 == 0 ) {
		p /= 2; a.push( 2 );
	}
	while( p >= 3 && p % 3 == 0 ) {
		p /= 3; a.push( 3 );
	}
	while( p >= 5 && p % 5 == 0 ) {
		p /= 5; a.push( 5 );
	}
	var	wheel = [4,2,4,2,4,6,2,6];
	for( var q = 7, i = 0; q * q <= p; q += wheel[i++], i %= 8 ) {
		while( p % q == 0 ) {
			p /= q; a.push( q );
		}
	}
	if( p > 1 ) a.push( p );
	return( a );
}

一応，計算結果を載せておく。

Wheel Factorization の比較詳細はコチラ
表３　Wheel Factorization の比較詳細
素数 時間（秒）
オリジナル 添削指導後 奇数サーチ Wheel
Factorization
2,500,000,000,009 6.050 0.180 0.117 0.117
5,000,000,000,053 9.460 0.233 0.145 0.143
10,000,000,000,037 14.930 0.305 0.183 0.179
20,000,000,000,021 22.790 0.408 0.236 0.234
40,000,000,000,001 36.020 0.554 0.311 0.307
80,000,000,000,027 56.620 0.764 0.420 0.411
160,000,000,000,069 91.540 1.061 0.570 0.562
320,000,000,000,029 146.120 1.480 0.785 0.775
640,000,000,000,033 234.660 2.075 1.094 1.078

JScript9 の結果

WSH の javascript エンジンを JScript9 に切り替えた場合の結果を示す。Wheel Factorization の効果が見えるようになった。なお，素数の大きさが $10^{13}$ を下回ると JIT コンパイラのコンパイル時間によるオーバーヘッドが見えてくる。

JScript9 の結果詳細はコチラ
表４　JScript9 の結果
素数 時間（秒）
オリジナル 添削指導後 奇数サーチ Wheel
Factorization
1,250,000,000,111 0.230 0.065 0.049 0.044
2,500,000,000,009 0.280 0.075 0.054 0.047
5,000,000,000,053 0.420 0.093 0.063 0.053
10,000,000,000,037 0.620 0.115 0.075 0.062
20,000,000,000,021 0.930 0.149 0.090 0.073
40,000,000,000,001 1.480 0.195 0.114 0.089
80,000,000,000,027 2.380 0.260 0.146 0.112
160,000,000,000,069 3.720 0.352 0.193 0.143
320,000,000,000,029 5.660 0.483 0.258 0.189
640,000,000,000,033 8.720 0.668 0.349 0.253

Chakra の結果

WSH の javascript エンジンを Chakra に切り替えた場合の結果を示す。Wheel Factorization の効果がよりはっきりと見えるようになった。なお，JIT コンパイラのオーバーヘッドは JScript9 よりも大きい。時間をかけて最適化を図っているものと考えられる。

Chakra の結果詳細はコチラ
表５　Chakra の結果
素数 時間（秒）
オリジナル 添削指導後 奇数サーチ Wheel
Factorization
1,250,000,000,111 0.180 0.122 0.079 0.054
2,500,000,000,009 0.210 0.157 0.097 0.060
5,000,000,000,053 0.310 0.208 0.122 0.069
10,000,000,000,037 0.480 0.278 0.157 0.084
20,000,000,000,021 0.720 0.376 0.208 0.102
40,000,000,000,001 1.140 0.517 0.279 0.133
80,000,000,000,027 1.770 0.715 0.376 0.172
160,000,000,000,069 2.760 0.995 0.516 0.237
320,000,000,000,029 4.210 1.390 0.712 0.316
640,000,000,000,033 6.890 1.950 0.994 0.416

備考

WSH の javascript エンジンを JScript9 や Chakra に切り替える方法は下記記事を参照のこと。

【京大2021】3ⁿ－2ⁿが素数ならnも素数になることを示せ - Qiita

ベンチマーク環境は下記の通り

  • OS: Windows10 Pro 2022H2
  • CPU: Core i5 10400 2.9GHz
  • RAM: 8GB

次回

次回に続きます。

謝辞

@fujitanozomu （藤田望）さんご指導ありがとうございました。

  1. 素数判定アルゴリズムいろいろ - Qiita ↩

  2. 素数定理 - Wikipedia ↩

  3. Wheel Factorizationについて - Qiita ↩

20
21
3

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
20
21

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?