Scheme中的Mandelbrot集合实现非常缓慢。

Question

我正在尝试学习Lisp/Scheme，并尝试在其中实现一个非常简单的曼德博集，以便练习。我遇到的问题是代码运行非常慢。起初，我以为是因为我使用了递归而不是命令式循环，但我尝试在Python中重新编写了几乎相同的代码（包括递归），它运行得非常顺畅。

所以我想知道我的代码中是否有明显的缺陷，以及如何使其运行更快。

以下是Scheme（Racket）中的代码片段。我也在SBCL中做了几乎相同的事情，速度是可比较的。

#lang racket

(define-syntax dotimes 
   (syntax-rules () 
     ((_ (var n res) . body) 
      (do ((limit n) 
           (var 0 (+ var 1))) 
          ((>= var limit) res) 
        . body)) 
     ((_ (var n) . body) 
      (do ((limit n) 
           (var 0 (+ var 1))) 
          ((>= var limit)) 
        . body))))

(define (print-brot zr zc)
  (if (< (+ (* zr zr) (* zc zc)) 2)
      (display "@")
      (display ".")))

(define (brot zr zc cr cc i)
  (if (= i 0)
      (print-brot zr zc)
      (let ((z2r (- (* zr zr) (* zc zc))) (z2c (* 2 zr zc)))
        (brot (+ z2r cr) (+ z2c cc) cr cc (- i 1)))))

(define (linspace i w)
  (/ (- i (/ w 2)) (/ w 4)))

(define (brot-grid w h n)
  (dotimes (i w)
           (dotimes (j h)
                    (let ((x (linspace i w)) (y (linspace j h)))
                      (brot 0 0 x y n)))
           (newline)))

(brot-grid 40 80 20)

我希望代码块不要太杂乱，将其简化是很困难的。

此外，我知道Scheme和Common Lisp内置了复数，但我想使用普通实数进行测试，我不认为这是运行缓慢的原因。

brot函数的参数“i”是迭代次数，“brot-grid”的参数“n”也是每个点要使用的迭代次数。当我将它设置为超过10时，代码运行时间变得异常长，这似乎不正常。时间花费的增加也似乎不是线性的，例如，在我的机器上，n = 10时只需要大约一秒钟，但n = 15时需要几分钟，而n = 20时甚至无法完成。

那么，是什么让这段代码运行如此缓慢呢？

提前感谢您。

- Etienne Boulais

---

“zr”和“zc”应该很大吗？我曾经在某一点暂停它，而“zr”有超过4000个数字。由于Scheme具有大整数，直到程序内存全部消耗完毕，它才会对此大小进行抱怨。 - Sylwester

“实数”？浮点数？如果你想要使用实数/浮点数进行计算，那么你应该确保你实际上使用了它们，并且你的操作也使用了它们。我看到很多整数和有理数的操作，这可能会潜在地变得很慢，例如当使用大整数或大有理数时。只需跟踪或逐步执行函数，你就可以看到函数使用的数字。 - Rainer Joswig

谢谢，问题似乎确实与数字的表示有关。将其更改为在Z的模大于2时停止迭代，可以使其运行速度更快。确保使用正确的数字表示应该解决剩下的问题。我想我还不习惯scheme似乎会根据对它们执行的操作神奇地改变数字的表示方式。 - Etienne Boulais

所以基本上只需将linspace中的常数2更改为2.0，整个过程就可以在十分之一秒内完成。@EtienneBoulais 那是一个错误。它只是对平方进行求和而没有开根号，因此实际值应该是4。 - Sylwester

实际上是@RainerJoswig解决了这个问题。我没有在4000个数字中发现“/” :-/ - Sylwester

Answer 1

看了你的代码，我认为你正在使用有理数进行测试。这意味着使用相当准确的算术，缺点是很快就会使用具有巨大bignums的有理数作为分子和分母。

确保您使用浮点数（我建议双倍浮点数）的一种方法是使用中间函数将所有输入值转换为双倍浮点数，以便更容易地键入（比如说）0而不是0d0。

一旦确定使用双精度浮点数可以使其更快，您就可以开始在整个过程中添加类型声明，以便编译器为您生成更好的代码。

- Vatine

---

你可以将有理数限制在任何精度内，例如小数点后100位。这比双精度慢，但比精确有理数快得多。 - Will Ness

我认为Common Lisp除了精确的有理数（除非作为实现扩展），没有其他东西，所以如果你想要非精确的有理数，你可能需要自己实现它们。 - Vatine

我的意思是，简单地乘以（10^n），向下截取到最近的整数，并用10^n分母形成一个新的比率。这不是正确的方法，但它可以做到，如果不行，就再增加n。 :) - Will Ness

@WillNess 那样做是可行的，但在那个点上，使用双精度浮点数可能更容易（而且在接近0.0d0时可能更精确...）。 :) - Vatine

嗯，在这个方案中提取尾数应该很容易，就像双精度浮点数一样。所以，有100（或1000）个有意义的数字，然后。 :) 当然，这是一个真正粗糙的“解决方案”（因为它无法控制累积误差）。我相信还有真正的自适应精度方案。 - Will Ness

获取浮点数计算的简单方法是使用带有小数点“0.0”或甚至“0.”的所有数字常量进行编写。 - soegaard

Answer 2

这里是一种通用Lisp变体：

(defun print-brot (zr zc)
  (write-char (if (< (+ (* zr zr)
                        (* zc zc))
                     2.0d0)
                  #\@
                #\.)))

(defun brot (zr zc cr cc i)
  (loop repeat i
        for z2r = (- (* zr zr) (* zc zc))
        for z2c = (* 2.0d0 zr zc)
        until (or (> (abs zr) 1.0d50)
                  (> (abs zc) 1.0d50))
        do (setf zr (+ z2r cr)
                 zc (+ z2c cc)))
  (print-brot zr zc))

(defun linspace (i w)
  (/ (- i (/ w 2.0d0)) (/ w 4.0d0)))

(defun brot-grid (w h n)
  (terpri)
  (loop for i from 0.0d0 by 1.0d0
        repeat w do
        (loop for j from 0.0d0 by 1.0d0
              repeat h do
              (brot 0.0d0 0.0d0 (linspace i w) (linspace j h) n))
    (terpri)))

请注意双精度浮点常量的使用。同时迭代双精度浮点数和整数。

在SBCL中运行，未经优化但编译为本地代码：

*  (time (brot-grid 20 40 20))

........................................
....................@...................
....................@...................
....................@...................
...................@@@..................
.................@@@@@@@................
...................@@@..................
................@@@@@@@@@...............
..............@@@@@@@@@@@@@.............
............@@@@@@@@@@@@@@@@@...........
..............@@@@@@@@@@@@@.............
...............@@@@@@@@@@@..............
..................@...@.................
........................................
........................................
........................................
........................................
........................................
........................................
........................................
Evaluation took:
  0.003 seconds of real time
  0.002577 seconds of total run time (0.001763 user, 0.000814 system)
  100.00% CPU
  6,691,716 processor cycles
  2,064,384 bytes consed

优化代码意味着：

• 更高的编译器优化设置
• 可能添加一些类型声明
• 消除浮点数临时分配

- Rainer Joswig