色々な素数判定とゴールドバッハ予想 with Python

# はじめに

小学生でも理解できる未解決問題として知られている**ゴールドバッハ予想**を、Python で確認していきます。

# ゴールドバッハ予想とは

$2$ より大きい偶数は、$2$ つの素数の和で表されるという予想です。
$4$ は $4=2+2$ の組み合わせなので、$6$ 以上の偶数は奇素数の和で表されるとも言い換えられます。
$$\begin{align*}
  4=2+2\\\\
  6=3+3\\\\
  8=3+5\\\\
  10=3+7\\\\
  10=5+5\\\\
  12=5+7\\\\
\end{align*}$$

コンピュータによって $4\times 10^{18}$ までの偶数で成立するそうで、統計学的な見地からおそらくこの予想は正しいとされていますが、ポリア予想のようなケースもあるので実際どんな数にも当てはまるかは分かっていません。

# ゴールドバッハ予想が成立するか確認してみよう

$4\times 10^{18}$ までみてられないので、$100000$ ぐらいまでを確認します。
素数の組み合わせが複数存在する場合もありますが、1 つでも発見されれば終了とします。

## 偶数判定

$100000$ までの $2$ を除く偶数を判定しますが、$2$ で割り切れるかの確認のみなので、難しいことはありません。
```python:偶数判定
for num in range(3, 1000001):
    if num % 2 == 0:
        pass  # 実際の処理
    else:
        continue
```
といいつつ、実際はループを $2$ 刻みにすれば解決です。
```python:こちらでよい
for num in range(4, 1000001, 2):
    pass  # 実際の処理
```

## 素数判定

今回のアルゴリズムで重要なのが、素数の判定です。
任意の自然数 $N$ において、 $2\leqq a<N$ を満足する約数 $a$ が存在するかを確認します。
```python:素数判定
N = 100  # 任意の自然数
for i in range(2, N):
    if N % i == 0:
        break  # 素数ではないのでループを抜ける
    else:
        continue
else:
    print("prime number")
```
`N` を素数に変更すると、`prime number` と表示されます。
ただしこの方法は数が大きくなるにつれて、効率が悪くなります。

### 素数判定を効率化

例えば $60$ の約数を記述すると、$1,2,3,4,5,6,10,12,15,20,30,60$ となります。
$$\begin{align*}
  1\times60\\\\
  2\times30\\\\
  3\times20\\\\
  4\times15\\\\
  5\times12\\\\
  6\times10\\\\
  10\times6\\\\
  12\times5\\\\
  15\times4\\\\
  20\times3\\\\
  30\times2\\\\
  60\times1\\\\
\end{align*}$$

$N$ の約数を小さい順から並べたものと、その逆順のものをそれぞれ乗算すると、$N$ になります。
この法則から $N$ の $2$ 番目に大きい約数は最大でも $\frac{N}{2}$ となるはずです。
よって素数判定のループ数を半分にしても問題ないことが分かります。

```python:効率化された素数判定
N = 100  # 任意の自然数
for i in range(2, int(N / 2) + 1):
    if N % i == 0:
        break  # 素数ではないのでループを抜ける
    else:
        continue
else:
    print("prime number")
```

### さらに素数判定を効率化

例えば $36$ の約数を記述すると、$1,2,3,4,6,9,12,18,36$ となります。
$$\begin{align*}
  1\times36\\\\
  2\times18\\\\
  3\times12\\\\
  4\times9\\\\
  6\times6\\\\
\end{align*}$$

乗算は順序を考慮しなくてよいので、同じ約数の組み合わせは省略しました。
この法則から考慮すべき $N$ の約数の最大値は $\sqrt{N}$ となるはずです。
よって素数判定のループ数は半分どころか平方根をとった数でもよいことが分かります。
```python:さらに効率化された素数判定
N = 100  # 任意の自然数
for i in range(2, int(N**0.5) + 1):
    if N % i == 0:
        break  # 素数ではないのでループを抜ける
    else:
        continue
else:
    print("prime number")
```

### もっと素数判定を効率化

