前回のあらすじ
- 自作の素因数分解プログラムを Qiita で公開した。
- プログラムに添削が入り，30倍以上速くなった。
- なぜ速くなったのか，逆になぜ自作プログラムは遅かったのか考えた。
- そのまま使うのもアレなので添削プログラムを改良して少しだけ速くした。
詳しくは下記の記事を参照されたい。
エラトステネスの篩（ふるい）
素数判定では片っ端から自然数で割って確認してみるのだが（試し割りという），除算は現代の最新プロセッサでも比較的コストのかかる計算であり，試し割りの回数はできるだけ抑制したい。
このため当初の自作プログラムでは，試し割りをする回数を抑えるために試し割りをする自然数の素数判定を行っていたのだが，素数判定ルーチンの中で試し割りをしていたので意味がなかった。というか逆に計算量が大幅に増えてしまっていた。親会社は経営合理化したつもりで，単に子会社に面倒を押し付けただけっていうアレよ。
だがしかし，世の中には「除算」を行わなくても素数の選別ができる方法がある。それがエラトステネスの篩（ふるい）だ。詳しくは脚注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 );
}
長いので全コードはコチラ
//------------------------------------------------------------------------------
// 与えられた自然数 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 よりも少し遅くなるからである。
測定時間の詳細はコチラ
| 素数 | 計算時間（秒） | ||
|---|---|---|---|
| オリジナル | 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 |
| 素数 | 計算時間（秒） | ||
|---|---|---|---|
| オリジナル | 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 |
| 素数 | 計算時間（秒） | ||
|---|---|---|---|
| オリジナル | 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 とは２倍近い差が開いている。
測定時間の詳細はコチラ
| 素数 | 計算時間（秒） | ||
|---|---|---|---|
| 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 |
| 素数 | 計算時間（秒） | ||
|---|---|---|---|
| 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 |
| 素数 | 計算時間（秒） | ||
|---|---|---|---|
| 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 さん，ご指摘ありがとうございました。
誤：エラストテネス
正：エラトステネス
次回
次回に続きます。
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(3) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(4) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(5) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(6) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(7) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(8) - Qiita
- Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(9) 最終回 - Qiita