Optimizing Numerical Code

From: David Steuber (david_at_david-steuber.com)
Date: 08/31/04 

Date: 30 Aug 2004 20:07:07 -0400

I need to speed up the following function:

(defun double-float-escape-fn (c-real c-imag)
  "Return iterations for c-real + c-imag i as floating-point numbers"
  (declare (type double-float c-real c-imag)
           (optimize (speed 3) (safety 0) (space 0)))
  (loop for iterations below *max-iterations*
     for period from 1
     for z-imag = 0.0D0 then (+ (* 2.0D0 z-real z-imag)
  c-imag)
     for z-real = 0.0D0 then (+ (- z-real-squared
  z-imag-squared) c-real)
     for z-real-squared = 0.0D0 then (* z-real z-real)
     for z-imag-squared = 0.0D0 then (* z-imag z-imag)
     with a-imag = 0.0D0
     with a-real = 0.0D0
     with check = 3
     unless (< (+ z-real-squared z-imag-squared) 4.0D0) return (the
  fixnum iterations)
     when (and (= a-imag z-imag) (= a-real z-real) (> iterations 1))
  return (the fixnum *max-iterations*)
     when (> period check) do (setf period 1 check (* 2 check) a-imag
  z-imag a-real z-real)
     finally (return (the fixnum iterations))))

The function calculates the escape value for the Mandelbrot set. I
use the returned value to decide what color to paint a pixel. I am
considering doing this function in assembler and using FFI. The only
problem with that is that it locks me into an architecture and
implementation. So before I go that route, I want to see what else
can be done in Lisp land.

Here is a sample timing of the function from SBCL on OS X running on a
PowerBook G4 866Mhz:

* (svref *real-array* 260)

-0.10488076384628028d0
* (svref *imag-array* 145)

0.9256099657646029d0
* *max-iterations*

2858197
* (time (double-float-escape-fn (svref *real-array* 260) (svref
*imag-array* 145)))

Evaluation took:
                 3.51 seconds of real time
                 2.35 seconds of user run time
                 1.07 seconds of system run time
                 0 page faults and
                 204,646,208 bytes consed.
2858197

This was one pixel in a 320x240 image.

Thanks. 

-- 
An ideal world is left as an excercise to the reader.
   --- Paul Graham, On Lisp 8.1