set - Common Lisp でのマンデルブロー集合の実装

Question

私はいくつかの異なる言語でのマンデルブロー集合の実装に取り​​組んできました。私は C++、C#、Java、および Python で動作する実装を持っていますが、Common Lisp の実装には、私には理解できないバグがいくつかあります。セットを生成しますが、パイプラインのどこかでセットが歪んでしまいます。ファイル I/O CLO が問題ではないことをテストし、ほぼ確実に知っています。可能性は低いですが、かなり徹底的にテストしました。

これらの実装の目的は、相互にベンチマークすることであることに注意してください。そのため、コードの実装を可能な限り類似させて、比較できるようにしています。

マンデルブロー集合 (ここでは Python 実装によって生成されます):

「マンデルブロー集合（Pythonで生成）」

しかし、私の Common Lisp プログラムはこれを生成します:

「Common 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 コードも投稿できます。

ありがとう！

編集: Edmund の応答のおかげで、私はバグを見つけました-移植時に亀裂をすり抜けたものです。変更されたコード:

(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: 宣言がなく、ループ内で 2 乗が 2 回計算されます

• mandelbrot: TMPREAL の計算では、古い値ではなく、TMPCLX の新しい値を使用しています。

• METHODS を使用してピクセルを設定する必要はありません。スロー。

• FLOORDIV は、Common Lisp の FLOOR または TRUNCATE (目的に応じて) のいずれかです。(FLOOR 10 3) を参照してください。

• 型宣言を使用する

• writexbm では、DATA と ARRSIZE を繰り返し呼び出さないでください。

• setpixel、unsetpixel は非常にコストがかかるように見えます。これも繰り返し構造を逆参照しています。

• draw-mandelbrot には、一度に実行できる反復計算が多数あります。

• Common Lisp には、コードを単純化する 2 次元配列があります。

• Common Lisp には複素数があり、コードも単純化されます

• 変数名 'self' は、Common Lisp では意味がありません。それが何であるかに名前を付けます。

一般的に、コードは無駄に満ちています。コードは Common Lisp で誰も使用しないようなスタイルで書かれているため、コードのベンチマークを行うことにはほとんど意味がありません。Common Lisp は、Macsyma のような大規模な数学ソフトウェアの経験に基づいて設計されており、数学コードを簡単な方法で記述できます (オブジェクトを使用せず、数値や配列に対する関数だけを使用するなど)。優れたコンパイラは、プリミティブ型、プリミティブ操作、および型宣言を利用できます。したがって、スタイルは、Python (通常はオブジェクト指向の Python または C コードの呼び出し) または Ruby で作成するものとは異なります。重い数値コードでは、通常、CLOS のように動的ディスパッチを使用することはお勧めできません。タイトなループでCLOS呼び出しを介してビットマップにピクセルを設定することは、実際には避けたいことです(最適化する方法を知っていない限り)。

より優れた Lisp コンパイラは、数値関数をコンパイルして直接のマシン コードにします。コンパイル中に、どの操作が一般的で最適化できないかについてのヒントが示されます (開発者が型情報を追加するまで)。開発者は、関数を「DISASSEMBLE」して、一般的なコードや不要な関数呼び出しを行うコードをチェックすることもできます。「TIME」はランタイム情報を提供し、開発者に「consed」のメモリ量を通知します。数値コードでは、「floats」の consing は通常のパフォーマンスの問題です。

要約すると、次のようになります。

• コードを書いて、異なる言語で同じことをすると思っても、コードが似ていたり構造が似ていたりしても、実際にはそうではないかもしれません - 両方の言語と両方の言語の実装をよく知っている場合を除きます。

• ある言語でコードを書き、それを同様のスタイルで別の言語に移植すると、この種の問題に対するソリューションを別の方法で作成するという既存の文化全体を見逃す可能性があります。たとえば、C++ でオブジェクト指向スタイルのコードを記述し、それを FORTRAN と同様の方法で移植できます。しかし、FORTRAN でそのようなコードを書く人はいません。FORTRAN スタイルで書かれていると、通常、コードが高速になります。特に、コンパイラは慣用的な FORTRAN コード用に大幅に最適化されているためです。

• 「ローマにいるときはローマ人のように話せ」

例：

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

上記の関数は、数値の 2 次元配列を返します。

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 ファイルとして書き込みます。1 つの数字「黒」は「黒」になり、他の数字は「白」になります

電話

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

Answer 2

この部分が正しいかどうかはわかりません：

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

tempcplx は最初の行で新しい値で上書きされていませんか?つまり、2 行目は元の値ではなく、新しい値を使用していますか?

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