Code Feedback Wanted (Generating more garbage)

From: David Steuber (david.steuber_at_verizon.net)
Date: 03/24/04 

Date: Wed, 24 Mar 2004 10:57:22 GMT

I've done a bit more hacking on my Lisp code for generating the
Mandelbrot Set. I've not been working on it constantly, that's why
I'm not further along than might be expected from the last time I
posted. I'm not using rationals anymore. I'm going with ints as a
way to do fixed radix math. I've not implemented certain features
yet like adjusting the number of bits used based on zoom factor. I'm
also doing a fairly brute-force approach to calculating the pixels.
I am considering using an array to hold the escape value data so that
I can do edge detection on the set so that I don't have to calculate
pixels in the interior. Anyway, here is the Lisp that I have:

;;;; fractal.lisp -- an exploration of the Mandelbrot Set using
;;;; arbitrary sized integer math.
;;;;
;;;; By David Steuber with feedback from comp.lang.lisp
;;;;
;;;; This software is put in the public domain by David Steuber
;;;;
;;;; Just for the record, I hate global variables because they
;;;; break modularity by being non-local.

;;; global variables related to the fractional bits in our fixed radix numbers
(defvar *fractional-bits* 16
  "Number of fractional bits in our number")
(defvar *escape-value* (ash 4 *fractional-bits*)
  "Numbers larger than this have escaped and are not part of the Mandelbrot Set")
(defvar *my-two* (ash 2 *fractional-bits*))

(defun set-fractional-bits (fb)
  "A nice way to set both *fractional-bits* and *escape-value* at the same time"
  (setf *escape-value* (ash 4 fb)
        *my-two* (ash 2 fb)
        *fractional-bits* fb))

(defun shifted-multiply (&rest args)
  "Multiplies the args as integers left-shifted by *fractional-bits* and right-shifts accordingly"
  (let ((n (apply #'* args)))
    (ash n (- (* (max 0 (- (list-length args) 1)) *fractional-bits*)))))

(defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations 10000))
  "Return iterations for c-real + c-imag i as integers left shifted by *fractional-bits*"
  (loop for iterations below max-iterations
     for z-imag = 0 then (+ (shifted-multiply z-real z-imag *my-two*) c-imag)
     for z-real = 0 then (+ (- z-real-squared z-imag-squared) c-real)
     for z-real-squared = 0 then (shifted-multiply z-real z-real)
     for z-imag-squared = 0 then (shifted-multiply z-imag z-imag)
     unless (< (+ z-real-squared z-imag-squared) *escape-value*) return iterations
     finally (return iterations)))

(defun mandelbrot-escape-fn-pc (c-real c-imag &key (max-iterations 10000))
  "Return iterations for c-real + c-imag i as integers left shifted by *fractional-bits* doing period checking"
  (loop for iterations below max-iterations
     for period from 1
     for z-imag = 0 then (+ (shifted-multiply z-real z-imag *my-two*) c-imag)
     for z-real = 0 then (+ (- z-real-squared z-imag-squared) c-real)
     for z-real-squared = 0 then (shifted-multiply z-real z-real)
     for z-imag-squared = 0 then (shifted-multiply z-imag z-imag)
     with a-imag = 0
     with a-real = 0
     with check = 3
     unless (< (+ z-real-squared z-imag-squared) *escape-value*) return iterations
     when (and (= a-imag z-imag) (= a-real z-real) (> iterations 1)) return max-iterations
     when (> period check) do (setf period 1 check (* 2 check) a-imag z-imag a-real z-real)
     finally (return iterations)))

(defun numbers-to-fixed-radix (&rest nums)
  "Return list of numbers converted to fixed-radix integers with *fractional-bits*"
  (mapcar #'(lambda (n)
              (round (* n (ash 1 *fractional-bits*))))
          nums))
  
(defun fixed-radix-to-numbers (&rest nums)
  "Almost the inverse of numbers-to-fixed-radix"
  (mapcar #'(lambda (n)
              (/ n (ash 1 *fractional-bits*)))
          nums))

