Implémentation de Mandelbrot Set en Common Lisp

Question

J'ai travaillé sur une implémentation du Mandelbrot Set en plusieurs langues. J'ai une implémentation fonctionnelle en C++, C#, Java et Python, mais l'implémentation de Common Lisp a quelques bugs que je n'arrive pas à comprendre. Il génère l'ensemble, mais quelque part dans le pipeline, l'ensemble est déformé. J'ai testé et je sais avec près de certitude que le fichier I/O CLO n'est pas le problème - il est peu probable mais possible, je l'ai testé assez à fond. Notez que le but de ces implémentations est de les comparer les unes aux autres. J'essaie donc de garder les implémentations de code aussi similaires que possible pour qu'elles soient comparables.

L'ensemble de Mandelbrot (ici généré par la mise en œuvre de Python):

http://www.freeimagehosting.net/uploads/65cb71a873.png http://www.freeimagehosting.net/uploads/65cb71a873.png "Mandelbrot (généré par Python)"

Mais mon programme Common Lisp génère ceci:

http://www.freeimagehosting.net/uploads/50bf29bcc9.png http://www.freeimagehosting.net/uploads/50bf29bcc9.png « Common Lisp l'ensemble de Mandelbrot déformé de la version "

Le bug est identique dans Clisp et 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) 

Et le plus proche de c'est 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(); 

Je peux poster le C++, C# ou code Java et besoin d'être.

Merci!

EDIT: Grâce à la réponse d'Edmund, j'ai trouvé le bug - Juste quelque chose qui a glissé à travers les fissures lors du portage. Code modifié:

(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) 

Bien que le code ne soit pas très LISP-ish (est-ce un mot?) Cela fonctionne. Merci à tous ceux qui a posté/a commenté/répondu :)

Answer 1

Je ne suis pas sûr que cette partie est correcte:

 (setq tmpcplx (+ (* (* tmpreal tmpcplx) 2) cplx)) 
     (setq tmpreal (+ (- (* tmpreal tmpreal) (* tmpcplx tmpcplx)) 
      real)) 

Isn » t tempcplx est écrasé avec sa nouvelle valeur sur la première ligne, ce qui signifie que la deuxième ligne utilise la nouvelle valeur, pas l'originale?

Dans la version Python vous éviter ce problème en utilisant TMPB:

tmpb = tmpcplx 
    tmpcplx = tmpreal*tmpcplx*2 
    tmpreal = tmpreal*tmpreal - tmpb*tmpb 
    tmpcplx += cplx 
    tmpreal += real 

Il me semble à la version Lisp devrait faire quelque chose de similaire, à savoir conserver la valeur d'origine de tmpcplx d'abord, et utiliser ce magasin pour le calcul de tmpreal:

 (setq tmpb cplx) 
     (setq tmpcplx (+ (* (* tmpreal tmpcplx) 2) cplx)) 
     (setq tmpreal (+ (- (* tmpreal tmpreal) (* tmpb tmpb)) 
      real)) 

Answer 2

Quelques remarques sur votre code:

• Mandelbrot: déclarations font défaut, les places sont calculées deux fois dans la boucle

• Mandelbrot: dans le calcul pour TMPREAL vous utilisez la nouvelle valeur de TMPCLX, et non l'ancien

• vous ne voulez pas utiliser des méthodes pour définir des pixels. LENT.

• FLOORDIV est l'un des FLOOR ou de tronquer (selon ce que vous voulez) en Common Lisp, voir (SOL 10 3)

• déclarations de type d'utilisation

• à writexbm n'appelle pas à plusieurs reprises DATA et ARRSIZE

• setPixel, unsetpixel semble très cher, encore une fois déréférencement à plusieurs reprises la structure

• tirage-Mandelbrot a beaucoup de r epeated calculs qui peuvent être faites une fois

• Common Lisp a des tableaux 2d qui visent à simplifier le code

• Common Lisp a des nombres complexes, ce qui simplifie également le code

• un nom de variable 'auto' ne fait aucune sens dans Lisp commun. Nomme-le à ce que c'est.

