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