0.7.13.7:
[sbcl.git] / 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       0 ; known property of IEEE floating point: 0.0 is represented as 0.
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       0 ; known property of IEEE floating point: 0.0d0 is represented as 0.
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       0
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   (if (zerop bits) ; IEEE float special case
214       0.0
215       (let ((sign (ecase (ldb (byte 1 31) bits)
216                     (0  1.0)
217                     (1 -1.0)))
218             (expt (- (ldb (byte 8 23) bits) 127))
219             (mant (* (logior (ldb (byte 23 0) bits)
220                              (ash 1 23))
221                      (expt 0.5 23))))
222         (* sign (kludge-opaque-expt 2.0 expt) mant))))
223
224 (defun make-double-float (hi lo)
225   (if (and (zerop hi) (zerop lo)) ; IEEE float special case
226       0.0d0
227       (let* ((bits (logior (ash hi 32) lo))
228              (sign (ecase (ldb (byte 1 63) bits)
229                      (0  1.0d0)
230                      (1 -1.0d0)))
231              (expt (- (ldb (byte 11 52) bits) 1023))
232              (mant (* (logior (ldb (byte 52 0) bits)
233                               (ash 1 52))
234                       (expt 0.5d0 52))))
235         (* sign (kludge-opaque-expt 2.0d0 expt) mant))))