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

Pythonで数論!未解決問題「コラッツ・角谷予想」|Pythonで数学を学ぼう! 第9回

コラッツ予想

前回はコラッツ・角谷予想を紹介し、予想を検証するプログラムをコーディングしてみました。

任意の自然数に対し、偶数なら2で割り、奇数なら3倍し1を足す、とい8うシンプルなアルゴリズムによって生成される数列──コラッツ・シークエンスは必ず1で終わるだろう。

これがコラッツ・角谷予想です。

1937年にローター・コラッツによって提起されて以来、今日まで証明されていない数論の未解決問題です。

前回のプログラム「col.gy」では次の4つを算出します。

  1. 与えられたnに対するコラッツ・シークエンス生成
  2. a≦n≦bに対するコラッツ・シークエンス生成
  3. a≦n≦b(数式入力可)に対するステップ数表示+時間計測
  4. コラッツ・アルゴリズムを変えた場合

いくつもの自然数に対してプログラムを実行させて実験を行ってみることで、予想の検討と同時にコラッツ・シークエンスの挙動について分かってくることがあります。

コラッツ・シークエンスのステップ数に注目する

nが27の場合に1に至るまで111ステップかかります。

そこで、nの値を大きくしていくとステップ数がどのように大きくなるのかが気になるでしょう。プログラム「col.gy」の3を選択することでステップ数だけを知ることができます。

次は計算結果の1例です。

予想に反し、nを大きくしてもステップ数はさほど大きくはなりません。11桁のn(10の10乗)に対してもステップ数は高々100台です。

そこで、a≦n≦bのnに対するコラッツ・シークエンスのステップ数の最大値を求めるプログラムをコーディングしてみます。

さらに、a≦n≦bのnに対するステップ数の頻度を調べるため、グラフを描くプログラムもコーディングしてみます。

プログラム「colg.py」

>>> # colg.py
>>> # ステップ数(最大ステップ数)と頻度のグラフ描画
>>> import numpy as np         # 配列を扱う数値計算ライブラリNumPy
>>> import matplotlib.pyplot as plt      # グラフ描画ライブラリmatplotlib
>>> import japanize_matplotlib       # matplotlibの日本語化
>>> import datetime as dt
>>> from decimal import Decimal
>>>
>>> while(1):
>>>  Model = input('1.a≦n≦bに対するコラッツ・シークエンスのステップ数の最大値とnを算出\r\n'
>>>         '2.a≦n≦bに対する横軸ステップ数、縦軸頻度の棒グラフ描画\r\n'
>>>         '1、2のどれかを入力>>> ')
>>> if Model.isdecimal():
>>>  Model = int(Model)
>>>  if 1 <= Model <= 2:
>>>   break
>>>
>>>
>>> # 横軸 ステップ数s s=0,1,...,s_max-1
>>> s_max = 1000
>>> s = np.empty(s_max)
>>>
>>> S = np.empty(s_max)
>>>
>>> for i in range(s_max):
>>>   s[i] = i
>>>   S[i] = 0
>>>
>>> # ステップ数が最大となるとき、nとステップ数を格納するリスト
>>> M = np.empty(2)
>>> M[0] = 0
>>> M[1] = 0
>>>
>>> def col(n):
>>>   N = n
>>>   c = 0
>>>   while True:
>>>     if n == 1:
>>>       break
>>>     if n % 2 == 0:
>>>       n = Decimal(n)/Decimal("2")
>>>       c = c + 1
>>>     else:
>>>      n = Decimal("3")*Decimal(n)+Decimal("1")
>>>      c = c + 1
>>>   if c < s_max :
>>>     S[c] = S[c] + 1
>>>   if c > M[1] :
>>>     M[0] = N
>>>     M[1] = c
>>>
>>> if Model == 1:
>>>   a = eval(input('nの始まりは 数式OK '))
>>>   b = eval(input('nの終わりは 数式OK '))
>>>
>>>   time = dt.datetime
>>>   start = time.now() # 時間計測開始
>>>   for i in range(a, b + 1):
>>>     col(i)
>>>   print('n=',int(M[0]),'max(Step)=',int(M[1]))
>>>   stop = time.now() # 時間計測終了
>>>   t = stop - start
>>>   print('計算時間',t)
>>>
>>> elif Model == 2:
>>>   a = eval(input('nの始まりは 数式OK '))
>>>   b = eval(input('nの終わりは 数式OK '))
>>>
>>>   #start = time.time() # 時間計測開始
>>>
>>>   #print('n'+'\t'+'Step')
>>>   for i in range(a, b + 1):
>>>     col(i)
>>>   print('n=',int(M[0]),'max(Step)=',int(M[1]))
>>>
>>>   #t = time.time() - start # 時間計測終了
>>>   #print('計算時間',t,'秒')
>>>   plt.bar(s, S, color = "blue")
>>>   plt.title('コラッツ予想の検証 コラッツ・シークエンス ステップ数の分布({}≦n≦{})'.format(a,b))
>>>   plt.xlabel('Step数')
>>>   plt.ylabel('頻度')
>>>   plt.xticks(np.arange(0,s_max+0.1,50))
>>>   plt.grid(True)
>>>   plt.show()

