Skip to content
Snippets Groups Projects
float.lisp 56.5 KiB
Newer Older
ram's avatar
ram committed
;;; -*- 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 $")
ram's avatar
ram committed
;;; **********************************************************************
;;;
;;;    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
dtc's avatar
dtc committed
;;; Long-float support by Douglas Crosher, 1998.
ram's avatar
ram committed
;;; 
(in-package "KERNEL")
(export '(%unary-truncate %unary-round %unary-ftruncate
	  %unary-ftruncate/single-float %unary-ftruncate/double-float))
ram's avatar
ram committed

#-x87
(export '(%unary-fround/single-float %unary-fround/double-float))

ram's avatar
ram committed
(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))
ram's avatar
ram committed

(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))
ram's avatar
ram committed

(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))
   (dpb exp vm:single-float-exponent-byte
	(dpb sig vm:single-float-significand-byte
	     (if (zerop sign) 0 -1)))))
ram's avatar
ram committed
;;;
(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)))))
dtc's avatar
dtc committed
;;;
#+(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)))
ram's avatar
ram committed
					

;;;; Float parameters:

(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))
dtc's avatar
dtc committed
#-long-float
(defconstant least-positive-long-float least-positive-double-float)
dtc's avatar
dtc committed
#+(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))
dtc's avatar
dtc committed
#-long-float
(defconstant least-negative-long-float least-negative-double-float)
dtc's avatar
dtc committed
#+(and long-float x86)
(defconstant least-negative-long-float (long-from-bits 1 0 1))
ram's avatar
ram committed

(defconstant least-positive-normalized-single-float
  (single-from-bits 0 vm:single-float-normal-exponent-min 0))
ram's avatar
ram committed
(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))
ram's avatar
ram committed
(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))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant least-positive-normalized-long-float
  least-positive-normalized-double-float)
dtc's avatar
dtc committed
#+(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)))
ram's avatar
ram committed
(defconstant least-negative-normalized-double-float
  (double-from-bits 1 vm:double-float-normal-exponent-min 0))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant least-negative-normalized-long-float
  least-negative-normalized-double-float)
dtc's avatar
dtc committed
#+(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)))
ram's avatar
ram committed

(defconstant most-positive-single-float
  (single-from-bits 0 vm:single-float-normal-exponent-max
		    (ldb vm:single-float-significand-byte -1)))
ram's avatar
ram committed
(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)))
ram's avatar
ram committed
(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)))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant most-positive-long-float most-positive-double-float)
dtc's avatar
dtc committed
#+(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)))
ram's avatar
ram committed
(defconstant most-negative-double-float
  (double-from-bits 1 vm:double-float-normal-exponent-max
		    (ldb (byte vm:double-float-digits 0) -1)))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant most-negative-long-float most-negative-double-float)
dtc's avatar
dtc committed
#+(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)))
ram's avatar
ram committed

(defconstant single-float-positive-infinity
  (single-from-bits 0 (1+ vm:single-float-normal-exponent-max) 0))
ram's avatar
ram committed
(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))
ram's avatar
ram committed
(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))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant long-float-positive-infinity double-float-positive-infinity)
dtc's avatar
dtc committed
#+(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)))
ram's avatar
ram committed
(defconstant double-float-negative-infinity
  (double-from-bits 1 (1+ vm:double-float-normal-exponent-max) 0))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant long-float-negative-infinity double-float-negative-infinity)
dtc's avatar
dtc committed
#+(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)))
ram's avatar
ram committed

(defconstant single-float-epsilon
  (single-from-bits 0 (- vm:single-float-bias (1- vm:single-float-digits)) 1))
ram's avatar
ram committed
(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))
ram's avatar
ram committed
(defconstant short-float-negative-epsilon single-float-negative-epsilon)
ram's avatar
ram committed
(defconstant double-float-epsilon
  (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)))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant long-float-epsilon double-float-epsilon)
dtc's avatar
dtc committed
#+(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))))
ram's avatar
ram committed
(defconstant double-float-negative-epsilon
  (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)))
dtc's avatar
dtc committed
#-long-float
ram's avatar
ram committed
(defconstant long-float-negative-epsilon double-float-negative-epsilon)
dtc's avatar
dtc committed
#+(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))))
ram's avatar
ram committed


;;;; Float predicates and environment query:

pw's avatar
pw committed
(declaim (maybe-inline float-denormalized-p float-infinity-p float-nan-p
ram's avatar
ram committed
			 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)))
ram's avatar
ram committed
	  (not (zerop x))))
    ((double-float)
     (and (zerop (ldb vm:double-float-exponent-byte
		      (double-float-high-bits x)))
dtc's avatar
dtc committed
	  (not (zerop x))))
    #+(and long-float x86)
    ((long-float)
     (and (zerop (ldb vm:long-float-exponent-byte (long-float-exp-bits x)))
ram's avatar
ram committed
	  (not (zerop x))))))

