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