※プログラムファイルはこちらからダウンロード(リンク先のファイルを別名でダウンロード)できます。

時分秒表示を行うdatetimeモジュール

前回のプログラム「col.gy」ではtimeモジュールを用いて秒の表示を行いました。

今回のプログラム「colg.gy」ではdatetimeモジュールを用いることで時分秒の表示を行うようにしました。

計算誤差を防ぐdecimalモジュール

pythonをノーマルで使用すると計算誤差が生じます。ある意味pythonの宿命的欠点です。

今回のプログラム「colg.gy」では大きな数を扱うので、計算誤差が問題になってきます。

そこで、任意精度の数値計算を行うdecimalモジュールの出番です。
>>> from decimal import Decimal
としてインポートします。

プログラム「col.gy」のn = n/2とn = 3*n + 1のコードを
n = n/2 → n = Decimal(n)/Decimal("2")
n = 3*n + 1 → n = Decimal("3")*Decimal(n)+Decimal("1")
のようにすれば正確に(デフォルトは28桁の精度)計算してくれます。

プログラム「colg.gy」の実行結果

$ python colg.py
1.a≦n≦bに対するコラッツ・シークエンスのステップ数の最大値とnを算出
2.a≦n≦bに対する横軸ステップ数、縦軸頻度の棒グラフ描画
1、2のどれかを入力>>> 1
nの始まりは 数式OK 1
nの終わりは 数式OK 10**7
n= 8400511 max(Step)= 685
計算時間 0:21:02.099879

$ python colg.py
1.a≦n≦bに対するコラッツ・シークエンスのステップ数の最大値とnを算出
2.a≦n≦bに対する横軸ステップ数、縦軸頻度の棒グラフ描画
1、2のどれかを入力>>> 1
nの始まりは 数式OK 1
nの終わりは 数式OK 10**8
n= 63728127 max(Step)= 949
計算時間 4:07:28.299884

1から10の7乗までのnの場合、ステップ数の最大値は685、計算時間は21分、1から10の8乗までのnの場合、ステップ数の最大値は949、計算時間は4時間7分でした。

私はpythonとNAS専用マシンにubuntu、Core i3マシンを使っています。先日、さらに記録更新を狙って1から10の13乗まで実行させました。1週間連続計算を続けても終了しませんでした。

確かに、10の8乗の場合の10の5乗倍(10万倍)だとすると4時間の10万倍で40万時間すなわち45年かかる計算になります。単純計算(アルゴリズム)の限界です。

記録を狙うには、ハードと理論武装の両面を充実させる必要があります。円周率や素数の計算と同じです。ここにコラッツ・シークエンスのステップ数の記録があります。

nが597京7996兆3043億4350万1855(19桁)のとき2389ステップが最高記録です。

このようにコラッツ・シークエンスのステップ数は高々数千です。そこで、a≦n≦bに対するステップ数の頻度を調べるプログラムをコーディングしてみます。

プログラム「colg.gy」の2を選ぶと、a≦n≦bに対するnのステップ数を計算することでステップ数毎の頻度を集計します。

次がその実行結果です。

$ python colg.py
1.a≦n≦bに対するコラッツ・シークエンスのステップ数の最大値とnを算出
2.a≦n≦bに対する横軸ステップ数、縦軸頻度の棒グラフ描画
1、2のどれかを入力>>> 2
nの始まりは 数式OK 1
nの終わりは 数式OK 10**8
n= 63728127 max(Step)= 949

1から10の6乗までのnについてステップ数を調べると、200あたりに集中していることが分かります。

2020年に2の68乗(21桁)まで検証はされています。予想は正しいように思われます。予想が成り立たないと考える数学者はいなのではないでしょうか。

しかし証明されていない以上、さらにその先の大きな数において反例が見つかる場合もあります。

コラッツ・角谷予想の解決者に賞金500ドルを与えると言ったのが、数学者ポール・エルデシュ(1913-1996)です。彼は予想について「数学はまだこの種の問題に対する用意ができていない」と語りました。

コラッツ・角谷予想は小学生でも分かる内容にもかかわらず、その正体が魑魅魍魎(ちみもうりょう)である未解決難問であるのが面白いところです。

Pythonを用いて、手計算では味わうことができない面白さを味わってみてください。

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

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

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

あわせて読みたい

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