setcommon-lispmandelbrot
Mandelbrot Set implementation in Common Lisp
I've been working on an implementation of the Mandelbrot Set in several different languages. I have a working implementation in C++, C#, Java, and Python, but the Common Lisp implementation has some bugs that I just can't figure out. It generates the set but somewhere in the pipeline the set gets distorted. I've tested and know with near certainty that the file I/O CLO isn't the problem - it's unlikely but possible, I've tested it pretty thoroughly.
Note that the intent of these implementations is to benchmark them against one another - so I'm trying to keep the code implementations as similar as possible so they are comparable.
The Mandelbrot set (here generated by the Python implementation):
"Mandelbrot Set (Generated by Python)"
But my Common Lisp program generates this:
"Common Lisp version's Distorted Mandelbrot Set"
The bug is identical in both Clisp and SBCL.
CODE:
Common Lisp:
(defun mandelbrot (real cplx num_iter)
(if (> (+ (* real real) (* cplx cplx)) 4)
1
(let ((tmpreal real) (tmpcplx cplx) (i 1))
(loop
(setq tmpcplx (+ (* (* tmpreal tmpcplx) 2) cplx))
(setq tmpreal (+ (- (* tmpreal tmpreal) (* tmpcplx tmpcplx))
real))
(setq i (+ i 1))
(cond
((> (+ (* tmpreal tmpreal)
(* tmpcplx tmpcplx)) 4) (return i))
((= i num_iter) (return 0)))))))
(defun floordiv (dend sor) (/ (- dend (mod dend sor)) sor))
(defclass xbm () (
(data :accessor data :initarg :data)
(dim :reader dim :initarg :dim)
(arrsize :reader arrsize :initarg :arrsize)))
(defmethod width ((self xbm)) (third (dim self)))
(defmethod height ((self xbm)) (second (dim self)))
(defun generate (width height)
(let ((dims (list 0 0 0)) (arrsize_tmp 0))
(setq dims (list 0 0 0))
(setf (second dims) height)
(setf (third dims) width)
(setf (first dims) (floordiv (third dims) 8))
(unless (= (mod width 8) 0) (setf (first dims) (+ (first dims) 1)))
(setq arrsize_tmp (* (first dims) (second dims)))
(make-instance 'xbm
:data (make-array arrsize_tmp :initial-element 0)
:dim dims
:arrsize arrsize_tmp)))
(defun writexbm (self f)
(with-open-file (stream f :direction :output :if-exists :supersede)
(let ((fout stream))
(format fout "#define mandelbrot_width ~d~&" (width self))
(format fout "#define mandelbrot_height ~d~&" (height self))
(format fout "#define mandelbrot_x_hot 1~&")
(format fout "#define mandelbrot_y_hot 1~&")
(format fout "static char mandelbrot_bits[] = {")
(let ((i 0))
(loop
(if (= (mod i 8) 0)
(format fout "~& ")
(format fout " "))
(format fout "0x~x," (svref (data self) i))
(unless (< (setf i (+ i 1)) (arrsize self))
(return t)))))))
(defmethod setpixel ((self xbm) (x integer) (y integer))
(if (and (< x (third (dim self))) (< y (second (dim self))))
(let ((val (+ (floordiv x 8) (* y (first (dim self))))))
(setf (svref (data self) val) (boole boole-ior (svref (data self) val) (ash 1 (mod x 8)))))))
(defmethod unsetpixel ((self xbm) (x integer) (y integer))
(if (and (< x (third (dim self))) (< y (second (dim self))))
(let ((val (+ (floordiv x 8) (* y (first (dim self))))))
(setf (svref (data self) val) (boole boole-xor (boole boole-ior
(svref (data self) val) (ash 1 (mod x 8))) (ash 1 (mod x 8)))))))
(defmethod draw_mandelbrot ((xbm xbm) (num_iter integer) (xmin number)
(xmax number) (ymin number) (ymax number))
(let ((img_width (width xbm)) (img_height (height xbm)) (xp 0))
(loop
(if (< xp img_width)
(let ((xcoord (+ (* (/ xp img_width) (- xmax xmin)) xmin)) (yp 0))
(loop
(if (< yp img_height)
(let (
(ycoord (+ (* (/ yp img_height) (- ymax ymin)) ymin)))
(let ((val (mandelbrot xcoord ycoord num_iter)))
(if (> val 0) (unsetpixel xbm xp yp) (setpixel xbm xp yp)))
(setq yp (+ yp 1)))
(return 0)))
(setq xp (+ xp 1)))
(return 0)))))
(defun main ()
(let ((maxiter 0) (xmin 0) (xmax 0) (ymin 0) (ymax 0) (file nil) (xsize 0) (ysize 0) (picture nil))
(format t "maxiter? ")
(setq maxiter (read))
(format t "xmin? ")
(setq xmin (read))
(format t "xmax? ")
(setq xmax (read))
(format t "ymin? ")
(setq ymin (read))
(format t "ymax? ")
(setq ymax (read))
(format t "file path: ")
(setq file (read-line))
(format t "picture width? ")
(setq xsize (read))
(format t "picture height? ")
(setq ysize (read))
(format t "~&")
(setq picture (generate xsize ysize))
(draw_mandelbrot picture maxiter xmin xmax ymin ymax)
(writexbm picture file)
(format t "File Written.")
0))
(main)
And the closest to it is Python:
from xbm import *
def mandelbrot(real_old,cplx_old,i):
real = float(real_old)
cplx = float(cplx_old)
if (real*real+cplx*cplx) > 4:
return 1
tmpreal = real
tmpcplx = cplx
for rep in range(1,i):
tmpb = tmpcplx
tmpcplx = tmpreal*tmpcplx*2
tmpreal = tmpreal*tmpreal - tmpb*tmpb
tmpcplx += cplx
tmpreal += real
tmpb = tmpcplx*tmpcplx + tmpreal*tmpreal
if tmpb > 4:
return rep+1
else:
return 0
def draw_mandelbrot(pic, num_iter, xmin, xmax, ymin, ymax):
img_width = pic.width()
img_height = pic.height()
for xp in range(img_width):
xcoord = (((float(xp)) / img_width) * (xmax - xmin)) + xmin
for yp in range(img_height):
ycoord = (((float(yp)) / img_height) * (ymax - ymin)) + ymin
val = mandelbrot(xcoord, ycoord, num_iter)
if (val):
pic.unsetpixel(xp, yp)
else:
pic.setpixel(xp, yp)
def main():
maxiter = int(raw_input("maxiter? "))
xmin = float(raw_input("xmin? "))
xmax = float(raw_input("xmax? "))
ymin = float(raw_input("ymin? "))
ymax = float(raw_input("ymax? "))
file = raw_input("file path: ")
xsize = int(raw_input("picture width? "))
ysize = int(raw_input("picture height? "))
print
picture = xbm(xsize, ysize)
draw_mandelbrot(picture, maxiter, xmin, xmax, ymin, ymax)
picture.writexbm(file)
print "File Written. "
return 0;
main()
[xbm.py]
from array import *
class xbm:
def __init__(self, width, height):
self.dim = [0, 0, 0]
self.dim[1] = height
self.dim[2] = width
self.dim[0] = self.dim[2] / 8
if width % 8 != 0:
self.dim[0] += 1
self.arrsize = self.dim[0] * self.dim[1]
self.data = array('B', (0 for x in range(self.arrsize)))
self.hex = ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'a', 'b', 'c', 'd', 'e', 'f']
def __nibbletochar__(self, a):
if a < 0 or a > 16:
return '0'
else:
return self.hex[a]
def setpixel(self, x, y):
if x < self.dim[2] and y < self.dim[1]:
self.data[(x / 8) + (y * self.dim[0])] |= 1 << (x % 8)
def unsetpixel(self, x, y):
if x < self.dim[2] and y < self.dim[1]:
self.data[(x / 8) + (y * self.dim[0])] |= 1 << (x % 8)
self.data[(x / 8) + (y * self.dim[0])] ^= 1 << (x % 8)
def width(self):
return self.dim[2]
def height(self):
return self.dim[1]
def writexbm(self, f):
fout = open(f, 'wt')
fout.write("#define mandelbrot_width ")
fout.write(str(self.dim[2]))
fout.write("\n#define mandelbrot_height ")
fout.write(str(self.dim[1]))
fout.write("\n#define mandelbrot_x_hot 1")
fout.write("\n#define mandelbrot_y_hot 1")
fout.write("\nstatic char mandelbrot_bits[] = {")
for i in range(self.arrsize):
if (i % 8 == 0): fout.write("\n\t")
else: fout.write(" ")
fout.write("0x")
fout.write(self.__nibbletochar__(((self.data[i] >> 4) & 0x0F)))
fout.write(self.__nibbletochar__((self.data[i] & 0x0F)))
fout.write(",")
fout.write("\n};\n")
fout.close();
I can post the C++, C#, or Java code as well need be.
Thanks!
EDIT: Thanks to Edmund's response i found the bug- Just something that slipped through the cracks on porting. Modified code:
(defun mandelbrot (real cplx num_iter)
(if (> (+ (* real real) (* cplx cplx)) 4)
1
(let ((tmpreal real) (tmpcplx cplx) (i 1) (tmpb cplx))
(loop
(setq tmpb tmpcplx)
(setq tmpcplx (+ (* (* tmpreal tmpcplx) 2) cplx))
(setq tmpreal (+ (- (* tmpreal tmpreal) (* tmpb tmpb))
real))
(setq i (+ i 1))
(cond
((> (+ (* tmpreal tmpreal)
(* tmpcplx tmpcplx)) 4) (return i))
((= i num_iter) (return 0)))))))
(defun floordiv (dend sor) (/ (- dend (mod dend sor)) sor))
(defclass xbm () (
(data :accessor data :initarg :data)
(dim :reader dim :initarg :dim)
(arrsize :reader arrsize :initarg :arrsize)))
(defun width (self) (third (dim self)))
(defun height (self) (second (dim self)))
(defun generate (width height)
(let ((dims (list 0 0 0)) (arrsize_tmp 0))
(setq dims (list 0 0 0))
(setf (second dims) height)
(setf (third dims) width)
(setf (first dims) (floordiv (third dims) 8))
(unless (= (mod width 8) 0) (setf (first dims) (+ (first dims) 1)))
(setq arrsize_tmp (* (first dims) (second dims)))
(make-instance 'xbm
:data (make-array arrsize_tmp :initial-element 0)
:dim dims
:arrsize arrsize_tmp)))
(defun writexbm (self f)
(with-open-file (stream f :direction :output :if-exists :supersede)
(let ((fout stream))
(format fout "#define mandelbrot_width ~d~&" (width self))
(format fout "#define mandelbrot_height ~d~&" (height self))
(format fout "#define mandelbrot_x_hot 1~&")
(format fout "#define mandelbrot_y_hot 1~&")
(format fout "static char mandelbrot_bits[] = {")
(let ((i 0))
(loop
(if (= (mod i 8) 0)
(format fout "~& ")
(format fout " "))
(format fout "0x~x," (svref (data self) i))
(unless (< (setf i (+ i 1)) (arrsize self))
(return t)))))))
(defun setpixel (self x y)
(if (and (< x (third (dim self))) (< y (second (dim self))))
(let ((val (+ (floordiv x 8) (* y (first (dim self))))))
(setf (svref (data self) val) (boole boole-ior (svref (data self) val) (ash 1 (mod x 8)))))))
(defun unsetpixel (self x y)
(if (and (< x (third (dim self))) (< y (second (dim self))))
(let ((val (+ (floordiv x 8) (* y (first (dim self))))))
(setf (svref (data self) val) (boole boole-xor (boole boole-ior
(svref (data self) val) (ash 1 (mod x 8))) (ash 1 (mod x 8)))))))
(defun draw_mandelbrot (xbm num_iter xmin xmax ymin ymax)
(let ((img_width (width xbm)) (img_height (height xbm)) (xp 0))
(loop
(if (< xp img_width)
(let ((xcoord (+ (* (/ xp img_width) (- xmax xmin)) xmin)) (yp 0))
(loop
(if (< yp img_height)
(let (
(ycoord (+ (* (/ yp img_height) (- ymax ymin)) ymin)))
(let ((val (mandelbrot xcoord ycoord num_iter)))
(if (> val 0) (unsetpixel xbm xp yp) (setpixel xbm xp yp)))
(setq yp (+ yp 1)))
(return 0)))
(setq xp (+ xp 1)))
(return 0)))))
(defun main ()
(let ((maxiter 0) (xmin 0) (xmax 0) (ymin 0) (ymax 0) (file nil) (xsize 0) (ysize 0) (picture nil))
(format t "maxiter? ")
(setq maxiter (read))
(format t "xmin? ")
(setq xmin (read))
(format t "xmax? ")
(setq xmax (read))
(format t "ymin? ")
(setq ymin (read))
(format t "ymax? ")
(setq ymax (read))
(format t "file path: ")
(setq file (read-line))
(format t "picture width? ")
(setq xsize (read))
(format t "picture height? ")
(setq ysize (read))
(format t "~&")
(setq picture (generate xsize ysize))
(draw_mandelbrot picture maxiter xmin xmax ymin ymax)
(writexbm picture file)
(format t "File Written.")
0))
(main)
Though the code isn't very LISP-ish (is that a word?) it works. Thanks to all who posted/commented/answered :)
Solution
I'm not sure this part is correct:
(setq tmpcplx (+ (* (* tmpreal tmpcplx) 2) cplx))
(setq tmpreal (+ (- (* tmpreal tmpreal) (* tmpcplx tmpcplx))
real))
Isn't tempcplx being overwritten with its new value on the first line, meaning that the second line is using the new value, not the original one?
In the Python version you're avoiding this problem by using tmpb:
tmpb = tmpcplx
tmpcplx = tmpreal*tmpcplx*2
tmpreal = tmpreal*tmpreal - tmpb*tmpb
tmpcplx += cplx
tmpreal += real
It seems to me the Lisp version should do something similar, i.e. store the original value of tmpcplx first, and use that store for the calculation of tmpreal:
(setq tmpb cplx)
(setq tmpcplx (+ (* (* tmpreal tmpcplx) 2) cplx))
(setq tmpreal (+ (- (* tmpreal tmpreal) (* tmpb tmpb))
real))
- Math.Sin() gives incorrect value
- How to run my python script when the sunOS is start booting
- Express-session: not resetting cookie expiration on each request
- Getting a stack overflow exception when normalizing a vector
- Edit default summary function in R gives error for multiple variables
- What was a For loop? Why isn't it needed in R?
- How to use download button in shiny and save results in various formats (csv, texte, pdf, spss...)?
- Why are there two assignment operators, `<-` and `->` in R?
- lm()$assign: what is it?
- How to get the value of list(...) in R and S functions
- Design matrix for MLM from library(lme4) with fixed and random effects
- how to generate elements not included in my sample
- Create a matrix with gradually changing values without a for loop
- Emacs ESS and S-plus ( S+ ) 8.1 compatability
- How to lag date-index in a time-series in R?
- Nonlinear regression in R / S
- Calling R from S-Plus?