(macrolet ((frob (name doc single double #+(and long-float x86) long
		       #+double-double double-double)
ram's avatar
ram committed
	     `(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)
ram's avatar
ram committed
			  ,single)))
		  ((double-float)
		   (let ((hi (double-float-high-bits x))
			 (lo (double-float-low-bits x)))
		     (and (> (ldb vm:double-float-exponent-byte hi)
			     vm:double-float-normal-exponent-max)
dtc's avatar
dtc committed
			  ,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)))))
ram's avatar
ram committed

  (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))
dtc's avatar
dtc committed
	 (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)))
ram's avatar
ram committed

  (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)))
dtc's avatar
dtc committed
	(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)))
ram's avatar
ram committed

  (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)
dtc's avatar
dtc committed
		   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))))
ram's avatar
ram committed


;;; 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.
;;;
pw's avatar
pw committed
(declaim (maybe-inline float-precision))
ram's avatar
ram committed
(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
ram's avatar
ram committed
	 integer-decode-single-denorm))
      ((double-float)
       (frob vm:double-float-digits vm:double-float-bias
dtc's avatar
dtc committed
	 integer-decode-double-denorm))
      #+long-float
      ((long-float)
       (frob vm:long-float-digits vm:long-float-bias
rtoy's avatar
rtoy committed
	     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)))))))
ram's avatar
ram committed


ram's avatar
ram committed
(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))
dtc's avatar
dtc committed
  (* (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))))
dtc's avatar
dtc committed
	 (float -1 float1)
	 (float 1 float1))
     (abs float2)))
ram's avatar
ram committed

(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)))

ram's avatar
ram committed
(defun float-format-digits (format)
  (ecase format
    ((short-float single-float) vm:single-float-digits)
dtc's avatar
dtc committed
    ((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)))
ram's avatar
ram committed

pw's avatar
pw committed
(declaim (inline float-digits float-radix))
ram's avatar
ram committed

(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)
dtc's avatar
dtc committed
    ((double-float) vm:double-float-digits)
    #+long-float
    ((long-float) vm:long-float-digits)
    #+double-double
    ((double-double-float) vm:double-double-float-digits)))
ram's avatar
ram committed

(defun float-radix (f)
  "Returns (as an integer) the radix b of its floating-point
   argument."
  (number-dispatch ((f float))
    ((float) 2)))
ram's avatar
ram committed



;;;; INTEGER-DECODE-FLOAT and DECODE-FLOAT:

pw's avatar
pw committed
(declaim (maybe-inline integer-decode-single-float
		       integer-decode-double-float))
ram's avatar
ram committed

;;; 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))
ram's avatar
ram committed
	 (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))
ram's avatar
ram committed
	(return))
      (setq sig (ash sig 1))
      (incf extra-bias))
    (values sig
	    (- (- vm:single-float-bias) vm:single-float-digits extra-bias)
ram's avatar
ram committed
	    (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))
ram's avatar
ram committed
	 (sign (if (minusp (float-sign x)) -1 1))
	 (biased (- exp vm:single-float-bias vm:single-float-digits)))
ram's avatar
ram committed
    (declare (fixnum biased))
    (unless (<= exp vm:single-float-normal-exponent-max)
      (error "Can't decode NAN or infinity: ~S." x))
ram's avatar
ram committed
    (cond ((and (zerop exp) (zerop sig))
	   (values 0 biased sign))
	  ((< exp vm:single-float-normal-exponent-min)
ram's avatar
ram committed
	   (integer-decode-single-denorm x))
	  (t
	   (values (logior sig vm:single-float-hidden-bit) biased sign)))))
ram's avatar
ram committed


;;; 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))
ram's avatar
ram committed
	 (low-bits (double-float-low-bits x))
	 (sign (if (minusp (float-sign x)) -1 1))
	 (biased (- (- vm:double-float-bias) vm:double-float-digits)))
ram's avatar
ram committed
    (if (zerop sig-high)
	(let ((sig low-bits)
	      (extra-bias (- vm:double-float-digits 33))
ram's avatar
ram committed
	      (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))
ram's avatar
ram committed
		  (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))
ram's avatar
ram committed
	      (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))
ram's avatar
ram committed
	 (sign (if (minusp (float-sign x)) -1 1))
	 (biased (- exp vm:double-float-bias vm:double-float-digits)))
ram's avatar
ram committed
    (declare (fixnum biased))
    (unless (<= exp vm:double-float-normal-exponent-max)
      (error "Can't decode NAN or infinity: ~S." x))
ram's avatar
ram committed
    (cond ((and (zerop exp) (zerop sig) (zerop lo))
	   (values 0 biased sign))
	  ((< exp vm:double-float-normal-exponent-min)
ram's avatar
ram committed
	   (integer-decode-double-denorm x))
	  (t
	   (values
	    (logior (ash (logior (ldb vm:double-float-significand-byte hi)
				 vm:double-float-hidden-bit)
ram's avatar
ram committed
			 32)
		    lo)
	    biased sign)))))


dtc's avatar
dtc committed
;;; 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)))))
dtc's avatar
dtc committed

;;; 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)))))


ram's avatar
ram committed
;;; 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)
dtc's avatar
dtc committed
     (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))))
ram's avatar
ram committed


pw's avatar
pw committed
(declaim (maybe-inline decode-single-float decode-double-float))
ram's avatar
ram committed

;;; 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))
ram's avatar
ram committed
	    (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))
ram's avatar
ram committed
	 (sign (float-sign x))
	 (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))
ram's avatar
ram committed
    (cond ((zerop x)
	   (values 0.0f0 biased sign))
	  ((< exp vm:single-float-normal-exponent-min)
ram's avatar
ram committed
	   (decode-single-denorm x))
	  (t
	   (values (make-single-float
		    (dpb vm:single-float-bias
			 vm:single-float-exponent-byte
			 bits))
ram's avatar
ram committed
		   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))
ram's avatar
ram committed
	     (ldb (byte 32 0) sig))
	    (truly-the fixnum (+ exp vm:double-float-digits))
ram's avatar
ram committed
	    (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))
ram's avatar
ram committed
	 (sign (float-sign x))
	 (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))
ram's avatar
ram committed
    (cond ((zerop x)
	   (values 0.0d0 biased sign))
	  ((< exp vm:double-float-normal-exponent-min)
ram's avatar
ram committed
	   (decode-double-denorm x))
	  (t
	   (values (make-double-float
		    (dpb vm:double-float-bias vm:double-float-exponent-byte hi)
ram's avatar
ram committed
		    lo)
		   biased sign)))))


dtc's avatar
dtc committed
;;; 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))))
dtc's avatar
dtc committed

ram's avatar
ram committed
;;; 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)
dtc's avatar
dtc committed
     (decode-double-float f))
    #+long-float
    ((long-float)
     (decode-long-float f))
    #+double-double
    ((double-double-float)
     (decode-double-double-float f))))
ram's avatar
ram committed


;;;; SCALE-FLOAT:

pw's avatar
pw committed
(declaim (maybe-inline scale-single-float scale-double-float))
ram's avatar
ram committed

;;; 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))))
ram's avatar
ram committed
	   (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
ram's avatar
ram committed
		 :operands (list x exp)))
	(when (vm:current-float-trap :underflow)
	  (error 'floating-point-underflow :operation 'scale-float
		 :operands (list x exp)))
ram's avatar
ram committed
	(let ((shift (1- new-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)))))))
ram's avatar
ram committed
       (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
ram's avatar
ram committed
	     :operands (list x exp)))
    x)
   (t
    (when (vm:current-float-trap :overflow)
      (error 'floating-point-overflow :operation 'scale-float
ram's avatar
ram committed
	     :operands (list x exp)))
    (when (vm:current-float-trap :inexact)
      (error 'floating-point-inexact :operation 'scale-float
ram's avatar
ram committed
	     :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))
ram's avatar
ram committed
	 (new-exp (+ old-exp exp)))
    (cond
     ((zerop x) x)
     ((or (< old-exp vm:single-float-normal-exponent-min)
	  (< new-exp vm:single-float-normal-exponent-min))
ram's avatar
ram committed
      (scale-float-maybe-underflow x exp))
     ((or (> old-exp vm:single-float-normal-exponent-max)
	  (> new-exp vm:single-float-normal-exponent-max))
ram's avatar
ram committed
      (scale-float-maybe-overflow x exp))
     (t
      (make-single-float (dpb new-exp vm:single-float-exponent-byte bits))))))
ram's avatar
ram committed
;;;
(declaim (maybe-inline scale-double-float))
ram's avatar
ram committed
(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))
ram's avatar
ram committed
	 (new-exp (+ old-exp exp)))
    (cond
     ((zerop x) x)
     ((or (< old-exp vm:double-float-normal-exponent-min)
	  (< new-exp vm:double-float-normal-exponent-min))
ram's avatar
ram committed
      (scale-float-maybe-underflow x exp))
     ((or (> old-exp vm:double-float-normal-exponent-max)
	  (> new-exp vm:double-float-normal-exponent-max))
ram's avatar
ram committed
      (scale-float-maybe-overflow x exp))
     (t
      (make-double-float (dpb new-exp vm:double-float-exponent-byte hi)
ram's avatar
ram committed
			 lo)))))

dtc's avatar
dtc committed
#+(and x86 long-float)
(defun scale-long-float (x exp)
  (declare (long-float x) (fixnum exp))
  (scale-float x exp))
ram's avatar
ram committed

#+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))))

ram's avatar
ram committed
;;; 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)
dtc's avatar
dtc committed
     (scale-double-float f ex))
    #+long-float
    ((long-float)
     (scale-long-float f ex))
    #+double-double
    ((double-double-float)
     (scale-double-double-float f ex))))
ram's avatar
ram committed


;;;; 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))
ram's avatar
ram committed
	 (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
dtc's avatar
dtc committed
			     fixnum))