LoginSignup
0
0
@TETSURO1999

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

Last updated at Posted at 2023-11-03

0. 前回までのあらすじ

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

1. 区分篩（ふるい）のさらなる高速化

前回の記事より，区分篩（ふるい）自然探索版のコア部分を抽出すると以下のようになる。素数リスト list1[] および 素数判定テーブル table[] は L1 データキャッシュ内に十分収まるはずである。

区分篩（ふるい）自然数探索版のコア部分
for(var i = 0; i < 1201; i++) {
	var	p = list1[i];
	for(var j = start[i]; j < 9740; j += p)
		table[j] = 1;
	start[i] = j - 9740;
}

一方，奇数探索版は以下のようになる。素数リスト list1[] は $2$ が抜けるのでサイズが一つ減る。また素数判定テーブル table[] のサイズが半分になることから，倍速化が期待できる。

区分篩（ふるい）奇数探索版のコア部分
for(var i = 0; i < 1200; i++) {
	var	p = list1[i];
	for(var j = start[i]; j < 4870; j += p)
		table[j] = 1;
	start[i] = j - 4870;
}

2. 奇数探索版の設計

2.1 素数リストの定義

$9739$ 以下の奇数の素数は $1200$ 個ある。素数 $p$ の素数判定テーブルおける位置を $q$ とする。

$i$ $0$ $1$ $2$ $3$ $4$ $5$ $6$ $7$ $\cdots$ $1199$
$q$ $1$ $2$ $3$ $5$ $6$ $8$ $9$ $11$ $\cdots$ $4869$
$p$ $3$ $5$ $7$ $11$ $13$ $17$ $19$ $23$ $\cdots$ $9739$

素数の値 $p$ と位置 $q$ は相互変換可能である。

$$
q = \frac{p - 1}{2} \iff p = 2 q + 1
$$

2.2 素数判定テーブル（奇数のみ）の定義

素数判定テーブル（奇数のみ）を定義する。位置 $j$ における値 $n$ とおく。

$j$ $0$ $1$ $2$ $3$ $4$ $5$ $6$ $7$ $\cdots$ $3246$ $4869$
$n$ $1$ $3$ $5$ $7$ $9$ $11$ $13$ $15$ $\cdots$ $9737$ $9739$

素数判定テーブルの値 $n$ と位置 $j$ も相互変換可能である。

j = \frac{n - 1}{2} \iff n = 2 j + 1

2.3 初回位置

通常のエラトステネスの篩（ふるい）では $p^2$ 以上の $p$ の倍数を合成数として消していく。すなわち初回位置は $p^2$ である。一方，今回の区分篩（ふるい）の場合，$9739$ 以下の素数リストは既に得られており，初回位置は $9741$ 以上の $p$ の倍数となる。すなわち $p^2$ と $9741$ の大小関係によって場合分けが必要である。

区分篩（ふるい）を行う最初の値 $n_0$，素数判定テーブルにおける位置 $j_0$ とおく。

  • $p^2 \ge 9741$ のとき（$\sqrt{9741} = 98.7 \cdots$ なので等号が成立することはない）

    \begin{aligned}
n_0 &= p^2 = (2q + 1)^2 = 4q^2 + 4q + 1 = 2 \cdot \left\lbrace 2q \cdot (q + 1) \right\rbrace + 1 \\
    &= 2j_0 + 1 \\
j_0 &= 2q \cdot (q + 1)
\end{aligned}

  • $p^2 < 9741$ のとき，$(2m + 1) \cdot p \ge 9741$ となる最小の自然数 $m$ とおく。

    \begin{aligned}
(2m + 1) \cdot p &= 2mp + p = 2mp + 2q + 1 \ge 9741 \\
mp + q &\ge 4870
\end{aligned}

    より

    m = \left\lceil \frac{4870 - q}{p} \right\rceil

    となる。なお $\lceil x \rceil $ は $x$ 以上の最小の整数という意味である。※いわゆる天井関数

    \begin{aligned}
n_0 &= (2m + 1) \cdot p \\
    &= (2m + 1) \cdot (2q + 1) \\
    &= 2\cdot(2mq + m + q) + 1 \\
    &= 2j_0 + 1 \\
j_0 &= 2mq + m + q \\
\end{aligned}
    • $(4870 - q) \bmod p = 0$ のとき，$(2m + 1) \cdot p = 9741$ となるので 
      \begin{aligned}
