src/code/float.lisp

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