;;;;;; Mandelbrot set generation                         -*- Scheme -*-

;;; Taylor Campbell wrote this code; he places it in the public domain.

(declare (usual-integrations))

(define (mandelbrot n)
  (let ((m (fix:- 8 (fix:remainder n 8)))
        (n-float (int:->flonum n))
        (out (current-output-port)))
    (display "P4" out)
    (newline out)
    (write n out)
    (display " " out)
    (write n out)
    (newline out)
    (let loop-y ((y 0) (byte 0) (bit 0))
      (let loop-x ((x 0) (byte byte) (bit bit))
        (let ((proceed (lambda (byte bit)
                         (cond ((fix:< x n)
                                (loop-x (fix:+ x 1) byte bit))
                               ((fix:< y n)
                                (loop-y (fix:+ y 1) byte bit)))))
              (byte (if (mandelbrot? x y n-float)
                        (+ (fix:lsh byte 1) 1)
                        (fix:lsh byte 1)))
              (bit (fix:+ bit 1)))
          (cond ((fix:= bit 8)
                 (write-char (integer->char byte) out)
                 (proceed 0 0))
                ((fix:= x (fix:- n 1))
                 (write-char (integer->char (fix:lsh byte m))
                             out)
                 (proceed 0 0))
                (else
                 (proceed byte bit))))))))

(define-integrable (mandelbrot? x y n)
  (let ((x (int:->flonum x))
        (y (int:->flonum y)))
    (let ((Cr (flo:- (flo:/ (flo:* x 2.0) n) 1.5))
          (Ci (flo:- (flo:/ (flo:* y 2.0) n) 1.0)))
      (let loop ((i 0) (Zr 0.0) (Zi 0.0))
        (let ((Zr (flo:+ (flo:* Zr Zr)
                         (flo:+ (flo:negate (flo:* Zi Zi))
                                Cr)))
              (Zi (flo:+ (flo:* 2.0 (flo:* Zr Zi))
                         Ci)))
          (cond ((flo:> (flo:+ (flo:* Zr Zr)
                               (flo:* Zi Zi))
                        4.0)
                 #f)
                ((fix:> i 50)
                 #t)
                (else
                 (loop (fix:+ i 1) Zr Zi))))))))