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

前回までのあらすじ

自作の素因数分解プログラムを Qiita で公開した。プログラムに添削が入り，30倍以上速くなった。ちなみに教えて頂いた方法は「試し割り法」と呼ばれている。教えられた方法をそのまま使うのもアレなので，もっと速い方法を編み出せないか考えてみた。

現代のプロセッサでも割り算は比較的重いとされるので，与えられた自然数 $n$ に対してエラトステネスの篩（ふるい）と呼ばれる手法を用いて $\sqrt{n}$ 以下の素数を列挙し，それらの素数だけで試し割りをすれば除算回数は最小化され，格段の高速化が期待できるはずだ！しかし，エラトステネスの篩（ふるい）は大量にメモリを使用するせいか，意図したほど速くならなかった。

一方，試し割り法は $\sqrt{n}$ 以下の自然数すべてで割るのではなく，奇数だけ，あるいは２・３の倍数を除いた場合，２・３・５の倍数を除いた場合など除算回数を減らす手法（Wheel Factorization）が知られている。Wheel Factorization は比較的に簡単に実装でき，実際に効果を確認できた。

同様にエラトステネスの篩（ふるい）でも，$\sqrt{n}$ 以下の自然数すべてを配列に確保するのではなく，奇数だけ，あるいは２・３の倍数を除く，さらに２・３・５の倍数を除けば使用メモリを低減でき，高速化が期待できる。しかし，Wheel Factorization と異なり，極めてコーディングが複雑になるため，２・３の倍数を除く場合のところでストップしていた。要するに２・３・５の倍数を除く場合のコードが複雑になり過ぎて手に負えなかった。

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

２・３の倍数を除く方法（改）

２・３・５の倍数を除く方法をどうやって実装したらよいのか考えている途中で思いついたアイディアをまず２・３の倍数を除く方法で試してみる。

非等間隔な数列

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

$\textsf{表１　２・３の倍数を除いた場合}$
$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$
$13$	$14$	$15$	$16$	$17$	$18$	$19$	$20$	$21$	$22$	$23$	$24$
$25$	$26$	$27$	$28$	$29$	$30$	$31$	$32$	$33$	$34$	$35$	$36$
$37$	$38$	$39$	$40$	$41$	$42$	$43$	$44$	$45$	$46$	$47$	$48$

表１を整理すると以下のようになる。この等間隔に並んでいない数列に対して，素数または合成数のマーキングを行っていかなくてはならない。

$\textsf{表２　２・３の倍数を除いた場合}$
$1$	$5$	$7$	$11$	$13$	$17$	$19$	$23$	$25$	$29$	$31$	$35$
$37$	$41$	$43$	$47$	$49$	$53$	$55$	$59$	$61$	$65$	$67$	$71$
$73$	$77$	$79$	$83$	$85$	$89$	$91$	$95$	$97$	$101$	$103$	$107$
$109$	$113$	$115$	$119$	$121$	$125$	$127$	$131$	$133$	$137$	$139$	$143$

エラトステネスの篩（ふるい）では，たとえば最初に素数 $5$ を見つけたとして，$5^2$ 以上の $5$ の倍数，すなわちこのケースでは $5\cdot 5, 5\cdot 7, 5\cdot 11, \cdots$（すなわち $5$ と $5$ 以上の要素の積）を合成数としてマーキングしていくが，$5$ の倍数の位置は非等間隔に並んでおり，非常に取り扱いにくい。同様に $7, 11, 13$ の倍数の位置を以下に示す。

表３　$5, 7, 11, 13$ の倍数の位置（1-origin）

$5$ の倍数			$7$ の倍数			$11$ の倍数			$13$ の倍数
合成数	位置		合成数	位置		合成数	位置		合成数	位置
$5 \cdot 5$	$9$		$7 \cdot 7$	$17$		$11 \cdot 11$	$41$		$13 \cdot 13$	$57$
$5 \cdot 7$	$12$		$7 \cdot 11$	$26$		$11 \cdot 13$	$48$		$13 \cdot 17$	$74$
$5 \cdot 11$	$19$		$7 \cdot 13$	$31$		$11 \cdot 17$	$63$		$13 \cdot 19$	$83$
$5 \cdot 13$	$22$		$7 \cdot 17$	$40$		$11 \cdot 19$	$70$		$13 \cdot 23$	$100$
$5 \cdot 17$	$29$		$7 \cdot 19$	$45$		$11 \cdot 23$	$85$		$13 \cdot 25$	$109$
$5 \cdot 19$	$32$		$7 \cdot 23$	$54$		$11 \cdot 25$	$92$		$13 \cdot 29$	$126$
$5 \cdot 23$	$39$		$7 \cdot 25$	$59$		$11 \cdot 29$	$107$		$13 \cdot 31$	$135$
$5 \cdot 25$	$42$		$7 \cdot 29$	$68$		$11 \cdot 31$	$114$		$13 \cdot 35$	$152$

