1.0.27.41: floating point implementation smoothing
authorChristophe Rhodes <csr21@cantab.net>
Fri, 24 Apr 2009 14:41:54 +0000 (14:41 +0000)
committerChristophe Rhodes <csr21@cantab.net>
Fri, 24 Apr 2009 14:41:54 +0000 (14:41 +0000)
To get floating point stuff exactly right, we should build a complete
IEEE float implementation to do calculations in for the cross-compiler.
Since that's not going to happen this millennium, instead try to be
careful when writing code that looks constant-foldable.  Some other
fixups on the way...

6 messages follow:

fix load-time tests in src/code/pred

It turns out that #c(1.1 0) is not portable: it's a REAL in clisp and a
COMPLEX in sbcl.

begin work on floats

Floats Are Hard.  The issue is that the host's float implementation,
even if it agrees with SBCL that SINGLE-FLOAT is IEEE single and
DOUBLE-FLOAT is IEEE double, may not match sbcl idiosyncracy for
idiosyncracy.  For example, clisp doesn't support denormals, so its
LEAST-FOOATIVE-QUUXLE-FLOAT constants are very different from sbcl's:
and sbcl's can't even be represented within the host.  Ugh.

Defining the print-related MIN-E constants is, however, easy enough.

comment (well, #!+long-float) out some floating point constants

The clauses in question were never taken absent #!+long-float anyway.

-0.0 is not portable: many lisps don't respect negative zeros

Use make-unportable-float instead, and hope that this doesn't matter
during cross-compilation...

host floating point differences

Not all lisps think (log 2d0 10d0) is the same.  Compute it accurately
and use LOAD-TIME-VALUE.

tentative attempt at smoothing over host floating point differences

Compute all the necessary constants as double-float bit patterns using
LOAD-TIME-VALUE.

src/code/irrat.lisp
src/code/pred.lisp
src/code/print.lisp
src/compiler/srctran.lisp
src/compiler/x86/float.lisp
version.lisp-expr

index 5e2c8d0..55dbb83 100644 (file)
              (values
               (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0)
               0))
