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