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.)
7 ;;;; This software is part of the SBCL system. See the README file for
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.
16 (in-package "SB!KERNEL")
20 (eval-when (:compile-toplevel :load-toplevel :execute)
22 ;;; These functions let us create floats from bits with the significand
23 ;;; uniformly represented as an integer. This is less efficient for double
24 ;;; floats, but is more convenient when making special values, etc.
25 (defun single-from-bits (sign exp sig)
26 (declare (type bit sign) (type (unsigned-byte 24) sig)
27 (type (unsigned-byte 8) exp))
29 (dpb exp sb!vm:single-float-exponent-byte
30 (dpb sig sb!vm:single-float-significand-byte
31 (if (zerop sign) 0 -1)))))
32 (defun double-from-bits (sign exp sig)
33 (declare (type bit sign) (type (unsigned-byte 53) sig)
34 (type (unsigned-byte 11) exp))
35 (make-double-float (dpb exp sb!vm:double-float-exponent-byte
36 (dpb (ash sig -32) sb!vm:double-float-significand-byte
37 (if (zerop sign) 0 -1)))
38 (ldb (byte 32 0) sig)))
39 #!+(and long-float x86)
40 (defun long-from-bits (sign exp sig)
41 (declare (type bit sign) (type (unsigned-byte 64) sig)
42 (type (unsigned-byte 15) exp))
43 (make-long-float (logior (ash sign 15) exp)
44 (ldb (byte 32 32) sig)
45 (ldb (byte 32 0) sig)))
51 (defconstant least-positive-single-float (single-from-bits 0 0 1))
52 (defconstant least-positive-short-float least-positive-single-float)
53 (defconstant least-negative-single-float (single-from-bits 1 0 1))
54 (defconstant least-negative-short-float least-negative-single-float)
55 (defconstant least-positive-double-float (double-from-bits 0 0 1))
57 (defconstant least-positive-long-float least-positive-double-float)
58 #!+(and long-float x86)
59 (defconstant least-positive-long-float (long-from-bits 0 0 1))
60 (defconstant least-negative-double-float (double-from-bits 1 0 1))
62 (defconstant least-negative-long-float least-negative-double-float)
63 #!+(and long-float x86)
64 (defconstant least-negative-long-float (long-from-bits 1 0 1))
66 (defconstant least-positive-normalized-single-float
67 (single-from-bits 0 sb!vm:single-float-normal-exponent-min 0))
68 (defconstant least-positive-normalized-short-float
69 least-positive-normalized-single-float)
70 (defconstant least-negative-normalized-single-float
71 (single-from-bits 1 sb!vm:single-float-normal-exponent-min 0))
72 (defconstant least-negative-normalized-short-float
73 least-negative-normalized-single-float)
74 (defconstant least-positive-normalized-double-float
75 (double-from-bits 0 sb!vm:double-float-normal-exponent-min 0))
77 (defconstant least-positive-normalized-long-float
78 least-positive-normalized-double-float)
79 #!+(and long-float x86)
80 (defconstant least-positive-normalized-long-float
81 (long-from-bits 0 sb!vm:long-float-normal-exponent-min
82 (ash sb!vm:long-float-hidden-bit 32)))
83 (defconstant least-negative-normalized-double-float
84 (double-from-bits 1 sb!vm:double-float-normal-exponent-min 0))
86 (defconstant least-negative-normalized-long-float
87 least-negative-normalized-double-float)
88 #!+(and long-float x86)
89 (defconstant least-negative-normalized-long-float
90 (long-from-bits 1 sb!vm:long-float-normal-exponent-min
91 (ash sb!vm:long-float-hidden-bit 32)))
93 (defconstant most-positive-single-float
94 (single-from-bits 0 sb!vm:single-float-normal-exponent-max
95 (ldb sb!vm:single-float-significand-byte -1)))
96 (defconstant most-positive-short-float most-positive-single-float)
97 (defconstant most-negative-single-float
98 (single-from-bits 1 sb!vm:single-float-normal-exponent-max
99 (ldb sb!vm:single-float-significand-byte -1)))
100 (defconstant most-negative-short-float most-negative-single-float)
101 (defconstant most-positive-double-float
102 (double-from-bits 0 sb!vm:double-float-normal-exponent-max
103 (ldb (byte sb!vm:double-float-digits 0) -1)))
105 (defconstant most-positive-long-float most-positive-double-float)
106 #!+(and long-float x86)
107 (defconstant most-positive-long-float
108 (long-from-bits 0 sb!vm:long-float-normal-exponent-max
109 (ldb (byte sb!vm:long-float-digits 0) -1)))
110 (defconstant most-negative-double-float
111 (double-from-bits 1 sb!vm:double-float-normal-exponent-max
112 (ldb (byte sb!vm:double-float-digits 0) -1)))
114 (defconstant most-negative-long-float most-negative-double-float)
115 #!+(and long-float x86)
116 (defconstant most-negative-long-float
117 (long-from-bits 1 sb!vm:long-float-normal-exponent-max
118 (ldb (byte sb!vm:long-float-digits 0) -1)))
120 ;;; We don't want to do these DEFCONSTANTs at cross-compilation time,
121 ;;; because the cross-compilation host might not support floating
122 ;;; point infinities.
123 (eval-when (:load-toplevel :execute)
124 (defconstant single-float-positive-infinity
125 (single-from-bits 0 (1+ sb!vm:single-float-normal-exponent-max) 0))
126 (defconstant short-float-positive-infinity single-float-positive-infinity)
127 (defconstant single-float-negative-infinity
128 (single-from-bits 1 (1+ sb!vm:single-float-normal-exponent-max) 0))
129 (defconstant short-float-negative-infinity single-float-negative-infinity)
130 (defconstant double-float-positive-infinity
131 (double-from-bits 0 (1+ sb!vm:double-float-normal-exponent-max) 0))
133 (defconstant long-float-positive-infinity double-float-positive-infinity)
134 #!+(and long-float x86)
135 (defconstant long-float-positive-infinity
136 (long-from-bits 0 (1+ sb!vm:long-float-normal-exponent-max)
137 (ash sb!vm:long-float-hidden-bit 32)))
138 (defconstant double-float-negative-infinity
139 (double-from-bits 1 (1+ sb!vm:double-float-normal-exponent-max) 0))
141 (defconstant long-float-negative-infinity double-float-negative-infinity)
142 #!+(and long-float x86)
143 (defconstant long-float-negative-infinity
144 (long-from-bits 1 (1+ sb!vm:long-float-normal-exponent-max)
145 (ash sb!vm:long-float-hidden-bit 32)))
148 (defconstant single-float-epsilon
149 (single-from-bits 0 (- sb!vm:single-float-bias
150 (1- sb!vm:single-float-digits)) 1))
151 (defconstant short-float-epsilon single-float-epsilon)
152 (defconstant single-float-negative-epsilon
153 (single-from-bits 0 (- sb!vm:single-float-bias sb!vm:single-float-digits) 1))
154 (defconstant short-float-negative-epsilon single-float-negative-epsilon)
155 (defconstant double-float-epsilon
156 (double-from-bits 0 (- sb!vm:double-float-bias
157 (1- sb!vm:double-float-digits)) 1))
159 (defconstant long-float-epsilon double-float-epsilon)
160 #!+(and long-float x86)
161 (defconstant long-float-epsilon
162 (long-from-bits 0 (- sb!vm:long-float-bias (1- sb!vm:long-float-digits))
163 (+ 1 (ash sb!vm:long-float-hidden-bit 32))))
164 (defconstant double-float-negative-epsilon
165 (double-from-bits 0 (- sb!vm:double-float-bias sb!vm:double-float-digits) 1))
167 (defconstant long-float-negative-epsilon double-float-negative-epsilon)
168 #!+(and long-float x86)
169 (defconstant long-float-negative-epsilon
170 (long-from-bits 0 (- sb!vm:long-float-bias sb!vm:long-float-digits)
171 (+ 1 (ash sb!vm:long-float-hidden-bit 32))))
173 ;;;; float predicates and environment query
176 (declaim (maybe-inline float-denormalized-p float-infinity-p float-nan-p
177 float-trapping-nan-p))
179 (defun float-denormalized-p (x)
181 "Return true if the float X is denormalized."
182 (number-dispatch ((x float))
184 (and (zerop (ldb sb!vm:single-float-exponent-byte (single-float-bits x)))
187 (and (zerop (ldb sb!vm:double-float-exponent-byte
188 (double-float-high-bits x)))
190 #!+(and long-float x86)
192 (and (zerop (ldb sb!vm:long-float-exponent-byte (long-float-exp-bits x)))
195 (macrolet ((def-frob (name doc single double #!+(and long-float x86) long)
198 (number-dispatch ((x float))
200 (let ((bits (single-float-bits x)))
201 (and (> (ldb sb!vm:single-float-exponent-byte bits)
202 sb!vm:single-float-normal-exponent-max)
205 (let ((hi (double-float-high-bits x))
206 (lo (double-float-low-bits x)))
207 (declare (ignorable lo))
208 (and (> (ldb sb!vm:double-float-exponent-byte hi)
209 sb!vm:double-float-normal-exponent-max)
211 #!+(and long-float x86)
213 (let ((exp (long-float-exp-bits x))
214 (hi (long-float-high-bits x))
215 (lo (long-float-low-bits x)))
216 (declare (ignorable lo))
217 (and (> (ldb sb!vm:long-float-exponent-byte exp)
218 sb!vm:long-float-normal-exponent-max)
221 (def-frob float-infinity-p
222 "Return true if the float X is an infinity (+ or -)."
223 (zerop (ldb sb!vm:single-float-significand-byte bits))
224 (and (zerop (ldb sb!vm:double-float-significand-byte hi))
226 #!+(and long-float x86)
227 (and (zerop (ldb sb!vm:long-float-significand-byte hi))
230 (def-frob float-nan-p
231 "Return true if the float X is a NaN (Not a Number)."
232 (not (zerop (ldb sb!vm:single-float-significand-byte bits)))
233 (or (not (zerop (ldb sb!vm:double-float-significand-byte hi)))
235 #!+(and long-float x86)
236 (or (not (zerop (ldb sb!vm:long-float-significand-byte hi)))
239 (def-frob float-trapping-nan-p
240 "Return true if the float X is a trapping NaN (Not a Number)."
241 (zerop (logand (ldb sb!vm:single-float-significand-byte bits)
242 sb!vm:single-float-trapping-nan-bit))
243 (zerop (logand (ldb sb!vm:double-float-significand-byte hi)
244 sb!vm:double-float-trapping-nan-bit))
245 #!+(and long-float x86)
246 (zerop (logand (ldb sb!vm:long-float-significand-byte hi)
247 sb!vm:long-float-trapping-nan-bit))))
249 ;;; If denormalized, use a subfunction from INTEGER-DECODE-FLOAT to find the
250 ;;; actual exponent (and hence how denormalized it is), otherwise we just
251 ;;; return the number of digits or 0.
252 #!-sb-fluid (declaim (maybe-inline float-precision))
253 (defun float-precision (f)
255 "Returns a non-negative number of significant digits in its float argument.
256 Will be less than FLOAT-DIGITS if denormalized or zero."
257 (macrolet ((frob (digits bias decode)
259 ((float-denormalized-p f)
260 (multiple-value-bind (ignore exp) (,decode f)
261 (declare (ignore ignore))
263 (+ ,digits (1- ,digits) ,bias exp))))
266 (number-dispatch ((f float))
268 (frob sb!vm:single-float-digits sb!vm:single-float-bias
269 integer-decode-single-denorm))
271 (frob sb!vm:double-float-digits sb!vm:double-float-bias
272 integer-decode-double-denorm))
275 (frob sb!vm:long-float-digits sb!vm:long-float-bias
276 integer-decode-long-denorm)))))
278 (defun float-sign (float1 &optional (float2 (float 1 float1)))
280 "Returns a floating-point number that has the same sign as
281 float1 and, if float2 is given, has the same absolute value
283 (declare (float float1 float2))
284 (* (if (etypecase float1
285 (single-float (minusp (single-float-bits float1)))
286 (double-float (minusp (double-float-high-bits float1)))
288 (long-float (minusp (long-float-exp-bits float1))))
293 (defun float-format-digits (format)
295 ((short-float single-float) sb!vm:single-float-digits)
296 ((double-float #!-long-float long-float) sb!vm:double-float-digits)
298 (long-float sb!vm:long-float-digits)))
300 #!-sb-fluid (declaim (inline float-digits float-radix))
302 (defun float-digits (f)
303 (number-dispatch ((f float))
304 ((single-float) sb!vm:single-float-digits)
305 ((double-float) sb!vm:double-float-digits)
307 ((long-float) sb!vm:long-float-digits)))
309 (defun float-radix (x)
311 "Returns (as an integer) the radix b of its floating-point
313 (declare (type float x) (ignore x))
316 ;;;; INTEGER-DECODE-FLOAT and DECODE-FLOAT
319 (declaim (maybe-inline integer-decode-single-float
320 integer-decode-double-float))
322 ;;; Handle the denormalized case of INTEGER-DECODE-FLOAT for SINGLE-FLOAT.
323 (defun integer-decode-single-denorm (x)
324 (declare (type single-float x))
325 (let* ((bits (single-float-bits (abs x)))
326 (sig (ash (ldb sb!vm:single-float-significand-byte bits) 1))
328 (declare (type (unsigned-byte 24) sig)
329 (type (integer 0 23) extra-bias))
331 (unless (zerop (logand sig sb!vm:single-float-hidden-bit))
333 (setq sig (ash sig 1))
336 (- (- sb!vm:single-float-bias)
337 sb!vm:single-float-digits
339 (if (minusp (float-sign x)) -1 1))))
341 ;;; Handle the single-float case of INTEGER-DECODE-FLOAT. If an infinity or
342 ;;; NaN, error. If a denorm, call i-d-s-DENORM to handle it.
343 (defun integer-decode-single-float (x)
344 (declare (single-float x))
345 (let* ((bits (single-float-bits (abs x)))
346 (exp (ldb sb!vm:single-float-exponent-byte bits))
347 (sig (ldb sb!vm:single-float-significand-byte bits))
348 (sign (if (minusp (float-sign x)) -1 1))
349 (biased (- exp sb!vm:single-float-bias sb!vm:single-float-digits)))
350 (declare (fixnum biased))
351 (unless (<= exp sb!vm:single-float-normal-exponent-max)
352 (error "can't decode NaN or infinity: ~S" x))
353 (cond ((and (zerop exp) (zerop sig))
354 (values 0 biased sign))
355 ((< exp sb!vm:single-float-normal-exponent-min)
356 (integer-decode-single-denorm x))
358 (values (logior sig sb!vm:single-float-hidden-bit) biased sign)))))
360 ;;; Like INTEGER-DECODE-SINGLE-DENORM, only doubly so.
361 (defun integer-decode-double-denorm (x)
362 (declare (type double-float x))
363 (let* ((high-bits (double-float-high-bits (abs x)))
364 (sig-high (ldb sb!vm:double-float-significand-byte high-bits))
365 (low-bits (double-float-low-bits x))
366 (sign (if (minusp (float-sign x)) -1 1))
367 (biased (- (- sb!vm:double-float-bias) sb!vm:double-float-digits)))
370 (extra-bias (- sb!vm:double-float-digits 33))
372 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
374 (unless (zerop (logand sig bit)) (return))
375 (setq sig (ash sig 1))
377 (values (ash sig (- sb!vm:double-float-digits 32))
378 (truly-the fixnum (- biased extra-bias))
380 (let ((sig (ash sig-high 1))
382 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
384 (unless (zerop (logand sig sb!vm:double-float-hidden-bit))
386 (setq sig (ash sig 1))
388 (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
389 (truly-the fixnum (- biased extra-bias))
392 ;;; Like INTEGER-DECODE-SINGLE-FLOAT, only doubly so.
393 (defun integer-decode-double-float (x)
394 (declare (double-float x))
396 (hi (double-float-high-bits abs))
397 (lo (double-float-low-bits abs))
398 (exp (ldb sb!vm:double-float-exponent-byte hi))
399 (sig (ldb sb!vm:double-float-significand-byte hi))
400 (sign (if (minusp (float-sign x)) -1 1))
401 (biased (- exp sb!vm:double-float-bias sb!vm:double-float-digits)))
402 (declare (fixnum biased))
403 (unless (<= exp sb!vm:double-float-normal-exponent-max)
404 (error "Can't decode NaN or infinity: ~S." x))
405 (cond ((and (zerop exp) (zerop sig) (zerop lo))
406 (values 0 biased sign))
407 ((< exp sb!vm:double-float-normal-exponent-min)
408 (integer-decode-double-denorm x))
411 (logior (ash (logior (ldb sb!vm:double-float-significand-byte hi)
412 sb!vm:double-float-hidden-bit)
417 #!+(and long-float x86)
418 (defun integer-decode-long-denorm (x)
419 (declare (type long-float x))
420 (let* ((high-bits (long-float-high-bits (abs x)))
421 (sig-high (ldb sb!vm:long-float-significand-byte high-bits))
422 (low-bits (long-float-low-bits x))
423 (sign (if (minusp (float-sign x)) -1 1))
424 (biased (- (- sb!vm:long-float-bias) sb!vm:long-float-digits)))
427 (extra-bias (- sb!vm:long-float-digits 33))
429 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
431 (unless (zerop (logand sig bit)) (return))
432 (setq sig (ash sig 1))
434 (values (ash sig (- sb!vm:long-float-digits 32))
435 (truly-the fixnum (- biased extra-bias))
437 (let ((sig (ash sig-high 1))
439 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
441 (unless (zerop (logand sig sb!vm:long-float-hidden-bit))
443 (setq sig (ash sig 1))
445 (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
446 (truly-the fixnum (- biased extra-bias))
449 #!+(and long-float x86)
450 (defun integer-decode-long-float (x)
451 (declare (long-float x))
452 (let* ((hi (long-float-high-bits x))
453 (lo (long-float-low-bits x))
454 (exp-bits (long-float-exp-bits x))
455 (exp (ldb sb!vm:long-float-exponent-byte exp-bits))
456 (sign (if (minusp exp-bits) -1 1))
457 (biased (- exp sb!vm:long-float-bias sb!vm:long-float-digits)))
458 (declare (fixnum biased))
459 (unless (<= exp sb!vm:long-float-normal-exponent-max)
460 (error "can't decode NaN or infinity: ~S" x))
461 (cond ((and (zerop exp) (zerop hi) (zerop lo))
462 (values 0 biased sign))
463 ((< exp sb!vm:long-float-normal-exponent-min)
464 (integer-decode-long-denorm x))
466 (values (logior (ash hi 32) lo) biased sign)))))
468 ;;; Dispatch to the correct type-specific i-d-f function.
469 (defun integer-decode-float (x)
471 "Returns three values:
472 1) an integer representation of the significand.
473 2) the exponent for the power of 2 that the significand must be multiplied
474 by to get the actual value. This differs from the DECODE-FLOAT exponent
475 by FLOAT-DIGITS, since the significand has been scaled to have all its
476 digits before the radix point.
477 3) -1 or 1 (i.e. the sign of the argument.)"
478 (number-dispatch ((x float))
480 (integer-decode-single-float x))
482 (integer-decode-double-float x))
485 (integer-decode-long-float x))))
487 #!-sb-fluid (declaim (maybe-inline decode-single-float decode-double-float))
489 ;;; Handle the denormalized case of DECODE-SINGLE-FLOAT. We call
490 ;;; INTEGER-DECODE-SINGLE-DENORM and then make the result into a float.
491 (defun decode-single-denorm (x)
492 (declare (type single-float x))
493 (multiple-value-bind (sig exp sign) (integer-decode-single-denorm x)
494 (values (make-single-float
495 (dpb sig sb!vm:single-float-significand-byte
496 (dpb sb!vm:single-float-bias
497 sb!vm:single-float-exponent-byte
499 (truly-the fixnum (+ exp sb!vm:single-float-digits))
502 ;;; Handle the single-float case of DECODE-FLOAT. If an infinity or NaN,
503 ;;; error. If a denorm, call d-s-DENORM to handle it.
504 (defun decode-single-float (x)
505 (declare (single-float x))
506 (let* ((bits (single-float-bits (abs x)))
507 (exp (ldb sb!vm:single-float-exponent-byte bits))
508 (sign (float-sign x))
509 (biased (truly-the single-float-exponent
510 (- exp sb!vm:single-float-bias))))
511 (unless (<= exp sb!vm:single-float-normal-exponent-max)
512 (error "can't decode NaN or infinity: ~S" x))
514 (values 0.0f0 biased sign))
515 ((< exp sb!vm:single-float-normal-exponent-min)
516 (decode-single-denorm x))
518 (values (make-single-float
519 (dpb sb!vm:single-float-bias
520 sb!vm:single-float-exponent-byte
524 ;;; Like DECODE-SINGLE-DENORM, only doubly so.
525 (defun decode-double-denorm (x)
526 (declare (double-float x))
527 (multiple-value-bind (sig exp sign) (integer-decode-double-denorm x)
528 (values (make-double-float
529 (dpb (logand (ash sig -32) (lognot sb!vm:double-float-hidden-bit))
530 sb!vm:double-float-significand-byte
531 (dpb sb!vm:double-float-bias
532 sb!vm:double-float-exponent-byte 0))
533 (ldb (byte 32 0) sig))
534 (truly-the fixnum (+ exp sb!vm:double-float-digits))
537 ;;; Like DECODE-SINGLE-FLOAT, only doubly so.
538 (defun decode-double-float (x)
539 (declare (double-float x))
541 (hi (double-float-high-bits abs))
542 (lo (double-float-low-bits abs))
543 (exp (ldb sb!vm:double-float-exponent-byte hi))
544 (sign (float-sign x))
545 (biased (truly-the double-float-exponent
546 (- exp sb!vm:double-float-bias))))
547 (unless (<= exp sb!vm:double-float-normal-exponent-max)
548 (error "can't decode NaN or infinity: ~S" x))
550 (values 0.0d0 biased sign))
551 ((< exp sb!vm:double-float-normal-exponent-min)
552 (decode-double-denorm x))
554 (values (make-double-float
555 (dpb sb!vm:double-float-bias
556 sb!vm:double-float-exponent-byte hi)
560 #!+(and long-float x86)
561 (defun decode-long-denorm (x)
562 (declare (long-float x))
563 (multiple-value-bind (sig exp sign) (integer-decode-long-denorm x)
564 (values (make-long-float sb!vm:long-float-bias (ash sig -32)
565 (ldb (byte 32 0) sig))
566 (truly-the fixnum (+ exp sb!vm:long-float-digits))
569 #!+(and long-float x86)
570 (defun decode-long-float (x)
571 (declare (long-float x))
572 (let* ((hi (long-float-high-bits x))
573 (lo (long-float-low-bits x))
574 (exp-bits (long-float-exp-bits x))
575 (exp (ldb sb!vm:long-float-exponent-byte exp-bits))
576 (sign (if (minusp exp-bits) -1l0 1l0))
577 (biased (truly-the long-float-exponent
578 (- exp sb!vm:long-float-bias))))
579 (unless (<= exp sb!vm:long-float-normal-exponent-max)
580 (error "can't decode NaN or infinity: ~S" x))
582 (values 0.0l0 biased sign))
583 ((< exp sb!vm:long-float-normal-exponent-min)
584 (decode-long-denorm x))
586 (values (make-long-float
587 (dpb sb!vm:long-float-bias sb!vm:long-float-exponent-byte
593 ;;; Dispatch to the appropriate type-specific function.
594 (defun decode-float (f)
596 "Returns three values:
597 1) a floating-point number representing the significand. This is always
598 between 0.5 (inclusive) and 1.0 (exclusive).
599 2) an integer representing the exponent.
600 3) -1.0 or 1.0 (i.e. the sign of the argument.)"
601 (number-dispatch ((f float))
603 (decode-single-float f))
605 (decode-double-float f))
608 (decode-long-float f))))
612 #!-sb-fluid (declaim (maybe-inline scale-single-float scale-double-float))
614 ;;; Handle float scaling where the X is denormalized or the result is
615 ;;; denormalized or underflows to 0.
616 (defun scale-float-maybe-underflow (x exp)
617 (multiple-value-bind (sig old-exp) (integer-decode-float x)
618 (let* ((digits (float-digits x))
619 (new-exp (+ exp old-exp digits
621 (single-float sb!vm:single-float-bias)
622 (double-float sb!vm:double-float-bias))))
623 (sign (if (minusp (float-sign x)) 1 0)))
627 (single-float sb!vm:single-float-normal-exponent-min)
628 (double-float sb!vm:double-float-normal-exponent-min)))
629 (when (sb!vm:current-float-trap :inexact)
630 (error 'floating-point-inexact :operation 'scale-float
631 :operands (list x exp)))
632 (when (sb!vm:current-float-trap :underflow)
633 (error 'floating-point-underflow :operation 'scale-float
634 :operands (list x exp)))
635 (let ((shift (1- new-exp)))
636 (if (< shift (- (1- digits)))
639 (single-float (single-from-bits sign 0 (ash sig shift)))
640 (double-float (double-from-bits sign 0 (ash sig shift)))))))
643 (single-float (single-from-bits sign new-exp sig))
644 (double-float (double-from-bits sign new-exp sig))))))))
646 ;;; Called when scaling a float overflows, or the original float was a
647 ;;; NaN or infinity. If overflow errors are trapped, then error,
648 ;;; otherwise return the appropriate infinity. If a NaN, signal or not
650 (defun scale-float-maybe-overflow (x exp)
652 ((float-infinity-p x)
653 ;; Infinity is infinity, no matter how small...
656 (when (and (float-trapping-nan-p x)
657 (sb!vm:current-float-trap :invalid))
658 (error 'floating-point-invalid-operation :operation 'scale-float
659 :operands (list x exp)))
662 (when (sb!vm:current-float-trap :overflow)
663 (error 'floating-point-overflow :operation 'scale-float
664 :operands (list x exp)))
665 (when (sb!vm:current-float-trap :inexact)
666 (error 'floating-point-inexact :operation 'scale-float
667 :operands (list x exp)))
670 (single-float single-float-positive-infinity)
671 (double-float double-float-positive-infinity))))))
673 ;;; Scale a single or double float, calling the correct over/underflow
675 (defun scale-single-float (x exp)
676 (declare (single-float x) (fixnum exp))
677 (let* ((bits (single-float-bits x))
678 (old-exp (ldb sb!vm:single-float-exponent-byte bits))
679 (new-exp (+ old-exp exp)))
682 ((or (< old-exp sb!vm:single-float-normal-exponent-min)
683 (< new-exp sb!vm:single-float-normal-exponent-min))
684 (scale-float-maybe-underflow x exp))
685 ((or (> old-exp sb!vm:single-float-normal-exponent-max)
686 (> new-exp sb!vm:single-float-normal-exponent-max))
687 (scale-float-maybe-overflow x exp))
689 (make-single-float (dpb new-exp
690 sb!vm:single-float-exponent-byte
692 (defun scale-double-float (x exp)
693 (declare (double-float x) (fixnum exp))
694 (let* ((hi (double-float-high-bits x))
695 (lo (double-float-low-bits x))
696 (old-exp (ldb sb!vm:double-float-exponent-byte hi))
697 (new-exp (+ old-exp exp)))
700 ((or (< old-exp sb!vm:double-float-normal-exponent-min)
701 (< new-exp sb!vm:double-float-normal-exponent-min))
702 (scale-float-maybe-underflow x exp))
703 ((or (> old-exp sb!vm:double-float-normal-exponent-max)
704 (> new-exp sb!vm:double-float-normal-exponent-max))
705 (scale-float-maybe-overflow x exp))
707 (make-double-float (dpb new-exp sb!vm:double-float-exponent-byte hi)
710 #!+(and x86 long-float)
711 (defun scale-long-float (x exp)
712 (declare (long-float x) (fixnum exp))
715 ;;; Dispatch to the correct type-specific scale-float function.
716 (defun scale-float (f ex)
718 "Returns the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
719 of precision or overflow."
720 (number-dispatch ((f float))
722 (scale-single-float f ex))
724 (scale-double-float f ex))
727 (scale-long-float f ex))))
729 ;;;; converting to/from floats
731 (defun float (number &optional (other () otherp))
733 "Converts any REAL to a float. If OTHER is not provided, it returns a
734 SINGLE-FLOAT if NUMBER is not already a FLOAT. If OTHER is provided, the
735 result is the same float format as OTHER."
737 (number-dispatch ((number real) (other float))
738 (((foreach rational single-float double-float #!+long-float long-float)
739 (foreach single-float double-float #!+long-float long-float))
740 (coerce number '(dispatch-type other))))
743 (coerce number 'single-float))))
745 (macrolet ((frob (name type)
747 (number-dispatch ((x real))
748 (((foreach single-float double-float #!+long-float long-float
752 (bignum-to-float x ',type))
754 (float-ratio x ',type))))))
755 (frob %single-float single-float)
756 (frob %double-float double-float)
758 (frob %long-float long-float))
760 ;;; Convert a ratio to a float. We avoid any rounding error by doing an
761 ;;; integer division. Accuracy is important to preserve read/print
762 ;;; consistency, since this is ultimately how the reader reads a float. We
763 ;;; scale the numerator by a power of two until the division results in the
764 ;;; desired number of fraction bits, then do round-to-nearest.
765 (defun float-ratio (x format)
766 (let* ((signed-num (numerator x))
767 (plusp (plusp signed-num))
768 (num (if plusp signed-num (- signed-num)))
769 (den (denominator x))
770 (digits (float-format-digits format))
772 (declare (fixnum digits scale))
773 ;; Strip any trailing zeros from the denominator and move it into the scale
774 ;; factor (to minimize the size of the operands.)
775 (let ((den-twos (1- (integer-length (logxor den (1- den))))))
776 (declare (fixnum den-twos))
777 (decf scale den-twos)
778 (setq den (ash den (- den-twos))))
779 ;; Guess how much we need to scale by from the magnitudes of the numerator
780 ;; and denominator. We want one extra bit for a guard bit.
781 (let* ((num-len (integer-length num))
782 (den-len (integer-length den))
783 (delta (- den-len num-len))
784 (shift (1+ (the fixnum (+ delta digits))))
785 (shifted-num (ash num shift)))
786 (declare (fixnum delta shift))
788 (labels ((float-and-scale (bits)
789 (let* ((bits (ash bits -1))
790 (len (integer-length bits)))
791 (cond ((> len digits)
792 (aver (= len (the fixnum (1+ digits))))
793 (scale-float (floatit (ash bits -1)) (1+ scale)))
795 (scale-float (floatit bits) scale)))))
797 (let ((sign (if plusp 0 1)))
800 (single-from-bits sign sb!vm:single-float-bias bits))
802 (double-from-bits sign sb!vm:double-float-bias bits))
805 (long-from-bits sign sb!vm:long-float-bias bits))))))
807 (multiple-value-bind (fraction-and-guard rem)
808 (truncate shifted-num den)
809 (let ((extra (- (integer-length fraction-and-guard) digits)))
810 (declare (fixnum extra))
813 ((oddp fraction-and-guard)
817 (if (zerop (logand fraction-and-guard 2))
819 (1+ fraction-and-guard)))
820 (float-and-scale (1+ fraction-and-guard)))))
822 (return (float-and-scale fraction-and-guard)))))
823 (setq shifted-num (ash shifted-num -1))
827 These might be useful if we ever have a machine w/o float/integer conversion
828 hardware. For now, we'll use special ops that uninterruptibly frob the
829 rounding modes & do ieee round-to-integer.
831 ;;; The compiler compiles a call to this when we are doing %UNARY-TRUNCATE
832 ;;; and the result is known to be a fixnum. We can avoid some generic
833 ;;; arithmetic in this case.
834 (defun %unary-truncate-single-float/fixnum (x)
835 (declare (single-float x) (values fixnum))
836 (locally (declare (optimize (speed 3) (safety 0)))
837 (let* ((bits (single-float-bits x))
838 (exp (ldb sb!vm:single-float-exponent-byte bits))
839 (frac (logior (ldb sb!vm:single-float-significand-byte bits)
840 sb!vm:single-float-hidden-bit))
841 (shift (- exp sb!vm:single-float-digits sb!vm:single-float-bias)))
842 (when (> exp sb!vm:single-float-normal-exponent-max)
843 (error 'floating-point-invalid-operation :operator 'truncate
845 (if (<= shift (- sb!vm:single-float-digits))
847 (let ((res (ash frac shift)))
848 (declare (type (unsigned-byte 31) res))
853 ;;; Double-float version of this operation (see above single op).
854 (defun %unary-truncate-double-float/fixnum (x)
855 (declare (double-float x) (values fixnum))
856 (locally (declare (optimize (speed 3) (safety 0)))
857 (let* ((hi-bits (double-float-high-bits x))
858 (exp (ldb sb!vm:double-float-exponent-byte hi-bits))
859 (frac (logior (ldb sb!vm:double-float-significand-byte hi-bits)
860 sb!vm:double-float-hidden-bit))
861 (shift (- exp (- sb!vm:double-float-digits sb!vm:word-bits)
862 sb!vm:double-float-bias)))
863 (when (> exp sb!vm:double-float-normal-exponent-max)
864 (error 'floating-point-invalid-operation :operator 'truncate
866 (if (<= shift (- sb!vm:word-bits sb!vm:double-float-digits))
868 (let* ((res-hi (ash frac shift))
869 (res (if (plusp shift)
872 (ash (double-float-low-bits x)
873 (- shift sb!vm:word-bits))))
875 (declare (type (unsigned-byte 31) res-hi res))
881 ;;; This function is called when we are doing a truncate without any funky
882 ;;; divisor, i.e. converting a float or ratio to an integer. Note that we do
883 ;;; *not* return the second value of truncate, so it must be computed by the
884 ;;; caller if needed.
886 ;;; In the float case, we pick off small arguments so that compiler can use
887 ;;; special-case operations. We use an exclusive test, since (due to round-off
888 ;;; error), (float most-positive-fixnum) may be greater than
889 ;;; most-positive-fixnum.
890 (defun %unary-truncate (number)
891 (number-dispatch ((number real))
893 ((ratio) (values (truncate (numerator number) (denominator number))))
894 (((foreach single-float double-float #!+long-float long-float))
895 (if (< (float most-negative-fixnum number)
897 (float most-positive-fixnum number))
898 (truly-the fixnum (%unary-truncate number))
899 (multiple-value-bind (bits exp) (integer-decode-float number)
900 (let ((res (ash bits exp)))
905 ;;; Similar to %UNARY-TRUNCATE, but rounds to the nearest integer. If we
906 ;;; can't use the round primitive, then we do our own round-to-nearest on the
907 ;;; result of i-d-f. [Note that this rounding will really only happen with
908 ;;; double floats, since the whole single-float fraction will fit in a fixnum,
909 ;;; so all single-floats larger than most-positive-fixnum can be precisely
910 ;;; represented by an integer.]
911 (defun %unary-round (number)
912 (number-dispatch ((number real))
914 ((ratio) (values (round (numerator number) (denominator number))))
915 (((foreach single-float double-float #!+long-float long-float))
916 (if (< (float most-negative-fixnum number)
918 (float most-positive-fixnum number))
919 (truly-the fixnum (%unary-round number))
920 (multiple-value-bind (bits exp) (integer-decode-float number)
921 (let* ((shifted (ash bits exp))
922 (rounded (if (and (minusp exp)
925 (lognot (ash -1 (- exp))))
935 "RATIONAL produces a rational number for any real numeric argument. This is
936 more efficient than RATIONALIZE, but it assumes that floating-point is
937 completely accurate, giving a result that isn't as pretty."
938 (number-dispatch ((x real))
939 (((foreach single-float double-float #!+long-float long-float))
940 (multiple-value-bind (bits exp) (integer-decode-float x)
943 (let* ((int (if (minusp x) (- bits) bits))
944 (digits (float-digits x))
947 (integer-/-integer int (ash 1 (+ digits (- ex))))
948 (integer-/-integer (ash int ex) (ash 1 digits)))))))
951 (defun rationalize (x)
953 "Converts any REAL to a RATIONAL. Floats are converted to a simple rational
954 representation exploiting the assumption that floats are only accurate to
955 their precision. RATIONALIZE (and also RATIONAL) preserve the invariant:
956 (= x (float (rationalize x) x))"
957 (number-dispatch ((x real))
958 (((foreach single-float double-float #!+long-float long-float))
959 ;; Thanks to Kim Fateman, who stole this function rationalize-float from
960 ;; macsyma's rational. Macsyma'a rationalize was written by the legendary
961 ;; Gosper (rwg). Guy Steele said about Gosper, "He has been called the
962 ;; only living 17th century mathematician and is also the best pdp-10
963 ;; hacker I know." So, if you can understand or debug this code you win
965 (cond ((minusp x) (- (rationalize (- x))))
968 (let ((eps (etypecase x
969 (single-float single-float-epsilon)
970 (double-float double-float-epsilon)
972 (long-float long-float-epsilon)))
975 (do ((xx x (setq y (/ (float 1.0 x) (- xx (float a x)))))
976 (num (setq a (truncate x))
977 (+ (* (setq a (truncate y)) num) onum))
978 (den 1 (+ (* a den) oden))
981 ((and (not (zerop den))
982 (not (> (abs (/ (- x (/ (float num x)
986 (integer-/-integer num den))
987 (declare ((dispatch-type x) xx)))))))