ここで位置を一つおきに見ると $5$ の倍数は $10$ 間隔，$7$ の倍数は $14$ 間隔，$11$ の倍数は $22$ 間隔，$13$ の倍数は $26$ 間隔となってることが分かる。すなわち $p$ の倍数は $2p$ 間隔になりそうである。

定式化

２，３の倍数を除く自然数の数列を $m(i) = 6j + w(k)$ とおく。
$i = 2j + k$，$j \ge 0$，$k = \lbrace 0, 1 \rbrace$，$w(k) = \lbrace -1, +1 \rbrace$ とする。

素数 $m(i)$ に対して，$m(i)^2 = m(I_0)$ となる初期位置 $I_0$ を求める。

\begin{align}
m(i)^2 &= (6j + w(k))^2 \\
&= 36j^2 + 12j\cdot w(k) + w(k)^2 \\
&= 6 \left( 6j^2 + 2j\cdot w(k) \right) + 1 \\
&= 6 J_0 + w(K_0) \\
J_0 &= 6j^2 + 2j\cdot w(k) \\
K_0 &= 1 \\
I_0 &= 2J_0 + K_0 = 12j^2 + 4j\cdot w(k) + 1
\end{align}

一つ置いた位置 $I_2$ を求める。

\begin{align}
m(i)\cdot m(i+2) &= (6j + w(k))(6(j + 1) + w(k)) \\
&= 36j\cdot(j+1) + 6j\cdot w(k) + 6(j+1) \cdot w(k) + w(k)^2 \\
&= 6 \left( 6j^2 + 6j + 2j\cdot w(k) + w(k) \right) + 1 \\
&= 6 J_2 + w(K_2) \\
J_2 &= 6j^2 + 6j + 2j\cdot w(k) + w(k) \\
K_2 &= 1 \\
I_2 &= 2J_2 + K_2 = 12j^2 + 12j + 4j\cdot w(k) + 2w(k) + 1

\end{align}

$I_2 - I_0 = 12j + 2w(k) = 2(6j + w(k)) = 2m(i)$ となる。

調べてみると非負の整数 $x \ge 0$ に対して

$I_{x + 2} - I_x =2m(i)$

が成立する。すなわち，初期配置から $I_0 \rightarrow I_1 \rightarrow I_2 \cdots $ と一つずつ移動しながらマーキングするのではなく，

\begin{alignat}{3}
I_0 &\rightarrow I_2 &\rightarrow I_4 \cdots \\
I_1 &\rightarrow I_3 &\rightarrow I_5 \cdots
\end{alignat}

と一つおきにマーキングしていけば良い。

参考までに位置 $I_1$ を求める。

$k = 0$ のとき

\begin{align}
m(i)\cdot m(i+1) &= (6j + w(0))(6j + w(1)) \\
&= 36j^2 + 6j \cdot \lbrace w(0) + w(1) \rbrace + w(0)\cdot w(1) \\
&= 6 \cdot 6j^2 - 1 \\
&= 6 J_1 + w(K_1) \\
J_1 &= 6j^2  \\
K_1 &= 0 \\
I_1 &= 2J_1 + K_1 = 12j^2
\end{align}

$k = 1$ のとき

\begin{align}
m(i)\cdot m(i+1) &= (6j + w(1))\lbrace6(j + 1) + w(0)\rbrace \\
&= 36j\cdot(j+1) + 6j \cdot w(0) + 6(j+1)\cdot w(1) + w(0)\cdot w(1) \\
&= 6 \cdot (6j^2 + 6j + 1 ) - 1 \\
&= 6 J_1 + w(K_1) \\
J_1 &= 6j^2 + 6j + 1 \\
K_1 &= 0 \\
I_1 &= 2J_1 + K_1 = 12j^2 + 12j + 2
\end{align}

となる。これらを統合すると

I_1 = 12j^2 + 2\cdot(6j + 1)\cdot k

