set - Common Lisp 中的 Mandelbrot Set 实现

Question

我一直致力于用几种不同的语言实现Mandelbrot 集。我在 C++、C#、Java 和 Python 中有一个可用的实现，但是 Common Lisp 实现有一些我无法弄清楚的错误。它生成集合，但在管道中的某个地方，集合被扭曲了。我已经测试并且几乎可以肯定地知道文件 I/O CLO 不是问题 - 这不太可能但可能，我已经对其进行了非常彻底的测试。

请注意，这些实现的目的是将它们相互进行基准测试 - 所以我试图使代码实现尽可能相似，以便它们具有可比性。

Mandelbrot 集（此处由 Python 实现生成）：

“Mandelbrot 集（由 Python 生成）”

但是我的 Common Lisp 程序会生成这个：

“普通 Lisp 版本的扭曲曼德布罗集”

该错误在 Clisp 和 SBCL 中是相同的。

代码：

普通 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)

最接近它的是 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：缺少声明，平方在循环中计算两次

• mandelbrot：在计算 TMPREAL 时，您使用的是 TMPCLX 的新值，而不是旧值

• 您不想使用 METHODS 设置像素。慢。

• FLOORDIV 是 Common Lisp 中的 FLOOR 或 TRUNCATE 之一（取决于您想要什么），请参阅 (FLOOR 10 3)

• 使用类型声明

• 在 writexbm 中不要重复调用 DATA 和 ARRSIZE

• setpixel, unsetpixel 看起来很贵，再次重复解引用结构

• draw-mandelbrot 有很多可以一次完成的重复计算

• Common Lisp 具有简化代码的二维数组

• Common Lisp 有复数，这也简化了代码

• 在 Common Lisp 中，变量名 'self' 没有意义。把它命名为它是什么。

通常代码充满了浪费。对代码进行基准测试没有什么意义，因为它是用一种希望没有人在 Common Lisp 中使用的风格编写的。Common Lisp 是根据 Macsyma 等大型数学软件的经验设计的，允许以直接的方式编写数学代码（没有对象，只是对数字、数组等的函数）。更好的编译器可以利用原始类型、原始操作和类型声明。因此，这种风格不同于用 Python（通常是面向对象的 Python 或调用某些 C 代码）或 Ruby 编写的风格。在繁重的数字代码中，像 CLOS 那样进行动态调度通常不是一个好主意。在紧密的循环中通过 CLOS 调用在位图中设置像素确实是人们想要避免的事情（除非您知道如何优化它）。

更好的 Lisp 编译器会将数字函数编译为直接机器代码。在编译期间，它们会提示哪些操作是通用的并且无法优化（直到开发人员添加更多类型信息）。开发人员还可以“分解”函数并检查通用代码或执行不必要的函数调用。“TIME”提供运行时信息，并告知开发人员“consed”的内存量。在数字代码中，“floats”的consing是一个常见的性能问题。

所以，总结一下：

• 如果您编写代码并认为它​​在不同语言中的作用相同，当代码看起来相似或具有相似结构时，情况可能并非如此——除非您真的了解两种语言和两种语言的实现。

• 如果你用一种语言编写代码并以类似的风格将其移植到另一种语言，你可能会错过整个现有的文化，以不同的方式为这类问题编写解决方案。例如，可以以面向对象的风格用 C++ 编写代码，并以类似于 FORTRAN 的方式将其移植。但是没有人在 FORTRAN 中编写这样的代码。以 FORTRAN 风格编写，通常会产生更快的代码 - 特别是因为编译器针对惯用的 FORTRAN 代码进行了高度优化。

• “在罗马时，像罗马人一样说话”

例子：

在 SETPIXEL 中有一个调用 (first (dim self))。为什么不首先将该值作为结构中的一个槽，而不是一直进行列表访问呢？但是在计算过程中该值是恒定的。仍然传递结构，并且一直检索值。为什么不直接在主循环之外获取值并直接传递呢？而不是对其进行多次计算？

为了让您了解如何编写代码（使用类型声明、循环、复数......），这里有一个稍微不同的 mandelbrot 计算版本。

核心算法：

(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

在我看来，Lisp 版本应该做类似的事情，即首先存储 tmpcplx 的原始值，然后使用该存储来计算 tmpreal：

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