ラマヌジャン

haskell

わたしはやはり小学時代にラマヌジャンの逸事を学んだ。ラマヌジャンは嘗て病院の中に一人の数学者と一しょになった。「1729という数字って何の変哲もないよね。」というふりに、「いやいや、二つの正整数の三乗の和で二通りに書ける最小の数だよ。」と答えた。−−こう云う逸事を学んだのである。
当時のわたしはこの逸事の中にラマヌジャンの鬼才ぶりを発見した。少くとも発見する為に努力したことは事実である。しかし今は不幸にも寸毫の教訓さえ発見出来ない。この逸事の今のわたしにも多少の興味を与えるは僅かに下のように考えるからである。


数学者のふりは如何に無茶を極めていたか！

できの悪いパロディー*1は放置して、ラマヌジャンを超えようという企画が、http://d.hatena.ne.jp/qqqlxl/20070604/p3 あたりであるようです。

つまり、二つの正整数の4乗和で二通りに書ける数を出力してみよう、ということ。

まあ Haskell 使ってみますよね。

手始めに
http://d.hatena.ne.jp/qqqlxl/20070601/p2
を移植したっぽい感じなのを。

import List
frommathematica = show$filter(\((x,_,_),(y,_,_))->x==y)$zip al$tail al
  where al = sort [(x^4+y^4,x,y)|x<-[1..200],y<-[1..x]]

まあ、これはこれでよい。ただ、これは(x,y)が200以下の範囲に存在していることを「知っている」から書けるのであるから、これはルール違反であろう。

もちろん、可算の範囲から探すのですから、ある程度\lambda x y . x^4+y^4に関する知識を使わなくてはなりません。
今回は、これがx,yの両方向に単調増加であることを仮定しました。


そりゃ無限リスト使いますよね。

http://d.hatena.ne.jp/qqqlxl/20070604/p3

main = print $ show $ take 10 $ filter (\((x,_,_),(y,_,_))->x==y) $ zip ml $tail ml
 where ml = merge [[(x^4+y^4,x,y)|y<-[1..x]]|x<-[1..]]
merge ([]:ys) = merge ys
merge (xx@(tx@(_,_,b):xs):ys)
  | b == 1   = tx:merge(xs:ys)
  |otherwise = mer (xx,yy)
    where yy = merge ys
	  mer ([],wy)=wy
          mer (wx@(ux@(ua,_,_):uxs),wy@(uy@(ub,_,_):uys))
	   | ub < ua  = uy : mer(wx,uys)
	   |otherwise = ux : mer(uxs,wy)

実行すると小さい順に10個出力します。

やっていることは、[ [(1,1)],[(2,1),(2,2)],[(3,1),(3,2),(3,3)]....] というものに、(1^4+1^4,1,1) というように比べる数を付加して、それによって、並べ替えて ml [(2,1,1),(17,2,1),(32,2,2),(82,3,1),(97,3,2),(162,3,3),(257,4,1),(272,4,2)...]を作ります。
あとは、これとこれをひとつずらしたものを並べて、値が同じものがあったらそれを出力するという寸法です。

普通に、可算個のリスト*2をマージソートすると困ったこと(最小を探すために、可算個のリストを繰り始める)になりますが、それを防ぐのは

  | b == 1   = tx:merge(xs:ys)

の行が担当しています。つまり、(_,1) というものがきたら、それ以降は繰らなくていいと教えているのです。



さて、

ml1 =  merge1 [[(x^4+y^4,x,y)|y<-[x..]]|x<-[1..]]
merge1 (xx@(tx@(_,a,b):xs):ys)
  | b == a   = tx:merge1(xs:ys)
  |otherwise = mer1 (xx,yy)
    where yy = merge1 ys
          mer1 (wx@(ux@(ua,_,_):uxs),wy@(uy@(ub,_,_):uys))
	   | ub < ua  = uy : mer1(wx,uys)
	   |otherwise = ux : mer1(uxs,wy)

で、置き換えたものも考えましたが、実はやってみるとこっちはだいぶ時間がかかるのです。


実はこれ、半径r、角度0-45度の円弧を0度方向に射影するか90度方向に射影するかを考えたときに0度方向に射影したほうが短いからのようなんです。

20個まで出力した結果は以下。

((635318657,59,158),(635318657,133,134)),
((3262811042,7,239),(3262811042,157,227)),
((8657437697,193,292),(8657437697,256,257)),
((10165098512,118,316),(10165098512,266,268)),
((51460811217,177,474),(51460811217,399,402)),
((52204976672,14,478),(52204976672,314,454)),
((68899596497,271,502),(68899596497,298,497)),
((86409838577,103,542),(86409838577,359,514)),
((138519003152,386,584),(138519003152,512,514)),
((160961094577,222,631),(160961094577,503,558)),
((162641576192,236,632),(162641576192,532,536)),
((264287694402,21,717),(264287694402,471,681)),
((397074160625,295,790),(397074160625,665,670)),
((701252453457,579,876),(701252453457,768,771)),
((823372979472,354,948),(823372979472,798,804)),
((835279626752,28,956),(835279626752,628,908)),
((1102393543952,542,1004),(1102393543952,596,994)),
((1382557417232,206,1084),(1382557417232,718,1028)),
((1525400095457,413,1106),(1525400095457,931,938)),
((2039256901250,35,1195),(2039256901250,785,1135))

それぞれにかかった時間は
15個
(36.44 secs, 2158601256 bytes)
(561.03 secs, 8606960412 bytes)
20個
(71.09 secs, 4102552440 bytes)
(-170.92 secs, 16840743808 bytes)
なんと時間が桁あふれしました(顔。

*1:芥川の侏儒の言葉を参照 http://www.aozora.gr.jp/cards/000879/files/158_15132.html

*2:無限リストではなくて、可算個のリスト