Généralement, le code est plein de déchets. Cela n'a pas de sens d'évaluer votre code, car il est écrit dans un style que, espérons-le, personne n'utilise dans Common Lisp. Common Lisp a été conçu avec l'expérience de grands logiciels mathématiques comme Macsyma et permet d'écrire du code mathématique de façon directe (pas d'objets, juste des fonctions sur les nombres, les tableaux, ...). Les meilleurs compilateurs peuvent tirer parti des types primitifs, des opérations primitives et des déclarations de types. Ainsi, le style est différent de ce que l'on pourrait écrire en Python (qui est généralement soit un Python orienté objet ou un appel à un code C), soit Ruby. En code numérique lourd, ce n'est généralement pas une bonne idée d'avoir une répartition dynamique comme avec CLOS. La définition de pixels dans des bitmaps via des appels CLOS dans une boucle LOOP est vraiment quelque chose que l'on veut éviter (sauf si vous savez comment l'optimiser).

Les meilleurs compilateurs Lisp compileront les fonctions numériques pour diriger le code machine. Pendant la compilation, ils donnent des indications sur les opérations qui sont génériques et qui ne peuvent pas être optimisées (jusqu'à ce que le développeur ajoute plus d'informations de type). Le développeur peut également 'DÉSASSEMBLER' les fonctions et vérifier le code générique ou faire des appels de fonction inutiles. «Dans consing de code numérique. « TIME » donne des informations d'exécution et informe également le développeur de la quantité de mémoire » CONSED flotte 'est un problème de performance habituelle

Donc, pour résumer:

• Si vous écrivez du code et pensez qu'il fait la même chose dans différentes langues, lorsque le code est similaire ou a une structure similaire, cela peut ne pas être le cas - sauf si vous connaissez vraiment les deux langages et les deux implémentations. Si vous écrivez du code dans une langue et que vous le portez dans un style similaire à celui d'une autre langue, vous risquez de manquer un ole culture existante pour écrire des solutions à ces types de problèmes d'une manière différente. Par exemple, on peut écrire du code en C++ dans un style orienté objet et le porter de la même façon que FORTRAN. Mais personne n'écrit un tel code en FORTRAN. Écrit dans le style FORTRAN, il se traduira généralement par un code plus rapide - d'autant plus que les compilateurs sont fortement optimisés pour le code FORTRAN idiomatique.

• "à Rome, parler comme les Romains"

Exemple:

dans setPixel il y a un appel à (premier (auto) faible). Pourquoi ne pas faire de cette valeur un slot dans la structure en premier lieu, au lieu de faire un accès à la liste tout le temps? Mais alors la valeur est constante pendant le calcul. La structure est toujours transmise et la valeur est récupérée tout le temps. Pourquoi ne pas simplement obtenir la valeur en dehors de la boucle principale et la passer directement? Au lieu d'en faire plusieurs calculs?

Pour vous donner une idée de la façon dont le code peut être écrit (avec des déclarations de type, des boucles, des nombres complexes, ...), voici une version légèrement différente du calcul de mandelbrot.

L'algorithme de base:

(defvar *num-x-cells* 1024) 
(defvar *num-y-cells* 1024) 
(defvar *depth* 60) 


(defun m (&key (left -1.5) (top -1.0) (right 0.5) (bottom 1.0) (depth *depth*)) 
    (declare (optimize (speed 3) (safety 0) (debug 0) (space 0))) 
    (loop with delta-x-cell float = (/ (- right left) *num-x-cells*) 
     and delta-y-cell float = (/ (- bottom top) *num-y-cells*) 
     and field = (make-array (list *num-x-cells* *num-y-cells*)) 
     for ix fixnum below *num-x-cells* 
     for x float = (+ (* (float ix) delta-x-cell) left) 
     do (loop for iy fixnum below *num-y-cells* 
       for y = (+ (* (float iy) delta-y-cell) top) 
       do (loop for i fixnum below depth 
          for z of-type complex = (complex x y) 
          then (+ (complex x y) (* z z)) 
          for exit = (> (+ (* (realpart z) (realpart z)) 
              (* (imagpart z) (imagpart z))) 
             4) 
          finally (setf (aref field ix iy) i) 
          until exit)) 
     finally (return field))) 

dessus fonction retourne un tableau 2D de nombres.

écriture d'un fichier XBM:

(defun writexbm (array pathname &key (black *depth*)) 
    (declare (fixnum black) 
      (optimize (speed 3) (safety 2) (debug 0) (space 0))) 
    (with-open-file (stream pathname :direction :output :if-exists :supersede) 
    (format stream "#define mandelbrot_width ~d~&" (array-dimension array 0)) 
    (format stream "#define mandelbrot_height ~d~&" (array-dimension array 1)) 
    (format stream "#define mandelbrot_x_hot 1~&") 
    (format stream "#define mandelbrot_y_hot 1~&") 
    (format stream "static char mandelbrot_bits[] = {") 
    (loop for j fixnum below (array-dimension array 1) do 
      (loop for i fixnum below (truncate (array-dimension array 0) 8) 
       for m fixnum = 0 then (mod (1+ m) 8) do 
       (when (zerop m) (terpri stream)) 
       (format stream "0x~2,'0x, " 
         (let ((v 0)) 
          (declare (fixnum v)) 
          (dotimes (k 8 v) 
          (declare (fixnum k)) 
          (setf v (logxor (ash (if (= (aref array 
                   (+ (* i 8) k) j) 
                 black) 
                1 0) 
               k) 
              v))))))) 
    (format stream "~&}~&"))) 

Au-dessus de la fonction prend un tableau et un chemin et écrit le tableau en tant que fichier XBM. Un numéro « noir » sera « noir » et les autres chiffres sont « blancs »

Appel

(writexbm (m) "/tmp/m.xbm")