Pythonでゴールドバッハ予想｜Pythonで数学を学ぼう！ 第10回

はじめに

1742年の手紙

1742年にプロイセンの数学者クリスティアン・ゴールドバッハ（1690-1764）がレオンハルト・オイラーへあてた手紙の中に

5より大きな任意の自然数は、三つの素数の和で表せる。…（☆）

と書きました。

6=2+2+2　偶数＝偶数の素数＋偶数の素数＋偶数の素数
7=2+2+3　奇数＝偶数の素数＋偶数の素数＋奇数の素数
8=2+3+3　偶数＝偶数の素数＋奇数の素数＋奇数の素数
9=3+3+3　奇数＝奇数の素数＋奇数の素数＋奇数の素数
10=2+3+5　偶数＝偶数の素数＋奇数の素数＋奇数の素数

このように、偶数を3つの素数で表す場合にはその1つが2であることがわかります。

したがって、

6=2+2+2　→　4=2+2
8=2+3+3　→　6=3+3
10=2+3+5　→　8=3+5

（☆）は

4以上の任意の偶数は、2つの素数の和で表せる。

または

6以上の任意の偶数は、2つの奇素数の和で表せる。…（★）

と言い換えることができます。ただし、同じ素数を2度使ってもよいとします。

前回のコラッツ予想同様にゴールドバッハ予想も正しいと考えられていますが、未解決の問題です。

Pythonで検証プログラムを作り挑戦してみましょう。

ゴールドバッハ予想を検証するプログラム　フロチャート

そこで（★）を検証するプログラムをPythonでコーディングしてみます。

以下のようなブロックになります。

[1] 検証する偶数の範囲を決める m以上n以下
[2] 3以上max以下の奇素数リストを生成
[3] 偶数even = 奇素数odd_prime + 奇素数even - odd_prime を調べる

[2]において、奇素数リストodd_prime=[3,5,7,11,…]を生成します。

最も単純なアルゴリズムによる「goldbach1.py」「goldbach2.py」と、連載（素数生成プログラム（エラトステネスの篩）を実装｜Pythonで数学を学ぼう！第7回）で紹介した、より高速なエラトステネスの篩による「goldbach3.py」の2通りとして速度（計算時間）を比較できるようにしました。

[3]において、偶数＝奇素数＋奇素数となる2つの奇素数を見つけだします。

コーディングのポイントは以下です。

>>>　if even - odd_prime[i] in odd_prime:

最初に選んだ偶数evenと次に奇素数リストodd_prime[]から選んだ奇素数odd_prime[i]に対して、差even - odd_prime[i] が奇素数リストodd_prime[]に含まれているか否かを判定します。

含まれていれば（★）を満たす2つの奇素数が見つかったことになります。含まれていなければ見つかるまで探します。

ゴールドバッハ予想検証プログラム「goldbach1.py」

[3]において（★）を満たす2つの奇素数の組をすべて見つけだし表示します。

>>>　#[1]検証する偶数の範囲を決める m以上n以下
>>>　min = int(input('ゴールドバッハ予想（6以上の偶数=2つの奇素数の和）をmin(≧6)からmaxの範囲の偶数で調べます。\nmin >> '))
>>>　if min % 2 == 0:
>>>　　　min = min
>>>　else:
>>>　　　min = min + 1
>>>
>>>　max = int(input('max >> '))
>>>

>>>　#[2] 3以上max以下の奇素数リストを生成
>>>　odd_prime = []
>>>　for i in range (3,max+1,2):
>>>　　　for j in range(3, i):
>>>　　　　　if i % j ==0:
>>>　　　　　　　break
>>>　　　else:
>>>　　　　　odd_prime.append(i)
>>>　print('奇素数リスト\n', odd_prime)
>>>
>>>
>>>　#[3]偶数even = 奇素数odd_prime + 奇素数even - odd_prime を調べる
>>>　# すべての組合せを表示
>>>　for even in range(min, max+1 , 2):
>>>　print(even, end=' ')
>>>　i = 0
>>>　while i <= len(odd_prime) - 1:
>>>　if even - odd_prime[i] in odd_prime:
>>>　print('=', odd_prime[i], '+', even - odd_prime[i], end=' ')
i += 1
>>>　print()

コードファイルのダウンロードはコチラ（https://drive.google.com/file/d/1I8bxuV6MfHSBzoJNRW2c1gTX1RUn_jh3/view?usp=sharing）から

次は6以上30以下の偶数についての結果です。

ゴールドバッハ予想検証プログラム「goldbach2.py」

素数生成アルゴリズムはgoldbach1.pyと同じものを使用します。

goldbach1.pyでは大きい偶数evenに対しては組合せも多く見つかるので、大きな偶数に対しては時間がかかりすぎます。

そこで、最初の1組の奇素数が見つかればOKとして、大きい偶数に対しても確認がスムーズにできるようにします。ついでに計算時間を計測してgoldbach3.pyと比較してみます。

>>>　import datetime as dt
>>>
>>>　#[1]検証する偶数の範囲を決める m以上n以下
min = int(input('ゴールドバッハ予想（6以上の偶数=2つの奇素数の和）をmin(≧6)からmaxの範囲の偶数で調べます。\nmin >> '))
>>>　if min % 2 == 0:
>>>　　　min = min
>>>　else:
>>>　　　min = min + 1
>>>
>>>　max = int(input('max >> '))

