Gauche練習帳複素数を使ってみる

Gaucheは有理数や複素数といった数値クラスの種類が豊富ですね。Lisp系言語の伝統なのかもしれません。いまどきのプログラミング言語だとComplexクラスのようなライブラリを持っていることも多いですが、(define z 1+2i)のように直接記述できるのはめずらしいように思います。
そんなわけで今回は複素数を使ってみます。複素数といえばマンデルブロ集合です。

複素数関係の主な関数には以下のようなものがあります。

(make-rectangular 1 2)  ; 複素数の作成
;=> 1.0+2.0i

(define z 3+4i)
;=> z
(real-part z)  ; 実数部の取得
;=> 3.0
(imag-part z)  ; 虚数部の取得
;=> 4.0
(magnitude z)  ; 絶対値の計算
;=> 5.0

マンデルブロ集合は、
　 $Z_{0} \quad = \quad 0$
　 $Z_{n+1} \quad = \quad Z_{n}^2 \quad + \quad C \quad (n=0,1,2, \dots)$
の計算を複素数Cに対してmax回数繰り返しても $\small |Z_{n}| \quad > \quad 2$ にならないCの集合、とのことです。また $\small |Z_{n}| \quad > \quad 2$ になるCも、そのときのnの値によって色分けすることが多いですね。Gaucheの複素数クラスは、この計算を実数部、虚数部を意識せずに自然に記述できるので良い感じです。

;; マンデルブロ集合計算本体
(define (mandelbrot C)
  (let loop ((Zn 0) (n 0))
    (cond ((>= n calc-max) 0)
          ((> (magnitude Zn) 2.0) n)
          (else (loop (+ (* Zn Zn) C) (+ n 1))))))

これを横軸に実数部、縦軸に虚数部をとった複素平面上に描画するとこうなります。

(use srfi-1)
(use sdl)

;; 描画領域定義
(define screen #f)
(define scrx 400)
(define scry 400)

;; SDLユーティリティ
(define (init-screen)
  (sdl-init SDL_INIT_EVERYTHING)
  (set! screen (sdl-set-video-mode scrx scry 0 SDL_SWSURFACE)))
(define (refresh-screen)
  (sdl-update-rect screen 0 0 0 0))
(define (wait-for-quit)
  (let ((e (sdl-make-event)))
    (while (not (= (sdl-event-type e) SDL_QUIT))
           (sdl-poll-event e))))

;; パレット初期化
(define calc-max 200)
(define palette-table
  (list->vector
   (iota calc-max 0 (truncate (/ #xFFFFFF calc-max)))))

;; マンデルブロ集合計算本体
(define (mandelbrot C)
  (let loop ((Zn 0) (n 0))
    (cond ((>= n calc-max) 0)
          ((> (magnitude Zn) 2.0) n)
          (else (loop (+ (* Zn Zn) C) (+ n 1))))))

;; 描画ループ
(define (draw-mandelbrot r_start i_start step)
  (define (plot x y n)
    (put-pixel screen x (- scry y 1) (vector-ref palette-table n)))
  (define (draw x y)
    (let ((z (make-rectangular
              (+ r_start (* x step))
              (+ i_start (* y step)))))
      (plot x y (mandelbrot z))
      '()))
  (let ((x_list (iota scrx 0))
        (y_list (iota scry 0)))
    (map (lambda (y)
           (map (lambda (x) (draw x y)) x_list)
           (refresh-screen)
           '())
         y_list)))

;; 実行
(init-screen)
(draw-mandelbrot -2.0 -2.0 0.01)
(wait-for-quit)

今回「プログラミングGauche」にあった「Lisp脳」の謎に迫る - Schemeプログラマの発想をまねて、ループ部分にdoではなくiotaを使ってリスト処理にしてみましたが、手続き的な書き方に慣れてるとやっぱりメモリ効率とか不安ですね。まだまだ練習が必要だ。