ラマヌジャンのタクシー数 -Elixirの無限リストを使って-

はじめに

Elixirのストリームの習熟のためにラマヌジャンのタクシー数を題材にして無限リスト生成に取り組みました。

天才数学者ラマヌジャン

インド生まれの数学者ラマヌジャンには天才的な才能がありました。夢の中で女神が教えてくれたという不思議な数式は証明は与えられてないものの、その多くは正しく不思議な結果であったといいます。そんなラマヌジャンが入院中に見舞に訪れた数学者ハーディーとの会話で次のような有名な逸話があります。

徳島大学のページから引用させていただきました。
http://www-math.ias.tokushima-u.ac.jp/~katayama/suori/taxi.pdf

1918 年の冬には，ラマヌジャンは既にイギリスで療養所に入っており，そこへ見舞いに来たハーディー が，ラマヌジャンに次のように言ったそうです． 「乗ってきたタクシーのナンバーを見たが 1729 だった．特に特徴のないつまらない数字だったよ．」 これを聞いたラマヌジャンは，即座に次のように答えたという. 「いいえ，それはとても興味深い数字です．なぜならそれは 2 通りの立方数の 2 つの和で表せる最小の 数だから．」 即ち 1729 は， 1729 = 12^3 + 1^3 = 10^3 + 9^3で 1729 より小さな数 n で n = a^3 + b^3 = c^3 + d^3 (a,b) /= (c,d) 1 ≤ a ≤ b, 1 ≤ c ≤ d となることはないと即答したというのです。

1729(1^3+12^3 = 9^3+10^3)はラマヌジャンのタクシー数と呼ばれています。

素朴なコード

Elixirでこのタクシー数を求めることを考えます。有限なリストを使えば例えば次のようになります。

defmodule Ram1 do
  def taxi(n) do
    for p <- 1..n,
        q <- 1..n,
        r <- 1..n,
        s <- 1..n,
        p<q,r<s,p<r,
        p*p*p+q*q*q == r*r*r+s*s*s do {{p,q},{r,s}} end 
    |> Enum.take(1)
  end
end


iex(1)> Ram1.taxi(20)
[{{1, 12}, {9, 10}}]
iex(2)>

変数p,q,r,sについて１からnまでの数を割り当てます。そして条件に合うものをリストとして生成します。一番小さいものですからEnum.take/1を使って1番目のものを求めます。これがタクシー数です。

無限リストを使って

さて、ここからがいよいよ本番です。調べる変数の範囲を限定せず、無限リストを生成することによってタクシー数を求められないだろうか？という試みです。

Elixirには多様なストリームがライブラリとして収録されています。ストリームはその都度、値が計算されます。０から無限大のストリームであれば０，１，２・・・と必要に応じて生成します。コンピューターのメモリは有限だからです。Streamモジュールの中のunfoldというのが今回のケースで使えそうです。

unfold

次の例はDrew Olsonさんが説明に使った例です。（参考文献）

fibs = Stream.unfold({1, 1}, fn {a, b} ->  
  {a, {b, a + b}}
end)

fibs \
|> Enum.take(10) \
|> IO.inspect

# => [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]

ストリームの初期値と後続の値の計算方法を記述した関数を与えることによりフィボナッチ数列を生成するストリームを記述しています。

話は単純ではなかった

変数p,q,r,sに１から無限大まで数を割り当てていくストリームを作ればいいのだから、簡単なように思えました。初期値を1とし、それに１を加算することにより次々と値を生成すれば良さそうにも思えます。しかし、それでは最初の変数pを調べるのにひたすら無限大を追跡してしまい、いつまでもq,r,s変数の吟味には取り掛かることができません。

これを解決するには p,q,r,sを　(1,1,1,1) (2,1,1,1) のようにまとめて小さな数から生成するストリームが必要でした。各変数の合計数が４、５、６．・・・になるように値を生成してくれるものが必要です。さて、これをどうやって実装したらいいでしょうか？

格子点の列挙問題

このストリームがやろうとしていることは格子点を求めることでした。x,yの二次元のグラフですと(1,0)　(2,1)などが格子点です。その格子点を原点に近いところから順次取り出していくものです。これの4次元バージョンを作れば良さそうです。

