この記事はりひにーさんのLisp Advent Calendar 2014 25日目の記事です。
期間を過ぎていますが余っていたので埋めました。
今までProcessingから始まりジャバ、C++、CUDAで並列演算、とやってきたクライン群の極限集合の描画をCommon Lispでやったので記事にします。
アルゴリズムとか数学的に詳しい話はインドラの真珠という本を読んでみてください。
インドラの真珠｜日本評論社
今年はこの本に出会えたのが色々と大きかったです。
何をやっているか
何をやっているかものすごくざっくりと説明すると、
メビウス変換を表す2x2複素数行列、a、bとその逆行列A、Bを様々な組み合わせで用いて無限に変換を行った後の極限点を求めて描画しているという感じでしょうか。
その極限点の集合をどう求めるのかというと、このようなグラフを探索するのです。
各ノードのラベルを有限長なので有限語と呼びます。今回は深さ優先で探索していきます。たどり着いたそれぞれのノードの有限語に無限長の語、aaa...やabABabAB...等の無限語を適用することで極限を求めます。無限語は変換の固定点を用います。例えばaaa...という無限語であれば、aの固定点、abABabAB...という無限語であればabABの固定点を用いる、といった具合です。
あらかじめ決めた深さに達するか、描画するのに十分な深さまで探索されたと判断されたならば、その点を描画し、次の枝に移ります。
実行結果
(defparameter *magnification* 200) (get-limiting-set 2 2 15 0.01) (draw-limiting-set)
(defparameter *magnification* 200) (get-limiting-set #c(1.95 0.02) 3 25 0.01) (draw-limiting-set)
(defparameter *magnification* 8) (get-limiting-set #c(1.91 0.05) #c(1.91 0.05) 20 0.1) (draw-limiting-set)
コードについて
相当酷いコードだと思いますがgistにおいてあります。処理系はclisp-2.49です。
色々な人のlispコードを見てちゃんとしたlispコードの書き方を勉強しなければなりません。
クライン群の極限集合を描画する
注意点ですが、配列のインデックスは１始まりになっています。インドラの真珠に載っていたコードの配列インデックスが１始まりだったからです。
無限大
Common Lispに無限大は規定されていないらしいです。処理系によってはサポートされているらしいのですがCLispにはないみたいです。とりあえずmost-positive-fixnumを無限としておきました。
(defparameter *inf* (* most-positive-fixnum #c(1 1)))
行列クラス
*演算子をオーバーロードできれば良いのですけど、できるかわかりませんでした。mat*という名前で定義しておきましたが気持ち悪い名前だと思います。
print-objectメソッドによってオブジェクトがプリントされる内容を制御できます。
固定点の計算
fix-plusという名前で定義してあります。√の前の符号で二つ固定点が出る時があります。今回は+の方でとってあります。僕が試した範囲では得られる図形に変わりはなかったように思います。
変換の取得
set-gensでは二つのパラメータからメビウス変換の行列を求めています。その中でt_abという変数を求めているのですが、ここも√の前の符号で二つの値がでます。この場合では+をとるか-をとるかで得られる図形が結構変わります。今回は-を取ってあります。
数式の入力
これまで極限集合の描画プログラムを書いてきて結構きつかったのが数式の入力です。複雑で読みにくく、バグの温床です。
lispでは改行を入れることでだいぶ見やすくなると感じました。これはなかなか面白いですね。マクロを使えばべき乗等も含め楽になったりするんですかねぇ。
(setf (z0 (/ (* (- t_ab 2) t_b) (+ (- (* t_b t_ab) (* t_a 2)) (* t_ab #c(0 2))))))
ジャバだとこうでした。めちゃくちゃつらい。
Complex z0 = t_ab.sub(2.0).mult(t_b).div(t_b.mult(t_ab).sub(t_a.mult(2.0)).add(t_ab.mult(new Complex(0, 2.0))));
complex<double> z0 = (t_ab - 2.0) * t_b / (t_b * t_ab - 2.0 * t_a + t_ab * complex<double>(0, 2.0));
末尾再帰の最適化について
極限集合の探索を深さ優先で行う関数を再帰で書いたのですが、ある一定の深さ以上を探索するとデバッガが起動せず、Lisp connection closed unexpectedly: connection broken by remote peerなどといったメッセージとともにSLIMEが落ちてしまいました。きちんとしたエラーメッセージが全く出なかったのですが、Twitterでデバッガに落ちなかったエラーメッセージは*inferior-lisp*バッファに表示されることもあるとのリプライをもらったので確認してみるとスタックオーバーフローのエラーでした。どうやら末尾再帰が最適化されていなかったようです。
末尾再帰の最適化はCommon Lispの規格で決まっているわけではないので処理系に依存しますが、CLispではコンパイル時に最適化されるそうです。compileコマンドで関数のみをコンパイルすることができるようなので、定義後にコンパイルし、最適化してやりました。
こんな感じ
(compile 'explore-limiting-set-with-dfs)
ちなみに最適化する方法がわからなかった時にループで探索を行うコードも書いたのでgistにはコメントアウトして置いてあります。
遅い
遅いです。以前ジャバで組んだものと比べてものすごく遅いです。パラメータによっては全く計算が終わりません。僕の組み方が非効率的なのでしょうけど型や高速化のオプションをつけての高速化でどのくらい速くなるんでしょうか。これは今後の課題ですね。
おわり
Lisperへの道は長く険しい。
