コラム 人と星とともにある数学 数学

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乗までという記録があります。

記事の評価をお願いします!

  • この記事を書いた人
  • 最新記事

桜井進(さくらいすすむ)様

1968年山形県生まれ。 サイエンスナビゲーター®。株式会社sakurAi Science Factory 代表取締役CEO。 (略歴) 東京工業大学理学部数学科卒、同大学大学院院社会理工学研究科博士課程中退。 東京理科大学大学院非常勤講師。 理数教育研究所Rimse「算数・数学の自由研究」中央審査委員。 高校数学教科書「数学活用」(啓林館)著者。 公益財団法人 中央教育研究所 理事。 国土地理院研究評価委員会委員。 2000年にサイエンスナビゲーターを名乗り、数学の驚きと感動を伝える講演活動をスタート。東京工業大学世界文明センターフェローを経て現在に至る。 子どもから大人までを対象とした講演会は年間70回以上。 全国で反響を呼び、テレビ・新聞・雑誌など様々なメディアに出演。 著書に『感動する!数学』『わくわく数の世界の大冒険』『面白くて眠れなくなる数学』など50冊以上。 サイエンスナビゲーターは株式会社sakurAi Science Factoryの登録商標です。

-コラム, 人と星とともにある数学, 数学
-

© 2021 株式会社インフォマティクス