前回までのあらすじ

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

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

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

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

だがしかし，ようやくアイディアが閃いた！というのが今回である。

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

２・３・５の倍数を除く場合

非等間隔な数列

２・３・５の倍数を除いた自然数 $1, 7, 11, 13, 17, 19, 23, 29, \cdots$ を考える。表１に示すように対象の密度は $8/30 = 0.267\cdots$ になる。自然数を全て探索する場合を $1$ とすると，２の倍数を除く場合（奇数サーチ）が $1/2$，２・３の倍数を除く場合が $1/3$ であるから，２・３の倍数を除いた場合に対して $(8/30)/(1/3) = 4/5$ 倍と探索範囲は僅か２割しか減らないのに対し，コーディングは極めて複雑になる。

$\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$	$49$	$50$	$51$	$52$	$53$	$54$	$55$	$56$	$57$	$58$	$59$	$60$

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

$\textsf{表２　２・３・５の倍数を除く自然数}$
$1$	$7$	$11$	$13$	$17$	$19$	$23$	$29$
$31$	$37$	$41$	$43$	$47$	$49$	$53$	$59$
$61$	$67$	$71$	$73$	$77$	$79$	$83$	$89$
$91$	$97$	$101$	$103$	$107$	$109$	$113$	$119$

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

表３　$7^2$ 以上の $7$ の倍数の合成数の位置（1-origin）
合成数	位置	合成数	位置	合成数	位置	合成数	位置
$7 \cdot 7$	$14$	$7 \cdot 37$	$70$	$7 \cdot 67$	$126$	$7 \cdot 97$	$182$
$7 \cdot 11$	$21$	$7 \cdot 41$	$77$	$7 \cdot 71$	$133$	$7 \cdot 101$	$189$
$7 \cdot 13$	$25$	$7 \cdot 43$	$81$	$7 \cdot 73$	$137$	$7 \cdot 103$	$193$
$7 \cdot 17$	$32$	$7 \cdot 47$	$88$	$7 \cdot 77$	$144$	$7 \cdot 107$	$200$
$7 \cdot 19$	$36$	$7 \cdot 49$	$92$	$7 \cdot 79$	$148$	$7 \cdot 109$	$204$
$7 \cdot 23$	$43$	$7 \cdot 53$	$99$	$7 \cdot 83$	$155$	$7 \cdot 113$	$211$
$7 \cdot 29$	$55$	$7 \cdot 59$	$111$	$7 \cdot 89$	$167$	$7 \cdot 119$	$213$
$7 \cdot 31$	$58$	$7 \cdot 61$	$114$	$7 \cdot 91$	$170$	$7 \cdot 121$	$226$

ところが表３を横方向に（すなわち８個周期で）見てみると，すべての行で位置の差分は $56$ であり，等間隔に並んでいることが分かる。

横方向への移動距離

数列 $m(i) = 30j + w(k)$ とおく。
$i = 8j + k$，$j \ge 0$，$0 \le k \le 7$ として，$w(k) = \lbrace 1, 7, 11, 13, 17, 19, 23, 29 \rbrace$ とする。

$m(i)$ の二乗の初期位置 $I_0$ を求める。

\begin{align}
m(i)^2 &= (30j + w(k))^2\\
&= 900j^2 + 60j\cdot w(k) +w(k)^2 \\
&= 30\left(30j^2 + 2j\cdot w(k) + \left\lfloor\frac{w(k)^2}{30}\right\rfloor\right) + (w(k)^2 \bmod 30) \\
&= 30 J_0 + w(K_0) \\
&= m(I_0) \\
J_0 &= 30j^2 + 2j\cdot w(k) + \left\lfloor\frac{w(k)^2}{30}\right\rfloor\\
K_0 &= w^{-1}(w(k)^2 \bmod 30) \\
I_0 &= 8 J_0 + K_0
\end{align}

初期位置から８個目の位置 $I_8$ を求める。

\begin{align}
m(i) \cdot m(i + 8) &= (30j + w(k))(30(j+1) + w(k)) \\
&= 900j(j+1) + 30j\cdot w(k)+30(j+1) \cdot w(k) + w(k)^2 \\
&= 30\left(30j^2 + 30j + 2j\cdot w(k) + w(k) + \left\lfloor\frac{w(k)^2}{30}\right\rfloor\right) \\
&+ (w(k)^2 \bmod 30) \\
&= 30 J_8 + w(K_8) \\
&= m(I_8) \\
J_8 &= 30j^2 + 30j + 2j\cdot w(k) + w(k) + \left\lfloor\frac{w(k)^2}{30}\right\rfloor\\
K_8 &= w^{-1}(w(k)^2 \bmod 30) = K_0\\
I_8 &= 8J_8 + K_8 = 8J_8 + K_0
\end{align}

