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 ;;; Don't need to define it in the host in both passes
  20 (eval-when (#-sb-xc :compile-toplevel :load-toplevel :execute)
  21 ;;; These functions let us create floats from bits with the
  22 ;;; significand uniformly represented as an integer. This is less
  23 ;;; efficient for double floats, but is more convenient when making
  24 ;;; special values, etc.
  25 (declaim (inline single-from-bits double-from-bits))
  26
  27 (defun single-from-bits (sign exp sig)
  28   (declare (type bit sign) (type (unsigned-byte 24) sig)
  29            (type (unsigned-byte 8) exp))
  30   (make-single-float
  31    (dpb exp sb-vm:single-float-exponent-byte
  32         (dpb sig sb-vm:single-float-significand-byte
  33              (if (zerop sign) 0 -1)))))
  34 (defun double-from-bits (sign exp sig)
  35   (declare (type bit sign) (type (unsigned-byte 53) sig)
  36            (type (unsigned-byte 11) exp))
  37   #-64-bit
  38   (make-double-float (dpb exp sb-vm:double-float-hi-exponent-byte
  39                           (dpb (ash sig -32)
  40                                sb-vm:double-float-hi-significand-byte
  41                                (if (zerop sign) 0 -1)))
  42                      (ldb (byte 32 0) sig))
  43   #+64-bit
  44   (%make-double-float
  45    (dpb exp sb-vm:double-float-exponent-byte
  46         (dpb sig sb-vm:double-float-significand-byte
  47              (if (zerop sign) 0 -1)))))
  48 #+(and long-float x86)
  49 (defun long-from-bits (sign exp sig)
  50   (declare (type bit sign) (type (unsigned-byte 64) sig)
  51            (type (unsigned-byte 15) exp))
  52   (make-long-float (logior (ash sign 15) exp)
  53                    (ldb (byte 32 32) sig)
  54                    (ldb (byte 32 0) sig)))
  55
  56 ) ; EVAL-WHEN
  57
  58 \f
  59 ;;;; float parameters
  60
  61 (defconstant least-positive-single-float (single-from-bits 0 0 1))
  62 (defconstant least-positive-short-float (single-from-bits 0 0 1))
  63 (defconstant least-negative-single-float (single-from-bits 1 0 1))
  64 (defconstant least-negative-short-float (single-from-bits 1 0 1))
  65 (defconstant least-positive-double-float (double-from-bits 0 0 1))
  66 #-long-float
  67 (defconstant least-positive-long-float (double-from-bits 0 0 1))
  68 #+(and long-float x86)
  69 (defconstant least-positive-long-float (long-from-bits 0 0 1))
  70 (defconstant least-negative-double-float (double-from-bits 1 0 1))
  71 #-long-float
  72 (defconstant least-negative-long-float (double-from-bits 1 0 1))
  73 #+(and long-float x86)
  74 (defconstant least-negative-long-float (long-from-bits 1 0 1))
  75
  76 (defconstant least-positive-normalized-single-float
  77   (single-from-bits 0 sb-vm:single-float-normal-exponent-min 0))
  78 (defconstant least-positive-normalized-short-float
  79   least-positive-normalized-single-float)
  80 (defconstant least-negative-normalized-single-float
  81   (single-from-bits 1 sb-vm:single-float-normal-exponent-min 0))
  82 (defconstant least-negative-normalized-short-float
  83   least-negative-normalized-single-float)
  84 (defconstant least-positive-normalized-double-float
  85   (double-from-bits 0 sb-vm:double-float-normal-exponent-min 0))
  86 #-long-float
  87 (defconstant least-positive-normalized-long-float
  88   least-positive-normalized-double-float)
  89 #+(and long-float x86)
  90 (defconstant least-positive-normalized-long-float
  91   (long-from-bits 0 sb-vm:long-float-normal-exponent-min
  92                   (ash sb-vm:long-float-hidden-bit 32)))
  93 (defconstant least-negative-normalized-double-float
  94   (double-from-bits 1 sb-vm:double-float-normal-exponent-min 0))
  95 #-long-float
  96 (defconstant least-negative-normalized-long-float
  97   least-negative-normalized-double-float)
  98 #+(and long-float x86)
  99 (defconstant least-negative-normalized-long-float
 100   (long-from-bits 1 sb-vm:long-float-normal-exponent-min
 101                   (ash sb-vm:long-float-hidden-bit 32)))
 102
 103 (defconstant most-positive-single-float
 104   (single-from-bits 0 sb-vm:single-float-normal-exponent-max
 105                     (ldb sb-vm:single-float-significand-byte -1)))
 106 (defconstant most-positive-short-float most-positive-single-float)
 107 (defconstant most-negative-single-float
 108   (single-from-bits 1 sb-vm:single-float-normal-exponent-max
 109                     (ldb sb-vm:single-float-significand-byte -1)))
 110 (defconstant most-negative-short-float most-negative-single-float)
 111 (defconstant most-positive-double-float
 112   (double-from-bits 0 sb-vm:double-float-normal-exponent-max
 113                     (ldb (byte sb-vm:double-float-digits 0) -1)))
 114 #-long-float
 115 (defconstant most-positive-long-float most-positive-double-float)
 116 #+(and long-float x86)
 117 (defconstant most-positive-long-float
 118   (long-from-bits 0 sb-vm:long-float-normal-exponent-max
 119                   (ldb (byte sb-vm:long-float-digits 0) -1)))
 120 (defconstant most-negative-double-float
 121   (double-from-bits 1 sb-vm:double-float-normal-exponent-max
 122                     (ldb (byte sb-vm:double-float-digits 0) -1)))
 123 #-long-float
 124 (defconstant most-negative-long-float most-negative-double-float)
 125 #+(and long-float x86)
 126 (defconstant most-negative-long-float
 127   (long-from-bits 1 sb-vm:long-float-normal-exponent-max
 128                   (ldb (byte sb-vm:long-float-digits 0) -1)))
 129
 130 (defconstant single-float-positive-infinity
 131   (single-from-bits 0 (1+ sb-vm:single-float-normal-exponent-max) 0))
 132 (defconstant short-float-positive-infinity single-float-positive-infinity)
 133 (defconstant single-float-negative-infinity
 134   (single-from-bits 1 (1+ sb-vm:single-float-normal-exponent-max) 0))
 135 (defconstant short-float-negative-infinity single-float-negative-infinity)
 136 (defconstant double-float-positive-infinity
 137   (double-from-bits 0 (1+ sb-vm:double-float-normal-exponent-max) 0))
 138 #+(not long-float)
 139 (defconstant long-float-positive-infinity double-float-positive-infinity)
 140 #+(and long-float x86)
 141 (defconstant long-float-positive-infinity
 142   (long-from-bits 0 (1+ sb-vm:long-float-normal-exponent-max)
 143                   (ash sb-vm:long-float-hidden-bit 32)))
 144 (defconstant double-float-negative-infinity
 145   (double-from-bits 1 (1+ sb-vm:double-float-normal-exponent-max) 0))
 146 #+(not long-float)
 147 (defconstant long-float-negative-infinity double-float-negative-infinity)
 148 #+(and long-float x86)
 149 (defconstant long-float-negative-infinity
 150   (long-from-bits 1 (1+ sb-vm:long-float-normal-exponent-max)
 151                   (ash sb-vm:long-float-hidden-bit 32)))
 152
 153 (defconstant single-float-epsilon
 154   (single-from-bits 0 (- sb-vm:single-float-bias
 155                          (1- sb-vm:single-float-digits)) 1))
 156 (defconstant short-float-epsilon single-float-epsilon)
 157 (defconstant single-float-negative-epsilon
 158   (single-from-bits 0 (- sb-vm:single-float-bias sb-vm:single-float-digits) 1))
 159 (defconstant short-float-negative-epsilon single-float-negative-epsilon)
 160 (defconstant double-float-epsilon
 161   (double-from-bits 0 (- sb-vm:double-float-bias
 162                          (1- sb-vm:double-float-digits)) 1))
 163 #-long-float
 164 (defconstant long-float-epsilon double-float-epsilon)
 165 #+(and long-float x86)
 166 (defconstant long-float-epsilon
 167   (long-from-bits 0 (- sb-vm:long-float-bias (1- sb-vm:long-float-digits))
 168                   (+ 1 (ash sb-vm:long-float-hidden-bit 32))))
 169 (defconstant double-float-negative-epsilon
 170   (double-from-bits 0 (- sb-vm:double-float-bias sb-vm:double-float-digits) 1))
 171 #-long-float
 172 (defconstant long-float-negative-epsilon double-float-negative-epsilon)
 173 #+(and long-float x86)
 174 (defconstant long-float-negative-epsilon
 175   (long-from-bits 0 (- sb-vm:long-float-bias sb-vm:long-float-digits)
 176                   (+ 1 (ash sb-vm:long-float-hidden-bit 32))))
 177
 178 ;;; Limits for floats that can be truncated into a fixnum
 179 ;;; with no loss of precision.
 180 ;;; (We don't have constants for "most-fooative-fixnum as a mumble-float")
 181 (defconstant most-positive-fixnum-single-float
 182   (single-from-bits 0 (+ sb-vm:n-fixnum-bits sb-vm:single-float-bias -1)
 183                                (ldb (byte (1- sb-vm:single-float-digits) 0) -1)))
 184
 185 (defconstant most-negative-fixnum-single-float
 186   (single-from-bits 1 (+ sb-vm:n-fixnum-bits sb-vm:single-float-bias) 0))
 187
 188 (defconstant most-positive-fixnum-double-float
 189   (double-from-bits 0 (+ sb-vm:n-fixnum-bits sb-vm:double-float-bias -1)
 190                                (ldb (byte (1- sb-vm:double-float-digits) 0) -1)))
 191
 192 (defconstant most-negative-fixnum-double-float
 193   (double-from-bits 1 (+ sb-vm:n-fixnum-bits sb-vm:double-float-bias) 0))
 194
 195 \f
 196 ;;;; float predicates and environment query
 197
 198 (declaim (maybe-inline float-denormalized-p float-infinity-p float-nan-p
 199                        float-trapping-nan-p))
 200
 201 (defmacro sfloat-bits-subnormalp (bits)
 202   `(zerop (ldb sb-vm:single-float-exponent-byte ,bits)))
 203 #-64-bit
 204 (defmacro dfloat-high-bits-subnormalp (bits)
 205   `(zerop (ldb sb-vm:double-float-hi-exponent-byte ,bits)))
 206 #+64-bit
 207 (defmacro dfloat-bits-subnormalp (bits)
 208   `(zerop (ldb sb-vm:double-float-exponent-byte ,bits)))
 209
 210 (defun float-denormalized-p (x)
 211   "Return true if the float X is denormalized."
 212   (declare (explicit-check))
 213   (number-dispatch ((x float))
 214     ((single-float)
 215      #+64-bit
 216      (let ((bits (single-float-bits x)))
 217        (and (ldb-test (byte 31 0) bits) ; is nonzero (disregard the sign bit)
 218             (sfloat-bits-subnormalp bits)))
 219      #-64-bit
 220      (and (zerop (ldb sb-vm:single-float-exponent-byte (single-float-bits x)))
 221           (not (zerop x))))
 222     ((double-float)
 223      #+64-bit
 224      (let ((bits (double-float-bits x)))
 225        ;; is nonzero after shifting out the sign bit
 226        (and (not (zerop (logand (ash bits 1) most-positive-word)))
 227             (dfloat-bits-subnormalp bits)))
 228      #-64-bit
 229      (and (zerop (ldb sb-vm:double-float-hi-exponent-byte
 230                       (double-float-high-bits x)))
 231           (not (zerop x))))
 232     #+(and long-float x86)
 233     ((long-float)
 234      (and (zerop (ldb sb-vm:long-float-exponent-byte (long-float-exp-bits x)))
 235           (not (zerop x))))))
 236
 237 (defmacro float-inf-or-nan-test (var single double #+(and long-float x86) long)
 238   `(number-dispatch ((,var float))
 239      ((single-float)
 240       (let ((bits (single-float-bits ,var)))
 241         (and (> (ldb sb-vm:single-float-exponent-byte bits)
 242                 sb-vm:single-float-normal-exponent-max)
 243              ,single)))
 244      ((double-float)
 245       #+64-bit
 246       (let ((bits (double-float-bits ,var)))
 247         (and (> (ldb sb-vm:double-float-exponent-byte bits)
 248                 sb-vm:double-float-normal-exponent-max)
 249              ,double))
 250       #-64-bit
 251       (let ((hi (double-float-high-bits ,var))
 252             (lo (double-float-low-bits ,var)))
 253         (declare (ignorable lo))
 254         (and (> (ldb sb-vm:double-float-hi-exponent-byte hi)
 255                 sb-vm:double-float-normal-exponent-max)
 256              ,double)))
 257      #+(and long-float x86)
 258      ((long-float)
 259       (let ((exp (long-float-exp-bits ,var))
 260             (hi (long-float-high-bits ,var))
 261             (lo (long-float-low-bits ,var)))
 262         (declare (ignorable lo))
 263         (and (> (ldb sb-vm:long-float-exponent-byte exp)
 264                 sb-vm:long-float-normal-exponent-max)
 265              ,long)))))
 266
 267 ;; Infinities and NANs have the maximum exponent
 268 (defun float-infinity-or-nan-p (x)
 269   (float-inf-or-nan-test x t t #+(and long-float x86) t))
 270
 271 ;; Infinity has 0 for the significand
 272 (defun float-infinity-p (x)
 273   "Return true if the float X is an infinity (+ or -)."
 274   (float-inf-or-nan-test
 275    x
 276    (zerop (ldb sb-vm:single-float-significand-byte bits))
 277
 278    #+64-bit (zerop (ldb sb-vm:double-float-significand-byte bits))
 279    #-64-bit (zerop (logior (ldb sb-vm:double-float-hi-significand-byte hi) lo))
 280
 281    #+(and long-float x86)
 282    (and (zerop (ldb sb-vm:long-float-significand-byte hi))
 283         (zerop lo))))
 284
 285 ;; NaNs have nonzero for the significand
 286 (defun float-nan-p (x)
 287   "Return true if the float X is a NaN (Not a Number)."
 288   (float-inf-or-nan-test
 289    x
 290    (not (zerop (ldb sb-vm:single-float-significand-byte bits)))
 291
 292    #+64-bit (not (zerop (ldb sb-vm:double-float-significand-byte bits)))
 293    #-64-bit (not (zerop (logior (ldb sb-vm:double-float-hi-significand-byte hi) lo)))
 294
 295    #+(and long-float x86)
 296    (or (not (zerop (ldb sb-vm:long-float-significand-byte hi)))
 297        (not (zerop lo)))))
 298
 299 (defmacro with-float-inf-or-nan-test (float infinity nan normal)
 300   `(block nil
 301      ,(if (equal infinity nan)
 302           `(float-inf-or-nan-test
 303             ,float
 304             (return ,nan)
 305             (return ,nan))
 306           `(float-inf-or-nan-test
 307             ,float
 308             (if (zerop (ldb sb-vm:single-float-significand-byte bits))
 309                 (return ,infinity)
 310                 (return ,nan))
 311             (if #+64-bit (zerop (ldb sb-vm:double-float-significand-byte bits))
 312                 #-64-bit (zerop (logior (ldb sb-vm:double-float-hi-significand-byte hi) lo))
 313                 (return ,infinity)
 314                 (return ,nan))))
 315      ,normal))
 316
 317 (defun float-trapping-nan-p (x)
 318   "Return true if the float X is a trapping NaN (Not a Number)."
 319   ;; MIPS has trapping NaNs (SNaNs) with the trapping-nan-bit SET.
 320   ;; All the others have trapping NaNs (SNaNs) with the
 321   ;; trapping-nan-bit CLEAR.  Note that the given implementation
 322   ;; considers infinities to be FLOAT-TRAPPING-NAN-P on most
 323   ;; architectures.
 324   (float-inf-or-nan-test
 325    x
 326    ;; SINGLE-FLOAT
 327    #+mips (logbitp 22 bits)
 328    #-mips (not (logbitp 22 bits))
 329
 330    ;; DOUBLE-FLOAT
 331    #+mips (logbitp 19 hi)
 332    #+(and (not mips) 64-bit) (not (logbitp 51 bits))
 333    #+(and (not mips) (not 64-bit)) (not (logbitp 19 hi))
 334
 335    ;; LONG-FLOAT (this code is dead anyway)
 336    #+(and long-float x86)
 337    (zerop (logand (ldb sb-vm:long-float-significand-byte hi)
 338                   (ash 1 30)))))
 339
 340 ;;; If denormalized, use a subfunction from INTEGER-DECODE-FLOAT to find the
 341 ;;; actual exponent (and hence how denormalized it is), otherwise we just
 342 ;;; return the number of digits or 0.
 343 (declaim (maybe-inline float-precision))
 344 (defun float-precision (f)
 345   "Return a non-negative number of significant digits in its float argument.
 346   Will be less than FLOAT-DIGITS if denormalized or zero."
 347   (declare (explicit-check))
 348   (integer-length
 349    (number-dispatch ((f float))
 350      ((single-float)
 351       (let ((bits (single-float-bits f)))
 352         (if (sfloat-bits-subnormalp bits)
 353             (ldb sb-vm:single-float-significand-byte bits)
 354             (return-from float-precision sb-vm:single-float-digits))))
 355      ((double-float)
 356       #+64-bit
 357       (let ((bits (double-float-bits f)))
 358         (if (dfloat-bits-subnormalp bits)
 359             (ldb sb-vm:double-float-significand-byte bits)
 360             (return-from float-precision sb-vm:double-float-digits)))
 361       #-64-bit
 362       (let ((high (double-float-high-bits f)))
 363         (if (not (dfloat-high-bits-subnormalp high))
 364             (return-from float-precision sb-vm:double-float-digits)
 365             (let ((n (integer-length (ldb sb-vm:double-float-hi-significand-byte high))))
 366               (if (/= 0 n)
 367                   (return-from float-precision (+ n 32))
 368                   (double-float-low-bits f)))))))))
 369
 370 (defun float-sign (float1 &optional (float2 (float 1 float1)))
 371   "Return a floating-point number that has the same sign as
 372    FLOAT1 and, if FLOAT2 is given, has the same absolute value
 373    as FLOAT2."
 374   (declare (float float1 float2) (explicit-check))
 375   (sb-xc:*
 376    (if (etypecase float1
 377          (single-float (minusp (single-float-bits float1)))
 378          ;; If 64-bits words, use all the bits. No need to right-shift them.
 379          (double-float (minusp #+64-bit (double-float-bits float1)
 380                                #-64-bit (double-float-high-bits float1)))
 381          #+long-float
 382          (long-float (minusp (long-float-exp-bits float1))))
 383        (float -1 float1)
 384        (float 1 float1))
 385    (abs float2)))
 386
 387 ;;; When all we want is the sign bit, there is a simpler way to extract it
 388 ;;; than via either integer-decode-float or float-sign. Just shift the msb
 389 ;;; over to the lsb position. FLOAT-SIGN produces some pretty horrific code
 390 ;;; if the specific subtype of float is unnown:
 391 ;;;  (minusp (float-sign x)) becomes (< (float-sign x) (float 0 x))
 392 ;;; which ends up calling not only FLOAT-SIGN, but also FLOAT merely to cast
 393 ;;; the integer 0 into a float of whatever type X is.
 394 #-sb-xc-host
 395 (defun float-sign-bit (x)      ; return 1 or 0, literally the sign bit
 396   (declare (explicit-check))
 397   (number-dispatch ((x float))
 398     ((single-float) (float-sign-bit x))
 399     ((double-float) (float-sign-bit x))))
 400 #-sb-xc-host
 401 (defun float-sign-bit-set-p (x)
 402   (declare (explicit-check))
 403   (number-dispatch ((x float))
 404     ((single-float) (float-sign-bit-set-p x))
 405     ((double-float) (float-sign-bit-set-p x))))
 406
 407 (declaim (inline float-digits float-radix))
 408
 409 (defun float-digits (f)
 410   "Return a non-negative number of radix-b digits used in the
 411    representation of its argument."
 412   (declare (explicit-check))
 413   (number-dispatch ((f float))
 414     ((single-float) sb-vm:single-float-digits)
 415     ((double-float) sb-vm:double-float-digits)
 416     #+long-float
 417     ((long-float) sb-vm:long-float-digits)))
 418
 419 (defun float-radix (x)
 420   "Return (as an integer) the radix b of its floating-point argument."
 421   (declare (ignore x) (type float x))
 422   2)
 423
 424 \f
 425 ;;;; INTEGER-DECODE-FLOAT and DECODE-FLOAT
 426
 427 (defconstant-eqx float-decoding-error "Can't decode NaN or infinity: ~S."
 428   #'string=)
 429
 430 (declaim (maybe-inline integer-decode-single-float
 431                        integer-decode-double-float))
 432
 433 ;;; binary point is to the left of the 23 represented mantissa bits,
 434 ;;; and normal exponent min is -126, so -149 is the effective exponent
 435 ;;; of a subnormal in common-lisp terms.
 436 (defconstant subnormal-sfloat-exponent -149)
 437 ;; binary point is to the left of the 52 represented mantissa bits,
 438 ;; and normal exponent min is -1022, so -1074 is the effective exponent.
 439 (defconstant subnormal-dfloat-exponent -1074)
 440
 441 ;;; Handle the single-float case of INTEGER-DECODE-FLOAT. If an infinity or
 442 ;;; NaN, error.
 443 (defun integer-decode-single-float (x)
 444   (declare (single-float x))
 445   (let* ((bits (single-float-bits x))
 446          (frac (ldb sb-vm:single-float-significand-byte bits))
 447          (sign (if (minusp bits) -1 1))
 448          (exp (ldb sb-vm:single-float-exponent-byte bits)))
 449     (cond ((= exp 0)
 450            (values frac (if (= frac 0) 0 subnormal-sfloat-exponent) sign))
 451           ((> exp sb-vm:single-float-normal-exponent-max)
 452            (error float-decoding-error x))
 453           (t
 454            (values (logior sb-vm:single-float-hidden-bit frac)
 455                    (- exp sb-vm:single-float-bias sb-vm:single-float-digits)
 456                    sign)))))
 457
 458 ;;; like INTEGER-DECODE-SINGLE-FLOAT, only doubly so
 459 (defun integer-decode-double-float (x)
 460   (declare (double-float x))
 461   #-64-bit ; treat high and low bits separately until the end
 462   (let* ((hi (double-float-high-bits x))
 463          (sign (if (minusp hi) -1 1))
 464          (lo (double-float-low-bits x))
 465          (mantissa (logior (ash (ldb sb-vm:double-float-hi-significand-byte hi) 32) lo))
 466          (exp (ldb sb-vm:double-float-hi-exponent-byte hi)))
 467     (cond ((zerop (logior (ldb (byte 31 0) hi) lo))
 468            (values 0 0 sign))
 469           ((< exp sb-vm:double-float-normal-exponent-min)
 470            (values mantissa subnormal-dfloat-exponent sign))
 471           ((> exp sb-vm:double-float-normal-exponent-max)
 472            (error float-decoding-error x))
 473           (t
 474            (values (logior sb-vm:double-float-hidden-bit mantissa)
 475                    (- exp sb-vm:double-float-bias sb-vm:double-float-digits)
 476                    sign))))
 477   #+64-bit ; don't split the high and low bits
 478   (let* ((bits (double-float-bits x))
 479          (frac (ldb sb-vm:double-float-significand-byte bits))
 480          (sign (if (minusp bits) -1 1))
 481          (exp (ldb sb-vm:double-float-exponent-byte bits)))
 482     (cond ((= exp 0)
 483            (values frac (if (= frac 0) 0 subnormal-dfloat-exponent) sign))
 484           ((> exp sb-vm:double-float-normal-exponent-max)
 485            (error float-decoding-error x))
 486           (t
 487            (values (logior sb-vm:double-float-hidden-bit frac)
 488                    (- exp sb-vm:double-float-bias sb-vm:double-float-digits)
 489                    sign)))))
 490
 491 ;;; Dispatch to the correct type-specific i-d-f function.
 492 (defun integer-decode-float (x)
 493   "Return three values:
 494    1) an integer representation of the significand.
 495    2) the exponent for the power of 2 that the significand must be multiplied
 496       by to get the actual value. This differs from the DECODE-FLOAT exponent
 497       by FLOAT-DIGITS, since the significand has been scaled to have all its
 498       digits before the radix point.
 499    3) -1 or 1 (i.e. the sign of the argument.)"
 500   (declare (explicit-check))
 501   (number-dispatch ((x float))
 502     ((single-float)
 503      (integer-decode-single-float x))
 504     ((double-float)
 505      (integer-decode-double-float x))))
 506
 507 ;;; Handle the single-float case of DECODE-FLOAT. If an infinity or NaN,
 508 ;;; error. For subnormals, we left-align the significant bits into a field
 509 ;;; that is FLOAT-DIGITS wide, and decrease the exponent.
 510 (defun decode-single-float (x)
 511   (declare (single-float x))
 512   (let* ((bits (single-float-bits x))
 513          (biased-exp (ldb sb-vm:single-float-exponent-byte bits)))
 514     (if (> biased-exp sb-vm:single-float-normal-exponent-max)
 515         (error float-decoding-error x)
 516         (let ((frac (ldb sb-vm:single-float-significand-byte bits)))
 517           (multiple-value-bind (new-exp new-frac lisp-exponent)
 518               (cond ((/= biased-exp 0) ; normal
 519                      ;; SINGLE-FLOAT-BIAS as the stored exponent yields
 520                      ;; an effective exponent of -1.
 521                      (values sb-vm:single-float-bias frac
 522                              (- biased-exp sb-vm:single-float-bias)))
 523                     ((= frac 0) (values 0 0 0))
 524                     (t ; subnormal. Normalize it and unset the implied 1 bit
 525                      (let ((prec (integer-length frac)))
 526                        (values sb-vm:single-float-bias
 527                                (ldb (byte (1- sb-vm:single-float-digits) 0)
 528                                     (ash frac (- sb-vm:single-float-digits prec)))
 529                                (+ subnormal-sfloat-exponent prec)))))
 530             (values (make-single-float (dpb new-exp sb-vm:single-float-exponent-byte new-frac))
 531                     lisp-exponent
 532                     (float-sign x)))))))
 533
 534 ;;; The double-float logic mostly follows the skeleton of the above code,
 535 ;;; but there is a consed bignum or two on 32-bit architectures.
 536 ;;; Consing for the sake of code clarity is worth it as far as I'm concerned.
 537 (defun decode-double-float (x)
 538   (declare (double-float x))
 539   (let* #+64-bit ((bits (double-float-bits x))
 540                   (biased-exp (ldb sb-vm:double-float-exponent-byte bits)))
 541         #-64-bit ((high (double-float-high-bits x))
 542                   (biased-exp (ldb sb-vm:double-float-hi-exponent-byte high)))
 543     (if (> biased-exp sb-vm:double-float-normal-exponent-max)
 544         (error float-decoding-error x)
 545         (let ((frac #+64-bit (ldb sb-vm:double-float-significand-byte bits)
 546                     #-64-bit (logior (ash (ldb sb-vm:double-float-hi-significand-byte high) 32)
 547                                      (double-float-low-bits x))))
 548           (multiple-value-bind (new-exp new-frac lisp-exponent)
 549               (cond ((/= biased-exp 0)  ; normal
 550                      ;; DOUBLE-FLOAT-BIAS as the stored exponent yields
 551                      ;; an effective exponent of -1.
 552                      (values sb-vm:double-float-bias frac
 553                              (- biased-exp sb-vm:double-float-bias)))
 554                     ((= frac 0) (values 0 0 0))
 555                     (t ; subnormal. Normalize it and unset the implied 1 bit
 556                      (let ((prec (integer-length frac)))
 557                        (values sb-vm:double-float-bias
 558                                (ldb (byte (1- sb-vm:double-float-digits) 0)
 559                                     (ash frac (- sb-vm:double-float-digits prec)))
 560                                (+ subnormal-dfloat-exponent prec)))))
 561             (values #-64-bit
 562                     (make-double-float (dpb new-exp sb-vm:double-float-hi-exponent-byte
 563                                             (ldb (byte 32 32) new-frac))
 564                                        (ldb (byte 32 0) new-frac))
 565                     #+64-bit
 566                     (%make-double-float (dpb new-exp sb-vm:double-float-exponent-byte new-frac))
 567                     lisp-exponent
 568                     (float-sign x)))))))
 569
 570 ;;; Dispatch to the appropriate type-specific function.
 571 (defun decode-float (f)
 572   "Return three values:
 573    1) a floating-point number representing the significand. This is always
 574       between 0.5 (inclusive) and 1.0 (exclusive).
 575    2) an integer representing the exponent.
 576    3) -1.0 or 1.0 (i.e. the sign of the argument.)"
 577   (declare (explicit-check))
 578   (number-dispatch ((f float))
 579     ((single-float)
 580      (decode-single-float f))
 581     ((double-float)
 582      (decode-double-float f))))
 583
 584 \f
 585 ;;;; SCALE-FLOAT
 586
 587 (declaim (maybe-inline scale-single-float scale-double-float))
 588
 589 ;;; Scale a single or double float, calling the correct over/underflow
 590 ;;; functions.
 591 (defun scale-single-float (x exp)
 592   (declare (single-float x) (integer exp))
 593   (etypecase exp
 594     (fixnum
 595      (let* ((bits (single-float-bits x))
 596             (old-exp (ldb sb-vm:single-float-exponent-byte bits))
 597             (new-exp (+ old-exp exp)))
 598        (cond
 599          ((zerop x) x)
 600          ((or (< old-exp sb-vm:single-float-normal-exponent-min)
 601               (< new-exp sb-vm:single-float-normal-exponent-min))
 602           (scale-single-float-maybe-underflow x exp))
 603          ((or (> old-exp sb-vm:single-float-normal-exponent-max)
 604               (> new-exp sb-vm:single-float-normal-exponent-max))
 605           (scale-single-float-maybe-overflow x exp))
 606          (t
 607           (make-single-float (dpb new-exp
 608                                   sb-vm:single-float-exponent-byte
 609                                   bits))))))
 610     (unsigned-byte (scale-single-float-maybe-overflow x exp))
 611     ((integer * 0) (scale-single-float-maybe-underflow x exp))))
 612
 613 (defun scale-double-float (x exp)
 614   (etypecase exp
 615     (fixnum
 616      #+64-bit
 617      (let* ((bits (double-float-bits x))
 618             (old-exp (ldb sb-vm:double-float-exponent-byte bits))
 619             (new-exp (+ old-exp exp)))
 620        (cond
 621          ((zerop x) x)
 622          ((or (< old-exp sb-vm:double-float-normal-exponent-min)
 623               (< new-exp sb-vm:double-float-normal-exponent-min))
 624           (scale-double-float-maybe-underflow x exp))
 625          ((or (> old-exp sb-vm:double-float-normal-exponent-max)
 626               (> new-exp sb-vm:double-float-normal-exponent-max))
 627           (scale-double-float-maybe-overflow x exp))
 628          (t
 629           (%make-double-float (dpb new-exp sb-vm:double-float-exponent-byte bits)))))
 630      #-64-bit
 631      (let* ((hi (double-float-high-bits x))
 632             (lo (double-float-low-bits x))
 633             (old-exp (ldb sb-vm:double-float-hi-exponent-byte hi))
 634             (new-exp (+ old-exp exp)))
 635        (cond
 636          ((zerop x) x)
 637          ((or (< old-exp sb-vm:double-float-normal-exponent-min)
 638               (< new-exp sb-vm:double-float-normal-exponent-min))
 639           (scale-double-float-maybe-underflow x exp))
 640          ((or (> old-exp sb-vm:double-float-normal-exponent-max)
 641               (> new-exp sb-vm:double-float-normal-exponent-max))
 642           (scale-double-float-maybe-overflow x exp))
 643          (t
 644           (make-double-float (dpb new-exp sb-vm:double-float-hi-exponent-byte hi)
 645                              lo)))))
 646     (unsigned-byte (scale-double-float-maybe-overflow x exp))
 647     ((integer * 0) (scale-double-float-maybe-underflow x exp))))
 648
 649 ;;; Dispatch to the correct type-specific scale-float function.
 650 (defun scale-float (f ex)
 651   "Return the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
 652   of precision or overflow."
 653   (declare (explicit-check f))
 654   (number-dispatch ((f float))
 655     ((single-float)
 656      (scale-single-float f ex))
 657     ((double-float)
 658      (scale-double-float f ex))))
 659 \f
 660 ;;;; converting to/from floats
 661
 662 (defun float (number &optional (other () otherp))
 663   "Converts any REAL to a float. If OTHER is not provided, it returns a
 664   SINGLE-FLOAT if NUMBER is not already a FLOAT. If OTHER is provided, the
 665   result is the same float format as OTHER."
 666   (declare (explicit-check))
 667   (if otherp
 668       (number-dispatch ((number real) (other float))
 669         (((foreach rational single-float double-float #+long-float long-float)
 670           (foreach single-float double-float #+long-float long-float))
 671          (coerce number '(dispatch-type other))))
 672       (if (floatp number)
 673           number
 674           (coerce number 'single-float))))
 675
 676 (macrolet ((frob (name type)
 677              `(defun ,name (x)
 678                 (declare (explicit-check x))
 679                 (number-dispatch ((x real))
 680                   (((foreach single-float double-float #+long-float long-float
 681                      sb-vm:signed-word
 682                      #-sb-xc-host
 683                      ,@(and (sb-c::template-translates-arg-p '%double-float 0 'word)
 684                             '(word))))
 685                    (coerce x ',type))
 686                   ((ratio)
 687                    (,(symbolicate type '-ratio) x))
 688                   ((bignum)
 689                    (,(symbolicate 'bignum-to- type) x))))))
 690   (frob %single-float single-float)
 691   (frob %double-float double-float)
 692   #+long-float
 693   (frob %long-float long-float))
 694
 695 ;;; Convert a ratio to a float. We avoid any rounding error by doing an
 696 ;;; integer division. Accuracy is important to preserve print-read
 697 ;;; consistency, since this is ultimately how the reader reads a float. We
 698 ;;; scale the numerator by a power of two until the division results in the
 699 ;;; desired number of fraction bits, then do round-to-nearest.
 700 (macrolet ((def (format)
 701              `(defun ,(symbolicate format '-ratio) (x)
 702                 (let* ((signed-num (numerator x))
 703                        (plusp (plusp signed-num))
 704                        (num (if plusp signed-num (- signed-num)))
 705                        (den (denominator x))
 706                        (digits ,(package-symbolicate :sb-vm format '-digits))
 707                        (scale 0))
 708                   (declare (fixnum digits scale))
 709                   ;; Strip any trailing zeros from the denominator and move it into the scale
 710                   ;; factor (to minimize the size of the operands.)
 711                   (let ((den-twos (1- (integer-length (logxor den (1- den))))))
 712                     (declare (fixnum den-twos))
 713                     (decf scale den-twos)
 714                     (setq den (ash den (- den-twos))))
 715                   ;; Guess how much we need to scale by from the magnitudes of the numerator
 716                   ;; and denominator. We want one extra bit for a guard bit.
 717                   (let* ((num-len (integer-length num))
 718                          (den-len (integer-length den))
 719                          (delta (- den-len num-len))
 720                          (shift (1+ (the fixnum (+ delta digits))))
 721                          (shifted-num (ash num shift)))
 722                     (declare (fixnum delta shift))
 723                     (decf scale delta)
 724                     (labels ((float-and-scale (bits)
 725                                (let* ((bits (ash bits -1))
 726                                       (len (integer-length bits)))
 727                                  (cond ((> len digits)
 728                                         (aver (= len (the fixnum (1+ digits))))
 729                                         (scale-float (floatit (ash bits -1)) (1+ scale)))
 730                                        (t
 731                                         (scale-float (floatit bits) scale)))))
 732                              (floatit (bits)
 733                                (let ((sign (if plusp 0 1)))
 734                                  ,(case format
 735                                     (single-float
 736                                      `(single-from-bits sign sb-vm:single-float-bias bits))
 737                                     (double-float
 738                                      `(double-from-bits sign sb-vm:double-float-bias bits))
 739                                     #+long-float
 740                                     (long-float
 741                                      `(long-from-bits sign sb-vm:long-float-bias bits))))))
 742                       (declare (inline floatit))
 743                       (loop
 744                        (multiple-value-bind (fraction-and-guard rem)
 745                            (truncate shifted-num den)
 746                          (let ((extra (- (integer-length fraction-and-guard) digits)))
 747                            (declare (fixnum extra))
 748                            (cond ((/= extra 1)
 749                                   (aver (> extra 1)))
 750                                  ((oddp fraction-and-guard)
 751                                   (return
 752                                     (if (zerop rem)
 753                                         (float-and-scale
 754                                          (if (zerop (logand fraction-and-guard 2))
 755                                              fraction-and-guard
 756                                              (1+ fraction-and-guard)))
 757                                         (float-and-scale (1+ fraction-and-guard)))))
 758                                  (t
 759                                   (return (float-and-scale fraction-and-guard)))))
 760                          (setq shifted-num (ash shifted-num -1))
 761                          (incf scale)))))))))
 762   (def double-float)
 763   (def single-float))
 764
 765 ;;; This function is called when we are doing a truncate without any funky
 766 ;;; divisor, i.e. converting a float or ratio to an integer. Note that we do
 767 ;;; *not* return the second value of truncate, so it must be computed by the
 768 ;;; caller if needed.
 769 ;;;
 770 ;;; In the float case, we pick off small arguments so that compiler
 771 ;;; can use special-case operations.
 772 (defun %unary-truncate (number)
 773   (declare (explicit-check number))
 774   (macrolet ((fits-fixnum (type)
 775                `(sb-xc:<= ,(symbol-value (symbolicate 'most-negative-fixnum- type))
 776                           number
 777                           ,(symbol-value (symbolicate 'most-positive-fixnum- type))))
 778              (shift (type integer count)
 779                `(,(case type
 780                     #-64-bit
 781                     (double-float 'ash)
 782                     (t 'bignum-ashift-left-fixnum))
 783                  ,integer ,count)))
 784     (number-dispatch ((number real))
 785       ((integer) number)
 786       ((ratio) (values (truncate (numerator number) (denominator number))))
 787       (((foreach single-float double-float #+long-float long-float))
 788        (if (fits-fixnum (dispatch-type number))
 789            (truly-the fixnum (%unary-truncate number))
 790            (multiple-value-bind (bits exp sign) (integer-decode-float number)
 791              (shift (dispatch-type number)
 792                     (if (minusp sign)
 793                         (- bits)
 794                         bits)
 795                     exp)))))))
 796
 797 ;;; Produce both values, unlike %unary-truncate
 798 (defun unary-truncate (number)
 799   (declare (explicit-check number))
 800   (macrolet ((fits-fixnum (type)
 801                `(sb-xc:<= ,(symbol-value (symbolicate 'most-negative-fixnum- type))
 802                           number
 803                           ,(symbol-value (symbolicate 'most-positive-fixnum- type))))
 804              (shift (type)
 805                (case type
 806                  #-64-bit
 807                  (double-float
 808                   ;; Shifting negatives right is different
 809                   `(let ((truncated (ash bits exp)))
 810                      (if (minusp sign)
 811                          (- truncated)
 812                          truncated)))
 813                  (t
 814                   `(bignum-ashift-left-fixnum (if (minusp sign)
 815                                                   (- bits)
 816                                                   bits)
 817                                               exp)))))
 818     (number-dispatch ((number real))
 819       ((integer) (values number 0))
 820       ((ratio)
 821        (let ((truncated (truncate (numerator number) (denominator number))))
 822          (values truncated
 823                  (- number truncated))))
 824       (((foreach single-float double-float #+long-float long-float))
 825        (if (fits-fixnum (dispatch-type number))
 826            (let* ((truncated (truly-the fixnum (%unary-truncate number))))
 827              (values truncated
 828                      (sb-xc:- number
 829                               (coerce truncated '(dispatch-type number)))))
 830            (multiple-value-bind (bits exp sign) (integer-decode-float number)
 831              (let ((truncated (shift (dispatch-type number))))
 832                (values
 833                 truncated
 834                 #+64-bit
 835                 (coerce 0 '(dispatch-type number))
 836                 #-64-bit
 837                 (if (eq '(dispatch-type number) 'single-float)
 838                     (coerce 0 '(dispatch-type number))
 839                     (sb-xc:- number (coerce truncated '(dispatch-type number))))))))))))
 840
 841 ;;; Specialized versions for floats.
 842 (macrolet ((def (type name)
 843              `(defun ,name (number)
 844                 (if (sb-xc:<= ,(symbol-value (symbolicate 'most-negative-fixnum- type))
 845                               number
 846                               ,(symbol-value (symbolicate 'most-positive-fixnum- type)))
 847                     (truly-the fixnum (,name number))
 848                     (error "bad float arg")))))
 849   (def single-float %unary-truncate/single-float)
 850   (def double-float %unary-truncate/double-float)
 851   #+long-float
 852   (def double-float %unary-truncate/long-float))
 853
 854 (defun %unary-ceiling (m)
 855   (values (ceiling m)))
 856 (defun %unary-floor (m)
 857   (values (floor m)))
 858
 859 ;;; Similar to %UNARY-TRUNCATE, but rounds to the nearest integer. If we
 860 ;;; can't use the round primitive, then we do our own round-to-nearest on the
 861 ;;; result of i-d-f. [Note that this rounding will really only happen
 862 ;;; with double floats on 32-bit architectures, where there are
 863 ;;; fractional floats past most-x-fixnum]
 864 (defun %unary-round (number)
 865   (declare (explicit-check))
 866   (macrolet ((fits-fixnum (type)
 867                `(sb-xc:<= ,(symbol-value (symbolicate 'most-negative-fixnum- type))
 868                           number
 869                           ,(symbol-value (symbolicate 'most-positive-fixnum- type)))))
 870     (number-dispatch ((number real))
 871       ((integer) number)
 872       ((ratio) (values (round (numerator number) (denominator number))))
 873       (((foreach single-float double-float #+long-float long-float))
 874        (if (fits-fixnum (dispatch-type number))
 875            (truly-the fixnum (%unary-round number))
 876            #+64-bit
 877            (multiple-value-bind (bits exp sign) (integer-decode-float number)
 878              (bignum-ashift-left-fixnum
 879               (if (minusp sign)
 880                   (- bits)
 881                   bits)
 882               exp))
 883            #-64-bit
 884            (multiple-value-bind (bits exp) (integer-decode-float number)
 885              (let* ((shifted (ash bits exp))
 886                     (rounded (if (minusp exp)
 887                                  (let ((fractional-bits (logand bits (lognot (ash -1 (- exp)))))
 888                                        (0.5bits (ash 1 (- -1 exp))))
 889                                    (cond
 890                                      ((> fractional-bits 0.5bits) (1+ shifted))
 891                                      ((< fractional-bits 0.5bits) shifted)
 892                                      (t (if (oddp shifted) (1+ shifted) shifted))))
 893                                  shifted)))
 894                (if (minusp number)
 895                    (- rounded)
 896                    rounded))))))))
 897
 898 #-round-float
 899 (defun %unary-fround (number)
 900   (number-dispatch ((number real))
 901     ((integer) (float number))
 902     ((ratio) (float (round (numerator number) (denominator number))))
 903     (((foreach single-float double-float)) (%unary-fround number))))
 904
 905 (defun rational (x)
 906   "RATIONAL produces a rational number for any real numeric argument. This is
 907   more efficient than RATIONALIZE, but it assumes that floating-point is
 908   completely accurate, giving a result that isn't as pretty."
 909   (declare (explicit-check))
 910   (number-dispatch ((x real))
 911     (((foreach single-float double-float #+long-float long-float))
 912      (multiple-value-bind (bits exp sign) (integer-decode-float x)
 913        (if (eql bits 0)
 914            0
 915            (let ((int (if (minusp sign) (- bits) bits)))
 916              (cond ((minusp exp)
 917                     ;; Instead of division (which also involves GCD)
 918                     ;; find the first set bit of the numerator and shift accordingly,
 919                     ;; as the denominator is a power of two.
 920                     (let* ((pexp (- exp))
 921                            (set (first-bit-set bits)))
 922                       (if (> pexp set)
 923                           (%make-ratio (ash int (- set))
 924                                        (let ((shift (- pexp set)))
 925                                          (if (< shift sb-vm:n-fixnum-bits)
 926                                              (ash 1 shift)
 927                                              (bignum-ashift-left-fixnum 1 shift))))
 928                           (ash int exp))))
 929                    (t
 930                     (ash int exp)))))))
 931     ((rational) x)))