色々な素数判定とゴールドバッハ予想 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$ で割り切れるかの確認のみなので、難しいことはありません。

偶数判定

for num in range(3, 1000001):
    if num % 2 == 0:
        pass  # 実際の処理
    else:
        continue

といいつつ、実際はループを $2$ 刻みにすれば解決です。

こちらでよい

for num in range(4, 1000001, 2):
    pass  # 実際の処理

素数判定

今回のアルゴリズムで重要なのが、素数の判定です。
任意の自然数 $N$ において、 $2\leqq a<N$ を満足する約数 $a$ が存在するかを確認します。

素数判定

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}$ となるはずです。
よって素数判定のループ数を半分にしても問題ないことが分かります。

効率化された素数判定

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}$ となるはずです。
よって素数判定のループ数は半分どころか平方根をとった数でもよいことが分かります。

さらに効率化された素数判定

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}$ 以下の素数を知っておく必要があります。

もっと効率化された素数判定

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$ として、素数と判定されたものはリストに格納され、次の判定に利用されます。

素数のリスト

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 は用意されているものとします。

ゴールドバッハ予想

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$ 分かかりましたが、問題なく最後まで計算されました。

その他素数判定

エラトステネスの篩

高速な素数判定のアルゴリズムとして知られており、素数をフラグで管理するのが特徴です。
最初にそこそこのメモリを食いますが、サイズ変更がなく効率的です。

エラトステネスの篩

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

numpy のインストールが必要です

numpy を利用すると、さらに高速かつ容易に記述できます。

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 のネストがなくなるだけで、可読性が段違いです。

素数判定の関数

sympy のインストールが必要です

あまりにあっけないですが、sympy には素数を判定するための関数があります。

sympy での素数判定

import sympy

sympy.isprime(99991)

--> True

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

おわりに

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