src/code/float.lisp

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