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