src/code/cross-float.lisp

   1 ;;;; portable implementations or stubs for nonportable floating point
   2 ;;;; things, useful for building Python as a cross-compiler when
   3 ;;;; running under an ordinary ANSI Common Lisp implementation
   4
   5 ;;;; This software is part of the SBCL system. See the README file for
   6 ;;;; more information.
   7 ;;;;
   8 ;;;; This software is derived from the CMU CL system, which was
   9 ;;;; written at Carnegie Mellon University and released into the
  10 ;;;; public domain. The software is in the public domain and is
  11 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
  12 ;;;; files for more information.
  13
  14 (in-package "SB!IMPL")
  15
  16 ;;; There seems to be no portable way to mask float traps, but we
  17 ;;; shouldn't encounter any float traps when cross-compiling SBCL
  18 ;;; itself, anyway, so we just make this a no-op.
  19 (defmacro sb!vm::with-float-traps-masked (traps &body body)
  20   (declare (ignore traps))
  21   ;; FIXME: should become STYLE-WARNING?
  22   (format *error-output*
  23           "~&(can't portably mask float traps, proceeding anyway)~%")
  24   `(progn ,@body))
  25
  26 ;;; a helper function for DOUBLE-FLOAT-FOO-BITS functions
  27 ;;;
  28 ;;; Return the low N bits of X as a signed N-bit value.
  29 (defun mask-and-sign-extend (x n)
  30   (assert (plusp n))
  31   (let* ((high-bit (ash 1 (1- n)))
  32          (mask (1- (ash high-bit 1)))
  33          (uresult (logand mask x)))
  34     (if (zerop (logand uresult high-bit))
  35         uresult
  36         (logior uresult
  37                 (logand -1 (lognot mask))))))
  38
  39 ;;; portable implementations of SINGLE-FLOAT-BITS,
  40 ;;; DOUBLE-FLOAT-LOW-BITS, and DOUBLE-FLOAT-HIGH-BITS
  41 ;;;
  42 ;;; KLUDGE: These will fail if the target's floating point isn't IEEE,
  43 ;;; and so I'd be more comfortable if there were an assertion
  44 ;;; "target's floating point is IEEE" in the code, but I can't see how
  45 ;;; to express that.
  46 ;;;
  47 ;;; KLUDGE: It's sort of weird that these functions return signed
  48 ;;; 32-bit values instead of unsigned 32-bit values. This is the way
  49 ;;; that the CMU CL machine-dependent functions behaved, and I've
  50 ;;; copied that behavior, but it seems to me that it'd be more
  51 ;;; idiomatic to return unsigned 32-bit values. Maybe someday the
  52 ;;; machine-dependent functions could be tweaked to return unsigned
  53 ;;; 32-bit values?
  54 (defun single-float-bits (x)
  55   (declare (type single-float x))
  56   (assert (= (float-radix x) 2))
  57   (if (zerop x)
  58       (if (eql x 0.0f0) 0 #x-80000000)
  59       (multiple-value-bind (lisp-significand lisp-exponent lisp-sign)
  60           (integer-decode-float x)
  61         (assert (plusp lisp-significand))
  62         ;; Calculate IEEE-style fields from Common-Lisp-style fields.
  63         ;;
  64         ;; KLUDGE: This code was written from my foggy memory of what IEEE
  65         ;; format looks like, augmented by some experiments with
  66         ;; the existing implementation of SINGLE-FLOAT-BITS, and what
  67         ;; I found floating around on the net at
  68         ;;   <http://www.scri.fsu.edu/~jac/MAD3401/Backgrnd/ieee.html>,
  69         ;;   <http://rodin.cs.uh.edu/~johnson2/ieee.html>,
  70         ;; and
  71         ;;   <http://www.ttu.ee/sidu/cas/IEEE_Floating.htm>.
  72         ;; And beyond the probable sheer flakiness of the code, all the bare
  73         ;; numbers floating around here are sort of ugly, too. -- WHN 19990711
  74         (let* ((significand lisp-significand)
  75                (exponent (+ lisp-exponent 23 127))
  76                (unsigned-result
  77                 (if (plusp exponent)    ; if not obviously denormalized
  78                     (do ()
  79                         (nil)
  80                       (cond (;; special termination case, denormalized
  81                              ;; float number
  82                              (zerop exponent)
  83                              ;; Denormalized numbers have exponent one
  84                              ;; greater than the exponent field.
  85                              (return (ash significand -1)))
  86                             (;; ordinary termination case
  87                              (>= significand (expt 2 23))
  88                              (assert (< 0 significand (expt 2 24)))
  89                              ;; Exponent 0 is reserved for
  90                              ;; denormalized numbers, and 255 is
  91                              ;; reserved for specials like NaN.
  92                              (assert (< 0 exponent 255))
  93                              (return (logior (ash exponent 23)
  94                                              (logand significand
  95                                                      (1- (ash 1 23))))))
  96
  97                             (t
  98                              ;; Shift as necessary to set bit 24 of
  99                              ;; significand.
 100                              (setf significand (ash significand 1)
 101                                    exponent (1- exponent)))))
 102                     (do ()
 103                         ((zerop exponent)
 104                          ;; Denormalized numbers have exponent one
 105                          ;; greater than the exponent field.
 106                          (ash significand -1))
 107                       (unless (zerop (logand significand 1))
 108                         (warn "denormalized SINGLE-FLOAT-BITS ~S losing bits"
 109                               x))
 110                       (setf significand (ash significand -1)
 111                             exponent (1+ exponent))))))
 112           (ecase lisp-sign
 113             (1 unsigned-result)
 114             (-1 (logior unsigned-result (- (expt 2 31)))))))))
 115
 116 (defun double-float-bits (x)
 117   (declare (type double-float x))
 118   (assert (= (float-radix x) 2))
 119   (if (zerop x)
 120       (if (eql x 0.0d0) 0 #x-8000000000000000)
 121       ;; KLUDGE: As per comments in SINGLE-FLOAT-BITS, above.
 122       (multiple-value-bind (lisp-significand lisp-exponent lisp-sign)
 123           (integer-decode-float x)
 124         (assert (plusp lisp-significand))
 125         (let* ((significand lisp-significand)
 126                (exponent (+ lisp-exponent 52 1023))
 127                (unsigned-result
 128                 (if (plusp exponent)    ; if not obviously denormalized
 129                     (do ()
 130                         (nil)
 131                       (cond (;; special termination case, denormalized
 132                              ;; float number
 133                              (zerop exponent)
 134                              ;; Denormalized numbers have exponent one
 135                              ;; greater than the exponent field.
 136                              (return (ash significand -1)))
 137                             (;; ordinary termination case
 138                              (>= significand (expt 2 52))
 139                              (assert (< 0 significand (expt 2 53)))
 140                              ;; Exponent 0 is reserved for
 141                              ;; denormalized numbers, and 2047 is
 142                              ;; reserved for specials like NaN.
 143                              (assert (< 0 exponent 2047))
 144                              (return (logior (ash exponent 52)
 145                                              (logand significand
 146                                                      (1- (ash 1 52))))))
 147                             (t
 148                              ;; Shift as necessary to set bit 53 of
 149                              ;; significand.
 150                              (setf significand (ash significand 1)
 151                                    exponent (1- exponent)))))
 152                     (do ()
 153                         ((zerop exponent)
 154                          ;; Denormalized numbers have exponent one
 155                          ;; greater than the exponent field.
 156                          (ash significand -1))
 157                       (unless (zerop (logand significand 1))
 158                         (warn "denormalized SINGLE-FLOAT-BITS ~S losing bits"
 159                               x))
 160                       (setf significand (ash significand -1)
 161                             exponent (1+ exponent))))))
 162           (ecase lisp-sign
 163             (1 unsigned-result)
 164             (-1 (logior unsigned-result (- (expt 2 63)))))))))
 165
 166 (defun double-float-low-bits (x)
 167   (declare (type double-float x))
 168   (if (zerop x)
 169       0
 170       ;; FIXME: Unlike DOUBLE-FLOAT-HIGH-BITS or SINGLE-FLOAT-BITS,
 171       ;; the CMU CL DOUBLE-FLOAT-LOW-BITS seemed to return a unsigned
 172       ;; value, not a signed value, so we've done the same. But it
 173       ;; would be nice to make the family of functions have a more
 174       ;; consistent return convention.
 175       (logand #xffffffff (double-float-bits x))))
 176
 177 (defun double-float-high-bits (x)
 178   (declare (type double-float x))
 179   (if (zerop x)
 180       (if (eql x 0.0d0) 0 #x-80000000)
 181       (mask-and-sign-extend (ash (double-float-bits x) -32) 32)))
 182
 183 ;;; KLUDGE: This is a hack to work around a bug in CMU CL 18c which
 184 ;;; causes the 18c compiler to die with a floating point exception
 185 ;;; when trying to optimize the EXPT forms in the MAKE-FOO-FLOAT
 186 ;;; functions below. See the message
 187 ;;;   Subject: Re: Compiler bug?
 188 ;;;   From: Raymond Toy
 189 ;;;   Date: 28 Mar 2001 08:19:59 -0500
 190 ;;; on the cmucl-imp mailing list. Once the CMU CL folks
 191 ;;; make a bug-fix release, we can get rid of this and go back to
 192 ;;; calling EXPT directly. -- WHN 2001-04-05
 193 (defun kludge-opaque-expt (x y)
 194   (expt x y))
 195
 196 ;;; KLUDGE: These functions will blow up on any cross-compilation
 197 ;;; host Lisp which has less floating point precision than the target
 198 ;;; Lisp. In practice, this may not be a major problem: IEEE
 199 ;;; floating point arithmetic is so common these days that most
 200 ;;; cross-compilation host Lisps are likely to have exactly the same
 201 ;;; floating point precision as the target Lisp. If it turns out to be
 202 ;;; a problem, there are possible workarounds involving portable
 203 ;;; representations for target floating point numbers, like
 204 ;;;   (DEFSTRUCT TARGET-SINGLE-FLOAT
 205 ;;;     (SIGN (MISSING-ARG) :TYPE BIT)
 206 ;;;     (EXPONENT (MISSING-ARG) :TYPE UNSIGNED-BYTE)
 207 ;;;     (MANTISSA (MISSING-ARG) :TYPE UNSIGNED-BYTE))
 208 ;;; with some sort of MAKE-LOAD-FORM-ish magic to cause them to be
 209 ;;; written out in the appropriate target format. (And yes, those
 210 ;;; workarounds *do* look messy to me, which is why I just went
 211 ;;; with this quick kludge instead.) -- WHN 19990711
 212 (defun make-single-float (bits)
 213   (cond
 214     ;; IEEE float special cases
 215     ((zerop bits) 0.0)
 216     ((= bits #x-80000000) -0.0)
 217     (t (let* ((sign (ecase (ldb (byte 1 31) bits)
 218                       (0  1.0)
 219                       (1 -1.0)))
 220               (iexpt (ldb (byte 8 23) bits))
 221               (expt (if (zerop iexpt) ; denormalized
 222                         -126
 223                         (- iexpt 127)))
 224               (mant (* (logior (ldb (byte 23 0) bits)
 225                                (if (zerop iexpt)
 226                                    0
 227                                    (ash 1 23)))
 228                        (expt 0.5 23))))
 229          (* sign (kludge-opaque-expt 2.0 expt) mant)))))
 230
 231 (defun make-double-float (hi lo)
 232   (cond
 233     ;; IEEE float special cases
 234     ((and (zerop hi) (zerop lo)) 0.0d0)
 235     ((and (= hi #x-80000000) (zerop lo)) -0.0d0)
 236     (t (let* ((bits (logior (ash hi 32) lo))
 237               (sign (ecase (ldb (byte 1 63) bits)
 238                       (0  1.0d0)
 239                       (1 -1.0d0)))
 240               (iexpt (ldb (byte 11 52) bits))
 241               (expt (if (zerop iexpt) ; denormalized
 242                         -1022
 243                         (- iexpt 1023)))
 244               (mant (* (logior (ldb (byte 52 0) bits)
 245                                (if (zerop iexpt)
 246                                    0
 247                                    (ash 1 52)))
 248                        (expt 0.5d0 52))))
 249          (* sign (kludge-opaque-expt 2.0d0 expt) mant)))))
 250