Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!uunet!mcsun!ukc!edcastle!aipna!cstr!tim
From: tim@cstr.ed.ac.uk (Tim Bradshaw)
Newsgroups: comp.lang.lisp
Subject: Floating point in CL?
Message-ID: <TIM.90Sep27140541@watt.cstr.ed.ac.uk>
Date: 27 Sep 90 14:05:41 GMT
Sender: news@aipna.ed.ac.uk
Distribution: comp
Organization: CSTR, University of Edinburgh
Lines: 124

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.

Here is what I used to test this.  The Lisp is written deliberately to
look as close as I could make it like the C that follows, so it's a
bit nasty:

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

(defun make-mandel-10 ()
  ;; Brute-force Mandelbrot thing.  No complex numbers!
  ;; This version has hard-wired assumptions, and looks like the C.
  (declare (optimize (speed 3)
		     (safety 0)))
 (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 float x y))
      (dotimes (xindex 10)
	(dotimes (yindex 10)
	  (setf x (the float (to-unity xindex 10))
		y (the float (to-unity yindex 10)))
	  (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 (the float (+ x (* x x) (- (* y y))))
		  y (the float (+ y (* 2.0 x y))))
	    (if (or (>= (abs x) 2.0)
		    (>= (abs y) 2.0))
		(progn 
		  (setf (aref image xindex yindex) 1)
		  (return nil))))))
      t))


And the C:

main ()
{
    int image[10][10];
    float x, y;
    int xindex, yindex, i;

    for (xindex = 0; xindex < 10; xindex++)
	for (yindex = 0; yindex < 10; yindex++)
	    image[xindex][yindex] = 0;

    for (xindex = 0; xindex < 10; xindex++) {
	for (yindex = 0; yindex < 10; yindex++) {
	    x = ((float) xindex)/10.0 - 0.5;
	    y = ((float) yindex)/10.0 - 0.5;
	    
	    for (i = 0; i < 100; i++) {
		x = x + x * x - y * y;
		y = y + 2.0 * x * y;

		if (x >= 2.0 || x <= -2.0 || y >= 2.0 || y <= -2.0) {
		    image[xindex][yindex] = 1;
		    break;
		}
	    }
	}
    }
}
		    
On A Sun4/65, with Sun cc, no optimization, and Allegro CL "3.1.13.1
[Sun4] (0/0/0)"  I get:

    kahlo* time ./mandel 
            0.1 real	0.0 user	0.0 sys

and

    <cl> (time (make-mandel-10))
    cpu time (non-gc) 1683 msec user, 0 msec system
    cpu time (gc)     250 msec user, 0 msec system
    cpu time (total)  1933 msec user, 0 msec system
    real time  2610 msec

Well I should use a bigger test so the C timings show up, but this
doesn't look good.

So the question: are there general guidelines on writing
floating-point code which will perform well in Lisp, given a good
compiler?  I'm pretty naive about numerical programming, but I tried
the following, as can be seen in the above code.

	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.

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?  Am I missing something
else?  Is Allegro just no good at floating point?  If so are there
good floating-pouint Lisp systems for Sun4s?

Tim Bradshaw.  Internet: tim%ed.cstr@nsfnet-relay.ac.uk
UUCP: ...!uunet!mcvax!ukc!cstr!tim  JANET: tim@uk.ac.ed.cstr
"Quis custodiet ipsos custodes?"