-            ((let ((threshold #.(/ least-positive-double-float
-                                   double-float-epsilon))
+            ((let ((threshold
+                    ;; (/ least-positive-double-float double-float-epsilon)
+                    (load-time-value
+                     #!-long-float
+                     (sb!kernel:make-double-float #x1fffff #xfffffffe)
+                     #!+long-float
+                     (error "(/ least-positive-long-float long-float-epsilon)")))
                    (traps (ldb sb!vm::float-sticky-bits
                                (sb!vm:floating-point-modes))))
                 ;; Overflow raised or (underflow raised and rho <
   ;; influences the choices of these constants but doesn't say how to
   ;; choose them.  We'll just assume his choices matches our
   ;; implementation of log1p.
-  (let ((t0 #.(/ 1 (sqrt 2.0d0)))
+  (let ((t0 (load-time-value
+             #!-long-float
+             (sb!kernel:make-double-float #x3fe6a09e #x667f3bcd)
+             #!+long-float
+             (error "(/ (sqrt 2l0))")))
+        ;; KLUDGE: if repeatable fasls start failing under some weird
+        ;; xc host, this 1.2d0 might be a good place to examine: while
+        ;; it _should_ be the same in all vaguely-IEEE754 hosts, 1.2
+        ;; is not exactly representable, so something could go wrong.
         (t1 1.2d0)
         (t2 3d0)
-        (ln2 #.(log 2d0))
+        (ln2 (load-time-value
+              #!-long-float
+              (sb!kernel:make-double-float #x3fe62e42 #xfefa39ef)
+              #!+long-float
+              (error "(log 2l0)")))
         (x (float (realpart z) 1.0d0))
         (y (float (imagpart z) 1.0d0)))
     (multiple-value-bind (rho k)
       ;; space 0 to get maybe-inline functions inlined
       (declare (optimize (speed 3) (space 0)))
     (cond ((> (abs x)
-              ;; FIXME: this form is hideously broken wrt
-              ;; cross-compilation portability.  Much else in this
-              ;; file is too, of course, sometimes hidden by
-              ;; constant-folding, but this one in particular clearly
-              ;; depends on host and target
-              ;; MOST-POSITIVE-DOUBLE-FLOATs being equal.  -- CSR,
-              ;; 2003-04-20
-              #.(/ (+ (log 2.0d0)
-                      (log most-positive-double-float))
-                   4d0))
+              (load-time-value
+               #!-long-float
+               (sb!kernel:make-double-float #x406633ce #x8fb9f87e)
+               #!+long-float
+               (error "(/ (+ (log 2l0) (log most-positive-long-float)) 4l0)")))
            (coerce-to-complex-type (float-sign x)
                                    (float-sign y) z))
           (t
index ee5c9e8..1bd09c7 100644 (file)
@@ -361,8 +361,9 @@ length and have identical components. Other arrays must be EQ to be EQUAL."
 #!+sb-test
 (let ((test-cases `((0.0 ,(load-time-value (make-unportable-float :single-float-negative-zero)) t)
                     (0.0 1.0 nil)
-                    (#c(1 0) #c(1.0 0) t)
-                    (#c(1.1 0) #c(11/10 0) nil) ; due to roundoff error
+                    (#c(1 0) #c(1.0 0.0) t)
+                    (#c(0 1) #c(0.0 1.0) t)
+                    (#c(1.1 0.0) #c(11/10 0) nil) ; due to roundoff error
                     ("Hello" "hello" t)
                     ("Hello" #(#\h #\E #\l #\l #\o) t)
                     ("Hello" "goodbye" nil))))
index 84f557a..9e302d9 100644 (file)
 ;;; possible extension for the enthusiastic: printing floats in bases
 ;;; other than base 10.
 (defconstant single-float-min-e
-  (nth-value 1 (decode-float least-positive-single-float)))
+  (- 2 sb!vm:single-float-bias sb!vm:single-float-digits))
 (defconstant double-float-min-e
-  (nth-value 1 (decode-float least-positive-double-float)))
+  (- 2 sb!vm:double-float-bias sb!vm:double-float-digits))
 #!+long-float
 (defconstant long-float-min-e
   (nth-value 1 (decode-float least-positive-long-float)))
           (values (float 0.0e0 original-x) 1)
           (let* ((ex (locally (declare (optimize (safety 0)))
                        (the fixnum
-                         (round (* exponent (log 2e0 10))))))
+                         (round (* exponent
+                                   ;; this is the closest double float
+                                   ;; to (log 2 10), but expressed so
+                                   ;; that we're not vulnerable to the
+                                   ;; host lisp's interpretation of
+                                   ;; arithmetic.  (FIXME: it turns
+                                   ;; out that sbcl itself is off by 1
+                                   ;; ulp in this value, which is a
+                                   ;; little unfortunate.)
+                                   (load-time-value
+                                    #!-long-float
+                                    (sb!kernel:make-double-float 1070810131 1352628735)
+                                    #!+long-float
+                                    (error "(log 2 10) not computed")))))))
                  (x (if (minusp ex)
                         (if (float-denormalized-p x)
                             #!-long-float
index 10987c2..fc0b639 100644 (file)
              (t
               ;; (float x (+0.0)) => (or (member -0.0) (float x (0.0)))
               ;; (float x -0.0) => (or (member -0.0) (float x (0.0)))
-              (list (make-member-type :members (list (float -0.0 hi-val)))
+              (list (make-member-type :members (list (float (load-time-value (make-unportable-float :single-float-negative-zero)) hi-val)))
                     (make-numeric-type :class (numeric-type-class type)
                                        :format (numeric-type-format type)
                                        :complexp :real
index 9e10b3c..183f979 100644 (file)
   ((fp-constant) (single-reg double-reg #!+long-float long-reg))
   (let ((value (sb!c::constant-value (sb!c::tn-leaf x))))
     (with-empty-tn@fp-top(y)
-      (cond ((zerop value)
+      (cond ((or (eql value 0f0) (eql value 0d0) #!+long-float (eql value 0l0))
              (inst fldz))
             ((= value 1e0)
              (inst fld1))
+            #!+long-float
             ((= value (coerce pi *read-default-float-format*))
              (inst fldpi))
+            #!+long-float
             ((= value (log 10e0 2e0))
              (inst fldl2t))
+            #!+long-float
             ((= value (log 2.718281828459045235360287471352662e0 2e0))
              (inst fldl2e))
+            #!+long-float
             ((= value (log 2e0 10e0))
              (inst fldlg2))
+            #!+long-float
             ((= value (log 2e0 2.718281828459045235360287471352662e0))
              (inst fldln2))
             (t (warn "ignoring bogus i387 constant ~A" value))))))
index 9d4c47f..a9d632f 100644 (file)
@@ -17,4 +17,4 @@
 ;;; checkins which aren't released. (And occasionally for internal
 ;;; versions, especially for internal versions off the main CVS
 ;;; branch, it gets hairier, e.g. "0.pre7.14.flaky4.13".)
-"1.0.27.40"
+"1.0.27.41"