すなわち $I_8 - I_0 = 8(J_8 - J_0) = 240j + 8w(k) = 8(30j + w(k)) = 8 m(i)$ となる。
調べると非負の整数 $x \ge 0$ に対して $I_{x + 8} - I_x = 8m(i)$ が成立する。

以上より，初期位置から $I_0 \rightarrow I_1 \rightarrow I_2 \rightarrow I_3 \cdots$ と一つずつ移動しながらマーキングするのは移動量の計算が大変であるのに対し，

\begin{alignat}{3}
I_0 &\rightarrow I_8 &\rightarrow I_{16}\cdots \\
I_1 &\rightarrow I_9 &\rightarrow I_{17} \cdots \\
I_2 &\rightarrow I_{10} &\rightarrow I_{18} \cdots \\
 & \cdots\ & \\
I_7 &\rightarrow I_{15} &\rightarrow I_{23} \cdots \\
\end{alignat}

のようにマーキングしていけば，それぞれの移動量はすべて $8m(i)$ となるので極めて簡素に計算できる。

各スタート位置

初期位置から $x$ 個目の位置 $I_x$ を求める。$0 \le x \le 7$ とする。
$k' = k + x$ として，まず $0 \le k' \le 7$ のときを考える。

\begin{align}
m(i) \cdot m(i + x) &= (30j + w(k))(30j + w(k')) \\
&= 900j^2 + 30j \cdot \lbrace w(k) + w(k') \rbrace +  w(k) \cdot w(k') \\
&= 30\left( 30j^2 + \cdot \lbrace w(k) + w(k') \rbrace + \left\lfloor\frac{w(k) \cdot w(k')}{30}\right\rfloor \right)  \\
&+ (w(k) \cdot w(k') \bmod 30) \\
&= 30 J_x + w(K_x) \\
&= m(I_x) \\
J_x &= 30j^2 + j\cdot (w(k) + w(k')) + \left\lfloor\frac{w(k) \cdot w(k')}{30}\right\rfloor\\
K_x &= w^{-1}(w(k) \cdot w(k') \bmod 30)\\
I_x &= 8J_x + K_x
\end{align}

次に $k'' = k + x - 8$として，$0 \le k'' \le 7$ のときを考える。

\begin{align}
m(i) \cdot m(i + x) &= (30j + w(k))\lbrace 30(j + 1) + w(k'')\rbrace \\
&= 900j\cdot(j + 1) + 30j \cdot w(k'') + 30(j + 1) \cdot w(k) + w(k) \cdot w(k'') \\
&= 30\left( 30j^2 + j \cdot \lbrace w(k) + w(k'') + 30\rbrace + w(k) + \left\lfloor\frac{w(k) \cdot w(k'')}{30}\right\rfloor \right)  \\
&+ (w(k) \cdot w(k'') \bmod 30) \\
&= 30 J_x + w(K_x) \\
&= m(I_x) \\
J_x &= 30j^2 + j\cdot (w(k) + w(k'') + 30) + w(k) + \left\lfloor\frac{w(k) \cdot w(k'')}{30}\right\rfloor\\
K_x &= w^{-1}(w(k) \cdot w(k'') \bmod 30)\\
I_x &= 8J_x + K_x
\end{align}

ここで $w(k) = \lbrace 1,7,11,13,17,19,23,29,31,37,41,43,47,49,53,59 \rbrace$ と拡張すれば，両者を統合して

\begin{align}
J_x &= 30j^2 + j\cdot (w(k) + w(k + x)) + \left\lfloor\frac{w(k) \cdot w(k + x)}{30}\right\rfloor\\
K_x &= w^{-1}(w(k) \cdot w(k + x) \bmod 30)\\
I_x &= 8J_x + K_x
\end{align}

とおける。

下記の要素は事前に計算しておき，テーブル化しておく。

$w(k) + w(k + x)$
$\lfloor w(k) \cdot w(k + x) / 30 \rfloor$
$w^{-1}(w(k) \cdot w(k+x) \bmod 30)$

$0 \le k \le 7$ かつ $0 \le x \le 7$ なので $8 \times 8$ の二次元テーブルとなる。

表４ $w(k) + w(k+x)$
		$x$
		$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$
$k$	$0$	$2$	$8$	$12$	$14$	$18$	$20$	$24$	$30$

	$1$	$14$	$18$	$20$	$24$	$26$	$30$	$36$	$38$

	$2$	$22$	$24$	$28$	$30$	$34$	$40$	$42$	$48$

	$3$	$26$	$30$	$32$	$36$	$42$	$44$	$50$	$54$

	$4$	$34$	$36$	$40$	$46$	$48$	$54$	$58$	$60$

	$5$	$38$	$42$	$48$	$50$	$56$	$60$	$62$	$66$

	$6$	$46$	$52$	$54$	$60$	$64$	$66$	$70$	$72$

	$7$	$58$	$60$	$66$	$70$	$72$	$76$	$78$	$82$

表５　$\lfloor w(k)\cdot w(k+x)/30 \rfloor$
		$x$
		$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$
$k$	$0$	$0$	$0$	$0$	$0$	$0$	$0$	$0$	$0$

	$1$	$1$	$2$	$3$	$3$	$4$	$5$	$6$	$0$

	$2$	$4$	$4$	$6$	$6$	$8$	$10$	$0$	$2$

	$3$	$5$	$7$	$8$	$9$	$12$	$0$	$3$	$4$

	$4$	$9$	$10$	$13$	$16$	$0$	$3$	$6$	$7$

	$5$	$12$	$14$	$18$	$0$	$4$	$6$	$8$	$10$

	$6$	$17$	$22$	$0$	$5$	$8$	$9$	$13$	$14$

	$7$	$28$	$0$	$6$	$10$	$12$	$16$	$18$	$22$

表６　$w^{-1}(w(k) \cdot w(k+x) \bmod 30)$
		$x$
		$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$
$k$	$0$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$

	$1$	$5$	$4$	$0$	$7$	$3$	$2$	$6$	$1$

	$2$	$0$	$6$	$1$	$7$	$3$	$5$	$2$	$4$

	$3$	$5$	$2$	$1$	$7$	$4$	$3$	$0$	$6$

	$4$	$5$	$6$	$0$	$3$	$4$	$7$	$1$	$2$

	$5$	$0$	$4$	$2$	$5$	$3$	$7$	$1$	$6$

	$6$	$5$	$1$	$6$	$2$	$3$	$7$	$0$	$4$

	$7$	$0$	$7$	$6$	$5$	$4$	$3$	$2$	$1$

実装コード

JavaScript による実装コードを以下に示す。
配列 a[k][x] は $w(k) + w(k+x)$，
配列 b[k][x] は $8 \lfloor w(k) \cdot w(k+x)/30 \rfloor + w^{-1}(w(k) \cdot w(k+x) \bmod 30)$ に相当する。

エラトステネスの篩（ふるい）２・３・５の倍数を除く場合

function create_list(n) {
	var	n1 = Math.floor(Math.sqrt(n));
	var	rem = [0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 4, 4,
			   4, 4, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8];
	var	table = new Uint8Array(8 * Math.floor(n1 / 30) + rem[n1 % 30]);
	var	size = table.length;
	var	w = [1, 7, 11, 13, 17, 19, 23, 29];
	var	a = [[  2,  8, 12, 14, 18, 20, 24, 30],[ 14, 18, 20, 24, 26, 30, 36, 38],
			 [ 22, 24, 28, 30, 34, 40, 42, 48],[ 26, 30, 32, 36, 42, 44, 50, 54],
			 [ 34, 36, 40, 46, 48, 54, 58, 60],[ 38, 42, 48, 50, 56, 60, 62, 66],
			 [ 46, 52, 54, 60, 64, 66, 70, 72],[ 58, 60, 66, 70, 72, 76, 78, 82]];
	var b = [[  0,  1,  2,  3,  4,  5,  6,  7],[ 13, 20, 24, 31, 35, 42, 54, 57],
			 [ 32, 38, 49, 55, 67, 85, 90,108],[ 45, 58, 65, 79,100,107,128,142],
			 [ 77, 86,104,131,140,167,185,194],[ 96,116,146,157,187,207,217,238],
			 [141,177,190,226,251,263,288,300],[224,239,286,317,332,363,378,409]];
	var	list = [];
	list.push(2);
	list.push(3);
	list.push(5);
	for(var i = 1; i < size; i++) {
		if(table[i] != 0) continue;
		var	j = i >>> 3;
		var	k = i & 7;
		var	m = 30 * j + w[k];
		list.push(m);
		for(var i0 = 8 * (30 * j * j + j * a[k][0]) + b[k][0],
				i1 = 8 * (30 * j * j + j * a[k][1]) + b[k][1],
				i2 = 8 * (30 * j * j + j * a[k][2]) + b[k][2],
				i3 = 8 * (30 * j * j + j * a[k][3]) + b[k][3],
				i4 = 8 * (30 * j * j + j * a[k][4]) + b[k][4],
				i5 = 8 * (30 * j * j + j * a[k][5]) + b[k][5],
				i6 = 8 * (30 * j * j + j * a[k][6]) + b[k][6],
				i7 = 8 * (30 * j * j + j * a[k][7]) + b[k][7];
			i7 < size;
			i0 += 8 * m, i1 += 8 * m, i2 += 8 * m, i3 += 8 * m,
			i4 += 8 * m, i5 += 8 * m, i6 += 8 * m, i7 += 8 * m)
		{
			table[i0] = 1; table[i1] = 1; table[i2] = 1; table[i3] = 1;
			table[i4] = 1; table[i5] = 1; table[i6] = 1; table[i7] = 1;
		}
		if(i0 >= size) continue; table[i0] = 1;
		if(i1 >= size) continue; table[i1] = 1;
		if(i2 >= size) continue; table[i2] = 1;
		if(i3 >= size) continue; table[i3] = 1;
		if(i4 >= size) continue; table[i4] = 1;
		if(i5 >= size) continue; table[i5] = 1;
		if(i6 >= size) continue; table[i6] = 1;
	}
	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 世代）で測定できる機会を得たので比較した結果を以下に示す。今回の結果は，エラトステネスの篩（ふるい）２・３・５の倍数除去にあたる。２・３の倍数除去に対して使用メモリ容量を２割低減できたにも関わらず，第 11 世代を除いて高速化できていないどころか，むしろ遅くなってしまっている。

もちろんライバルの試し割り法にも追いつけていない。

結論

インテル第 13 世代 Raptor Lake はめっちゃ速い！
試し割り法はプロセッサの世代進化に応じて順当に速くなっている。
エラトステネスの篩（ふるい）はアルゴリズムの改善（使用メモリ容量の低減）ともプロセッサの世代進化とも綺麗な相関関係は見られない。

参考

参考までに以下のコードを以下に示す。

エラトステネスの篩（ふるい）
- 自然数探索
- 奇数探索
- ２・３の倍数除去
試し割り法
- ２・３・５の倍数除去（Wheel Factorization）手動ループアンローリング版

エラトステネスの篩（ふるい）の自然数探索版から奇数版，２・３の倍数を除く場合，そして今回の２・３・５の倍数を除く場合と進化するにつれて急激にプログラムが複雑になっていることが分かる。

エラトステネスの篩（ふるい）自然数探索の場合

function create_list(n) {
	var	n1 = Math.floor(Math.sqrt(n));
	var	table = new Uint8Array(1 + n1);
	var	size = table.length;
	var	list = [];
	for(var i = 2; i < size; i++) {
		if(table[i] != 0) continue;
		list.push(i);
		for(var j = i * i; j < size; j += i)
			table[j] = 1;
	}
	return list;
}

エラトステネスの篩（ふるい）奇数探索の場合

function create_list(n) {
	var	n1 = Math.floor(Math.sqrt(n));
	var	n2 = (n1 & 1 == 0) ? n1 - 1 : n1;
	var	table = new Uint8Array(1 + (n2 - 1) / 2);
	var	size = table.length;
	var	list = [];
	list.push(2);
	for(var i = 1; i < size; i++) {
		if(table[i] != 0) continue;
		var	m = 2 * i + 1;
		list.push(m);
		for(var j = 2 * i * (i + 1); j < size; j += m)
			table[j] = 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;
}

エラトステネスの篩（ふるい）共通素因数分解

function prime_factorize(n, list) {
	var	a = [];
	for(var i = 0, p = n; i < list.length && (q = list[i]) * q <= p; i++) {
		while(p % q == 0) {
			p /= q; a.push(q);
		}
	}
	if(p > 1) a.push(p);
	return a;
}

試し割り法（２・３・５の倍数を除く）Wheel Factorization による素因数分解

function prime_factorize(p) {
	var	a = [];
	var	r;
	while((r = p / 2) == Math.round(r)) { p = r; a.push(2); }
	while((r = p / 3) == Math.round(r)) { p = r; a.push(3); }
	while((r = p / 5) == Math.round(r)) { p = r; a.push(5); }
	for(var q = 7; q * q <= p;) {
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 4;
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 2;
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 4;
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 2;
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 4;
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 6;
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 2;
		while((r = p / q) == Math.round(r)) { p = r; a.push(q); } q += 6;
	}
	if(p > 1) a.push(p);
	return a;
}

次回

次回に続きます。まだだ・・・まだ終わらんよ！