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

# 0. 前回までのあらすじ
詳しくは下記の記事を参照されたい。
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件 - Qiita](https://qiita.com/TETSURO1999/items/ffc1cae68c16a5c26fde)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(2) - Qiita](https://qiita.com/TETSURO1999/items/28ca00f958628423d64f?1)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(3) - Qiita](https://qiita.com/TETSURO1999/items/38ddcc0f4f3411924fe9)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(4) - Qiita](https://qiita.com/TETSURO1999/items/a8321c80405a70755c41)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(5) - Qiita](https://qiita.com/TETSURO1999/items/4957acd1614aa319f787)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(6) - Qiita](https://qiita.com/TETSURO1999/items/a51747cfac9f89b51a11)
[Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(7) - Qiita](https://qiita.com/TETSURO1999/items/392fd3a8f999a018edfd)

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

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

```javascript:区分篩（ふるい）自然数探索版のコア部分
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[]` のサイズが半分になることから，倍速化が期待できる。

```javascript:区分篩（ふるい）奇数探索版のコア部分
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$ も相互変換可能である。

```math
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$ なので等号が成立することはない）
  ```math
  \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$ とおく。
  ```math
  \begin{aligned}
  (2m + 1) \cdot p &= 2mp + p = 2mp + 2q + 1 \ge 9741 \\
  mp + q &\ge 4870
  \end{aligned}
  ```
  より
  ```math
  m = \left\lceil \frac{4870 - q}{p} \right\rceil
  ```
  となる。なお $\lceil x \rceil $ は $x$ 以上の最小の整数という意味である。※いわゆる天井関数
  ```math
  \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$ となるので
    ```math
    \begin{aligned}
    n_0 &= (2m + 1) \cdot p \\
        &= 9741 \\
        &= 2j_0 + 1 \\
    j_0 &= 4870
    \end{aligned}
    ```
  - $(4870 - q) \bmod p \ne 0$ のとき
    ```math
    \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$ を用いて
```math
\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}$ は
```math
\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$ に依存せず
```math
\Delta j = j_{k+1} - j_k = 2q + 1 = p
```
一定となる。

# 3. 実装コード
実装コードを以下に示す。
```javascript
//------------------------------------------------------------------------------
// 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 エンジンでは動作しない。

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

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
  fontcolor black
  fontsize 16
  fontstyle normal
}
</style>
title
<font color=gray>■</font>第8世代　<font color=blue>■</font>第10世代　<font color=red>■</font>第11世代　<font color=orange>■</font>第13世代
<img https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/ccbdd4ec-25d6-f77d-2802-8f9706e6b602.png>

図１　区分篩（ふるい）奇数探索版の結果
end title
@enduml
```

# 5. 次回
次回に続きます・・・

- [Qiitaで自作の素因数分解プログラムを公開したら添削されて30倍速くなった件(9) 最終回 - Qiita](https://qiita.com/TETSURO1999/items/5d69e183af27487de7bc)