Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!sun-barr!apple!bionet!agate!saab.stanford.edu!neon!lucid.com!jwz From: jwz@lucid.com (Jamie Zawinski) Newsgroups: comp.lang.lisp Subject: Re: Floating point in CL? Message-ID: Date: 28 Sep 90 20:22:18 GMT References: Sender: jwz@lucid.com Organization: Lucid, Inc., Menlo Park, CA Lines: 153 In-reply-to: tim@cstr.ed.ac.uk's message of 27 Sep 90 14:05:41 GMT In article Tim Bradshaw writes: > > I wrote a trivial function (loop from 1 to n, incrementing a float); > and sure enough, Allegro comes quite close to C. Good. > Alas, when I write more complex code, it becomes appallingly slow. [...] > (defmacro to-unity (var maxvar) > ;; Convert var to be from -1 to 1 > ;; ?? Typing. Is coerce the right thing? Don't want a rational here. > `(- (/ (coerce ,var 'float) > (coerce ,maxvar 'float)) > 0.5)) COERCE is almost never what you want, certainly not in tight code; what you should say is something like (defmacro to-unity (var maxvar) "Convert var to be from -1 to 1" ;; ?? Typing. Is coerce the right thing? Don't want a rational here. `(- (/ (float ,var) (float ,maxvar)) 0.5)) The compiler will probably be smart enough to realize that the result of FLOAT is always a float, and will emit code for "/" that does a fp-divide, instead of code that does a type-dispatch on the arguments to decide what kind of division should take place. (I suppose it's possible for the compiler to make similar assumptions about the case where the second argument to COERCE is a constant, but that would be effectively the same as the compiler rewriting (COERCE FOO 'FLOAT) to be (FLOAT FOO).) > I assume it's crucially important that floats should be represented > directly rather than as pointers? How do you drop unsubtle enough > hints to the compiler that it listens to this? Things of type short-float will be immediate objects; the other float types will be represented as multi-word objects to which you have a pointer. Math on short-floats will probably be faster, but with a loss of precision. (defmacro to-unity (var maxvar) `(- (/ (float ,var 1.0s0) (float ,maxvar 1.0s0)) 0.5s0)) The compiler should realize that the result of FLOAT will always be a short-float, and will emit faster code for "/" and "-". But you may need to add calls to THE if your compiler isn't clever enough; your mileage may vary. > So the question: are there general guidelines on writing > floating-point code which will perform well in Lisp, [...] > Declaring types > Compiling with SPEED=3, SAFETY=0 > Making sure that results of intermediate parts of the > calculation are going to be FLOATs -- i.e. ensuring that / > doesn't get a chance to produce RATIONALs. > Using THE to allow the compiler to make assumptions about > types that it might not otherwise. These are the big ones. You probably shouldn't have to use THE very much, unless your compiler is failing to notice something. Another thing to consider is avoiding the overhead of a function call. When writing code that will turn into relatively few machine instructions (like numerical code often does) the overhead of building a new stack frame, saving registers, etc can be excessive. You can do this with macros, but I think it's cleaner to use functions declared or proclaimed INLINE when possible. This can make debugging easier, too - just turn off the inline declarations and you've got more info on the stack. > (defun make-mandel-10 () [...] Here is the inner loop of a Mandelbrot generator I implemented a while back for TI Explorer Lisp Machines. On a Lispm, type-checking isn't as important, but it does help. This function takes a coordinate in mandelbrot space, and some parameters, and returns the number of iterations before stability. It does the fp math in the precision of its arguments - if you pass in short floats, only short float arithmetic will be used; likewise for long-floats, etc. I structured it this way because the floating-point precision to be used is a runtime option of the program. But since this function it is proclaimed inline, if the compiler realizes that the arguments that its caller is passing in are all short-floats, it will emit short-float-only code for the divisions, subtractions, etc. (proclaim '(inline mandelbrot-one-point)) (defun mandelbrot-one-point (x y w h range depth real-origin imaginary-origin drop-out-distance &optional seed-x seed-y) "Computes the Mandelbrot mapping of a pixel XY. Returns a fixnum, the number of iterations to stability, or NIL if DROP-OUT-DISTANCE is reached. X,Y: The pixel we are computing (fixnums). W, H: The pixel size of the target area (fixnums). RANGE: The size of the source area in m-space (float). DEPTH: Maximum number of iterations (fixnum). REAL, IMAG: The origin of the source area in m-space (floats). DROP-OUT: float. SEED-X,Y: The image seed; if these are NIL, then the seed is the initial point (correct for Mandelbrot; for Julia, these are constant). All floating-point computations take place in the precision of the numbers passed in. So, if you want to use only short-float math, then all floating point values supplied should be short-floats. X,Y are really the complex number (REAL+X) + (IMAG+Y)i." ; (declare (values . iterations-or-nil)) (declare (optimize (speed 3) (safety 0)) (fixnum x y w h depth drop-out-distance) (float range real-origin imaginary-origin)) ;; ;; z0 = x + yi ;; z1 = z0^2 + mu ;; For Julia, mu is a constant; for Mandelbrot, mu is z0. ;; (let* ((s (max w h)) ;; real-z: real-origin + xr/s with same precision as r. ;; imaginary-z: imag-origin + yr/s with same precision as r. (real-z (+ real-origin (float (/ (* x range) s) range))) (imaginary-z (- imaginary-origin (float (/ (* y range) s) range))) (real-mu (or seed-x real-z)) (imaginary-mu (or seed-y imaginary-z)) (value nil)) (declare (float real-mu imaginary-mu real-z imaginary-z) (fixnum s) (type (or null fixnum) value)) (dotimes (iteration depth) (declare (fixnum iteration)) (if (>= (+ (* real-z real-z) (* imaginary-z imaginary-z)) drop-out-distance) (return) (setq value iteration)) (let ((new-real-z (- (- (* real-z real-z) (* imaginary-z imaginary-z)) real-mu))) ; (real-z^2 - imag-z^2) - real-mu (declare (float new-real-z)) ;; imaginary-z = (2 * real-z * imag-z) - imag-mu (setq imaginary-z (- (* (* real-z 2) imaginary-z) imaginary-mu)) (setq real-z new-real-z) )) value)) If you want to look at the rest of the code, you can anonymous FTP it from spice.cs.cmu.edu (128.32.150.27) in /usr/jwz/public/mandelbrot.lisp. (There's lots of other lisp code in that directory too - see _readme.text.) > "Quis custodiet ipsos custodes?" Hurm. -- Jamie