書棚にあった古いHaskellの入門書が目に留まりました。何気なくページをめくったところが、正にこの問題の無限ストリームをHaskellでどうするか？という説明でした。Haskellは遅延評価の仕組みの言語なのでうまいやり方があります。しかし、Elixirは基本的には正格言語であり部分的にストリームに遅延評価を採用しています。Haskellの方法とは違う方法を考える必要がありました。

本にあったヒントをもとに考えてみました。

変数がs１つの一次元なら簡単です。０，１，２、と生成すればいいです。
変数がr,sの２つの二次元で、値の合計数がｎの場合を考えてみます。変数rにまずｎを与えると残りのsはn-1の場合を考えればいいことになります。さらにrがn-1、n-2の場合を考えます。
以下同様に変数３つの場合、４つの場合を再帰的に考えていきます。

具体例

二次元の場合

(1,0)->(0,1)->(2,0)->(1,1)->(0,2)->(3,0)->(2,1)->(1,2)->(0,3)-> ...

三次元の場合

(1,0,0)->(0,?,?){二次元に任せる} -> (2,0,0)->(1,?,?){二次元に任せる} ...

四次元の場合

(0,0,0,0)->(1,0,0,0)->(0,?,?,?)｛三次元に任せる｝->(2,0,0,0)->(1,?,?,?)｛三次元に任せる｝

こんな風にして合計数がｎで四次元の無限リストを生成することを考えました。

試行錯誤の末に完成したコードは下記の通りです。

defmodule Lat do
  def lattice4() do
    Stream.unfold([0,0,0,0], fn([p,q,r,s]) -> {[p,q,r,s],next4([p,q,r,s])} end)
  end

  def next2([r,0]) do [r-1,1] end
  def next2([r,s]) do [r-1,s+1] end

  def next3([q,0,0]) do [q-1,1,0] end
  def next3({q,0,s}) do [q-1,s+1,0] end
  def next3([q,r,s]) do [q|next2([r,s])] end

  def next4([0,0,0,0]) do [1,0,0,0] end
  def next4([0,0,0,s]) do [s+1,0,0,0] end
  def next4([p,0,0,0]) do [p-1,1,0,0] end
  def next4([p,0,0,s]) do [p-1,s+1,0,0] end
  def next4([p,q,0,s]) do [p,q-1,s+1,0] end
  def next4([p,q,r,s]) do [p|next3([q,r,s])] end

end

iex(3)> Lat.lattice4 |> Enum.take(10)
[

  [0, 0, 0, 0],
  [1, 0, 0, 0],
  [0, 1, 0, 0],
  [0, 0, 1, 0],
  [0, 0, 0, 1],
  [2, 0, 0, 0],
  [1, 1, 0, 0],
  [1, 0, 1, 0],
  [1, 0, 0, 1],
  [0, 2, 0, 0]
]

無限リストバージョンのタクシー数

さあ、準備はできました。無限リストを使ってラマヌジャンのタクシー数を求めるコードは次のようになります。

defmodule Ram do
  def taxi(n) do
    Lat.lattice4
    |> Stream.filter(fn([p,q,r,s]) -> taxi_number(p,q,r,s) end)
    |> Enum.take(n)
  end

  def taxi_number(p,q,r,s) do
    if p<q && r<s && p<r && p*p*p+q*q*q == r*r*r+s*s*s do
         true
    else
        false
    end
  end
end

taxi(n) の引数ｎは小さいほうからｎ番目のタクシー数を意味しています。1番小さなのが1729になっているばずです。

iex(46)> Ram.taxi(1)
[[1, 12, 9, 10]]

1^3+12^3 = 9^3+10^3 = 1729 です。うまくいきました。

その次はどうでしょう？

iex(47)> Ram.taxi(2)
[[1, 12, 9, 10], [2, 16, 9, 15]]

2^3+16^3 = 9^3+15^3 = 4104

どうやら、うまくいっているようです。

参考文献

[翻訳] Elixirのストリーム
https://qiita.com/HirofumiTamori/items/abf9a9478bfc1161000c

「入門Haskell　はじめて学ぶ関数型言語」　向井　淳　著