n_0 &= (2m + 1) \cdot p \\
    &= 9741 \\
    &= 2j_0 + 1 \\
j_0 &= 4870
\end{aligned}
    • $(4870 - q) \bmod p \ne 0$ のとき 
      \begin{aligned}
n_0 &= (2m + 1 ) \cdot p \\
    &= 2mp + 2q + 1 \\
    &= 2 \left\lbrace \left\lceil \frac{4870 - q}{p} \right\rceil \cdot p + q \right\rbrace + 1 \\
    &= 2j_0 + 1 \\
j_0 &= \left\lceil \frac{4870 - q}{p} \right\rceil \cdot p + q \\
    &= \left\lbrace \left\lfloor \frac{4870 - q}{p} \right\rfloor + 1 \right\rbrace \cdot p + q \\
    &= (4870 - q) - \lbrace (4870 - q) \bmod p \rbrace + p + q \\
    &= p + 4870 + \lbrace (4870 - q) \bmod p \rbrace \\
\end{aligned}
      なお，$\lfloor x \rfloor$ は $x$ 以下の最大の整数という意味である。※いわゆる床関数

2.4 移動幅

区分篩（ふるい）を行う $k$ 番目の値 $n_k$ および位置 $j_k$ とする。$n_k$ は素数 $p$ の奇数倍であるから，適当な自然数 $x$ を用いて

\begin{aligned}
n_k &= (2x + 1) \cdot p \\
    &= (2x + 1) \cdot (2q + 1) \\
    &= 2\cdot(2xq + x + q) + 1 \\
    &= 2j_k + 1 \\
j_k &= 2xq + x + q
\end{aligned}

と表せる。一方，$k + 1$ 番目の値 $n_{k+1}$ および位置 $j_{k+1}$ は

\begin{aligned}
n_{k+1} &= \left\lbrace 2(x + 1) + 1\right\rbrace \cdot p \\
        &= (2x + 3) \cdot (2q + 1) \\
        &= 2\cdot(2xq + x + 3q + 1) + 1 \\
        &= 2j_{k+1} + 1 \\
j_{k+1} &= 2xq + x + 3q + 1 \\
\end{aligned}

となる。これより素数判定テーブルにおける位置の移動幅 $\Delta j$ は $k$ に依存せず

\Delta j = j_{k+1} - j_k = 2q + 1 = p

一定となる。

3. 実装コード

実装コードを以下に示す。

//------------------------------------------------------------------------------
// 3 以上 9,740 以下の素数リスト ※要素数 1,200個
//------------------------------------------------------------------------------
var	list1 = new Uint16Array([
	3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,

	/* 中略 */
 
	9697,9719,9721,9733,9739]);
//------------------------------------------------------------------------------
// 与えられた自然数 n に対して √n 以下の素数リストを作成する
//------------------------------------------------------------------------------
function create_list(n) {
	var	nsqrt = Math.floor(Math.sqrt(n));
	var	list2 = [];							// 得られた素数を格納するリスト
	var	start = new Uint32Array(1200);		// 初期位置
	var	table = new Uint8Array(4870);		// 区分篩（ふるい）を行う素数判定テーブル
 	//--------------------------------------------------------------------------
	// 初期位置の計算
  	//--------------------------------------------------------------------------
	list2.push(2);
	for(var i = 0; i < 1200; i++) {
		var	p = list1[i];
		list2.push(p);
		var	q = (p - 1) >>> 1;
		var	r, s;
		if((r = 2 * q * (q + 1) - 4870) >= 0)
			start[i] = r;
		else if((s = (4870 - q) % p) != 0)
			start[i] = p - s;
	}
  	//--------------------------------------------------------------------------
	// 区分篩（ふるい）の実行
  	//--------------------------------------------------------------------------
	for(var offset = 9740; offset < nsqrt; offset += 9740) {
		for(var i = 0; i < 1200; i++) {
			var	p = list1[i];
			for(var j = start[i]; j < 4870; j += p)
				table[j] = 1;
			start[i] = j - 4870;
		}
		for(var j = 0; j < 4870; j++) {
			if(table[j] == 0) list2.push(offset + 2 * j + 1);
			table[j] = 0;
		}
	}
	return list2;
}

4. 結果

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

残念ながら自然数探索版の２倍速とはならなかったが，どの世代のプロセッサでも著しい改善効果が見られた。試し割り法に迫る勢いである。あとちょいでアレである。

5. 次回

次回に続きます・・・

0
0
0

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
0
0