(defun mandelbrot-set (&key (x 0) (y 0) (zoom 1) (aspect 1)
                       (bitmap-width 100) (bitmap-height 100)
                       (data-file-name "mandelbrot-data.txt")
                       (max-iterations 1000))
  "Calculate bitmap data for the Mandelbrot Set"
  (let ((real-array (make-array bitmap-width)) ; to contain the real values left to right
        (imag-array (make-array bitmap-height)) ; to contain the imag values top to bottom
        (escape-fn #'mandelbrot-escape-fn) ; escape function to use
        (data-file (make-pathname :name data-file-name))
        (iterations 0) ; just another local variable
        (total-iterations 0)) ; useful to know?
    ;; initialize real-array
    (loop for i below bitmap-width
         with size = (/ (* 4 aspect) zoom)
         with pixelsize = (/ size bitmap-width)
         with start = (- x (* (/ bitmap-width 2) pixelsize))
         do (setf (svref real-array i) (car (numbers-to-fixed-radix (+ start (* i pixelsize))))))
    ;; initialize imag-array
    (loop for i below bitmap-height
         with size = (/ 4 zoom)
         with pixelsize = (/ size bitmap-height)
         with start = (- y (* (/ bitmap-height 2) pixelsize))
         do (setf (svref imag-array i) (car (numbers-to-fixed-radix (+ start (* i pixelsize)))))
         finally (nreverse imag-array)) ; this is so the first element is on the top row
    ;; we will start by doing this the brute force way and calculate every damn pixel
    (with-open-file (stream data-file :direction :output
                            :if-exists :supersede)
      (format stream "/* FRACTAL MANDELBROT */~%/* left top right bottom */~%") ; old format for now
      (format stream "~A ~A ~A ~A~%"
              (svref real-array 0) (svref imag-array 0)
              (svref real-array (1- bitmap-width)) (svref imag-array (1- bitmap-height)))
      (format stream "/* center_x center_y aspect zoom */~%")
      (format stream "~A ~A ~A ~A~%" x y aspect zoom)
      (format stream "/* width height maxitt */~%~A ~A ~A~%/* data */~%" bitmap-width bitmap-height max-iterations)
      (dotimes (j bitmap-height)
        (dotimes (i bitmap-width)
          (setf iterations (funcall escape-fn (svref real-array i) (svref imag-array j) :max-iterations max-iterations))
          (if (= iterations max-iterations)
              (setf escape-fn #'mandelbrot-escape-fn-pc)
              (setf escape-fn #'mandelbrot-escape-fn))
          (format stream "~A " iterations)
          (incf total-iterations iterations))
        (format stream "~%")))
    total-iterations))

Here is a test run:

CL-USER> (time (mandelbrot-set :x -3/4 :zoom 2 :bitmap-width 640 :bitmap-height 640))
(MANDELBROT-SET :X -3/4 :ZOOM 2 :BITMAP-WIDTH 640 :BITMAP-HEIGHT 640) took 126,626 milliseconds (126.626 seconds) to run.
Of that, 87,830 milliseconds (87.830 seconds) were spent in user mode
         32,520 milliseconds (32.520 seconds) were spent in system mode
         6,276 milliseconds (6.276 seconds) were spent executing other OS processes.
24,051 milliseconds (24.051 seconds) was spent in GC.
 2,252,959,120 bytes of memory allocated.
145962887

By running this Perl script on the data written out

  http://www.david-steuber.com/~david/files/fractal2image.txt

I generate an XPM file that I converted to this GIF file with Graphic
Converter:

  http://www.david-steuber.com/~david/files/mandelbrot3.gif

This was all done on a PowerBook G4 12" with the 866Mhz G4 and 640MB
RAM. 

-- 
It would not be too unfair to any language to refer to Java as a
stripped down Lisp or Smalltalk with a C syntax.
--- Ken Anderson
    http://openmap.bbn.com/~kanderso/performance/java/index.html

Relevant Pages

  • Re: Douady-Hubbard Potential
    ... I want to draw external rays of Mandelbrot set. ... Mandelbrot set), phican be aproximated by c. ... This last one needs to know about derivatives to get it. ... escape radious big or use many iterations at all. ...
    (sci.fractals)
  • Re: Why does the Mandelbrot set work?
    ... pixels according to how many iterations it takes to bailout. ... Mandelbrot set) was given the problem of iterating a formula like ... would be a series of colored concentric circles. ... the result is most astonishing and unexpected. ...
    (sci.math)
  • Re: Why does the Mandelbrot set work?
    ... The Mandelbrot set has always fascinated me for one ... pixels according to how many iterations it takes to ... the circle and that's it. ... result is most astonishing and unexpected. ...
    (sci.math)
  • Why does the Mandelbrot set work?
    ... pixels according to how many iterations it takes to bailout. ... Mandelbrot set) was given the problem of iterating a formula like ... would be a series of colored concentric circles. ... the result is most astonishing and unexpected. ...
    (sci.math)