5

I am trying to learn Lisp/Scheme and I tried implementing a very simple version of the mandelbrot set in it to get practice. The problem I ran into is that the code runs very, very slow. At first I thought it was because I was using recursion instead of imperative loops, but I tried re-writing more or less the same code (recursion included) in python (which doesn't even have tail-call optimisation), and it ran very smooth

So I am wondering if there is something obvious I am missing in my code and what I could do to make it run faster.

Here is the code snippet in Scheme (racket). I also did pretty much the same thing in SBCL and the speed was comparable

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

(I hope the code block is not too clustery, it was hard to strip it to something more simple)

Also, I know Scheme and Common Lisp have complex numbers built in but I wanted to test it using regular real numbers, I don't think this is the reason it runs so slow.

The parameter "i" of the brot function is the number of iterations, and the parameter "n" of brot-grid is also the number of iterations to use for each point. When I set it to more than like 10, the code takes forever to run, which doesn't seem normal. The increase in time taken also doesn't seem to be linear, for instance it only takes about a second on my machine with n = 10 but takes several minutes with n = 15 and doesn't even finish with n = 20

So, what is it that makes this code run so slow?

Thanks in advance

Etienne Boulais
  • 83
  • 6
  • Is `zr` and `zc` supposed to be very big? I paused it at one point and `zr` had over 4000 digits. Since Scheme has bigintegers it will not complain about the size until all the program memory is consumed. – Sylwester Aug 10 '15 at 17:07
  • 1
    'real numbers'? Floats? If you want to compute with real/float numbers, then you should make sure you actually use them and that your operations also use them. I see lots of integer and rational operations, which can be potentially slow, for example when using bignums or big rationals. Just trace or step the functions and you see what numbers the functions use. – Rainer Joswig Aug 10 '15 at 17:44
  • More than 10 iterations? Pshaw, imagine 1000! Even if you can implement *really snappy* ways to iterate, producing an Mset will be slow. If you just code it in basic ways, it will be *very* slow. It is computationally intensive, so maybe you need a better choice of language. And, you don't need complex number functions, you just need very quick ways to manipulate either fixed-point integers, or FPU real numbers. – Weather Vane Aug 10 '15 at 17:47
  • Thank you, the problem seems to be indeed with the representation of the numbers. Changing it to stop iterating when the module of Z is larger than 2 already makes it much faster. Making sure to use the right numbers representation should take care of the rest. I guess I'm not use to the way scheme seems to magically change the representation of the number depending on the operation you do on them – Etienne Boulais Aug 10 '15 at 18:31
  • I don't read Lisp, but you will always be quicker checking the magnitude's square against 4, than checking its square root against 2. – Weather Vane Aug 10 '15 at 19:01
  • 1
    So basically just changing the constant `2` in `linspace` to `2.0` makes the whole thing finish in a tenth of a second. @EtienneBoulais It's a bug there. It just sums the squares without square root so indeed the value should be `4` – Sylwester Aug 10 '15 at 19:04
  • @Sylwester Yeah I just did that before reading your answer, and indeed, changing the 2 for 2.0 in linspace makes it so the whole program runs using inexact floats rather than using ints and representing the result of divisions as exact fractions. Doing this, the program now runs at the same speed it would on another language that doesn't have exact numbers representation. Thank you all for your help! – Etienne Boulais Aug 10 '15 at 19:16
  • It's actually @RainerJoswig who solved this. I didn't spot the `/` amoungst the 4000 digits :-/ – Sylwester Aug 10 '15 at 19:19

2 Answers2

4

Looking at your code, I think you're testing it using rational numbers. This means pretty accurate arithmetics, with the drawback being atht you quickly end up using rationals with huge bignums as both numerator and denominator.

One way to ensure you're using floats (I'd suggest double-floats) is to have an intermediate function that converts all inputs to doubles, to make it easier to just type (say) 0 instead of 0d0.

Once you've established that using doubles makes it faster, you can start sprinkling type declarations throughout, to make it possible for the compiler to generate better code for you.

Vatine
  • 20,782
  • 4
  • 54
  • 70
  • 1
    you could limit your rationals to any precision, like e.g. 100 digits after the decimal point. that would be slower than doubles, but significantly faster than precise rationals. – Will Ness Aug 11 '15 at 10:53
  • I don't think Common Lisp has anything but precise rationals (unless as an implementation extension), so if you want non-precise ones, you'd probably have to implement them yourself. – Vatine Aug 11 '15 at 19:27
  • 1
    I meant, simply multiply by (10^n), truncate down to nearest integer and form a new ratio with 10^n denominator. that's not the Right Thing, but it'll do, and if not, increase `n` some more. :) – Will Ness Aug 11 '15 at 19:30
  • @WillNess That'd work, but at that point, it's probably easier to simply use double-floats (and it'll probably be more precise close to 0.0d0 anyway...). :) – Vatine Aug 11 '15 at 19:43
  • well, it should be easy to extract the mantissa in this scheme just as the doubles do it. So, 100 (or 1000) *meaningful digits*, then. :) that's a real coarse "solution" of course (because it has no control over the accumulating error). there are true adaptive precision schemes I'm sure. – Will Ness Aug 11 '15 at 20:08
  • An easy way to get floating point calculations is to write all number constants with decimal dot `0.0` or even `0.` – soegaard Aug 12 '15 at 14:19
2

Here is a Common Lisp variant:

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

Notice the use of double float constants. Also iterating both over double floats and integers.

Running it in SBCL unoptimized, but compiled to native code:

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

Optimizing the code then would mean:

  • higher compiler optimize setting
  • possibly adding some type declarations
  • getting rid of the float consing
Rainer Joswig
  • 136,269
  • 10
  • 221
  • 346