Path: utzoo!attcan!uunet!fernwood!franz!jim
From: jim@Franz.COM (Jim Veitch)
Newsgroups: comp.lang.lisp
Subject: Re: Floating point in CL?
Message-ID: <JIM.90Oct2093959@crisp.Franz.COM>
Date: 2 Oct 90 16:39:59 GMT
References: <TIM.90Sep27140541@watt.cstr.ed.ac.uk>
Sender: news@Franz.COM
Distribution: comp
Organization: Franz Inc., Berkeley, CA
Lines: 81

I thought I'd reply to your posting to comp.lang.lisp:

> For the first time in my life, I need to write some floating-point
> code.  Since I hate C, I want to do it in Lisp.  However my initial
> experiments indicate that our Lisp is *immeasurably* slower than C.

> CLtL states at several places that a good Common Lisp compiler
> can generate extremely good floating-point code.  Since we have at
> least a reasonable compiler (Franz Allegro, Sun4), I had some hope
> that it would generate at least reasonable floating-point code.

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

1. Allegro CL supports 2 types of floats (like C) -- double-floats and
single-floats.  Unlike C, where 'float' has a precise meaning, Common
Lisp 'float' can be a composite type.  In your code, simply declared
'float', the compiler doesn't know which one you mean, and generates
generic arithmetic instead.  If you declare all one type and if the
compiler trusts the declarations (in Allegro this is only true if
speed > safety), everything is in-lined.

It is true that there is no standard for which float types are
supported and which compiler optimizations are done.  I guess each
compiler uses differing sets and the only real ways to see which are to
read the documentation, use the metering tools (in the add-on Composer
product with Franz), and ask the vendor directly (which I would
encourage any Franz user to do!).

2. Where did all the time go in your example?

(a) mostly in generic floating point arithmetic (+, -, *)

(b) much of the rest in calling 'abs' and 'coerce' (which I removed in the
    example below).  'abs' is an unnecessary function call which entails
    boxing up the float and  handing a pointer to 'abs'.  You might argue
    the compiler should be smart enough to recognize a special 'abs' for
    declared floating point arguments, but Allegro does not do that at
    the current time.

(c) call to 'coerce' (which is a an extremely complex Common Lisp function)
    is dangerous unless the compiler is very clever.  Ours isn't clever
    enough to do the "right" thing.  In general, it's wisest to use the
    simpler thing, which the compiler is more likely to know about.  So
    the call to 'coerce' is replaced by a call to 'float', which the
    Allegro CL compiler does know about.

3. I rewrote your example to look exactly like the C code and ran it with
the following time (it is one clock-tick or less on a Sun4/110): 16 msecs.

(defun make-mandel-10 ()
  ;; Brute-force Mandelbrot thing.  No complex numbers!
  ;; This version has hard-wired assumptions, and looks like the C.
  (let ((image (make-array '(10 10)
			   :element-type
			   'bit :initial-element 0))
	;; image components set to 0 if in, 1 if not: different
	;; than make-mandel!
	(x 0.0)
	(y 0.0))
    (declare (type single-float x y)  		; originally declared 'float'
	     (optimize (speed 3) (safety 1)))
    (dotimes (xindex 10)
      (dotimes (yindex 10)
	(setf x (- (/ (float xindex) 10.0) 0.5)  ; 'float' replaces 'coerce'
	      y (- (/ (float yindex) 10.0) 0.5))
	(dotimes (i 100)
	  ;; z = x + iy => z * z = x*x +2ixy - y*y =>
	  ;; Re(z^2) =X^2 -y^2, Im(z^2) = 2xy
	  (setq x (+ x (* x x) (- (* y y)))
		y (+ y (* 2.0 x y)))
	  (if (or (>= x 2.0)	; 'abs' test expanded into 4 cases
		  (<= x -2.0)   ; just like the C code
		  (>= y 2.0)
		  (<= y 2.0))
	      (progn
		(setf (aref image xindex yindex) 1)
		(return nil))))))
    t))