Mandelbrot with complex numbers: Lisp vs C

cmandelbrot.c

#include <stdio.h>
#include <complex.h>
#include <math.h>

#define MSTEPS 50000
#define MRADIUS 2.0

#define IMG_SIZE 1024

/* Squared absolute value. */
static inline double cabs2(complex z) {
    double x = creal(z), y = cimag(z);
    return x*x + y*y;
}

/* Coloring according to
   http://fractals.nsu.ru/construction/mandelbrot_coloring.html
 */
static inline float mandelbrot(complex z0) {
  /* zp contains producnt of z_1 * z_2 * .. z_{n-1}, where z_n is
     first value out of MRADIUS circle. */
  complex z = z0, zp = 1;
  int i;

  for (i = 0; i < MSTEPS && (cabs2(z) < MRADIUS*MRADIUS); ++i) {
    zp *= z; // Accumulate zp
    z = z*z + z0; // Mandelbrot calculation
  }
  if (i == MSTEPS) {
    return 0; // Point inside
  } else {
    return carg(zp/z);
  }
}

int main(int argc, char *argv[]) {
  double delta = 4.0/(IMG_SIZE-1);
  FILE *outf = NULL;
  char buf;

  char name[]="test.pgm";
  outf = fopen(name, "w");
  fprintf(outf, "P5 %d %d 255\n", IMG_SIZE, IMG_SIZE/2);
 
  for (int ln = 0; ln < IMG_SIZE/2; ++ln) {
    double y0 = ln*delta;

    for (int p = 0; p < IMG_SIZE; ++p) {
      double x0 = (p - IMG_SIZE/2-IMG_SIZE/8)*delta;
      float b = mandelbrot(x0 + I*y0);
    
      if (b)
        buf = floor(255*(0.55 + 0.45*sin(b)));
      else
        buf = 0;
      fwrite(&buf, 1, 1, outf);
    }
  }

  fclose(outf);

  return 0;
}
/*
 * Compilation:
 * gcc -O3 -std=c99 -m cmandelbrot.c
 */

cmandelbrot.lisp

(declaim (optimize (speed 3) (safety 0) (debug 0)))


(deftype dcomplex () '(complex double-float))

(defconstant +nsteps+ 50000)
(defconstant +imgsize+ 1024)

(declaim (inline cabs2))

(defun cabs2 (z)
  (declare (type dcomplex z))
  (+ (expt (realpart z) 2)
     (expt (imagpart z) 2)))

(declaim (inline mandelbrot))
;; Mandelbrot algorithm with coloring from
;; http://fractals.nsu.ru/construction/mandelbrot_coloring.html
(defun mandelbrot (z0)
  (declare (type dcomplex z0))
  (let ((z z0)
        (zp (complex 1d0 0d0)))
    (declare (type dcomplex z zp))
    (dotimes (i +nsteps+ 0f0)
      (when (>= (cabs2 z) 4d0)
        (return (coerce (phase (/ zp z)) 'single-float)))
      (setf zp (* zp z))
      (setf z (+ z0 (* z z))))))


(defun calculate ()
  (let ((step (/ 4d0 (1- +imgsize+))))
    (with-open-file (img "test.pgm"
                         :element-type :default ; Bivalent stream
                         :if-exists :supersede
                         :direction :output)
      (format img "P5 ~A ~A 255~%" +imgsize+ (floor +imgsize+ 2))
        (dotimes (y (/ +imgsize+ 2))
          (dotimes (x +imgsize+)
            (let ((b (mandelbrot (complex (* step (- x
                                                     (floor +imgsize+ 2)
                                                     (floor +imgsize+ 8)))
                                          (* step y)))))
              (if (zerop b)
                  (write-byte 0 img)
                  (write-byte (floor (* 255 (+ 0.55 (* 0.45 (sin b)))))
                              img))))))))