Реализация Мандельброта в Common Lisp

Question

Я работал над реализацией набора Мандельброта на нескольких разных языках. У меня есть рабочая реализация на C ++, C #, Java и Python, но в реализации Common Lisp есть некоторые ошибки, которые я просто не могу понять. Он генерирует набор, но где-то в конвейере набор искажается. Я проверял и знаю с почти уверенностью, что файловый ввод-вывод CLO не проблема - это маловероятно, но возможно, я довольно тщательно его протестировал.

Обратите внимание, что цель этих реализаций состоит в том, чтобы сравнить их друг с другом - поэтому я стараюсь, чтобы реализации кода были как можно более похожими, чтобы они были сопоставимыми.

Набор Мандельброта (здесь сгенерированный реализацией Python):

«Набор Мандельброта (сгенерированный Python)»

Но моя программа Common Lisp генерирует это:

"Искаженный набор Мандельброта для версии Common Lisp"

Ошибка одинакова как в Clisp, так и в SBCL.

КОД:

Обычный Лисп:

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

И ближайший к нему 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();

Я также могу опубликовать код C ++, C # или Java.

Спасибо!

РЕДАКТИРОВАТЬ: Благодаря ответу Эдмунда я нашел ошибку - просто что-то, что ускользнуло от трещин при портировании. Модифицированный код:

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

Хотя код не очень LISP-иш (это слово?), Он работает. Спасибо всем, кто писал / комментировал / отвечал:)

Answer 1

Несколько замечаний по поводу вашего кода:

• Мандельброт: отсутствуют объявления, квадраты вычисляются дважды в цикле

• mandelbrot: в вычислениях для TMPREAL вы используете новое значение TMPCLX, а не старое

• Вы не хотите использовать МЕТОДЫ для установки пикселей. SLOW.

• FLOORDIV - это один из FLOOR или TRUNCATE (в зависимости от того, что вы хотите) в Common Lisp, см. (FLOOR 10 3)

• использовать объявления типов

• в writexbm повторно не вызывают DATA и ARRSIZE

• setpixel, unsetpixel выглядит очень дорого, снова многократно разыменовывая структуру

• draw-mandelbrot имеет много повторных вычислений, которые можно выполнить один раз

• Common Lisp имеет 2d массивы, которые упрощают код

• Common Lisp имеет комплексные числа, которые также упрощают код

• имя переменной 'self' не имеет смысла в Common Lisp. Назовите это тем, чем оно является.

Обычно код полон отходов. Нет смысла сравнивать ваш код, поскольку он написан в стиле, который, надеюсь, никто не использует в Common Lisp. Common Lisp был разработан с опытом работы с большим математическим программным обеспечением, таким как Macsyma, и позволяет писать математический код простым способом (без объектов, только функции над числами, массивами, ...). Лучшие компиляторы могут использовать преимущества примитивных типов, примитивных операций и объявлений типов. Таким образом, стиль отличается от того, что можно написать в Python (который обычно является либо объектно-ориентированным Python, либо вызовами некоторого кода на C), либо Ruby. В тяжелом числовом коде, как правило, не очень хорошая идея динамической отправки, как в случае с CLOS. Установка пикселей в растровых изображениях с помощью вызовов CLOS в тесной LOOP - это действительно то, чего нужно избегать (если вы не знаете, как его оптимизировать).

Лучшие компиляторы Lisp будут компилировать числовые функции для прямого машинного кода. Во время компиляции они дают подсказки, какие операции являются общими и не могут быть оптимизированы (пока разработчик не добавит больше информации о типе). Разработчик также может «РАЗБРАСИТЬ» функции и проверять код, который является общим или выполняет ненужные вызовы функций. «TIME» предоставляет информацию о времени выполнения, а также информирует разработчика об объеме памяти, «ограниченном». В числовом коде использование «float» является обычной проблемой производительности.

Итак, подведем итог:

• если вы пишете код и думаете, что он делает то же самое на разных языках, когда код выглядит одинаково или имеет похожую структуру, это может быть не так - если вы действительно не знаете оба языка и обе языковые реализации.

• если вы пишете код на одном языке и переносите его в похожем стиле на другой язык, вы можете пропустить всю существующую культуру, чтобы по-разному решать проблемы такого рода. Например, можно написать код на C ++ в объектно-ориентированном стиле и перенести его аналогично FORTRAN. Но никто не пишет такой код на Фортране. Написанный в стиле FORTRAN, обычно приводит к более быстрому коду - тем более что компиляторы сильно оптимизированы для идиоматического кода FORTRAN.

• "в Риме говорите как римляне"

Пример: * 1 069 *

в SETPIXEL есть вызов (first (dim self)). Почему бы не сделать это значение главным элементом слота в структуре, вместо того, чтобы делать доступ к списку постоянно? Но тогда значение является постоянным во время вычисления. Тем не менее структура передается, и значение извлекается все время. Почему бы просто не получить значение вне основного цикла и передать его напрямую? Вместо того, чтобы делать несколько вычислений этого?

Чтобы дать вам представление о том, как может быть написан код (с объявлениями типов, циклами, комплексными числами, ...), вот немного другая версия вычисления Мандельброта.

Основной алгоритм:

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

Функция выше возвращает двумерный массив чисел.

Запись файла 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 "~&}~&")))

Вышеупомянутая функция берет массив и путь и записывает массив в виде файла XBM. Одно число «черное» будет «черным», а другое число «белым»

Вызов

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

Answer 2

Я не уверен, что эта часть верна:

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

Разве tempcplx не перезаписывается новым значением в первой строке, что означает, что вторая строка использует новое значение, а не исходное?

В версии Python вы избегаете этой проблемы, используя tmpb:

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

Мне кажется, что версия на Лиспе должна делать что-то похожее, т.е. сначала сохранять первоначальное значение tmpcplx и использовать это хранилище для вычисления tmpreal:

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