src/code/float.lisp

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