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