Common Lispでマンデルブロ集合を計算してみた

Common Lispで複素数が扱えたので計算してProcessingで描画してみた。
1か月前くらいに下書きを書き始めたのだけれど、改めて記事にしてコードを眺めると粗が目立ったりでどうにも…
コードはここ
https://gist.github.com/soma-arc/9564293

マンデルブロ集合についてはマンデルブロ集合の不思議な世界さんがとても詳しく、わかりやすいです。初めて書いた時にとても参考になりました。

簡単に計算の仕方をまとめると、複素平面上に格子点を取ってその座標を複素数Cとし、以下の漸化式を繰り返し計算する。

が無限大に発散しないものがマンデルブロ集合となる。
また、となると無限大に発散することがわかっている。
今回のマンデルブロ集合の描画は、漸化式の計算回数で色分けを行っている。

(defparameter *width* 500)
(defparameter *height* 500)
(defparameter *reduction* 100.0)

(defun get-latice-points (width height reduction)
  (apply #'append (loop for x from (* -1 (/ width 2)) below (1+ (/ width 2))
                        collect (loop for y from (* -1 (/ height 2)) below (1+ (/ height 2))
                                 collect (complex (/ x reduction) (/ y reduction))))))

第一引数の関数にリストの要素すべてをその関数の引数として渡すapplyとリストの結合を行うappendを使って格子点のリストを得る。
格子点の間隔は１だと広すぎるので*reduction*で割って0.01刻みにしている。
また、0.01などの少数は、二進数にすると無限小数になってしまう場合があり、計算に誤差が出るので注意。0.01を直接かけると誤差が出るので100で割った

(defun calc-mandelbrot (c)
  (labels ((f (z c n)
    (cond ((= n 27) `(,c nil))
          ((< 2 (abs z)) `(,c ,n))
          (t (f (+ c (expt z 2)) c (1+ n))))))
    (f 0 c 0)))

漸化式の計算部分。
２７回計算しても発散しないとわかったら計算を打ち切ってマンデルブロ集合とする。
途中でzの絶対値が２より大きくなれば発散するものとして計算回数をくっつけてCと一緒に返す。

(defun get-mandelbrot ()
  (mapcar #'calc-mandelbrot (get-latice-points  *width* *height* *reduction*)))

(defun to-csv (coordinates-count)
  (let ((coordinates (car coordinates-count))
        (recursion-count (cadr coordinates-count)))
    (concatenate 'string 
                 (write-to-string  (realpart coordinates) ) "\," 
                 (write-to-string  (imagpart coordinates) ) "\," 
                 (write-to-string recursion-count))))

(defun to-csv-file (calculated-latice-points)
  (with-open-file (*standard-output* "mandelbrot.csv"
                                     :direction :output
                                     :if-exists :supersede)
    (mapcan (lambda (s)
              (format t (concatenate 'string s "~%")))
            (mapcar #'to-csv calculated-latice-points))))

残りの部分は格子点リストに対して計算を行う関数、座標と結果のリストをcsvにする関数、そしてそれらをファイルにする関数。
with-open-file中ではprincで出力するとうまく改行が行われなかったので、formatを使い、改行文字と結合してみた。メモ帳で開くと見た目は改行されていないが問題ない。
それとcuspを使ってファイル出力を行う時は絶対パスで指定する必要があるみたいだ。おそらくcuspのインタープリタの場所がカレントディレクトリになっているのだと思う。

(to-csv-file (get-mandelbrot))

あとは手に入ったcsvをごにょごにょして描画しよう。
Processingではこんな感じ。

void setup(){
  size(500, 500);
  translate(250, 250);
  String lines[] = loadStrings("mandelbrot.csv");
  for(String line : lines){
    String tmp[] = line.split(",");
    float x = (float(tmp[0]));
    float y = (float(tmp[1]));
    int count = 0;
    if(!tmp[2].equals("NIL")){
      count = int(tmp[2]);
    }
    stroke(0, count * 255 / 12, count * 255 / 12);
    point(round(x * 100), round(y * 100));
  }
}

縮小してあった100をかけてプロットする。しかし、ここでも誤差が出るらしく、roundで丸める必要がある。

今回は緑と青を混ぜた感じで描画してみた。
今のツイッターアイコンです

黒い部分がマンデルブロ集合。端の部分をうまく拡大できれば自己相似性が見える感じなのだけれど色々面倒くさそうなので僕はまだ試してません。

とりあえず自分なりに一つかけたので次はどうしようか。
Newmanアルゴリズムを一度組んであるのでlispでも組んでみようと思ってはいる。