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