例えば $37$ が素数か判定する際に、$2,3,4,5,6$ で除せるかを確認しますが、$6$ で除せるということは $2$ や $3$ でも除せるはずです。
つまり合成数である約数を判定に利用せず、既知の素数のみを判定に利用すればよいことになります。
ただしこの場合、事前に $\sqrt{N}$ 以下の素数を知っておく必要があります。
```python:もっと効率化された素数判定
primes = [2, 3, 5, 7]  # 素数
N = 100  # 任意の自然数
for prime in primes:
    if N % prime == 0:
        break  # 素数ではないのでループを抜ける
    else:
        continue
else:
    print("prime number")
```
$2$ 桁の数であれば `[2, 3, 5, 7]`、$3$ 桁であっても `[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]` で事足ります。

### 素数を格納するリストを生成

前述のロジックを利用して、素数のみを格納したリストを生成します。
初期値を $2,3$ として、素数と判定されたものはリストに格納され、次の判定に利用されます。
```python:素数のリスト
primes = [2, 3]  # 素数
M = 100000  # 上限
for i in range(4, M + 1):
    for prime in primes:
        if i % prime == 0:
            break  # 素数ではないのでループを抜ける
        else:
            if prime > i**0.5:
                primes.append(i)
                break
            continue
```
ちなみに $100000$ 以下で最大の素数は $99991$ です。

## ゴールドバッハ予想の実装

`primes` は用意されているものとします。
```python:ゴールドバッハ予想
M = 100000  # 上限
for num in range(4, M + 1, 2):
    for prime in primes:
        if num - prime in primes:
            print(f"{num}={prime}+{num - prime}")
            break
    else:
        print("ゴールドバッハ予想は否定された!")
        break
```
古い PC なので約 $1.5$ 分かかりましたが、問題なく最後まで計算されました。

# その他素数判定

## エラトステネスの篩

高速な素数判定のアルゴリズムとして知られており、素数をフラグで管理するのが特徴です。
最初にそこそこのメモリを食いますが、サイズ変更がなく効率的です。
```python:エラトステネスの篩
M = 100000  # 上限
primes = [True for _ in range(M + 1)]
primes[0:2] = [False, False]  # 0 と 1 は素数ではない

for num in range(2, M + 1):
    if not primes[num]:  # 合成数であればスキップ
        continue
    for comp in range(num**2, M + 1, num):  # 素数の倍数は合成数
        primes[comp] = False
primes
```
このアルゴリズムの手順は以下の通りです。
1. 上限 $+1$ のサイズで `True` のリストを定義 ($0$ を含むため $+1$)
2. $0,1$ は素数ではないので、フラグを `False` に変更
3. フラグが `False` であれば、素数ではないのでスキップ
4. 素数を発見したら、その倍数は素数ではないのでフラグを `False` に変更

4 の `for comp in range(num**2, M + 1, num)` が分かりにくいので、細かく説明します。
例えば $7$ を例にすると、これは素数でその倍数は合成数となります。
ただし既に $7$ よりも小さい素数の倍数に含まれる $2\times 7,\quad3\times 7,\quad5\times 7$ は、既にふるい落とされていることになります。
よってまだ判定がなされていない最小の合成数は自乗である $7\times 7$ であり、それ以降の倍数をふるい落とすために `range` のステップ数を `num` としています。

## エラトステネスの篩 with numpy

:::note warn
numpy のインストールが必要です
:::
numpy を利用すると、さらに高速かつ容易に記述できます。
```python:numpy でエラトステネスの篩
import numpy as np

M = 100000  # 上限
primes_array = np.ones(M + 1, dtype=bool)
primes_array[0:2] = False  # 0, 1 は素数ではない

for num in range(2, M + 1):
    if not primes_array[num]:  # 合成数であればスキップ
        continue
    primes_array[num**2 :: num] = False  # 素数の倍数は合成数
primes_array
```
`for` のネストがなくなるだけで、可読性が段違いです。

## 素数判定の関数

:::note warn
sympy のインストールが必要です
:::
あまりにあっけないですが、sympy には素数を判定するための関数があります。
```python:sympy での素数判定
import sympy

sympy.isprime(99991)

--> True
```
`isprime` 関数は引数が素数であれば `True`、合成数であれば `False` を返却します。
内部的には**ミラーラビン素数判定法**が利用されているので、かなり高速です。
他にも sympy には素数に関する関数があるので、興味あれば色々調べてみてください。

# おわりに

ゴールドバッハ予想の解決には $\yen120000000$ の懸賞金がかかっているそうです。
金欠の方は挑戦してみてはいかがでしょうか ??