はじめに
ラマヌジャンのタクシー数に関連して、とてもエレガントなPrologコードを発見しましたので、書き留めます。
ラマヌジャンのタクシー数とは
インド生まれの数学者ラマヌジャンには天才的な才能がありました。夢の中で女神が教えてくれたという不思議な数式は証明は与えられてないものの、その多くは正しく不思議な結果であったといいます。そんなラマヌジャンが入院中に見舞に訪れた数学者ハーディーとの会話で次のような有名な逸話があります。
徳島大学のページから引用させていただきました。
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)はラマヌジャンのタクシー数と呼ばれています。
解法
これをコンピューターで求めるとして、どのように計算したらいいでしょうか？ 素朴な方法として、解の範囲を推測してその範囲をしらみつぶしに調べるという方法があります。C言語など命令型の言語であれば for構文を利用する方法です。
for(i=1;i<100;i++)
for(j=1;j<100;j++)
for(k=1;k<100;k++)
for(l=1;l<100,l++)
if(.... )
上限を限定せずに探すにはどうしたらいいでしょうか？ ストリームが使える言語ならストリームを作る方法がよさそうです。Prologにはストリームはない代わりにバックトラックがあります。これをうまく使えないだろうか？
私の頭では思いつかなかったので、ググってみました。とてもエレガントなコードを発見しました。
作者は犬童健良さんという方でした。http://www.xkindo.net/cog_dec/author.html
犬童さんのコードを引用させていただきます。
% the Hardy–Ramanujan number
% r1729.pl
% 3 Jan 2011
% reference:
% Hardy, G. H.:
% Ramanujan: Twelve Lectures on Subjects Suggested by His Life and Work, 3rd ed.,
% New York: Chelsea, 1999.
% the Hardy–Ramanujan number (the Taxi-cab numbers)
% Z = X^3 + Y ^3 = P^3 + Q^3,
% where X, Y, P, Q > 0 are integers, X < Y, P < Q, and not (X = P).
% (i.e., sums of 2 cubes in more than 1 way.)
i( X, Y):- length(L, Y), nth0( X, L,_).
c( X, Z, Y):- i( Z, Y), X is Z ^3.
d( X, Y, [Z, W]):-
c( P, Z, Y),
c( Q, W, Y),
P < Q,
X is P + Q.
r(X, Y, [Z, W]):-
d(X, Y, Z),
d(X, Y, W),
Z @< W .
| ?- r(X,Y,Z).
X = 1729
Y = 12
Z = [[1,12],[9,10]]
お見事！としか言いようがありません。エレガントです。肝心な部分はi/2です。格子点を生成しています。
格子点
いろいろ調べたり、考えたるうちに、格子点の列挙問題だということがわかりました。二次元グラフの原点に近いところから順次、(1,1)(2,1)...と数の組を生成していきます。これはグラフ上の格子の部分を表しています。
| ?- i(X,Y).
X = 1
Y = 1;
X = 2
Y = 2;
X = 1
Y = 2;
...
ISO-Prolog規格でのlength/3は面白い定義になっていて、変数が与えられた場合には順次リストと整数を生成することとなっています。
O-Prolog Ver 1.62(Chika)
| ?- length(X,Y).
X = []
Y = 0;
X = [_]
Y = 1;
X = [_,_]
Y = 2;
X = [_,_,_]
Y = 3;
これをうまく応用して格子点を生成しているのでした。
ISO-Prologの範囲内で
オリジナルのコードではi/2を記述するのにnth0/2を利用しています。これはSWI-Prologの独自拡張のようです。ISO-Prologには含まれていません。そこで、この部分を次のように書き換えました。
i(X,Y) :- length(L,Y),count(L,X).
count(L,N) :-
length(L,N).
count([L|Ls],N) :-
count(Ls,N).
さあ、うまく動作するでしょうか？ 自作のO-Prologで試してみました。
O-Prolog Ver 1.62(Chika)
| ?- ['taxi.pl'].
yes
| ?- r(X,Y,Z).
X = 1729
Y = 12
Z = [[1,12],[9,10]]
yes
|
成功です。
終わりに
エレガントなPrologコードを公開されている犬童さんに大拍手をお送りします。ブラボー！