4f411f20c275999b44aa73b4e33fc3e53df6f87a
[sbcl.git] / src / code / float.lisp
1 ;;;; This file contains the definitions of float-specific number
2 ;;;; support (other than irrational stuff, which is in irrat.) There is
3 ;;;; code in here that assumes there are only two float formats: IEEE
4 ;;;; single and double. (LONG-FLOAT support has been added, but bugs
5 ;;;; may still remain due to old code which assumes this dichotomy.)
6
7 ;;;; This software is part of the SBCL system. See the README file for
8 ;;;; more information.
9 ;;;;
10 ;;;; This software is derived from the CMU CL system, which was
11 ;;;; written at Carnegie Mellon University and released into the
12 ;;;; public domain. The software is in the public domain and is
13 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
14 ;;;; files for more information.
15
16 (in-package "SB!KERNEL")
17 \f
18 ;;;; utilities
19
20 (eval-when (:compile-toplevel :load-toplevel :execute)
21
22 ;;; These functions let us create floats from bits with the
23 ;;; significand uniformly represented as an integer. This is less
24 ;;; efficient for double floats, but is more convenient when making
25 ;;; special values, etc.
26 (defun single-from-bits (sign exp sig)
27   (declare (type bit sign) (type (unsigned-byte 24) sig)
28            (type (unsigned-byte 8) exp))
29   (make-single-float
30    (dpb exp sb!vm:single-float-exponent-byte
31         (dpb sig sb!vm:single-float-significand-byte
32              (if (zerop sign) 0 -1)))))
33 (defun double-from-bits (sign exp sig)
34   (declare (type bit sign) (type (unsigned-byte 53) sig)
35            (type (unsigned-byte 11) exp))
36   (make-double-float (dpb exp sb!vm:double-float-exponent-byte
37                           (dpb (ash sig -32)
38                                sb!vm:double-float-significand-byte
39                                (if (zerop sign) 0 -1)))
40                      (ldb (byte 32 0) sig)))
41 #!+(and long-float x86)
42 (defun long-from-bits (sign exp sig)
43   (declare (type bit sign) (type (unsigned-byte 64) sig)
44            (type (unsigned-byte 15) exp))
45   (make-long-float (logior (ash sign 15) exp)
46                    (ldb (byte 32 32) sig)
47                    (ldb (byte 32 0) sig)))
48                                         
49 ) ; EVAL-WHEN
50 \f
51 ;;;; float parameters
52
53 (defconstant least-positive-single-float (single-from-bits 0 0 1))
54 (defconstant least-positive-short-float least-positive-single-float)
55 (defconstant least-negative-single-float (single-from-bits 1 0 1))
56 (defconstant least-negative-short-float least-negative-single-float)
57 (defconstant least-positive-double-float (double-from-bits 0 0 1))
58 #!-long-float
59 (defconstant least-positive-long-float least-positive-double-float)
60 #!+(and long-float x86)
61 (defconstant least-positive-long-float (long-from-bits 0 0 1))
62 (defconstant least-negative-double-float (double-from-bits 1 0 1))
63 #!-long-float
64 (defconstant least-negative-long-float least-negative-double-float)
65 #!+(and long-float x86)
66 (defconstant least-negative-long-float (long-from-bits 1 0 1))
67
68 (defconstant least-positive-normalized-single-float
69   (single-from-bits 0 sb!vm:single-float-normal-exponent-min 0))
70 (defconstant least-positive-normalized-short-float
71   least-positive-normalized-single-float)
72 (defconstant least-negative-normalized-single-float
73   (single-from-bits 1 sb!vm:single-float-normal-exponent-min 0))
74 (defconstant least-negative-normalized-short-float
75   least-negative-normalized-single-float)
76 (defconstant least-positive-normalized-double-float
77   (double-from-bits 0 sb!vm:double-float-normal-exponent-min 0))
78 #!-long-float
79 (defconstant least-positive-normalized-long-float
80   least-positive-normalized-double-float)
81 #!+(and long-float x86)
82 (defconstant least-positive-normalized-long-float
83   (long-from-bits 0 sb!vm:long-float-normal-exponent-min
84                   (ash sb!vm:long-float-hidden-bit 32)))
85 (defconstant least-negative-normalized-double-float
86   (double-from-bits 1 sb!vm:double-float-normal-exponent-min 0))
87 #!-long-float
88 (defconstant least-negative-normalized-long-float
89   least-negative-normalized-double-float)
90 #!+(and long-float x86)
91 (defconstant least-negative-normalized-long-float
92   (long-from-bits 1 sb!vm:long-float-normal-exponent-min
93                   (ash sb!vm:long-float-hidden-bit 32)))
94
95 (defconstant most-positive-single-float
96   (single-from-bits 0 sb!vm:single-float-normal-exponent-max
97                     (ldb sb!vm:single-float-significand-byte -1)))
98 (defconstant most-positive-short-float most-positive-single-float)
99 (defconstant most-negative-single-float
100   (single-from-bits 1 sb!vm:single-float-normal-exponent-max
101                     (ldb sb!vm:single-float-significand-byte -1)))
102 (defconstant most-negative-short-float most-negative-single-float)
103 (defconstant most-positive-double-float
104   (double-from-bits 0 sb!vm:double-float-normal-exponent-max
105                     (ldb (byte sb!vm:double-float-digits 0) -1)))
106 #!-long-float
107 (defconstant most-positive-long-float most-positive-double-float)
108 #!+(and long-float x86)
109 (defconstant most-positive-long-float
110   (long-from-bits 0 sb!vm:long-float-normal-exponent-max
111                   (ldb (byte sb!vm:long-float-digits 0) -1)))
112 (defconstant most-negative-double-float
113   (double-from-bits 1 sb!vm:double-float-normal-exponent-max
114                     (ldb (byte sb!vm:double-float-digits 0) -1)))
115 #!-long-float
116 (defconstant most-negative-long-float most-negative-double-float)
117 #!+(and long-float x86)
118 (defconstant most-negative-long-float
119   (long-from-bits 1 sb!vm:long-float-normal-exponent-max
120                   (ldb (byte sb!vm:long-float-digits 0) -1)))
121
122 ;;; We don't want to do these DEFCONSTANTs at cross-compilation time,
123 ;;; because the cross-compilation host might not support floating
124 ;;; point infinities. Putting them inside a LET removes
125 ;;; toplevel-formness, so that any EVAL-WHEN trickiness in the
126 ;;; DEFCONSTANT forms is suppressed.
127 (let ()
128 (defconstant single-float-positive-infinity
129   (single-from-bits 0 (1+ sb!vm:single-float-normal-exponent-max) 0))
130 (defconstant short-float-positive-infinity single-float-positive-infinity)
131 (defconstant single-float-negative-infinity
132   (single-from-bits 1 (1+ sb!vm:single-float-normal-exponent-max) 0))
133 (defconstant short-float-negative-infinity single-float-negative-infinity)
134 (defconstant double-float-positive-infinity
135   (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0))
136 #!+(not long-float)
137 (defconstant long-float-positive-infinity double-float-positive-infinity)
138 #!+(and long-float x86)
139 (defconstant long-float-positive-infinity
140   (long-from-bits 0 (1+ sb!vm:long-float-normal-exponent-max)
141                   (ash sb!vm:long-float-hidden-bit 32)))
142 (defconstant double-float-negative-infinity
143   (double-from-bits 1 (1+ sb!vm:double-float-normal-exponent-max) 0))
144 #!+(not long-float)
145 (defconstant long-float-negative-infinity double-float-negative-infinity)
146 #!+(and long-float x86)
147 (defconstant long-float-negative-infinity
148   (long-from-bits 1 (1+ sb!vm:long-float-normal-exponent-max)
149                   (ash sb!vm:long-float-hidden-bit 32)))
150 ) ; LET-to-suppress-possible-EVAL-WHENs
151
152 (defconstant single-float-epsilon
153   (single-from-bits 0 (- sb!vm:single-float-bias
154                          (1- sb!vm:single-float-digits)) 1))
155 (defconstant short-float-epsilon single-float-epsilon)
156 (defconstant single-float-negative-epsilon
157   (single-from-bits 0 (- sb!vm:single-float-bias sb!vm:single-float-digits) 1))
158 (defconstant short-float-negative-epsilon single-float-negative-epsilon)
159 (defconstant double-float-epsilon
160   (double-from-bits 0 (- sb!vm:double-float-bias
161                          (1- sb!vm:double-float-digits)) 1))
162 #!-long-float
163 (defconstant long-float-epsilon double-float-epsilon)
164 #!+(and long-float x86)
165 (defconstant long-float-epsilon
166   (long-from-bits 0 (- sb!vm:long-float-bias (1- sb!vm:long-float-digits))
167                   (+ 1 (ash sb!vm:long-float-hidden-bit 32))))
168 (defconstant double-float-negative-epsilon
169   (double-from-bits 0 (- sb!vm:double-float-bias sb!vm:double-float-digits) 1))
170 #!-long-float
171 (defconstant long-float-negative-epsilon double-float-negative-epsilon)
172 #!+(and long-float x86)
173 (defconstant long-float-negative-epsilon
174   (long-from-bits 0 (- sb!vm:long-float-bias sb!vm:long-float-digits)
175                   (+ 1 (ash sb!vm:long-float-hidden-bit 32))))
176 \f
177 ;;;; float predicates and environment query
178
179 #!-sb-fluid
180 (declaim (maybe-inline float-denormalized-p float-infinity-p float-nan-p
181                        float-trapping-nan-p))
182
183 (defun float-denormalized-p (x)
184   #!+sb-doc
185   "Return true if the float X is denormalized."
186   (number-dispatch ((x float))
187     ((single-float)
188      (and (zerop (ldb sb!vm:single-float-exponent-byte (single-float-bits x)))
189           (not (zerop x))))
190     ((double-float)
191      (and (zerop (ldb sb!vm:double-float-exponent-byte
192                       (double-float-high-bits x)))
193           (not (zerop x))))
194     #!+(and long-float x86)
195     ((long-float)
196      (and (zerop (ldb sb!vm:long-float-exponent-byte (long-float-exp-bits x)))
197           (not (zerop x))))))
198
199 (macrolet ((def (name doc single double #!+(and long-float x86) long)
200              `(defun ,name (x)
201                 ,doc
202                 (number-dispatch ((x float))
203                   ((single-float)
204                    (let ((bits (single-float-bits x)))
205                      (and (> (ldb sb!vm:single-float-exponent-byte bits)
206                              sb!vm:single-float-normal-exponent-max)
207                           ,single)))
208                   ((double-float)
209                    (let ((hi (double-float-high-bits x))
210                          (lo (double-float-low-bits x)))
211                      (declare (ignorable lo))
212                      (and (> (ldb sb!vm:double-float-exponent-byte hi)
213                              sb!vm:double-float-normal-exponent-max)
214                           ,double)))
215                   #!+(and long-float x86)
216                   ((long-float)
217                    (let ((exp (long-float-exp-bits x))
218                          (hi (long-float-high-bits x))
219                          (lo (long-float-low-bits x)))
220                      (declare (ignorable lo))
221                      (and (> (ldb sb!vm:long-float-exponent-byte exp)
222                              sb!vm:long-float-normal-exponent-max)
223                           ,long)))))))
224
225   (def float-infinity-p
226     "Return true if the float X is an infinity (+ or -)."
227     (zerop (ldb sb!vm:single-float-significand-byte bits))
228     (and (zerop (ldb sb!vm:double-float-significand-byte hi))
229          (zerop lo))
230     #!+(and long-float x86)
231     (and (zerop (ldb sb!vm:long-float-significand-byte hi))
232          (zerop lo)))
233
234   (def float-nan-p
235     "Return true if the float X is a NaN (Not a Number)."
236     (not (zerop (ldb sb!vm:single-float-significand-byte bits)))
237     (or (not (zerop (ldb sb!vm:double-float-significand-byte hi)))
238         (not (zerop lo)))
239     #!+(and long-float x86)
240     (or (not (zerop (ldb sb!vm:long-float-significand-byte hi)))
241         (not (zerop lo))))
242
243   (def float-trapping-nan-p
244     "Return true if the float X is a trapping NaN (Not a Number)."
245     (zerop (logand (ldb sb!vm:single-float-significand-byte bits)
246                    sb!vm:single-float-trapping-nan-bit))
247     (zerop (logand (ldb sb!vm:double-float-significand-byte hi)
248                    sb!vm:double-float-trapping-nan-bit))
249     #!+(and long-float x86)
250     (zerop (logand (ldb sb!vm:long-float-significand-byte hi)
251                    sb!vm:long-float-trapping-nan-bit))))
252
253 ;;; If denormalized, use a subfunction from INTEGER-DECODE-FLOAT to find the
254 ;;; actual exponent (and hence how denormalized it is), otherwise we just
255 ;;; return the number of digits or 0.
256 #!-sb-fluid (declaim (maybe-inline float-precision))
257 (defun float-precision (f)
258   #!+sb-doc
259   "Return a non-negative number of significant digits in its float argument.
260   Will be less than FLOAT-DIGITS if denormalized or zero."
261   (macrolet ((frob (digits bias decode)
262                `(cond ((zerop f) 0)
263                       ((float-denormalized-p f)
264                        (multiple-value-bind (ignore exp) (,decode f)
265                          (declare (ignore ignore))
266                          (truly-the fixnum
267                                     (+ ,digits (1- ,digits) ,bias exp))))
268                       (t
269                        ,digits))))
270     (number-dispatch ((f float))
271       ((single-float)
272        (frob sb!vm:single-float-digits sb!vm:single-float-bias
273          integer-decode-single-denorm))
274       ((double-float)
275        (frob sb!vm:double-float-digits sb!vm:double-float-bias
276          integer-decode-double-denorm))
277       #!+long-float
278       ((long-float)
279        (frob sb!vm:long-float-digits sb!vm:long-float-bias
280          integer-decode-long-denorm)))))
281
282 (defun float-sign (float1 &optional (float2 (float 1 float1)))
283   #!+sb-doc
284   "Return a floating-point number that has the same sign as
285    float1 and, if float2 is given, has the same absolute value
286    as float2."
287   (declare (float float1 float2))
288   (* (if (etypecase float1
289            (single-float (minusp (single-float-bits float1)))
290            (double-float (minusp (double-float-high-bits float1)))
291            #!+long-float
292            (long-float (minusp (long-float-exp-bits float1))))
293          (float -1 float1)
294          (float 1 float1))
295      (abs float2)))
296
297 (defun float-format-digits (format)
298   (ecase format
299     ((short-float single-float) sb!vm:single-float-digits)
300     ((double-float #!-long-float long-float) sb!vm:double-float-digits)
301     #!+long-float
302     (long-float sb!vm:long-float-digits)))
303
304 #!-sb-fluid (declaim (inline float-digits float-radix))
305
306 (defun float-digits (f)
307   (number-dispatch ((f float))
308     ((single-float) sb!vm:single-float-digits)
309     ((double-float) sb!vm:double-float-digits)
310     #!+long-float
311     ((long-float) sb!vm:long-float-digits)))
312
313 (setf (fdefinition 'float-radix)
314       ;; FIXME: Python flushes unused variable X in CLAMBDA, then
315       ;; flushes unused reference to X in XEP together with type
316       ;; check. When this is fixed, rewrite this definition in an
317       ;; ordinary form. -- APD, 2002-10-21
318       (lambda (x)
319         #!+sb-doc
320         "Return (as an integer) the radix b of its floating-point argument."
321         (unless (floatp x)
322           (error 'type-error :datum x :expected-type 'float))
323         2))
324 \f
325 ;;;; INTEGER-DECODE-FLOAT and DECODE-FLOAT
326
327 #!-sb-fluid
328 (declaim (maybe-inline integer-decode-single-float
329                        integer-decode-double-float))
330
331 ;;; Handle the denormalized case of INTEGER-DECODE-FLOAT for SINGLE-FLOAT.
332 (defun integer-decode-single-denorm (x)
333   (declare (type single-float x))
334   (let* ((bits (single-float-bits (abs x)))
335          (sig (ash (ldb sb!vm:single-float-significand-byte bits) 1))
336          (extra-bias 0))
337     (declare (type (unsigned-byte 24) sig)
338              (type (integer 0 23) extra-bias))
339     (loop
340       (unless (zerop (logand sig sb!vm:single-float-hidden-bit))
341         (return))
342       (setq sig (ash sig 1))
343       (incf extra-bias))
344     (values sig
345             (- (- sb!vm:single-float-bias)
346                sb!vm:single-float-digits
347                extra-bias)
348             (if (minusp (float-sign x)) -1 1))))
349
350 ;;; Handle the single-float case of INTEGER-DECODE-FLOAT. If an infinity or
351 ;;; NaN, error. If a denorm, call i-d-s-DENORM to handle it.
352 (defun integer-decode-single-float (x)
353   (declare (single-float x))
354   (let* ((bits (single-float-bits (abs x)))
355          (exp (ldb sb!vm:single-float-exponent-byte bits))
356          (sig (ldb sb!vm:single-float-significand-byte bits))
357          (sign (if (minusp (float-sign x)) -1 1))
358          (biased (- exp sb!vm:single-float-bias sb!vm:single-float-digits)))
359     (declare (fixnum biased))
360     (unless (<= exp sb!vm:single-float-normal-exponent-max)
361       (error "can't decode NaN or infinity: ~S" x))
362     (cond ((and (zerop exp) (zerop sig))
363            (values 0 biased sign))
364           ((< exp sb!vm:single-float-normal-exponent-min)
365            (integer-decode-single-denorm x))
366           (t
367            (values (logior sig sb!vm:single-float-hidden-bit) biased sign)))))
368
369 ;;; like INTEGER-DECODE-SINGLE-DENORM, only doubly so
370 (defun integer-decode-double-denorm (x)
371   (declare (type double-float x))
372   (let* ((high-bits (double-float-high-bits (abs x)))
373          (sig-high (ldb sb!vm:double-float-significand-byte high-bits))
374          (low-bits (double-float-low-bits x))
375          (sign (if (minusp (float-sign x)) -1 1))
376          (biased (- (- sb!vm:double-float-bias) sb!vm:double-float-digits)))
377     (if (zerop sig-high)
378         (let ((sig low-bits)
379               (extra-bias (- sb!vm:double-float-digits 33))
380               (bit (ash 1 31)))
381           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
382           (loop
383             (unless (zerop (logand sig bit)) (return))
384             (setq sig (ash sig 1))
385             (incf extra-bias))
386           (values (ash sig (- sb!vm:double-float-digits 32))
387                   (truly-the fixnum (- biased extra-bias))
388                   sign))
389         (let ((sig (ash sig-high 1))
390               (extra-bias 0))
391           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
392           (loop
393             (unless (zerop (logand sig sb!vm:double-float-hidden-bit))
394               (return))
395             (setq sig (ash sig 1))
396             (incf extra-bias))
397           (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
398                   (truly-the fixnum (- biased extra-bias))
399                   sign)))))
400
401 ;;; like INTEGER-DECODE-SINGLE-FLOAT, only doubly so
402 (defun integer-decode-double-float (x)
403   (declare (double-float x))
404   (let* ((abs (abs x))
405          (hi (double-float-high-bits abs))
406          (lo (double-float-low-bits abs))
407          (exp (ldb sb!vm:double-float-exponent-byte hi))
408          (sig (ldb sb!vm:double-float-significand-byte hi))
409          (sign (if (minusp (float-sign x)) -1 1))
410          (biased (- exp sb!vm:double-float-bias sb!vm:double-float-digits)))
411     (declare (fixnum biased))
412     (unless (<= exp sb!vm:double-float-normal-exponent-max)
413       (error "Can't decode NaN or infinity: ~S." x))
414     (cond ((and (zerop exp) (zerop sig) (zerop lo))
415            (values 0 biased sign))
416           ((< exp sb!vm:double-float-normal-exponent-min)
417            (integer-decode-double-denorm x))
418           (t
419            (values
420             (logior (ash (logior (ldb sb!vm:double-float-significand-byte hi)
421                                  sb!vm:double-float-hidden-bit)
422                          32)
423                     lo)
424             biased sign)))))
425
426 #!+(and long-float x86)
427 (defun integer-decode-long-denorm (x)
428   (declare (type long-float x))
429   (let* ((high-bits (long-float-high-bits (abs x)))
430          (sig-high (ldb sb!vm:long-float-significand-byte high-bits))
431          (low-bits (long-float-low-bits x))
432          (sign (if (minusp (float-sign x)) -1 1))
433          (biased (- (- sb!vm:long-float-bias) sb!vm:long-float-digits)))
434     (if (zerop sig-high)
435         (let ((sig low-bits)
436               (extra-bias (- sb!vm:long-float-digits 33))
437               (bit (ash 1 31)))
438           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
439           (loop
440             (unless (zerop (logand sig bit)) (return))
441             (setq sig (ash sig 1))
442             (incf extra-bias))
443           (values (ash sig (- sb!vm:long-float-digits 32))
444                   (truly-the fixnum (- biased extra-bias))
445                   sign))
446         (let ((sig (ash sig-high 1))
447               (extra-bias 0))
448           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
449           (loop
450             (unless (zerop (logand sig sb!vm:long-float-hidden-bit))
451               (return))
452             (setq sig (ash sig 1))
453             (incf extra-bias))
454           (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
455                   (truly-the fixnum (- biased extra-bias))
456                   sign)))))
457
458 #!+(and long-float x86)
459 (defun integer-decode-long-float (x)
460   (declare (long-float x))
461   (let* ((hi (long-float-high-bits x))
462          (lo (long-float-low-bits x))
463          (exp-bits (long-float-exp-bits x))
464          (exp (ldb sb!vm:long-float-exponent-byte exp-bits))
465          (sign (if (minusp exp-bits) -1 1))
466          (biased (- exp sb!vm:long-float-bias sb!vm:long-float-digits)))
467     (declare (fixnum biased))
468     (unless (<= exp sb!vm:long-float-normal-exponent-max)
469       (error "can't decode NaN or infinity: ~S" x))
470     (cond ((and (zerop exp) (zerop hi) (zerop lo))
471            (values 0 biased sign))
472           ((< exp sb!vm:long-float-normal-exponent-min)
473            (integer-decode-long-denorm x))
474           (t
475            (values (logior (ash hi 32) lo) biased sign)))))
476
477 ;;; Dispatch to the correct type-specific i-d-f function.
478 (defun integer-decode-float (x)
479   #!+sb-doc
480   "Return three values:
481    1) an integer representation of the significand.
482    2) the exponent for the power of 2 that the significand must be multiplied
483       by to get the actual value. This differs from the DECODE-FLOAT exponent
484       by FLOAT-DIGITS, since the significand has been scaled to have all its
485       digits before the radix point.
486    3) -1 or 1 (i.e. the sign of the argument.)"
487   (number-dispatch ((x float))
488     ((single-float)
489      (integer-decode-single-float x))
490     ((double-float)
491      (integer-decode-double-float x))
492     #!+long-float
493     ((long-float)
494      (integer-decode-long-float x))))
495
496 #!-sb-fluid (declaim (maybe-inline decode-single-float decode-double-float))
497
498 ;;; Handle the denormalized case of DECODE-SINGLE-FLOAT. We call
499 ;;; INTEGER-DECODE-SINGLE-DENORM and then make the result into a float.
500 (defun decode-single-denorm (x)
501   (declare (type single-float x))
502   (multiple-value-bind (sig exp sign) (integer-decode-single-denorm x)
503     (values (make-single-float
504              (dpb sig sb!vm:single-float-significand-byte
505                   (dpb sb!vm:single-float-bias
506                        sb!vm:single-float-exponent-byte
507                        0)))
508             (truly-the fixnum (+ exp sb!vm:single-float-digits))
509             (float sign x))))
510
511 ;;; Handle the single-float case of DECODE-FLOAT. If an infinity or NaN,
512 ;;; error. If a denorm, call d-s-DENORM to handle it.
513 (defun decode-single-float (x)
514   (declare (single-float x))
515   (let* ((bits (single-float-bits (abs x)))
516          (exp (ldb sb!vm:single-float-exponent-byte bits))
517          (sign (float-sign x))
518          (biased (truly-the single-float-exponent
519                             (- exp sb!vm:single-float-bias))))
520     (unless (<= exp sb!vm:single-float-normal-exponent-max)
521       (error "can't decode NaN or infinity: ~S" x))
522     (cond ((zerop x)
523            (values 0.0f0 biased sign))
524           ((< exp sb!vm:single-float-normal-exponent-min)
525            (decode-single-denorm x))
526           (t
527            (values (make-single-float
528                     (dpb sb!vm:single-float-bias
529                          sb!vm:single-float-exponent-byte
530                          bits))
531                    biased sign)))))
532
533 ;;; like DECODE-SINGLE-DENORM, only doubly so
534 (defun decode-double-denorm (x)
535   (declare (double-float x))
536   (multiple-value-bind (sig exp sign) (integer-decode-double-denorm x)
537     (values (make-double-float
538              (dpb (logand (ash sig -32) (lognot sb!vm:double-float-hidden-bit))
539                   sb!vm:double-float-significand-byte
540                   (dpb sb!vm:double-float-bias
541                        sb!vm:double-float-exponent-byte 0))
542              (ldb (byte 32 0) sig))
543             (truly-the fixnum (+ exp sb!vm:double-float-digits))
544             (float sign x))))
545
546 ;;; like DECODE-SINGLE-FLOAT, only doubly so
547 (defun decode-double-float (x)
548   (declare (double-float x))
549   (let* ((abs (abs x))
550          (hi (double-float-high-bits abs))
551          (lo (double-float-low-bits abs))
552          (exp (ldb sb!vm:double-float-exponent-byte hi))
553          (sign (float-sign x))
554          (biased (truly-the double-float-exponent
555                             (- exp sb!vm:double-float-bias))))
556     (unless (<= exp sb!vm:double-float-normal-exponent-max)
557       (error "can't decode NaN or infinity: ~S" x))
558     (cond ((zerop x)
559            (values 0.0d0 biased sign))
560           ((< exp sb!vm:double-float-normal-exponent-min)
561            (decode-double-denorm x))
562           (t
563            (values (make-double-float
564                     (dpb sb!vm:double-float-bias
565                          sb!vm:double-float-exponent-byte hi)
566                     lo)
567                    biased sign)))))
568
569 #!+(and long-float x86)
570 (defun decode-long-denorm (x)
571   (declare (long-float x))
572   (multiple-value-bind (sig exp sign) (integer-decode-long-denorm x)
573     (values (make-long-float sb!vm:long-float-bias (ash sig -32)
574                              (ldb (byte 32 0) sig))
575             (truly-the fixnum (+ exp sb!vm:long-float-digits))
576             (float sign x))))
577
578 #!+(and long-float x86)
579 (defun decode-long-float (x)
580   (declare (long-float x))
581   (let* ((hi (long-float-high-bits x))
582          (lo (long-float-low-bits x))
583          (exp-bits (long-float-exp-bits x))
584          (exp (ldb sb!vm:long-float-exponent-byte exp-bits))
585          (sign (if (minusp exp-bits) -1l0 1l0))
586          (biased (truly-the long-float-exponent
587                             (- exp sb!vm:long-float-bias))))
588     (unless (<= exp sb!vm:long-float-normal-exponent-max)
589       (error "can't decode NaN or infinity: ~S" x))
590     (cond ((zerop x)
591            (values 0.0l0 biased sign))
592           ((< exp sb!vm:long-float-normal-exponent-min)
593            (decode-long-denorm x))
594           (t
595            (values (make-long-float
596                     (dpb sb!vm:long-float-bias sb!vm:long-float-exponent-byte
597                          exp-bits)
598                     hi
599                     lo)
600                    biased sign)))))
601
602 ;;; Dispatch to the appropriate type-specific function.
603 (defun decode-float (f)
604   #!+sb-doc
605   "Return three values:
606    1) a floating-point number representing the significand. This is always
607       between 0.5 (inclusive) and 1.0 (exclusive).
608    2) an integer representing the exponent.
609    3) -1.0 or 1.0 (i.e. the sign of the argument.)"
610   (number-dispatch ((f float))
611     ((single-float)
612      (decode-single-float f))
613     ((double-float)
614      (decode-double-float f))
615     #!+long-float
616     ((long-float)
617      (decode-long-float f))))
618 \f
619 ;;;; SCALE-FLOAT
620
621 #!-sb-fluid (declaim (maybe-inline scale-single-float scale-double-float))
622
623 ;;; Handle float scaling where the X is denormalized or the result is
624 ;;; denormalized or underflows to 0.
625 (defun scale-float-maybe-underflow (x exp)
626   (multiple-value-bind (sig old-exp) (integer-decode-float x)
627     (let* ((digits (float-digits x))
628            (new-exp (+ exp old-exp digits
629                        (etypecase x
630                          (single-float sb!vm:single-float-bias)
631                          (double-float sb!vm:double-float-bias))))
632            (sign (if (minusp (float-sign x)) 1 0)))
633       (cond
634        ((< new-exp
635            (etypecase x
636              (single-float sb!vm:single-float-normal-exponent-min)
637              (double-float sb!vm:double-float-normal-exponent-min)))
638         (when (sb!vm:current-float-trap :inexact)
639           (error 'floating-point-inexact :operation 'scale-float
640                  :operands (list x exp)))
641         (when (sb!vm:current-float-trap :underflow)
642           (error 'floating-point-underflow :operation 'scale-float
643                  :operands (list x exp)))
644         (let ((shift (1- new-exp)))
645           (if (< shift (- (1- digits)))
646               (float-sign x 0.0)
647               (etypecase x
648                 (single-float (single-from-bits sign 0 (ash sig shift)))
649                 (double-float (double-from-bits sign 0 (ash sig shift)))))))
650        (t
651         (etypecase x
652           (single-float (single-from-bits sign new-exp sig))
653           (double-float (double-from-bits sign new-exp sig))))))))
654
655 ;;; Called when scaling a float overflows, or the original float was a
656 ;;; NaN or infinity. If overflow errors are trapped, then error,
657 ;;; otherwise return the appropriate infinity. If a NaN, signal or not
658 ;;; as appropriate.
659 (defun scale-float-maybe-overflow (x exp)
660   (cond
661    ((float-infinity-p x)
662     ;; Infinity is infinity, no matter how small...
663     x)
664    ((float-nan-p x)
665     (when (and (float-trapping-nan-p x)
666                (sb!vm:current-float-trap :invalid))
667       (error 'floating-point-invalid-operation :operation 'scale-float
668              :operands (list x exp)))
669     x)
670    (t
671     (when (sb!vm:current-float-trap :overflow)
672       (error 'floating-point-overflow :operation 'scale-float
673              :operands (list x exp)))
674     (when (sb!vm:current-float-trap :inexact)
675       (error 'floating-point-inexact :operation 'scale-float
676              :operands (list x exp)))
677     (* (float-sign x)
678        (etypecase x
679          (single-float single-float-positive-infinity)
680          (double-float double-float-positive-infinity))))))
681
682 ;;; Scale a single or double float, calling the correct over/underflow
683 ;;; functions.
684 (defun scale-single-float (x exp)
685   (declare (single-float x) (fixnum exp))
686   (let* ((bits (single-float-bits x))
687          (old-exp (ldb sb!vm:single-float-exponent-byte bits))
688          (new-exp (+ old-exp exp)))
689     (cond
690      ((zerop x) x)
691      ((or (< old-exp sb!vm:single-float-normal-exponent-min)
692           (< new-exp sb!vm:single-float-normal-exponent-min))
693       (scale-float-maybe-underflow x exp))
694      ((or (> old-exp sb!vm:single-float-normal-exponent-max)
695           (> new-exp sb!vm:single-float-normal-exponent-max))
696       (scale-float-maybe-overflow x exp))
697      (t
698       (make-single-float (dpb new-exp
699                               sb!vm:single-float-exponent-byte
700                               bits))))))
701 (defun scale-double-float (x exp)
702   (declare (double-float x) (fixnum exp))
703   (let* ((hi (double-float-high-bits x))
704          (lo (double-float-low-bits x))
705          (old-exp (ldb sb!vm:double-float-exponent-byte hi))
706          (new-exp (+ old-exp exp)))
707     (cond
708      ((zerop x) x)
709      ((or (< old-exp sb!vm:double-float-normal-exponent-min)
710           (< new-exp sb!vm:double-float-normal-exponent-min))
711       (scale-float-maybe-underflow x exp))
712      ((or (> old-exp sb!vm:double-float-normal-exponent-max)
713           (> new-exp sb!vm:double-float-normal-exponent-max))
714       (scale-float-maybe-overflow x exp))
715      (t
716       (make-double-float (dpb new-exp sb!vm:double-float-exponent-byte hi)
717                          lo)))))
718
719 #!+(and x86 long-float)
720 (defun scale-long-float (x exp)
721   (declare (long-float x) (fixnum exp))
722   (scale-float x exp))
723
724 ;;; Dispatch to the correct type-specific scale-float function.
725 (defun scale-float (f ex)
726   #!+sb-doc
727   "Return the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
728   of precision or overflow."
729   (number-dispatch ((f float))
730     ((single-float)
731      (scale-single-float f ex))
732     ((double-float)
733      (scale-double-float f ex))
734     #!+long-float
735     ((long-float)
736      (scale-long-float f ex))))
737 \f
738 ;;;; converting to/from floats
739
740 (defun float (number &optional (other () otherp))
741   #!+sb-doc
742   "Converts any REAL to a float. If OTHER is not provided, it returns a
743   SINGLE-FLOAT if NUMBER is not already a FLOAT. If OTHER is provided, the
744   result is the same float format as OTHER."
745   (if otherp
746       (number-dispatch ((number real) (other float))
747         (((foreach rational single-float double-float #!+long-float long-float)
748           (foreach single-float double-float #!+long-float long-float))
749          (coerce number '(dispatch-type other))))
750       (if (floatp number)
751           number
752           (coerce number 'single-float))))
753
754 (macrolet ((frob (name type)
755              `(defun ,name (x)
756                 (number-dispatch ((x real))
757                   (((foreach single-float double-float #!+long-float long-float
758                              fixnum))
759                    (coerce x ',type))
760                   ((bignum)
761                    (bignum-to-float x ',type))
762                   ((ratio)
763                    (float-ratio x ',type))))))
764   (frob %single-float single-float)
765   (frob %double-float double-float)
766   #!+long-float
767   (frob %long-float long-float))
768
769 ;;; Convert a ratio to a float. We avoid any rounding error by doing an
770 ;;; integer division. Accuracy is important to preserve read/print
771 ;;; consistency, since this is ultimately how the reader reads a float. We
772 ;;; scale the numerator by a power of two until the division results in the
773 ;;; desired number of fraction bits, then do round-to-nearest.
774 (defun float-ratio (x format)
775   (let* ((signed-num (numerator x))
776          (plusp (plusp signed-num))
777          (num (if plusp signed-num (- signed-num)))
778          (den (denominator x))
779          (digits (float-format-digits format))
780          (scale 0))
781     (declare (fixnum digits scale))
782     ;; Strip any trailing zeros from the denominator and move it into the scale
783     ;; factor (to minimize the size of the operands.)
784     (let ((den-twos (1- (integer-length (logxor den (1- den))))))
785       (declare (fixnum den-twos))
786       (decf scale den-twos)
787       (setq den (ash den (- den-twos))))
788     ;; Guess how much we need to scale by from the magnitudes of the numerator
789     ;; and denominator. We want one extra bit for a guard bit.
790     (let* ((num-len (integer-length num))
791            (den-len (integer-length den))
792            (delta (- den-len num-len))
793            (shift (1+ (the fixnum (+ delta digits))))
794            (shifted-num (ash num shift)))
795       (declare (fixnum delta shift))
796       (decf scale delta)
797       (labels ((float-and-scale (bits)
798                  (let* ((bits (ash bits -1))
799                         (len (integer-length bits)))
800                    (cond ((> len digits)
801                           (aver (= len (the fixnum (1+ digits))))
802                           (scale-float (floatit (ash bits -1)) (1+ scale)))
803                          (t
804                           (scale-float (floatit bits) scale)))))
805                (floatit (bits)
806                  (let ((sign (if plusp 0 1)))
807                    (case format
808                      (single-float
809                       (single-from-bits sign sb!vm:single-float-bias bits))
810                      (double-float
811                       (double-from-bits sign sb!vm:double-float-bias bits))
812                      #!+long-float
813                      (long-float
814                       (long-from-bits sign sb!vm:long-float-bias bits))))))
815         (loop
816           (multiple-value-bind (fraction-and-guard rem)
817               (truncate shifted-num den)
818             (let ((extra (- (integer-length fraction-and-guard) digits)))
819               (declare (fixnum extra))
820               (cond ((/= extra 1)
821                      (aver (> extra 1)))
822                     ((oddp fraction-and-guard)
823                      (return
824                       (if (zerop rem)
825                           (float-and-scale
826                            (if (zerop (logand fraction-and-guard 2))
827                                fraction-and-guard
828                                (1+ fraction-and-guard)))
829                           (float-and-scale (1+ fraction-and-guard)))))
830                     (t
831                      (return (float-and-scale fraction-and-guard)))))
832             (setq shifted-num (ash shifted-num -1))
833             (incf scale)))))))
834
835 #|
836 These might be useful if we ever have a machine without float/integer
837 conversion hardware. For now, we'll use special ops that
838 uninterruptibly frob the rounding modes & do ieee round-to-integer.
839
840 ;;; The compiler compiles a call to this when we are doing %UNARY-TRUNCATE
841 ;;; and the result is known to be a fixnum. We can avoid some generic
842 ;;; arithmetic in this case.
843 (defun %unary-truncate-single-float/fixnum (x)
844   (declare (single-float x) (values fixnum))
845   (locally (declare (optimize (speed 3) (safety 0)))
846     (let* ((bits (single-float-bits x))
847            (exp (ldb sb!vm:single-float-exponent-byte bits))
848            (frac (logior (ldb sb!vm:single-float-significand-byte bits)
849                          sb!vm:single-float-hidden-bit))
850            (shift (- exp sb!vm:single-float-digits sb!vm:single-float-bias)))
851       (when (> exp sb!vm:single-float-normal-exponent-max)
852         (error 'floating-point-invalid-operation :operator 'truncate
853                :operands (list x)))
854       (if (<= shift (- sb!vm:single-float-digits))
855           0
856           (let ((res (ash frac shift)))
857             (declare (type (unsigned-byte 31) res))
858             (if (minusp bits)
859                 (- res)
860                 res))))))
861
862 ;;; Double-float version of this operation (see above single op).
863 (defun %unary-truncate-double-float/fixnum (x)
864   (declare (double-float x) (values fixnum))
865   (locally (declare (optimize (speed 3) (safety 0)))
866     (let* ((hi-bits (double-float-high-bits x))
867            (exp (ldb sb!vm:double-float-exponent-byte hi-bits))
868            (frac (logior (ldb sb!vm:double-float-significand-byte hi-bits)
869                          sb!vm:double-float-hidden-bit))
870            (shift (- exp (- sb!vm:double-float-digits sb!vm:n-word-bits)
871                      sb!vm:double-float-bias)))
872       (when (> exp sb!vm:double-float-normal-exponent-max)
873         (error 'floating-point-invalid-operation :operator 'truncate
874                :operands (list x)))
875       (if (<= shift (- sb!vm:n-word-bits sb!vm:double-float-digits))
876           0
877           (let* ((res-hi (ash frac shift))
878                  (res (if (plusp shift)
879                           (logior res-hi
880                                   (the fixnum
881                                        (ash (double-float-low-bits x)
882                                             (- shift sb!vm:n-word-bits))))
883                           res-hi)))
884             (declare (type (unsigned-byte 31) res-hi res))
885             (if (minusp hi-bits)
886                 (- res)
887                 res))))))
888 |#
889
890 ;;; This function is called when we are doing a truncate without any funky
891 ;;; divisor, i.e. converting a float or ratio to an integer. Note that we do
892 ;;; *not* return the second value of truncate, so it must be computed by the
893 ;;; caller if needed.
894 ;;;
895 ;;; In the float case, we pick off small arguments so that compiler can use
896 ;;; special-case operations. We use an exclusive test, since (due to round-off
897 ;;; error), (float most-positive-fixnum) may be greater than
898 ;;; most-positive-fixnum.
899 (defun %unary-truncate (number)
900   (number-dispatch ((number real))
901     ((integer) number)
902     ((ratio) (values (truncate (numerator number) (denominator number))))
903     (((foreach single-float double-float #!+long-float long-float))
904      (if (< (float most-negative-fixnum number)
905             number
906             (float most-positive-fixnum number))
907          (truly-the fixnum (%unary-truncate number))
908          (multiple-value-bind (bits exp) (integer-decode-float number)
909            (let ((res (ash bits exp)))
910              (if (minusp number)
911                  (- res)
912                  res)))))))
913
914 ;;; Similar to %UNARY-TRUNCATE, but rounds to the nearest integer. If we
915 ;;; can't use the round primitive, then we do our own round-to-nearest on the
916 ;;; result of i-d-f. [Note that this rounding will really only happen with
917 ;;; double floats, since the whole single-float fraction will fit in a fixnum,
918 ;;; so all single-floats larger than most-positive-fixnum can be precisely
919 ;;; represented by an integer.]
920 (defun %unary-round (number)
921   (number-dispatch ((number real))
922     ((integer) number)
923     ((ratio) (values (round (numerator number) (denominator number))))
924     (((foreach single-float double-float #!+long-float long-float))
925      (if (< (float most-negative-fixnum number)
926             number
927             (float most-positive-fixnum number))
928          (truly-the fixnum (%unary-round number))
929          (multiple-value-bind (bits exp) (integer-decode-float number)
930            (let* ((shifted (ash bits exp))
931                   (rounded (if (and (minusp exp)
932                                     (oddp shifted)
933                                     (eql (logand bits
934                                                  (lognot (ash -1 (- exp))))
935                                          (ash 1 (- -1 exp))))
936                                (1+ shifted)
937                                shifted)))
938              (if (minusp number)
939                  (- rounded)
940                  rounded)))))))
941
942 (defun rational (x)
943   #!+sb-doc
944   "RATIONAL produces a rational number for any real numeric argument. This is
945   more efficient than RATIONALIZE, but it assumes that floating-point is
946   completely accurate, giving a result that isn't as pretty."
947   (number-dispatch ((x real))
948     (((foreach single-float double-float #!+long-float long-float))
949      (multiple-value-bind (bits exp) (integer-decode-float x)
950        (if (eql bits 0)
951            0
952            (let* ((int (if (minusp x) (- bits) bits))
953                   (digits (float-digits x))
954                   (ex (+ exp digits)))
955              (if (minusp ex)
956                  (integer-/-integer int (ash 1 (+ digits (- ex))))
957                  (integer-/-integer (ash int ex) (ash 1 digits)))))))
958     ((rational) x)))
959
960 (defun rationalize (x)
961   #!+sb-doc
962   "Converts any REAL to a RATIONAL. Floats are converted to a simple rational
963   representation exploiting the assumption that floats are only accurate to
964   their precision. RATIONALIZE (and also RATIONAL) preserve the invariant:
965       (= x (float (rationalize x) x))"
966   (number-dispatch ((x real))
967     (((foreach single-float double-float #!+long-float long-float))
968      ;; Thanks to Kim Fateman, who stole this function rationalize-float from
969      ;; macsyma's rational. Macsyma'a rationalize was written by the legendary
970      ;; Gosper (rwg). Guy Steele said about Gosper, "He has been called the
971      ;; only living 17th century mathematician and is also the best pdp-10
972      ;; hacker I know." So, if you can understand or debug this code you win
973      ;; big.
974      (cond ((minusp x) (- (rationalize (- x))))
975            ((zerop x) 0)
976            (t
977             (let ((eps (etypecase x
978                            (single-float single-float-epsilon)
979                            (double-float double-float-epsilon)
980                            #!+long-float
981                            (long-float long-float-epsilon)))
982                   (y ())
983                   (a ()))
984               (do ((xx x (setq y (/ (float 1.0 x) (- xx (float a x)))))
985                    (num (setq a (truncate x))
986                         (+ (* (setq a (truncate y)) num) onum))
987                    (den 1 (+ (* a den) oden))
988                    (onum 1 num)
989                    (oden 0 den))
990                   ((and (not (zerop den))
991                         (not (> (abs (/ (- x (/ (float num x)
992                                                 (float den x)))
993                                         x))
994                                 eps)))
995                    (integer-/-integer num den))
996                 (declare ((dispatch-type x) xx)))))))
997     ((rational) x)))