0.8.8.21:
[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   (not (zerop (ldb sb!vm:single-float-significand-byte bits)))
79   (or (not (zerop (ldb sb!vm:double-float-significand-byte hi)))
80       (not (zerop lo)))
81   #!+(and long-float x86)
82   (or (not (zerop (ldb sb!vm:long-float-significand-byte hi)))
83       (not (zerop lo))))
84
85 (!define-float-dispatching-function float-trapping-nan-p
86   "Return true if the float X is a trapping NaN (Not a Number)."
87   (zerop (logand (ldb sb!vm:single-float-significand-byte bits)
88                  sb!vm:single-float-trapping-nan-bit))
89   (zerop (logand (ldb sb!vm:double-float-significand-byte hi)
90                  sb!vm:double-float-trapping-nan-bit))
91   #!+(and long-float x86)
92   (zerop (logand (ldb sb!vm:long-float-significand-byte hi)
93                  sb!vm:long-float-trapping-nan-bit)))
94
95 ;;; If denormalized, use a subfunction from INTEGER-DECODE-FLOAT to find the
96 ;;; actual exponent (and hence how denormalized it is), otherwise we just
97 ;;; return the number of digits or 0.
98 #!-sb-fluid (declaim (maybe-inline float-precision))
99 (defun float-precision (f)
100   #!+sb-doc
101   "Return a non-negative number of significant digits in its float argument.
102   Will be less than FLOAT-DIGITS if denormalized or zero."
103   (macrolet ((frob (digits bias decode)
104                `(cond ((zerop f) 0)
105                       ((float-denormalized-p f)
106                        (multiple-value-bind (ignore exp) (,decode f)
107                          (declare (ignore ignore))
108                          (truly-the fixnum
109                                     (+ ,digits (1- ,digits) ,bias exp))))
110                       (t
111                        ,digits))))
112     (number-dispatch ((f float))
113       ((single-float)
114        (frob sb!vm:single-float-digits sb!vm:single-float-bias
115          integer-decode-single-denorm))
116       ((double-float)
117        (frob sb!vm:double-float-digits sb!vm:double-float-bias
118          integer-decode-double-denorm))
119       #!+long-float
120       ((long-float)
121        (frob sb!vm:long-float-digits sb!vm:long-float-bias
122          integer-decode-long-denorm)))))
123
124 (defun float-sign (float1 &optional (float2 (float 1 float1)))
125   #!+sb-doc
126   "Return a floating-point number that has the same sign as
127    FLOAT1 and, if FLOAT2 is given, has the same absolute value
128    as FLOAT2."
129   (declare (float float1 float2))
130   (* (if (etypecase float1
131            (single-float (minusp (single-float-bits float1)))
132            (double-float (minusp (double-float-high-bits float1)))
133            #!+long-float
134            (long-float (minusp (long-float-exp-bits float1))))
135          (float -1 float1)
136          (float 1 float1))
137      (abs float2)))
138
139 (defun float-format-digits (format)
140   (ecase format
141     ((short-float single-float) sb!vm:single-float-digits)
142     ((double-float #!-long-float long-float) sb!vm:double-float-digits)
143     #!+long-float
144     (long-float sb!vm:long-float-digits)))
145
146 #!-sb-fluid (declaim (inline float-digits float-radix))
147
148 (defun float-digits (f)
149   (number-dispatch ((f float))
150     ((single-float) sb!vm:single-float-digits)
151     ((double-float) sb!vm:double-float-digits)
152     #!+long-float
153     ((long-float) sb!vm:long-float-digits)))
154
155 (defun float-radix (x)
156   #!+sb-doc
157   "Return (as an integer) the radix b of its floating-point argument."
158   (declare (ignore x))
159   2)
160 \f
161 ;;;; INTEGER-DECODE-FLOAT and DECODE-FLOAT
162
163 #!-sb-fluid
164 (declaim (maybe-inline integer-decode-single-float
165                        integer-decode-double-float))
166
167 ;;; Handle the denormalized case of INTEGER-DECODE-FLOAT for SINGLE-FLOAT.
168 (defun integer-decode-single-denorm (x)
169   (declare (type single-float x))
170   (let* ((bits (single-float-bits (abs x)))
171          (sig (ash (ldb sb!vm:single-float-significand-byte bits) 1))
172          (extra-bias 0))
173     (declare (type (unsigned-byte 24) sig)
174              (type (integer 0 23) extra-bias))
175     (loop
176       (unless (zerop (logand sig sb!vm:single-float-hidden-bit))
177         (return))
178       (setq sig (ash sig 1))
179       (incf extra-bias))
180     (values sig
181             (- (- sb!vm:single-float-bias)
182                sb!vm:single-float-digits
183                extra-bias)
184             (if (minusp (float-sign x)) -1 1))))
185
186 ;;; Handle the single-float case of INTEGER-DECODE-FLOAT. If an infinity or
187 ;;; NaN, error. If a denorm, call i-d-s-DENORM to handle it.
188 (defun integer-decode-single-float (x)
189   (declare (single-float x))
190   (let* ((bits (single-float-bits (abs x)))
191          (exp (ldb sb!vm:single-float-exponent-byte bits))
192          (sig (ldb sb!vm:single-float-significand-byte bits))
193          (sign (if (minusp (float-sign x)) -1 1))
194          (biased (- exp sb!vm:single-float-bias sb!vm:single-float-digits)))
195     (declare (fixnum biased))
196     (unless (<= exp sb!vm:single-float-normal-exponent-max)
197       (error "can't decode NaN or infinity: ~S" x))
198     (cond ((and (zerop exp) (zerop sig))
199            (values 0 biased sign))
200           ((< exp sb!vm:single-float-normal-exponent-min)
201            (integer-decode-single-denorm x))
202           (t
203            (values (logior sig sb!vm:single-float-hidden-bit) biased sign)))))
204
205 ;;; like INTEGER-DECODE-SINGLE-DENORM, only doubly so
206 (defun integer-decode-double-denorm (x)
207   (declare (type double-float x))
208   (let* ((high-bits (double-float-high-bits (abs x)))
209          (sig-high (ldb sb!vm:double-float-significand-byte high-bits))
210          (low-bits (double-float-low-bits x))
211          (sign (if (minusp (float-sign x)) -1 1))
212          (biased (- (- sb!vm:double-float-bias) sb!vm:double-float-digits)))
213     (if (zerop sig-high)
214         (let ((sig low-bits)
215               (extra-bias (- sb!vm:double-float-digits 33))
216               (bit (ash 1 31)))
217           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
218           (loop
219             (unless (zerop (logand sig bit)) (return))
220             (setq sig (ash sig 1))
221             (incf extra-bias))
222           (values (ash sig (- sb!vm:double-float-digits 32))
223                   (truly-the fixnum (- biased extra-bias))
224                   sign))
225         (let ((sig (ash sig-high 1))
226               (extra-bias 0))
227           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
228           (loop
229             (unless (zerop (logand sig sb!vm:double-float-hidden-bit))
230               (return))
231             (setq sig (ash sig 1))
232             (incf extra-bias))
233           (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
234                   (truly-the fixnum (- biased extra-bias))
235                   sign)))))
236
237 ;;; like INTEGER-DECODE-SINGLE-FLOAT, only doubly so
238 (defun integer-decode-double-float (x)
239   (declare (double-float x))
240   (let* ((abs (abs x))
241          (hi (double-float-high-bits abs))
242          (lo (double-float-low-bits abs))
243          (exp (ldb sb!vm:double-float-exponent-byte hi))
244          (sig (ldb sb!vm:double-float-significand-byte hi))
245          (sign (if (minusp (float-sign x)) -1 1))
246          (biased (- exp sb!vm:double-float-bias sb!vm:double-float-digits)))
247     (declare (fixnum biased))
248     (unless (<= exp sb!vm:double-float-normal-exponent-max)
249       (error "Can't decode NaN or infinity: ~S." x))
250     (cond ((and (zerop exp) (zerop sig) (zerop lo))
251            (values 0 biased sign))
252           ((< exp sb!vm:double-float-normal-exponent-min)
253            (integer-decode-double-denorm x))
254           (t
255            (values
256             (logior (ash (logior (ldb sb!vm:double-float-significand-byte hi)
257                                  sb!vm:double-float-hidden-bit)
258                          32)
259                     lo)
260             biased sign)))))
261
262 #!+(and long-float x86)
263 (defun integer-decode-long-denorm (x)
264   (declare (type long-float x))
265   (let* ((high-bits (long-float-high-bits (abs x)))
266          (sig-high (ldb sb!vm:long-float-significand-byte high-bits))
267          (low-bits (long-float-low-bits x))
268          (sign (if (minusp (float-sign x)) -1 1))
269          (biased (- (- sb!vm:long-float-bias) sb!vm:long-float-digits)))
270     (if (zerop sig-high)
271         (let ((sig low-bits)
272               (extra-bias (- sb!vm:long-float-digits 33))
273               (bit (ash 1 31)))
274           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
275           (loop
276             (unless (zerop (logand sig bit)) (return))
277             (setq sig (ash sig 1))
278             (incf extra-bias))
279           (values (ash sig (- sb!vm:long-float-digits 32))
280                   (truly-the fixnum (- biased extra-bias))
281                   sign))
282         (let ((sig (ash sig-high 1))
283               (extra-bias 0))
284           (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
285           (loop
286             (unless (zerop (logand sig sb!vm:long-float-hidden-bit))
287               (return))
288             (setq sig (ash sig 1))
289             (incf extra-bias))
290           (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
291                   (truly-the fixnum (- biased extra-bias))
292                   sign)))))
293
294 #!+(and long-float x86)
295 (defun integer-decode-long-float (x)
296   (declare (long-float x))
297   (let* ((hi (long-float-high-bits x))
298          (lo (long-float-low-bits x))
299          (exp-bits (long-float-exp-bits x))
300          (exp (ldb sb!vm:long-float-exponent-byte exp-bits))
301          (sign (if (minusp exp-bits) -1 1))
302          (biased (- exp sb!vm:long-float-bias sb!vm:long-float-digits)))
303     (declare (fixnum biased))
304     (unless (<= exp sb!vm:long-float-normal-exponent-max)
305       (error "can't decode NaN or infinity: ~S" x))
306     (cond ((and (zerop exp) (zerop hi) (zerop lo))
307            (values 0 biased sign))
308           ((< exp sb!vm:long-float-normal-exponent-min)
309            (integer-decode-long-denorm x))
310           (t
311            (values (logior (ash hi 32) lo) biased sign)))))
312
313 ;;; Dispatch to the correct type-specific i-d-f function.
314 (defun integer-decode-float (x)
315   #!+sb-doc
316   "Return three values:
317    1) an integer representation of the significand.
318    2) the exponent for the power of 2 that the significand must be multiplied
319       by to get the actual value. This differs from the DECODE-FLOAT exponent
320       by FLOAT-DIGITS, since the significand has been scaled to have all its
321       digits before the radix point.
322    3) -1 or 1 (i.e. the sign of the argument.)"
323   (number-dispatch ((x float))
324     ((single-float)
325      (integer-decode-single-float x))
326     ((double-float)
327      (integer-decode-double-float x))
328     #!+long-float
329     ((long-float)
330      (integer-decode-long-float x))))
331
332 #!-sb-fluid (declaim (maybe-inline decode-single-float decode-double-float))
333
334 ;;; Handle the denormalized case of DECODE-SINGLE-FLOAT. We call
335 ;;; INTEGER-DECODE-SINGLE-DENORM and then make the result into a float.
336 (defun decode-single-denorm (x)
337   (declare (type single-float x))
338   (multiple-value-bind (sig exp sign) (integer-decode-single-denorm x)
339     (values (make-single-float
340              (dpb sig sb!vm:single-float-significand-byte
341                   (dpb sb!vm:single-float-bias
342                        sb!vm:single-float-exponent-byte
343                        0)))
344             (truly-the fixnum (+ exp sb!vm:single-float-digits))
345             (float sign x))))
346
347 ;;; Handle the single-float case of DECODE-FLOAT. If an infinity or NaN,
348 ;;; error. If a denorm, call d-s-DENORM to handle it.
349 (defun decode-single-float (x)
350   (declare (single-float x))
351   (let* ((bits (single-float-bits (abs x)))
352          (exp (ldb sb!vm:single-float-exponent-byte bits))
353          (sign (float-sign x))
354          (biased (truly-the single-float-exponent
355                             (- exp sb!vm:single-float-bias))))
356     (unless (<= exp sb!vm:single-float-normal-exponent-max)
357       (error "can't decode NaN or infinity: ~S" x))
358     (cond ((zerop x)
359            (values 0.0f0 biased sign))
360           ((< exp sb!vm:single-float-normal-exponent-min)
361            (decode-single-denorm x))
362           (t
363            (values (make-single-float
364                     (dpb sb!vm:single-float-bias
365                          sb!vm:single-float-exponent-byte
366                          bits))
367                    biased sign)))))
368
369 ;;; like DECODE-SINGLE-DENORM, only doubly so
370 (defun decode-double-denorm (x)
371   (declare (double-float x))
372   (multiple-value-bind (sig exp sign) (integer-decode-double-denorm x)
373     (values (make-double-float
374              (dpb (logand (ash sig -32) (lognot sb!vm:double-float-hidden-bit))
375                   sb!vm:double-float-significand-byte
376                   (dpb sb!vm:double-float-bias
377                        sb!vm:double-float-exponent-byte 0))
378              (ldb (byte 32 0) sig))
379             (truly-the fixnum (+ exp sb!vm:double-float-digits))
380             (float sign x))))
381
382 ;;; like DECODE-SINGLE-FLOAT, only doubly so
383 (defun decode-double-float (x)
384   (declare (double-float x))
385   (let* ((abs (abs x))
386          (hi (double-float-high-bits abs))
387          (lo (double-float-low-bits abs))
388          (exp (ldb sb!vm:double-float-exponent-byte hi))
389          (sign (float-sign x))
390          (biased (truly-the double-float-exponent
391                             (- exp sb!vm:double-float-bias))))
392     (unless (<= exp sb!vm:double-float-normal-exponent-max)
393       (error "can't decode NaN or infinity: ~S" x))
394     (cond ((zerop x)
395            (values 0.0d0 biased sign))
396           ((< exp sb!vm:double-float-normal-exponent-min)
397            (decode-double-denorm x))
398           (t
399            (values (make-double-float
400                     (dpb sb!vm:double-float-bias
401                          sb!vm:double-float-exponent-byte hi)
402                     lo)
403                    biased sign)))))
404
405 #!+(and long-float x86)
406 (defun decode-long-denorm (x)
407   (declare (long-float x))
408   (multiple-value-bind (sig exp sign) (integer-decode-long-denorm x)
409     (values (make-long-float sb!vm:long-float-bias (ash sig -32)
410                              (ldb (byte 32 0) sig))
411             (truly-the fixnum (+ exp sb!vm:long-float-digits))
412             (float sign x))))
413
414 #!+(and long-float x86)
415 (defun decode-long-float (x)
416   (declare (long-float x))
417   (let* ((hi (long-float-high-bits x))
418          (lo (long-float-low-bits x))
419          (exp-bits (long-float-exp-bits x))
420          (exp (ldb sb!vm:long-float-exponent-byte exp-bits))
421          (sign (if (minusp exp-bits) -1l0 1l0))
422          (biased (truly-the long-float-exponent
423                             (- exp sb!vm:long-float-bias))))
424     (unless (<= exp sb!vm:long-float-normal-exponent-max)
425       (error "can't decode NaN or infinity: ~S" x))
426     (cond ((zerop x)
427            (values 0.0l0 biased sign))
428           ((< exp sb!vm:long-float-normal-exponent-min)
429            (decode-long-denorm x))
430           (t
431            (values (make-long-float
432                     (dpb sb!vm:long-float-bias sb!vm:long-float-exponent-byte
433                          exp-bits)
434                     hi
435                     lo)
436                    biased sign)))))
437
438 ;;; Dispatch to the appropriate type-specific function.
439 (defun decode-float (f)
440   #!+sb-doc
441   "Return three values:
442    1) a floating-point number representing the significand. This is always
443       between 0.5 (inclusive) and 1.0 (exclusive).
444    2) an integer representing the exponent.
445    3) -1.0 or 1.0 (i.e. the sign of the argument.)"
446   (number-dispatch ((f float))
447     ((single-float)
448      (decode-single-float f))
449     ((double-float)
450      (decode-double-float f))
451     #!+long-float
452     ((long-float)
453      (decode-long-float f))))
454 \f
455 ;;;; SCALE-FLOAT
456
457 #!-sb-fluid (declaim (maybe-inline scale-single-float scale-double-float))
458
459 ;;; Handle float scaling where the X is denormalized or the result is
460 ;;; denormalized or underflows to 0.
461 (defun scale-float-maybe-underflow (x exp)
462   (multiple-value-bind (sig old-exp) (integer-decode-float x)
463     (let* ((digits (float-digits x))
464            (new-exp (+ exp old-exp digits
465                        (etypecase x
466                          (single-float sb!vm:single-float-bias)
467                          (double-float sb!vm:double-float-bias))))
468            (sign (if (minusp (float-sign x)) 1 0)))
469       (cond
470        ((< new-exp
471            (etypecase x
472              (single-float sb!vm:single-float-normal-exponent-min)
473              (double-float sb!vm:double-float-normal-exponent-min)))
474         (when (sb!vm:current-float-trap :inexact)
475           (error 'floating-point-inexact :operation 'scale-float
476                  :operands (list x exp)))
477         (when (sb!vm:current-float-trap :underflow)
478           (error 'floating-point-underflow :operation 'scale-float
479                  :operands (list x exp)))
480         (let ((shift (1- new-exp)))
481           (if (< shift (- (1- digits)))
482               (float-sign x 0.0)
483               (etypecase x
484                 (single-float (single-from-bits sign 0 (ash sig shift)))
485                 (double-float (double-from-bits sign 0 (ash sig shift)))))))
486        (t
487         (etypecase x
488           (single-float (single-from-bits sign new-exp sig))
489           (double-float (double-from-bits sign new-exp sig))))))))
490
491 ;;; Called when scaling a float overflows, or the original float was a
492 ;;; NaN or infinity. If overflow errors are trapped, then error,
493 ;;; otherwise return the appropriate infinity. If a NaN, signal or not
494 ;;; as appropriate.
495 (defun scale-float-maybe-overflow (x exp)
496   (cond
497    ((float-infinity-p x)
498     ;; Infinity is infinity, no matter how small...
499     x)
500    ((float-nan-p x)
501     (when (and (float-trapping-nan-p x)
502                (sb!vm:current-float-trap :invalid))
503       (error 'floating-point-invalid-operation :operation 'scale-float
504              :operands (list x exp)))
505     x)
506    (t
507     (when (sb!vm:current-float-trap :overflow)
508       (error 'floating-point-overflow :operation 'scale-float
509              :operands (list x exp)))
510     (when (sb!vm:current-float-trap :inexact)
511       (error 'floating-point-inexact :operation 'scale-float
512              :operands (list x exp)))
513     (* (float-sign x)
514        (etypecase x
515          (single-float
516           ;; SINGLE-FLOAT-POSITIVE-INFINITY
517           (single-from-bits 0 (1+ sb!vm:single-float-normal-exponent-max) 0))
518          (double-float
519           ;; DOUBLE-FLOAT-POSITIVE-INFINITY
520           (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0)))))))
521
522 ;;; Scale a single or double float, calling the correct over/underflow
523 ;;; functions.
524 (defun scale-single-float (x exp)
525   (declare (single-float x) (fixnum exp))
526   (let* ((bits (single-float-bits x))
527          (old-exp (ldb sb!vm:single-float-exponent-byte bits))
528          (new-exp (+ old-exp exp)))
529     (cond
530      ((zerop x) x)
531      ((or (< old-exp sb!vm:single-float-normal-exponent-min)
532           (< new-exp sb!vm:single-float-normal-exponent-min))
533       (scale-float-maybe-underflow x exp))
534      ((or (> old-exp sb!vm:single-float-normal-exponent-max)
535           (> new-exp sb!vm:single-float-normal-exponent-max))
536       (scale-float-maybe-overflow x exp))
537      (t
538       (make-single-float (dpb new-exp
539                               sb!vm:single-float-exponent-byte
540                               bits))))))
541 (defun scale-double-float (x exp)
542   (declare (double-float x) (fixnum exp))
543   (let* ((hi (double-float-high-bits x))
544          (lo (double-float-low-bits x))
545          (old-exp (ldb sb!vm:double-float-exponent-byte hi))
546          (new-exp (+ old-exp exp)))
547     (cond
548      ((zerop x) x)
549      ((or (< old-exp sb!vm:double-float-normal-exponent-min)
550           (< new-exp sb!vm:double-float-normal-exponent-min))
551       (scale-float-maybe-underflow x exp))
552      ((or (> old-exp sb!vm:double-float-normal-exponent-max)
553           (> new-exp sb!vm:double-float-normal-exponent-max))
554       (scale-float-maybe-overflow x exp))
555      (t
556       (make-double-float (dpb new-exp sb!vm:double-float-exponent-byte hi)
557                          lo)))))
558
559 #!+(and x86 long-float)
560 (defun scale-long-float (x exp)
561   (declare (long-float x) (fixnum exp))
562   (scale-float x exp))
563
564 ;;; Dispatch to the correct type-specific scale-float function.
565 (defun scale-float (f ex)
566   #!+sb-doc
567   "Return the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
568   of precision or overflow."
569   (number-dispatch ((f float))
570     ((single-float)
571      (scale-single-float f ex))
572     ((double-float)
573      (scale-double-float f ex))
574     #!+long-float
575     ((long-float)
576      (scale-long-float f ex))))
577 \f
578 ;;;; converting to/from floats
579
580 (defun float (number &optional (other () otherp))
581   #!+sb-doc
582   "Converts any REAL to a float. If OTHER is not provided, it returns a
583   SINGLE-FLOAT if NUMBER is not already a FLOAT. If OTHER is provided, the
584   result is the same float format as OTHER."
585   (if otherp
586       (number-dispatch ((number real) (other float))
587         (((foreach rational single-float double-float #!+long-float long-float)
588           (foreach single-float double-float #!+long-float long-float))
589          (coerce number '(dispatch-type other))))
590       (if (floatp number)
591           number
592           (coerce number 'single-float))))
593
594 (macrolet ((frob (name type)
595              `(defun ,name (x)
596                 (number-dispatch ((x real))
597                   (((foreach single-float double-float #!+long-float long-float
598                              fixnum))
599                    (coerce x ',type))
600                   ((bignum)
601                    (bignum-to-float x ',type))
602                   ((ratio)
603                    (float-ratio x ',type))))))
604   (frob %single-float single-float)
605   (frob %double-float double-float)
606   #!+long-float
607   (frob %long-float long-float))
608
609 ;;; Convert a ratio to a float. We avoid any rounding error by doing an
610 ;;; integer division. Accuracy is important to preserve read/print
611 ;;; consistency, since this is ultimately how the reader reads a float. We
612 ;;; scale the numerator by a power of two until the division results in the
613 ;;; desired number of fraction bits, then do round-to-nearest.
614 (defun float-ratio (x format)
615   (let* ((signed-num (numerator x))
616          (plusp (plusp signed-num))
617          (num (if plusp signed-num (- signed-num)))
618          (den (denominator x))
619          (digits (float-format-digits format))
620          (scale 0))
621     (declare (fixnum digits scale))
622     ;; Strip any trailing zeros from the denominator and move it into the scale
623     ;; factor (to minimize the size of the operands.)
624     (let ((den-twos (1- (integer-length (logxor den (1- den))))))
625       (declare (fixnum den-twos))
626       (decf scale den-twos)
627       (setq den (ash den (- den-twos))))
628     ;; Guess how much we need to scale by from the magnitudes of the numerator
629     ;; and denominator. We want one extra bit for a guard bit.
630     (let* ((num-len (integer-length num))
631            (den-len (integer-length den))
632            (delta (- den-len num-len))
633            (shift (1+ (the fixnum (+ delta digits))))
634            (shifted-num (ash num shift)))
635       (declare (fixnum delta shift))
636       (decf scale delta)
637       (labels ((float-and-scale (bits)
638                  (let* ((bits (ash bits -1))
639                         (len (integer-length bits)))
640                    (cond ((> len digits)
641                           (aver (= len (the fixnum (1+ digits))))
642                           (scale-float (floatit (ash bits -1)) (1+ scale)))
643                          (t
644                           (scale-float (floatit bits) scale)))))
645                (floatit (bits)
646                  (let ((sign (if plusp 0 1)))
647                    (case format
648                      (single-float
649                       (single-from-bits sign sb!vm:single-float-bias bits))
650                      (double-float
651                       (double-from-bits sign sb!vm:double-float-bias bits))
652                      #!+long-float
653                      (long-float
654                       (long-from-bits sign sb!vm:long-float-bias bits))))))
655         (loop
656           (multiple-value-bind (fraction-and-guard rem)
657               (truncate shifted-num den)
658             (let ((extra (- (integer-length fraction-and-guard) digits)))
659               (declare (fixnum extra))
660               (cond ((/= extra 1)
661                      (aver (> extra 1)))
662                     ((oddp fraction-and-guard)
663                      (return
664                       (if (zerop rem)
665                           (float-and-scale
666                            (if (zerop (logand fraction-and-guard 2))
667                                fraction-and-guard
668                                (1+ fraction-and-guard)))
669                           (float-and-scale (1+ fraction-and-guard)))))
670                     (t
671                      (return (float-and-scale fraction-and-guard)))))
672             (setq shifted-num (ash shifted-num -1))
673             (incf scale)))))))
674
675 #|
676 These might be useful if we ever have a machine without float/integer
677 conversion hardware. For now, we'll use special ops that
678 uninterruptibly frob the rounding modes & do ieee round-to-integer.
679
680 ;;; The compiler compiles a call to this when we are doing %UNARY-TRUNCATE
681 ;;; and the result is known to be a fixnum. We can avoid some generic
682 ;;; arithmetic in this case.
683 (defun %unary-truncate-single-float/fixnum (x)
684   (declare (single-float x) (values fixnum))
685   (locally (declare (optimize (speed 3) (safety 0)))
686     (let* ((bits (single-float-bits x))
687            (exp (ldb sb!vm:single-float-exponent-byte bits))
688            (frac (logior (ldb sb!vm:single-float-significand-byte bits)
689                          sb!vm:single-float-hidden-bit))
690            (shift (- exp sb!vm:single-float-digits sb!vm:single-float-bias)))
691       (when (> exp sb!vm:single-float-normal-exponent-max)
692         (error 'floating-point-invalid-operation :operator 'truncate
693                :operands (list x)))
694       (if (<= shift (- sb!vm:single-float-digits))
695           0
696           (let ((res (ash frac shift)))
697             (declare (type (unsigned-byte 31) res))
698             (if (minusp bits)
699                 (- res)
700                 res))))))
701
702 ;;; Double-float version of this operation (see above single op).
703 (defun %unary-truncate-double-float/fixnum (x)
704   (declare (double-float x) (values fixnum))
705   (locally (declare (optimize (speed 3) (safety 0)))
706     (let* ((hi-bits (double-float-high-bits x))
707            (exp (ldb sb!vm:double-float-exponent-byte hi-bits))
708            (frac (logior (ldb sb!vm:double-float-significand-byte hi-bits)
709                          sb!vm:double-float-hidden-bit))
710            (shift (- exp (- sb!vm:double-float-digits sb!vm:n-word-bits)
711                      sb!vm:double-float-bias)))
712       (when (> exp sb!vm:double-float-normal-exponent-max)
713         (error 'floating-point-invalid-operation :operator 'truncate
714                :operands (list x)))
715       (if (<= shift (- sb!vm:n-word-bits sb!vm:double-float-digits))
716           0
717           (let* ((res-hi (ash frac shift))
718                  (res (if (plusp shift)
719                           (logior res-hi
720                                   (the fixnum
721                                        (ash (double-float-low-bits x)
722                                             (- shift sb!vm:n-word-bits))))
723                           res-hi)))
724             (declare (type (unsigned-byte 31) res-hi res))
725             (if (minusp hi-bits)
726                 (- res)
727                 res))))))
728 |#
729
730 ;;; This function is called when we are doing a truncate without any funky
731 ;;; divisor, i.e. converting a float or ratio to an integer. Note that we do
732 ;;; *not* return the second value of truncate, so it must be computed by the
733 ;;; caller if needed.
734 ;;;
735 ;;; In the float case, we pick off small arguments so that compiler can use
736 ;;; special-case operations. We use an exclusive test, since (due to round-off
737 ;;; error), (float most-positive-fixnum) may be greater than
738 ;;; most-positive-fixnum.
739 (defun %unary-truncate (number)
740   (number-dispatch ((number real))
741     ((integer) number)
742     ((ratio) (values (truncate (numerator number) (denominator number))))
743     (((foreach single-float double-float #!+long-float long-float))
744      (if (< (float most-negative-fixnum number)
745             number
746             (float most-positive-fixnum number))
747          (truly-the fixnum (%unary-truncate number))
748          (multiple-value-bind (bits exp) (integer-decode-float number)
749            (let ((res (ash bits exp)))
750              (if (minusp number)
751                  (- res)
752                  res)))))))
753
754 ;;; Similar to %UNARY-TRUNCATE, but rounds to the nearest integer. If we
755 ;;; can't use the round primitive, then we do our own round-to-nearest on the
756 ;;; result of i-d-f. [Note that this rounding will really only happen with
757 ;;; double floats, since the whole single-float fraction will fit in a fixnum,
758 ;;; so all single-floats larger than most-positive-fixnum can be precisely
759 ;;; represented by an integer.]
760 (defun %unary-round (number)
761   (number-dispatch ((number real))
762     ((integer) number)
763     ((ratio) (values (round (numerator number) (denominator number))))
764     (((foreach single-float double-float #!+long-float long-float))
765      (if (< (float most-negative-fixnum number)
766             number
767             (float most-positive-fixnum number))
768          (truly-the fixnum (%unary-round number))
769          (multiple-value-bind (bits exp) (integer-decode-float number)
770            (let* ((shifted (ash bits exp))
771                   (rounded (if (and (minusp exp)
772                                     (oddp shifted)
773                                     (eql (logand bits
774                                                  (lognot (ash -1 (- exp))))
775                                          (ash 1 (- -1 exp))))
776                                (1+ shifted)
777                                shifted)))
778              (if (minusp number)
779                  (- rounded)
780                  rounded)))))))
781
782 (defun rational (x)
783   #!+sb-doc
784   "RATIONAL produces a rational number for any real numeric argument. This is
785   more efficient than RATIONALIZE, but it assumes that floating-point is
786   completely accurate, giving a result that isn't as pretty."
787   (number-dispatch ((x real))
788     (((foreach single-float double-float #!+long-float long-float))
789      (multiple-value-bind (bits exp) (integer-decode-float x)
790        (if (eql bits 0)
791            0
792            (let* ((int (if (minusp x) (- bits) bits))
793                   (digits (float-digits x))
794                   (ex (+ exp digits)))
795              (if (minusp ex)
796                  (integer-/-integer int (ash 1 (+ digits (- ex))))
797                  (integer-/-integer (ash int ex) (ash 1 digits)))))))
798     ((rational) x)))
799
800 (defun rationalize (x)
801   #!+sb-doc
802   "Converts any REAL to a RATIONAL. Floats are converted to a simple rational
803   representation exploiting the assumption that floats are only accurate to
804   their precision. RATIONALIZE (and also RATIONAL) preserve the invariant:
805       (= x (float (rationalize x) x))"
806   (number-dispatch ((x real))
807     (((foreach single-float double-float #!+long-float long-float))
808      ;; Thanks to Kim Fateman, who stole this function rationalize-float from
809      ;; macsyma's rational. Macsyma'a rationalize was written by the legendary
810      ;; Gosper (rwg). Guy Steele said about Gosper, "He has been called the
811      ;; only living 17th century mathematician and is also the best pdp-10
812      ;; hacker I know." So, if you can understand or debug this code you win
813      ;; big.
814      (cond ((minusp x) (- (rationalize (- x))))
815            ((zerop x) 0)
816            (t
817             (let ((eps (etypecase x
818                            (single-float single-float-epsilon)
819                            (double-float double-float-epsilon)
820                            #!+long-float
821                            (long-float long-float-epsilon)))
822                   (y ())
823                   (a ()))
824               (do ((xx x (setq y (/ (float 1.0 x) (- xx (float a x)))))
825                    (num (setq a (truncate x))
826                         (+ (* (setq a (truncate y)) num) onum))
827                    (den 1 (+ (* a den) oden))
828                    (onum 1 num)
829                    (oden 0 den))
830                   ((and (not (zerop den))
831                         (not (> (abs (/ (- x (/ (float num x)
832                                                 (float den x)))
833                                         x))
834                                 eps)))
835                    (integer-/-integer num den))
836                 (declare ((dispatch-type x) xx)))))))
837     ((rational) x)))