目次
コラッツ予想
前回はコラッツ・角谷予想を紹介し、予想を検証するプログラムをコーディングしてみました。
任意の自然数に対し、偶数なら2で割り、奇数なら3倍し1を足す、とい8うシンプルなアルゴリズムによって生成される数列──コラッツ・シークエンスは必ず1で終わるだろう。
これがコラッツ・角谷予想です。
1937年にローター・コラッツによって提起されて以来、今日まで証明されていない数論の未解決問題です。
前回のプログラム「col.gy」では次の4つを算出します。
- 与えられたnに対するコラッツ・シークエンス生成
- a≦n≦bに対するコラッツ・シークエンス生成
- a≦n≦b（数式入力可）に対するステップ数表示＋時間計測
- コラッツ・アルゴリズムを変えた場合
いくつもの自然数に対してプログラムを実行させて実験を行ってみることで、予想の検討と同時にコラッツ・シークエンスの挙動について分かってくることがあります。
コラッツ・シークエンスのステップ数に注目する
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を用いて、手計算では味わうことができない面白さを味わってみてください。