となる。

改善版のコード

改善版のコードを以下に示す。既存のコードと比較してみると，内側の for() ループの中身が極めて簡素化されていることが分かる。使用するメモリ容量は変わらないが，１バイト書き込みあたりの処理は激減しているはずだ。

function create_list(n) {
	var n1 = Math.floor(Math.sqrt(n));
	var n2 = (n1 % 2 == 0) ? n1 - 1 : n1;
	var	n3 = (n2 % 3 == 0) ? n2 - 2 : n2;
	var	table = new Uint8Array(2 + Math.floor(n3 / 3));
	var	size = table.length;
	var	list = [];
	list.push(2);
	list.push(3);
	for(var i = 2; i < size; i++) {
		if(table[i] != 0) continue;
		var	j = i >>> 1;
		var	k = i & 1;
		var	m = 6 * j + 2 * k - 1;
		list.push(m);
		for(var i0 = 12 * j * j + 4 * j * (2 * k - 1) + 1,
				i1 = 12 * j * j + 2 * (6 * j + 1) * k;
			i1 < size; i0 += 2 * m, i1 += 2 * m) {
			table[i0] = 1;
			table[i1] = 1;
		}
		if(i0 < size) table[i0] = 1;
	}
	return list;
}

従来版のコード

参考までに従来版のコードを示す。

function create_list(n) {
    var n1 = Math.floor(Math.sqrt(n));
    var n2 = (n1 % 2 == 0) ? n1 - 1 : n1;
	var	n3 = (n2 % 3 == 0) ? n2 - 2 : n2;
	var	table = new Uint8Array(2 + Math.floor(n3 / 3));
	var	size = table.length;
	var	list = [];
	list.push(2);
	list.push(3);
	for(var i = 2; i < size; i++) {
		if( table[i] != 0 ) continue;
		var	j = i >>> 1;
		var k = i & 1;
		var	m = 6 * j + 2 * k - 1;
		list.push(m);
		for(var jj = 6 * j * j + 2 * j * (2 * k - 1), kk = 1;
			(ii = 2 * jj + kk) < size; kk ^= 1) {
			table[ii] = 1;
			jj += (kk == k) ? 4 * j + (2 * k - 1) : 2 * j;
		}
	}
	return list;
}

測定結果

JavaScript において安全に使用できる整数最大値 Number.MAX_SAFE_INTEGER すなわち $9,007,199,254,740,991$ 以下の最大素数 $9,007,199,254,740,881$ を入力し，素数判定に要した時間を計測した。なお JavaScript エンジンは Chakra エンジンを使用した。TypedArray を使用してるので，Internet Explorer の JavaScript エンジンでは動作しない。

また，今回さまざまなプロセッサ（インテル第 8 ～ 13 世代）で測定できる機会を得たので比較した結果を以下に示す。

結論

結論，インテル第13世代はハンパなく速い・・・いや，これはアルゴリズムの比較だった。

残念ながら，ループ処理の軽減によるアルゴリズム改善は必ずしも高速化に寄与していない。ただし，第11世代では僅かに速くなっている。インテルコアプロセッサも第11世代以降はL2キャッシュ容量を急激に増大させており，その影響かとも思われたが，第13世代の結果を見ると説明がつかない。

表１　プロセッサ仕様の比較

世代	8th	10th	11th	13th
プロセッサ名	i7-8550U	i5-10400	i5-1135G7	i7-13700
コード名	Kaby Lake R	Comet Lake	Tiger Lake	Raptor Lake
コア/スレッド	4/8	6/12	4/8	(P)8/16
ベースクロック	1.8GHz	3.0GHz	2.4GHz	(P)2.1GHz
ブースト時クロック	4.0GHz	4.4GHz	4.2GHz	(P)5.2GHz
L1 命令キャッシュ容量	32kB/core	32kB/core	32kB/core	(P)32KB/core
L1 データキャッシュ容量	32kB/core	32kB/core	48kB/core	(P)48kB/core
L2 キャッシュ容量	256kB/core	256kB/core	1.25MB/core	(P)2MB/core
L3 キャッシュ容量	8MB	12MB	8MB	30MB

※第 13 世代は Performance(P) コアの仕様のみを表示しています。

次回

次回に続きます。次回はついに「２・３・５の倍数を除く」に挑戦します。

• Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(6) - Qiita
• Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(7) - Qiita
• Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(8) - Qiita
• Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(9) 最終回 - Qiita