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 shouldn't
  17 ;;; encounter any float traps when cross-compiling SBCL itself, anyway, so we
  18 ;;; 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, DOUBLE-FLOAT-LOW-BITS, and
  40 ;;; DOUBLE-FLOAT-HIGH-BITS
  41 ;;;
  42 ;;; KLUDGE: These will fail if the target's floating point isn't IEEE, and so
  43 ;;; I'd be more comfortable if there were an assertion "target's floating point
  44 ;;; is IEEE" in the code, but I can't see how to express that.
  45 ;;;
  46 ;;; KLUDGE: It's sort of weird that these functions return signed 32-bit values
  47 ;;; instead of unsigned 32-bit values. This is the way that the CMU CL
  48 ;;; machine-dependent functions behaved, and I've copied that behavior, but it
  49 ;;; seems to me that it'd be more idiomatic to return unsigned 32-bit values.
  50 ;;; Maybe someday the machine-dependent functions could be tweaked to return
  51 ;;; unsigned 32-bit values?
  52 (defun single-float-bits (x)
  53   (declare (type single-float x))
  54   (assert (= (float-radix x) 2))
  55   (if (zerop x)
  56     0 ; known property of IEEE floating point: 0.0 is represented as 0.
  57     (multiple-value-bind (lisp-significand lisp-exponent lisp-sign)
  58         (integer-decode-float x)
  59       (assert (plusp lisp-significand))
  60       ;; Calculate IEEE-style fields from Common-Lisp-style fields.
  61       ;;
  62       ;; KLUDGE: This code was written from my foggy memory of what IEEE
  63       ;; format looks like, augmented by some experiments with
  64       ;; the existing implementation of SINGLE-FLOAT-BITS, and what
  65       ;; I found floating around on the net at
  66       ;;   <http://www.scri.fsu.edu/~jac/MAD3401/Backgrnd/ieee.html>,
  67       ;;   <http://rodin.cs.uh.edu/~johnson2/ieee.html>,
  68       ;; and
  69       ;;   <http://www.ttu.ee/sidu/cas/IEEE_Floating.htm>.
  70       ;; And beyond the probable sheer flakiness of the code, all the bare
  71       ;; numbers floating around here are sort of ugly, too. -- WHN 19990711
  72       (let* ((significand lisp-significand)
  73              (exponent (+ lisp-exponent 23 127))
  74              (unsigned-result
  75               (if (plusp exponent) ; if not obviously denormalized
  76                 (do ()
  77                     (nil)
  78                   (cond (;; ordinary termination case
  79                          (>= significand (expt 2 23))
  80                          (assert (< 0 significand (expt 2 24)))
  81                          ;; Exponent 0 is reserved for denormalized numbers,
  82                          ;; and 255 is reserved for specials like NaN.
  83                          (assert (< 0 exponent 255))
  84                          (return (logior (ash exponent 23)
  85                                          (logand significand
  86                                                  (1- (ash 1 23))))))
  87                         (;; special termination case, denormalized float number
  88                          (zerop exponent)
  89                          ;; Denormalized numbers have exponent one greater than
  90                          ;; the exponent field.
  91                          (return (ash significand -1)))
  92                         (t
  93                          ;; Shift as necessary to set bit 24 of significand.
  94                          (setf significand (ash significand 1)
  95                                exponent (1- exponent)))))
  96                 (do ()
  97                     ((zerop exponent)
  98                      ;; Denormalized numbers have exponent one greater than the
  99                      ;; exponent field.
 100                      (ash significand -1))
 101                   (unless (zerop (logand significand 1))
 102                     (warn "denormalized SINGLE-FLOAT-BITS ~S losing bits" x))
 103                   (setf significand (ash significand -1)
 104                         exponent (1+ exponent))))))
 105         (ecase lisp-sign
 106           (1 unsigned-result)
 107           (-1 (logior unsigned-result (- (expt 2 31)))))))))
 108 (defun double-float-bits (x)
 109   (declare (type double-float x))
 110   (assert (= (float-radix x) 2))
 111   (if (zerop x)
 112     0 ; known property of IEEE floating point: 0.0d0 is represented as 0.
 113     ;; KLUDGE: As per comments in SINGLE-FLOAT-BITS, above.
 114     (multiple-value-bind (lisp-significand lisp-exponent lisp-sign)
 115         (integer-decode-float x)
 116       (assert (plusp lisp-significand))
 117       (let* ((significand lisp-significand)
 118              (exponent (+ lisp-exponent 52 1023))
 119              (unsigned-result
 120               (if (plusp exponent) ; if not obviously denormalized
 121                 (do ()
 122                     (nil)
 123                   (cond (;; ordinary termination case
 124                          (>= significand (expt 2 52))
 125                          (assert (< 0 significand (expt 2 53)))
 126                          ;; Exponent 0 is reserved for denormalized numbers,
 127                          ;; and 2047 is reserved for specials like NaN.
 128                          (assert (< 0 exponent 2047))
 129                          (return (logior (ash exponent 52)
 130                                          (logand significand
 131                                                  (1- (ash 1 52))))))
 132                         (;; special termination case, denormalized float number
 133                          (zerop exponent)
 134                          ;; Denormalized numbers have exponent one greater than
 135                          ;; the exponent field.
 136                          (return (ash significand -1)))
 137                         (t
 138                          ;; Shift as necessary to set bit 53 of significand.
 139                          (setf significand (ash significand 1)
 140                                exponent (1- exponent)))))
 141                 (do ()
 142                     ((zerop exponent)
 143                      ;; Denormalized numbers have exponent one greater than the
 144                      ;; exponent field.
 145                      (ash significand -1))
 146                   (unless (zerop (logand significand 1))
 147                     (warn "denormalized SINGLE-FLOAT-BITS ~S losing bits" x))
 148                   (setf significand (ash significand -1)
 149                         exponent (1+ exponent))))))
 150         (ecase lisp-sign
 151           (1 unsigned-result)
 152           (-1 (logior unsigned-result (- (expt 2 63)))))))))
 153 (defun double-float-low-bits (x)
 154   (declare (type double-float x))
 155   (if (zerop x)
 156     0
 157     ;; Unlike DOUBLE-FLOAT-HIGH-BITS or SINGLE-FLOAT-BITS, the CMU CL
 158     ;; DOUBLE-FLOAT-LOW-BITS seems to return a unsigned value, not a signed
 159     ;; value.
 160     (logand #xffffffff (double-float-bits x))))
 161 (defun double-float-high-bits (x)
 162   (declare (type double-float x))
 163   (if (zerop x)
 164     0
 165     (mask-and-sign-extend (ash (double-float-bits x) -32) 32)))
 166
 167 ;;; KLUDGE: These functions will blow up on any cross-compilation
 168 ;;; host Lisp which has less floating point precision than the target
 169 ;;; Lisp. In practice, this may not be a major problem: IEEE
 170 ;;; floating point arithmetic is so common these days that most
 171 ;;; cross-compilation host Lisps are likely to have exactly the same
 172 ;;; floating point precision as the target Lisp. If it turns out to be
 173 ;;; a problem, there are possible workarounds involving portable
 174 ;;; representations for target floating point numbers, like
 175 ;;;   (DEFSTRUCT TARGET-SINGLE-FLOAT
 176 ;;;     (SIGN (REQUIRED-ARGUMENT) :TYPE BIT)
 177 ;;;     (EXPONENT (REQUIRED-ARGUMENT) :TYPE UNSIGNED-BYTE)
 178 ;;;     (MANTISSA (REQUIRED-ARGUMENT) :TYPE UNSIGNED-BYTE))
 179 ;;; with some sort of MAKE-LOAD-FORM-ish magic to cause them to be
 180 ;;; written out in the appropriate target format. (And yes, those
 181 ;;; workarounds *do* look messy to me, which is why I just went
 182 ;;; with this quick kludge instead.) -- WHN 19990711
 183 (defun make-single-float (bits)
 184   (if (zerop bits) ; IEEE float special case
 185     0.0
 186     (let ((sign (ecase (ldb (byte 1 31) bits)
 187                   (0  1.0)
 188                   (1 -1.0)))
 189           (expt (- (ldb (byte 8 23) bits) 127))
 190           (mant (* (logior (ldb (byte 23 0) bits)
 191                            (ash 1 23))
 192                    (expt 0.5 23))))
 193       (* sign (expt 2.0 expt) mant))))
 194 (defun make-double-float (hi lo)
 195   (if (and (zerop hi) (zerop lo)) ; IEEE float special case
 196     0.0d0
 197     (let* ((bits (logior (ash hi 32) lo))
 198            (sign (ecase (ldb (byte 1 63) bits)
 199                    (0  1.0d0)
 200                    (1 -1.0d0)))
 201            (expt (- (ldb (byte 11 52) bits) 1023))
 202            (mant (* (logior (ldb (byte 52 0) bits)
 203                             (ash 1 52))
 204                     (expt 0.5d0 52))))
 205       (* sign (expt 2.0d0 expt) mant))))