コラッツ予想の計算用コードを遊びで書いていてPythonとRubyの比較が気になって興味がズレていった話

普段あまりYouTubeは見ないのですが（時間が溶けるので）、昨日、たまたま関係する動画がオススメされたので観てました。面白いなと思いWikipediaでちらってアルゴリズムをみて、こんなにシンプルなのに動作が荒ぶるんだな。。。と単純に興味がわいてPythonで簡単なコードを書き始めました。

a = []

for i in range(1, 10000):
    x = i
    c = 0
    while x != 1:
        if x % 2 == 1:
            x = x * 3 + 1
        else:
            x = int(x / 2)
        c += 1
    #print(c)
    a.append(c)

a.sort()

with open('a.txt', 'w') as f:
    f.writelines("\n".join(list(map(str, a))))

10⁶くらいになると単純な計算では時間がかかるようになり、メモ化するしかないな、と以下のようなコードを書きました。お、少し早くなった。おもろいなと。趣味の世界です。

a = []
memo = {}

for i in range(1, 100000000):
    x = i
    c = 0
    while x != 1:
        if x in memo:
            c += memo[x]
            break

        if x % 2 == 1:
            x = x * 3 + 1
        else:
            x = int(x / 2)
        c += 1
    #print(c)

    a.append(c)
    memo[i] = c

a.sort()

with open('a.txt', 'w') as f:
    f.writelines("\n".join(list(map(str, a))))

ちなみにこれをグラフ化するコードは以下です。最初、Excelでやってたのですがデータ数が多くなると表示できなくなったのでmatplotlibの力を借りました。簡単でいいですね〜Pythonは。ほんと素晴らしい。

import numpy as np
import matplotlib.pyplot as plt

with open('a.txt') as f:
    a = f.readlines()

a = np.array(a)

plt.hist(a, bins=100)
#plt.plot(a)

#ax = plt.gca()
#ax.spines['top'].set_color('none')
#ax.set_yscale('log')

plt.show()

繰り返し回数をグラフにしたり、分布を見てみたりしながら遊んでいました。

軸ラベルがないのでわかりにくいですが、以下は1から10⁸までのそれぞれの数値が1に至るまでの計算回数です。最初は100回くらいですが、10⁸あたりから500回くらいまで増えていきます。

以下は計算回数の分布です。Wikipediaにも同じようなグラフがありましたがほぼ同じですね。

偶数であれば2で割るので元より小さくなるけど、2で割った結果が奇数であれば3倍＋1になるので元の値より大きくなります。2で割ったときに偶数になるか奇数になるかは数によるので、素数と同様に数の分布の不思議ということなんでしょうか。素人なので何とも言えませんが、確かに法則らしいものはあるようでありません。面白い問題だなと思います。

で、PythonでこれならRubyだったらどうなるのかな？とちょっと邪？な考えが浮かびました。PythonのRubyも動的言語ですからそれほど大きな差はないけど。。。Pythonがこれだけ盛り上がっていますから改善も進んでいるのかな。。。と思いました。それでもRubyには負けてほしくはいですしね。

ということで、とりあえずPythonと等価なコードを書いてみました。Profileとりながら微妙な改善をしたあとがありますが、アルゴリズムは同じにしないと評価できないのでほぼ同じにしています。

a = []
memo = {}

#i = 0
#while i < 1000000
#    i += 1
#for i in 1..1000000-1 do
for i in 1..100000000-1 do
    x = i
    c = 0
    while x != 1
        if memo[x]
#        if memo.key?(x)
            c += memo[x]
            break
        end

        if x.odd?
#       if x % 2 == 1
            x = x * 3 + 1
        else
#            x = x >> 1
            x = x / 2
#            x = (x / 2).to_i #Pythonのように実数にはならない
        end

        c += 1
    end
    #print(c)

    a.push(c)
    memo[i] = c
end

a.sort!

IO.write("a.txt", a.join("\n"))

計測結果を以下に示します。ちなみにM1 Macでrbenvの環境で動作させています。試行回数は１回だけです。timeコマンドのuser値を取っています。数回やって平均を取るほど時間がなかったので。。。すんません。

Pythonはpipenv環境です。

今、授業では3.2.2を使用しているためMacのデフォルト(rbenv global)も3.2.2になっています。これだと、圧倒的にPythonに負けてしまいました。YJITってどうなんだろうとやってみたんですがそれほど効きませんでした。関数コールとかありませんし、こういう処理には効かないのかなと仕組みも調べず想像していました。

それでも型推論を動的に行わないことで早くなるという記事も見つけたのでこういう処理でも効くはずだ、と思い。3.3はさらに早くなっているという記事を見つけましたのでダメ元で3.3.3を入れて試したら。。。素晴らしい🎉 Pythonの半分以下です。Rubyのメンテナーの方々、お疲れ様です。ありがとうございました。どこから目線？

ということでコラッツ問題、時間があればまた遊んでみたいと思います。

memorandums

日々の生活で問題解決したこと、知ってよかったことなどを自分が思い出すために記録しています。

コラッツ予想の計算用コードを遊びで書いていてPythonとRubyの比較が気になって興味がズレていった話