Newer
Older
;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
;;;
;;; **********************************************************************
;;; This code was written as part of the CMU Common Lisp project at
;;; Carnegie Mellon University, and has been placed in the public domain.
;;;
(ext:file-comment
"$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/float.lisp,v 1.44 2010/02/05 18:10:58 rtoy Exp $")
;;; **********************************************************************
;;;
;;; This file contains the definitions of float specific number support
;;; (other than irrational stuff, which is in irrat.) There is code in here
;;; that assumes there are only two float formats: IEEE single and double.
;;;
;;; Author: Rob MacLachlan
(export '(%unary-truncate %unary-round %unary-ftruncate
%unary-ftruncate/single-float %unary-ftruncate/double-float))
#-x87
(export '(%unary-fround/single-float %unary-fround/double-float))
(in-package "LISP")
(export '(least-positive-normalized-short-float
least-positive-normalized-single-float
least-positive-normalized-double-float
least-positive-normalized-long-float
least-negative-normalized-short-float
least-negative-normalized-single-float
least-negative-normalized-double-float
least-negative-normalized-long-float
least-positive-single-float
least-positive-short-float
least-negative-single-float
least-negative-short-float
least-positive-double-float
least-positive-long-float
least-negative-double-float
least-negative-long-float
most-positive-single-float
most-positive-short-float
most-negative-single-float
most-negative-short-float
most-positive-double-float
most-positive-long-float
most-negative-double-float
most-negative-long-float))
(in-package "EXTENSIONS")
(export '(single-float-positive-infinity short-float-positive-infinity
double-float-positive-infinity long-float-positive-infinity
single-float-negative-infinity short-float-negative-infinity
double-float-negative-infinity long-float-negative-infinity
set-floating-point-modes float-denormalized-p float-nan-p
float-trapping-nan-p float-infinity-p))
(in-package "KERNEL")
;;;; Utilities:
;;; SINGLE-FROM-BITS, DOUBLE-FROM-BITS -- Internal
;;;
;;; These functions let us create floats from bits with the significand
;;; uniformly represented as an integer. This is less efficient for double
;;; floats, but is more convenient when making special values, etc.
;;;
(defun single-from-bits (sign exp sig)
(declare (type bit sign) (type (unsigned-byte 24) sig)
(type (unsigned-byte 8) exp))
(make-single-float
(dpb exp vm:single-float-exponent-byte
(dpb sig vm:single-float-significand-byte
(if (zerop sign) 0 -1)))))
;;;
(defun double-from-bits (sign exp sig)
(declare (type bit sign) (type (unsigned-byte 53) sig)
(type (unsigned-byte 11) exp))
(make-double-float (dpb exp vm:double-float-exponent-byte
(dpb (ash sig -32) vm:double-float-significand-byte
(if (zerop sign) 0 -1)))
(ldb (byte 32 0) sig)))
#+double-double
(defun double-double-from-bits (sign exp sig)
(declare (type bit sign) (type (unsigned-byte #.vm:double-double-float-digits) sig)
(type (unsigned-byte 11) exp)
(ignore exp))
(let ((lo (ldb (byte vm:double-float-digits 0) sig))
(hi (ldb (byte vm:double-float-digits vm:double-float-digits) sig))
(fsign (- 1 (* 2 sign))))
;; We can't return a double-double because exp might be very
;; small, and the low part of the double-double would get scaled
;; to 0. So we return both parts and expect the caller to figure
;; out what to do.
(values (scale-float (* fsign (float hi 1d0)) #.(- vm:double-float-digits))
(scale-float (* fsign (float lo 1d0)) #.(- vm:double-float-digits)))))
;;;
#+(and long-float x86)
(defun long-from-bits (sign exp sig)
(declare (type bit sign) (type (unsigned-byte 64) sig)
(type (unsigned-byte 15) exp))
(make-long-float (logior (ash sign 15) exp)
(ldb (byte 32 32) sig)
(ldb (byte 32 0) sig)))
(defconstant least-positive-single-float (single-from-bits 0 0 1))
(defconstant least-positive-short-float least-positive-single-float)
(defconstant least-negative-single-float (single-from-bits 1 0 1))
(defconstant least-negative-short-float least-negative-single-float)
(defconstant least-positive-double-float (double-from-bits 0 0 1))
(defconstant least-positive-long-float least-positive-double-float)
#+(and long-float x86)
(defconstant least-positive-long-float (long-from-bits 0 0 1))
(defconstant least-negative-double-float (double-from-bits 1 0 1))
(defconstant least-negative-long-float least-negative-double-float)
#+(and long-float x86)
(defconstant least-negative-long-float (long-from-bits 1 0 1))
(single-from-bits 0 vm:single-float-normal-exponent-min 0))
(defconstant least-positive-normalized-short-float
least-positive-normalized-single-float)
(defconstant least-negative-normalized-single-float
(single-from-bits 1 vm:single-float-normal-exponent-min 0))
(defconstant least-negative-normalized-short-float
least-negative-normalized-single-float)
(defconstant least-positive-normalized-double-float
(double-from-bits 0 vm:double-float-normal-exponent-min 0))
(defconstant least-positive-normalized-long-float
least-positive-normalized-double-float)
#+(and long-float x86)
(defconstant least-positive-normalized-long-float
(long-from-bits 0 vm:long-float-normal-exponent-min
(ash vm:long-float-hidden-bit 32)))
(double-from-bits 1 vm:double-float-normal-exponent-min 0))
(defconstant least-negative-normalized-long-float
least-negative-normalized-double-float)
#+(and long-float x86)
(defconstant least-negative-normalized-long-float
(long-from-bits 1 vm:long-float-normal-exponent-min
(ash vm:long-float-hidden-bit 32)))
(single-from-bits 0 vm:single-float-normal-exponent-max
(ldb vm:single-float-significand-byte -1)))
(defconstant most-positive-short-float most-positive-single-float)
(defconstant most-negative-single-float
(single-from-bits 1 vm:single-float-normal-exponent-max
(ldb vm:single-float-significand-byte -1)))
(defconstant most-negative-short-float most-negative-single-float)
(defconstant most-positive-double-float
(double-from-bits 0 vm:double-float-normal-exponent-max
(ldb (byte vm:double-float-digits 0) -1)))
(defconstant most-positive-long-float most-positive-double-float)
#+(and long-float x86)
(defconstant most-positive-long-float
(long-from-bits 0 vm:long-float-normal-exponent-max
(ldb (byte vm:long-float-digits 0) -1)))
(double-from-bits 1 vm:double-float-normal-exponent-max
(ldb (byte vm:double-float-digits 0) -1)))
(defconstant most-negative-long-float most-negative-double-float)
#+(and long-float x86)
(defconstant most-negative-long-float
(long-from-bits 1 vm:long-float-normal-exponent-max
(ldb (byte vm:long-float-digits 0) -1)))
(single-from-bits 0 (1+ vm:single-float-normal-exponent-max) 0))
(defconstant short-float-positive-infinity single-float-positive-infinity)
(defconstant single-float-negative-infinity
(single-from-bits 1 (1+ vm:single-float-normal-exponent-max) 0))
(defconstant short-float-negative-infinity single-float-negative-infinity)
(defconstant double-float-positive-infinity
(double-from-bits 0 (1+ vm:double-float-normal-exponent-max) 0))
(defconstant long-float-positive-infinity double-float-positive-infinity)
#+(and long-float x86)
(defconstant long-float-positive-infinity
(long-from-bits 0 (1+ vm:long-float-normal-exponent-max)
(ash vm:long-float-hidden-bit 32)))
(double-from-bits 1 (1+ vm:double-float-normal-exponent-max) 0))
(defconstant long-float-negative-infinity double-float-negative-infinity)
#+(and long-float x86)
(defconstant long-float-negative-infinity
(long-from-bits 1 (1+ vm:long-float-normal-exponent-max)
(ash vm:long-float-hidden-bit 32)))
(single-from-bits 0 (- vm:single-float-bias (1- vm:single-float-digits)) 1))
(defconstant short-float-epsilon single-float-epsilon)
(defconstant single-float-negative-epsilon
(single-from-bits 0 (- vm:single-float-bias vm:single-float-digits) 1))
(defconstant short-float-negative-epsilon single-float-negative-epsilon)
#-(and long-float x86)
(double-from-bits 0 (- vm:double-float-bias (1- vm:double-float-digits)) 1))
#+(and long-float x86)
(defconstant double-float-epsilon
(double-from-bits 0 (- vm:double-float-bias (1- vm:double-float-digits))
(expt 2 42)))
#+(and long-float x86)
(defconstant long-float-epsilon
(long-from-bits 0 (- vm:long-float-bias (1- vm:long-float-digits))
(+ 1 (ash vm:long-float-hidden-bit 32))))
#-(and long-float x86)
(double-from-bits 0 (- vm:double-float-bias vm:double-float-digits) 1))
#+(and long-float x86)
(defconstant double-float-negative-epsilon
(double-from-bits 0 (- vm:double-float-bias vm:double-float-digits)
(expt 2 42)))
(defconstant long-float-negative-epsilon double-float-negative-epsilon)
#+(and long-float x86)
(defconstant long-float-negative-epsilon
(long-from-bits 0 (- vm:long-float-bias vm:long-float-digits)
(+ 1 (ash vm:long-float-hidden-bit 32))))
(declaim (maybe-inline float-denormalized-p float-infinity-p float-nan-p
float-trapping-nan-p))
;;; FLOAT-DENORMALIZED-P -- Public
;;;
(defun float-denormalized-p (x)
"Return true if the float X is denormalized."
(number-dispatch ((x float))
((single-float)
(and (zerop (ldb vm:single-float-exponent-byte (single-float-bits x)))
(and (zerop (ldb vm:double-float-exponent-byte
(double-float-high-bits x)))
(not (zerop x))))
#+(and long-float x86)
((long-float)
(and (zerop (ldb vm:long-float-exponent-byte (long-float-exp-bits x)))
(macrolet ((frob (name doc single double #+(and long-float x86) long
#+double-double double-double)
`(defun ,name (x)
,doc
(number-dispatch ((x float))
((single-float)
(let ((bits (single-float-bits x)))
(and (> (ldb vm:single-float-exponent-byte bits)
vm:single-float-normal-exponent-max)
,single)))
((double-float)
(let ((hi (double-float-high-bits x))
(lo (double-float-low-bits x)))
(declare (ignorable lo))
(and (> (ldb vm:double-float-exponent-byte hi)
vm:double-float-normal-exponent-max)
,double)))
#+(and long-float x86)
((long-float)
(let ((exp (long-float-exp-bits x))
(hi (long-float-high-bits x))
(lo (long-float-low-bits x)))
(declare (ignorable lo))
(and (> (ldb vm:long-float-exponent-byte exp)
vm:long-float-normal-exponent-max)
,long)))
#+double-double
((double-double-float)
,double-double)))))
(frob float-infinity-p "Return true if the float X is an infinity (+ or -)."
(zerop (ldb vm:single-float-significand-byte bits))
(and (zerop (ldb vm:double-float-significand-byte hi))
(zerop lo))
#+(and long-float x86)
(and (zerop (ldb vm:long-float-significand-byte hi))
(zerop lo))
#+double-double
(float-infinity-p (double-double-hi x)))
(frob float-nan-p "Return true if the float X is a NaN (Not a Number)."
(not (zerop (ldb vm:single-float-significand-byte bits)))
(or (not (zerop (ldb vm:double-float-significand-byte hi)))
(not (zerop lo)))
#+(and long-float x86)
(or (not (zerop (ldb vm:long-float-significand-byte hi)))
(not (zerop lo)))
#+double-double
(float-nan-p (double-double-hi x)))
(frob float-trapping-nan-p
"Return true if the float X is a trapping NaN (Not a Number)."
(zerop (logand (ldb vm:single-float-significand-byte bits)
vm:single-float-trapping-nan-bit))
(zerop (logand (ldb vm:double-float-significand-byte hi)
vm:double-float-trapping-nan-bit))
#+(and long-float x86)
(zerop (logand (ldb vm:long-float-significand-byte hi)
vm:long-float-trapping-nan-bit))
#+double-double
(float-trapping-nan-p (double-double-hi x))))
;;; FLOAT-PRECISION -- Public
;;;
;;; If denormalized, use a subfunction from INTEGER-DECODE-FLOAT to find the
;;; actual exponent (and hence how denormalized it is), otherwise we just
;;; return the number of digits or 0.
;;;
(defun float-precision (f)
"Returns a non-negative number of significant digits in it's float argument.
Will be less than FLOAT-DIGITS if denormalized or zero."
(macrolet ((frob (digits bias decode)
`(cond ((zerop f) 0)
((float-denormalized-p f)
(multiple-value-bind (ignore exp)
(,decode f)
(declare (ignore ignore))
(truly-the fixnum
(+ ,digits (1- ,digits) ,bias exp))))
(t
,digits))))
(number-dispatch ((f float))
((single-float)
(frob vm:single-float-digits vm:single-float-bias
(frob vm:double-float-digits vm:double-float-bias
integer-decode-double-denorm))
#+long-float
((long-float)
(frob vm:long-float-digits vm:long-float-bias
integer-decode-long-denorm))
#+double-double
((double-double-float)
;; What exactly is the precision for a double-double? We make
;; it the sum of the precisions of the two components.
(let ((hi (double-double-hi f)))
(cond ((zerop hi)
;; ANSI CL says the precision is 0
0)
((< (abs hi) (scale-float least-positive-normalized-double-float 53))
;; For every power of 2 below this, we lose a bit of
;; precision. More or less.
(let ((cutoff
(nth-value
1
(decode-float
(scale-float least-positive-normalized-double-float 53)))))
(multiple-value-bind (f exp)
(decode-float hi)
(declare (ignore f))
(- 106 (- cutoff exp)))))
(t
;; Normally we have 106 bits of precision (twice the
;; double-float precision)
106)))))))
(defun float-sign (float1 &optional (float2 (float 1 float1)))
"Returns a floating-point number that has the same sign as
float1 and, if float2 is given, has the same absolute value
as float2."
(declare (float float1 float2))
(* (if (etypecase float1
(single-float (minusp (single-float-bits float1)))
(double-float (minusp (double-float-high-bits float1)))
#+long-float
(long-float (minusp (long-float-exp-bits float1)))
#+double-double
(double-double-float (minusp (double-double-hi float1))))
(defun float-sign (float1 &optional float2)
"Returns a floating-point number that has the same sign as
float1 and, if float2 is given, has the same absolute value
as float2."
(declare (float float1)
(type (or null float) float2))
(let ((f1-sign (if (etypecase float1
(single-float (minusp (single-float-bits float1)))
(double-float (minusp (double-float-high-bits float1)))
#+long-float
(long-float (minusp (long-float-exp-bits float1)))
#+double-double
(double-double-float (minusp (float-sign (double-double-hi float1)))))
(float -1 float1)
(float 1 float1))))
;; Multiplication of double-double-float's doesn't preserve the
;; sign of signed-zeroes, so we split the case of float2 this way
;; so we can get the right sign.
(if float2
(if (minusp f1-sign)
(- (abs float2))
(abs float2))
f1-sign)))
((short-float single-float) vm:single-float-digits)
((double-float #-long-float long-float) vm:double-float-digits)
#+long-float
(long-float vm:long-float-digits)
#+double-double
(double-double-float vm:double-double-float-digits)))
(defun float-digits (f)
"Returns a non-negative number of radix-b digits used in the
representation of it's argument. See Common Lisp: The Language
by Guy Steele for more details."
(number-dispatch ((f float))
((single-float) vm:single-float-digits)
((long-float) vm:long-float-digits)
#+double-double
((double-double-float) vm:double-double-float-digits)))
(defun float-radix (f)
"Returns (as an integer) the radix b of its floating-point
argument."
(number-dispatch ((f float))
((float) 2)))
(declaim (maybe-inline integer-decode-single-float
integer-decode-double-float))
;;; INTEGER-DECODE-SINGLE-DENORM -- Internal
;;;
;;; Handle the denormalized case of INTEGER-DECODE-FLOAT for SINGLE-FLOAT.
;;;
(defun integer-decode-single-denorm (x)
(declare (type single-float x))
(let* ((bits (single-float-bits (abs x)))
(sig (ash (ldb vm:single-float-significand-byte bits) 1))
(extra-bias 0))
(declare (type (unsigned-byte 24) sig)
(type (integer 0 23) extra-bias))
(loop
(unless (zerop (logand sig vm:single-float-hidden-bit))
(return))
(setq sig (ash sig 1))
(incf extra-bias))
(values sig
(- (- vm:single-float-bias) vm:single-float-digits extra-bias)
(if (minusp (float-sign x)) -1 1))))
;;; INTEGER-DECODE-SINGLE-FLOAT -- Internal
;;;
;;; Handle the single-float case of INTEGER-DECODE-FLOAT. If an infinity or
;;; NAN, error. If a denorm, call i-d-s-DENORM to handle it.
;;;
(defun integer-decode-single-float (x)
(declare (single-float x))
(let* ((bits (single-float-bits (abs x)))
(exp (ldb vm:single-float-exponent-byte bits))
(sig (ldb vm:single-float-significand-byte bits))
(biased (- exp vm:single-float-bias vm:single-float-digits)))
(unless (<= exp vm:single-float-normal-exponent-max)
(error "Can't decode NAN or infinity: ~S." x))
((< exp vm:single-float-normal-exponent-min)
(values (logior sig vm:single-float-hidden-bit) biased sign)))))
;;; INTEGER-DECODE-DOUBLE-DENORM -- Internal
;;;
;;; Like INTEGER-DECODE-SINGLE-DENORM, only doubly so.
;;;
(defun integer-decode-double-denorm (x)
(declare (type double-float x))
(let* ((high-bits (double-float-high-bits (abs x)))
(sig-high (ldb vm:double-float-significand-byte high-bits))
(low-bits (double-float-low-bits x))
(sign (if (minusp (float-sign x)) -1 1))
(biased (- (- vm:double-float-bias) vm:double-float-digits)))
(extra-bias (- vm:double-float-digits 33))
(bit (ash 1 31)))
(declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
(loop
(unless (zerop (logand sig bit)) (return))
(setq sig (ash sig 1))
(incf extra-bias))
(values (ash sig (- vm:double-float-digits 32))
(truly-the fixnum (- biased extra-bias))
sign))
(let ((sig (ash sig-high 1))
(extra-bias 0))
(declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
(loop
(unless (zerop (logand sig vm:double-float-hidden-bit))
(return))
(setq sig (ash sig 1))
(incf extra-bias))
(values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
(truly-the fixnum (- biased extra-bias))
sign)))))
;;; INTEGER-DECODE-DOUBLE-FLOAT -- Internal
;;;
;;; Like INTEGER-DECODE-SINGLE-FLOAT, only doubly so.
;;;
(defun integer-decode-double-float (x)
(declare (double-float x))
(let* ((abs (abs x))
(hi (double-float-high-bits abs))
(lo (double-float-low-bits abs))
(exp (ldb vm:double-float-exponent-byte hi))
(sig (ldb vm:double-float-significand-byte hi))
(biased (- exp vm:double-float-bias vm:double-float-digits)))
(unless (<= exp vm:double-float-normal-exponent-max)
(error "Can't decode NAN or infinity: ~S." x))
(cond ((and (zerop exp) (zerop sig) (zerop lo))
(values 0 biased sign))
((< exp vm:double-float-normal-exponent-min)
(logior (ash (logior (ldb vm:double-float-significand-byte hi)
vm:double-float-hidden-bit)
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
;;; INTEGER-DECODE-LONG-DENORM -- Internal
;;;
#+(and long-float x86)
(defun integer-decode-long-denorm (x)
(declare (type long-float x))
(let* ((high-bits (long-float-high-bits (abs x)))
(sig-high (ldb vm:long-float-significand-byte high-bits))
(low-bits (long-float-low-bits x))
(sign (if (minusp (float-sign x)) -1 1))
(biased (- (- vm:long-float-bias) vm:long-float-digits)))
(if (zerop sig-high)
(let ((sig low-bits)
(extra-bias (- vm:long-float-digits 33))
(bit (ash 1 31)))
(declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
(loop
(unless (zerop (logand sig bit)) (return))
(setq sig (ash sig 1))
(incf extra-bias))
(values (ash sig (- vm:long-float-digits 32))
(truly-the fixnum (- biased extra-bias))
sign))
(let ((sig (ash sig-high 1))
(extra-bias 0))
(declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
(loop
(unless (zerop (logand sig vm:long-float-hidden-bit))
(return))
(setq sig (ash sig 1))
(incf extra-bias))
(values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
(truly-the fixnum (- biased extra-bias))
sign)))))
#+double-double
(defun integer-decode-double-double-float (x)
(declare (type double-double-float x))
(multiple-value-bind (hi-int hi-exp sign)
(integer-decode-float (double-double-hi x))
(if (zerop (double-double-lo x))
(values (ash hi-int 53) (- hi-exp 53) sign)
(multiple-value-bind (lo-int lo-exp lo-sign)
(integer-decode-float (double-double-lo x))
;; We have x = 2^e1*i1 + 2^e2*i2
;; = 2^e2*(2^(e1-e2)*i1 + i2)
;;
;; NOTE: The hi and lo parts could actually have different
;; signs, so we need to add the two parts together with the
;; right sign!
(values (+ (* (* sign lo-sign) lo-int)
(ash hi-int (- hi-exp lo-exp)))
lo-exp
sign)))))
;;; INTEGER-DECODE-LONG-FLOAT -- Internal
;;;
#+(and long-float x86)
(defun integer-decode-long-float (x)
(declare (long-float x))
(let* ((hi (long-float-high-bits x))
(lo (long-float-low-bits x))
(exp-bits (long-float-exp-bits x))
(exp (ldb vm:long-float-exponent-byte exp-bits))
(sign (if (minusp exp-bits) -1 1))
(biased (- exp vm:long-float-bias vm:long-float-digits)))
(declare (fixnum biased))
(unless (<= exp vm:long-float-normal-exponent-max)
(error "Can't decode NAN or infinity: ~S." x))
(cond ((and (zerop exp) (zerop hi) (zerop lo))
(values 0 biased sign))
((< exp vm:long-float-normal-exponent-min)
(integer-decode-long-denorm x))
(t
(values (logior (ash hi 32) lo) biased sign)))))
;;; INTEGER-DECODE-FLOAT -- Public
;;;
;;; Dispatch to the correct type-specific i-d-f function.
;;;
(defun integer-decode-float (x)
"Returns three values:
1) an integer representation of the significand.
2) the exponent for the power of 2 that the significand must be multiplied
by to get the actual value. This differs from the DECODE-FLOAT exponent
by FLOAT-DIGITS, since the significand has been scaled to have all its
digits before the radix point.
3) -1 or 1 (i.e. the sign of the argument.)"
(number-dispatch ((x float))
((single-float)
(integer-decode-single-float x))
((double-float)
(integer-decode-double-float x))
#+long-float
((long-float)
(integer-decode-long-float x))
#+double-double
((double-double-float)
(integer-decode-double-double-float x))))
(declaim (maybe-inline decode-single-float decode-double-float))
;;; DECODE-SINGLE-DENORM -- Internal
;;;
;;; Handle the denormalized case of DECODE-SINGLE-FLOAT. We call
;;; INTEGER-DECODE-SINGLE-DENORM and then make the result into a float.
;;;
(defun decode-single-denorm (x)
(declare (type single-float x))
(multiple-value-bind (sig exp sign)
(integer-decode-single-denorm x)
(values (make-single-float
(dpb sig vm:single-float-significand-byte
(dpb vm:single-float-bias vm:single-float-exponent-byte 0)))
(truly-the fixnum (+ exp vm:single-float-digits))
(float sign x))))
;;; DECODE-SINGLE-FLOAT -- Internal
;;;
;;; Handle the single-float case of DECODE-FLOAT. If an infinity or NAN,
;;; error. If a denorm, call d-s-DENORM to handle it.
;;;
(defun decode-single-float (x)
(declare (single-float x))
(let* ((bits (single-float-bits (abs x)))
(exp (ldb vm:single-float-exponent-byte bits))
(biased (truly-the single-float-exponent
(- exp vm:single-float-bias))))
(unless (<= exp vm:single-float-normal-exponent-max)
(error "Can't decode NAN or infinity: ~S." x))
((< exp vm:single-float-normal-exponent-min)
(dpb vm:single-float-bias
vm:single-float-exponent-byte
bits))
biased sign)))))
;;; DECODE-DOUBLE-DENORM -- Internal
;;;
;;; Like DECODE-SINGLE-DENORM, only doubly so.
;;;
(defun decode-double-denorm (x)
(declare (double-float x))
(multiple-value-bind (sig exp sign)
(integer-decode-double-denorm x)
(values (make-double-float
(dpb (logand (ash sig -32) (lognot vm:double-float-hidden-bit))
vm:double-float-significand-byte
(dpb vm:double-float-bias vm:double-float-exponent-byte 0))
(truly-the fixnum (+ exp vm:double-float-digits))
(float sign x))))
;;; DECODE-DOUBLE-FLOAT -- Public
;;;
;;; Like DECODE-SINGLE-FLOAT, only doubly so.
;;;
(defun decode-double-float (x)
(declare (double-float x))
(let* ((abs (abs x))
(hi (double-float-high-bits abs))
(lo (double-float-low-bits abs))
(exp (ldb vm:double-float-exponent-byte hi))
(biased (truly-the double-float-exponent
(- exp vm:double-float-bias))))
(unless (<= exp vm:double-float-normal-exponent-max)
(error "Can't decode NAN or infinity: ~S." x))
((< exp vm:double-float-normal-exponent-min)
(dpb vm:double-float-bias vm:double-float-exponent-byte hi)
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
;;; DECODE-LONG-DENORM -- Internal
;;;
#+(and long-float x86)
(defun decode-long-denorm (x)
(declare (long-float x))
(multiple-value-bind (sig exp sign)
(integer-decode-long-denorm x)
(values (make-long-float vm:long-float-bias (ash sig -32)
(ldb (byte 32 0) sig))
(truly-the fixnum (+ exp vm:long-float-digits))
(float sign x))))
;;; DECODE-LONG-FLOAT -- Public
;;;
#+(and long-float x86)
(defun decode-long-float (x)
(declare (long-float x))
(let* ((hi (long-float-high-bits x))
(lo (long-float-low-bits x))
(exp-bits (long-float-exp-bits x))
(exp (ldb vm:long-float-exponent-byte exp-bits))
(sign (if (minusp exp-bits) -1l0 1l0))
(biased (truly-the long-float-exponent (- exp vm:long-float-bias))))
(unless (<= exp vm:long-float-normal-exponent-max)
(error "Can't decode NAN or infinity: ~S." x))
(cond ((zerop x)
(values 0.0l0 biased sign))
((< exp vm:long-float-normal-exponent-min)
(decode-long-denorm x))
(t
(values (make-long-float
(dpb vm:long-float-bias vm:long-float-exponent-byte
exp-bits)
hi
lo)
biased sign)))))
;;; DECODE-DOUBLE-DOUBLE-FLOAT -- Public
#+double-double
(defun decode-double-double-float (x)
(declare (type double-double-float x))
(multiple-value-bind (hi-frac hi-exp sign)
(decode-float (double-double-hi x))
(values (make-double-double-float hi-frac
(scale-float (double-double-lo x) (- hi-exp)))
hi-exp
(coerce sign 'double-double-float))))
;;; DECODE-FLOAT -- Public
;;;
;;; Dispatch to the appropriate type-specific function.
;;;
(defun decode-float (f)
"Returns three values:
1) a floating-point number representing the significand. This is always
between 0.5 (inclusive) and 1.0 (exclusive).
2) an integer representing the exponent.
3) -1.0 or 1.0 (i.e. the sign of the argument.)"
(number-dispatch ((f float))
((single-float)
(decode-single-float f))
((double-float)
(decode-long-float f))
#+double-double
((double-double-float)
(decode-double-double-float f))))
(declaim (maybe-inline scale-single-float scale-double-float))
;;; SCALE-FLOAT-MAYBE-UNDERFLOW -- Internal
;;;
;;; Handle float scaling where the X is denormalized or the result is
;;; denormalized or underflows to 0.
;;;
(defun scale-float-maybe-underflow (x exp)
(multiple-value-bind (sig old-exp)
(integer-decode-float x)
(let* ((digits (float-digits x))
(new-exp (+ exp old-exp digits
(etypecase x
(single-float vm:single-float-bias)
(double-float vm:double-float-bias))))
(sign (if (minusp (float-sign x)) 1 0)))
(cond
((< new-exp
(etypecase x
(single-float vm:single-float-normal-exponent-min)
(double-float vm:double-float-normal-exponent-min)))
(when (vm:current-float-trap :inexact)
(error 'floating-point-inexact :operation 'scale-float
(when (vm:current-float-trap :underflow)
(error 'floating-point-underflow :operation 'scale-float
:operands (list x exp)))
;; Is it necessary to have this IF here? Is there any case
;; where (ash sig shift) won't return 0 when
;; shift < -(digits-1)?
(if (< shift (- (1- digits)))
(single-float (single-from-bits sign 0 0))
(double-float (double-from-bits sign 0 0)))
(etypecase x
(single-float (single-from-bits sign 0 (ash sig shift)))
(double-float (double-from-bits sign 0 (ash sig shift)))))))
(t
(etypecase x
(single-float (single-from-bits sign new-exp sig))
(double-float (double-from-bits sign new-exp sig))))))))
;;; SCALE-FLOAT-MAYBE-OVERFLOW -- Internal
;;;
;;; Called when scaling a float overflows, or the oringinal float was a NaN
;;; or infinity. If overflow errors are trapped, then error, otherwise return
;;; the appropriate infinity. If a NaN, signal or not as appropriate.
;;;
(defun scale-float-maybe-overflow (x exp)
(cond
((float-infinity-p x)
;; Infinity is infinity, no matter how small...
x)
((float-nan-p x)
(when (and (float-trapping-nan-p x)
(vm:current-float-trap :invalid))
(error 'floating-point-invalid-operation :operation 'scale-float
(when (vm:current-float-trap :overflow)
(error 'floating-point-overflow :operation 'scale-float
(when (vm:current-float-trap :inexact)
(error 'floating-point-inexact :operation 'scale-float
:operands (list x exp)))
(* (float-sign x)
(etypecase x
(single-float single-float-positive-infinity)
(double-float double-float-positive-infinity))))))
;;; SCALE-SINGLE-FLOAT, SCALE-DOUBLE-FLOAT -- Internal
;;;
;;; Scale a single or double float, calling the correct over/underflow
;;; functions.
;;;
(defun scale-single-float (x exp)
(declare (single-float x) (fixnum exp))
(let* ((bits (single-float-bits x))
(old-exp (ldb vm:single-float-exponent-byte bits))
((or (< old-exp vm:single-float-normal-exponent-min)
(< new-exp vm:single-float-normal-exponent-min))
((or (> old-exp vm:single-float-normal-exponent-max)
(> new-exp vm:single-float-normal-exponent-max))
(make-single-float (dpb new-exp vm:single-float-exponent-byte bits))))))
(declaim (maybe-inline scale-double-float))
(defun scale-double-float (x exp)
(declare (double-float x) (fixnum exp))
(let* ((hi (double-float-high-bits x))
(lo (double-float-low-bits x))
(old-exp (ldb vm:double-float-exponent-byte hi))
((or (< old-exp vm:double-float-normal-exponent-min)
(< new-exp vm:double-float-normal-exponent-min))
((or (> old-exp vm:double-float-normal-exponent-max)
(> new-exp vm:double-float-normal-exponent-max))
(make-double-float (dpb new-exp vm:double-float-exponent-byte hi)
#+(and x86 long-float)
(defun scale-long-float (x exp)
(declare (long-float x) (fixnum exp))
(scale-float x exp))
#+double-double
(defun scale-double-double-float (x exp)
(declare (type double-double-float x) (fixnum exp))
(let ((hi (double-double-hi x))
(lo (double-double-lo x)))
(make-double-double-float (scale-double-float hi exp)
(scale-double-float lo exp))))
;;; SCALE-FLOAT -- Public
;;;
;;; Dispatch to the correct type-specific scale-float function.
;;;
(defun scale-float (f ex)
"Returns the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
of precision or overflow."
(number-dispatch ((f float))
((single-float)
(scale-single-float f ex))
((double-float)
(scale-long-float f ex))
#+double-double
((double-double-float)
(scale-double-double-float f ex))))
;;;; Converting to/from floats:
(defun float (number &optional (other () otherp))
"Converts any REAL to a float. If OTHER is not provided, it returns a
SINGLE-FLOAT if NUMBER is not already a FLOAT. If OTHER is provided, the
result is the same float format as OTHER."
(if otherp
(number-dispatch ((number real) (other float))
(((foreach rational single-float double-float #+long-float long-float
#+double-double double-double-float)
(foreach single-float double-float #+long-float long-float
#+double-double double-double-float))
(coerce number '(dispatch-type other))))
(if (floatp number)
number
(coerce number 'single-float))))
(macrolet ((frob (name type)
`(defun ,name (x)
(number-dispatch ((x real))
(((foreach single-float double-float
#+long-float long-float
#+double-double double-double-float