src/code/float.lisp

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