>>>　time = dt.datetime
>>>　start = time.now() 　　　　# 計算開始時刻マーク
>>>
>>>　#[2] 3以上max以下の奇素数リストを生成
>>>　odd_prime = []
>>>　for i in range (3, max+1, 2):
>>>　　　for j in range(3, i):
>>>　　　　　if i % j ==0:
>>>　　　　　　　break
>>>　　　else:
>>>　　　　　odd_prime.append(i)
>>>　print('奇素数リスト\n', odd_prime)
>>>
>>>
>>>　#[3]偶数even = 奇素数odd_prime + 奇素数even - odd_prime を調べる
>>>　# 1組だけを表示
>>>　for even in range(min, max + 1, 2):
>>>　　　i = 0
>>>　　　while i <= len(odd_prime) - 1:
>>>　　　　　if even - odd_prime[i] in odd_prime:
>>>　　　　　　　print(even, end=' ')
>>>　　　　　　　print('=', odd_prime[i], '+', even - odd_prime[i], end=' ')
>>>　　　　　　　i += 1
>>>　　　　　　　break
>>>　　　　　else:
>>>　　　　　　　i += 1
>>>　　　　　　　continue
>>>　　　　　else:
>>>　　　　　　　print(even, ' = NG') 　　　　# 予想の反例
>>>　　　　　print()
>>>
>>>
>>>　stop = time.now() 　　　　# 計算終了時刻マーク
>>>　t = stop - start
>>>　print('計算時間',t)

コードファイルのダウンロードはコチラ（https://drive.google.com/file/d/1sBVb-9EF2zLF91PIHepbo_f5CSoEpNlA/view?usp=sharing）から

以下は1000以上1030以下の偶数についての結果です。

ゴールドバッハ予想検証プログラム「goldbach3.py」

goldbach2.pyの[2]をエラトステネスのふるいに変えたプログラムです。

>>>　import math
>>>　import datetime as dt
>>>
>>>
>>>　#[1]検証する偶数の範囲を決める m以上n以下
>>>　min = int(input('ゴールドバッハ予想（6以上の偶数=2つの奇素数の和）をmin(≧6)からmaxの範囲の偶数で調べます。\nmin >> '))
>>>　if min % 2 == 0:
>>>　　　min = min
>>>　else:
>>>　　　min = min + 1
>>>
>>>　max = int(input('max >> '))
>>>
>>>　time = dt.datetime
>>>　start = time.now() 　# 計算開始時刻マーク
>>>
>>>　#[2] 3以上max以下の奇素数リストを生成 エラトステネスのふるい高速
>>>　p = [i for i in range(max + 1)]
>>>
>>>　for i in p[3:]:
>>>　　　if p[i] % 2 == 0:
>>>　　　　　p[i] = 0
>>>
>>>　sqrt_max = math.sqrt(max) 　　　　# √max の計算
>>>　for i in range(3, max):
>>>　　　if i > sqrt_max:
>>>　　　　　break
>>>　　　if p[i] != 0:
>>>　　　　　for j in range(i, max + 1, 2):
>>>　　　　　　　if i * j >= max + 1 :
>>>　　　　　　　　　break
>>>　　　　　　　p[i * j] = 0
>>>　odd_prime = sorted(list(set(p)))[2:]
>>>　print('奇素数リスト\n', odd_prime)
>>>
>>>
>>>　#[3]偶数even = 奇素数odd_prime + 奇素数even - odd_prime を調べる
>>>　# 1組だけを表示
>>>　for even in range(min, max + 1 , 2):
>>>　　　i = 0
>>>　　　while i <= len(odd_prime) - 1:
>>>　　　　　if even - odd_prime[i] in odd_prime:
>>>　　　　　　　print(even, end=' ')
>>>　　　　　　　print('=', odd_prime[i], '+', even - odd_prime[i], end=' ')
>>>　　　　　　　i += 1
>>>　　　　　　　break
>>>　　　　　else:
>>>　　　　　　　i += 1
>>>　　　　　　　continue
>>>　　　else:
>>>　　　　　print(num, ' = NG') 　　　　 # 予想の反例
>>>　　　print()
>>>
>>>　stop = time.now()　　　　# 計算終了時刻マーク
>>>　t = stop - start
>>>　print('計算時間',t)

コードファイルのダウンロードはコチラ（https://drive.google.com/file/d/1B5Iol8hs3V7x67PLh7R2YuwsoUnhnCAC/view?usp=sharing）から

次は1000以上1030以下の偶数についての結果です。

6以上100万以下についての検証結果

maxが10万オーダーであれば計算時間1分ほどです。そこでmaxを100万（10の6乗）としてgoldbach2.pyとgoldbach3.pyを走らせてみました。

goldbach2.py　計算時間 1:05:27.289655
goldbach3.py　計算時間 0:50:40.932465

1時間ほどかかりました。この差は[2]の素数リスト生成の差です。

コンピューターによるゴールドバッハ予想の検証は、2012年に4×10の18乗までという記録があります。