floating-point.lisp 5.68 KB
 lhealy committed Jul 24, 2008 1 ``````;; Comparison of floating point numbers using sequence integers `````` liam committed Feb 18, 2008 2 ``````;; Liam Healy 2008-01-22 19:00:17EST floating-point.lisp `````` lhealy committed Apr 04, 2008 3 ``````;; Time-stamp: <2008-04-04 11:35:20EDT floating-point.lisp> `````` lhealy committed Jul 24, 2008 4 ``````;; \$Id\$ `````` liam committed Feb 18, 2008 5 `````` `````` lhealy committed Jul 24, 2008 6 7 8 9 10 11 12 13 ``````;;; All floats can be represented by integers. There are two slightly ;;; different ways of doing this. The sequence integer ;;; has the same sign as the float; the IEEE754 integer is the ;;; unsigned version of the sequence integer, with the sign marked in ;;; the high order bit. The sequence integer is used for comparing ;;; numbers, where order preservation is important; the IEEE754 integer ;;; is used when reading a float storage location as an integer. `````` liam committed Feb 18, 2008 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 ``````(in-package :gsl) (export '(float-as-integer integer-as-float format-IEEE754-bits)) (defun decode-IEEE754 (float) "The significand (mantissa), exponent, and sign of the IEEE 754 representation of a floating point number, given as integers. It does not matter whether the actual representation follows IEEE. Values returned are significand, exponent, sign, bits in significand, bits in exponent." (flet ((esize (sigbits) (if (< sigbits 30) 8 11))) (if (zerop float) ;; zero must be treated specially; the whole word is 0 (let ((bits (1- (float-precision (float 1 float))))) (values 0 0 0 bits (esize bits))) (let* ((bits (1- (float-digits float))) (expbits (esize bits))) (multiple-value-bind (significand exponent sign) (integer-decode-float float) (if (< (float-precision float) (float-digits float)) ;; subnormalized number (values (ash significand ; shift right (- (float-precision float) (float-digits float))) 0 ; exponent is 0 for subnormalized numbers sign bits expbits) ;; normalized number (values (mask-field (byte bits 0) significand) ; leading "1" implicit (+ (1- (expt 2 (1- expbits))) bits exponent) ; bias exponent sign bits expbits))))))) (defun format-IEEE754-bits (float &optional (stream t)) "Format as binary each of the three pieces that make the IEEE 754 floating point representation for a float." (multiple-value-bind (mant exp sign sigbits expbits) (decode-IEEE754 float) (format stream "~1d ~v,'0b ~v,'0b" (if (minusp sign) 1 0) expbits exp sigbits mant))) `````` lhealy committed Jul 24, 2008 57 58 59 60 61 62 63 ``````(defun ieee754-sign-bit (float-type) "The bytespec to access the sign bit in the IEEE754 specification." (byte 1 (if (eq float-type 'double-float) 63 31))) (defun float-as-integer (float &optional ieee754) "The sequence integer corresponding to the float which satisfies three properties: `````` liam committed Feb 18, 2008 64 65 66 67 68 69 70 `````` 1) For two floats (< a b), then (< (float-as-integer a) (float-as-integer b)). 2) If two floats (< a b) are adjacent, then (= (1+ (float-as-integer a)) (float-as-integer b)). 3) (zerop (float-as-integer 0.0)) The absolute value of the integer is the integer of the IEEE754 representation without the sign bit, and the `````` lhealy committed Jul 24, 2008 71 72 73 `````` sign of the integer agrees with the sign of the float. To get the IEEE754 integer, specify 'single-float or 'double-float to ieee754." `````` liam committed Feb 18, 2008 74 75 76 `````` (multiple-value-bind (mant exp sign sigbits expbits) (decode-ieee754 float) (declare (ignore expbits)) `````` lhealy committed Apr 04, 2008 77 `````` (when (eq ieee754 t) (setf ieee754 (type-of float))) `````` lhealy committed Jul 24, 2008 78 79 80 `````` (if ieee754 (* (dpb 1 (ieee754-sign-bit ieee754) (+ (ash exp sigbits) mant))) (* sign (+ (ash exp sigbits) mant))))) `````` liam committed Feb 18, 2008 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 `````` ;; (float-as-integer most-negative-single-float) ;; -2139095039 ;; (float-as-integer least-negative-single-float) ;; -1 ;; (float-as-integer 0.0f0) ;; 0 ;; (float-as-integer least-positive-single-float) ;; 1 ;; (float-as-integer (- 1.0f0 single-float-negative-epsilon)) ;; 1065353215 ;; (float-as-integer 1.0f0) ;; 1065353216 ;; (float-as-integer (+ 1.0f0 single-float-epsilon)) ;; 1065353217 ;; (float-as-integer most-positive-single-float) ;; 2139095039 (defun integer-as-float (integer float-type) `````` lhealy committed Jul 24, 2008 100 101 102 103 104 105 `````` "Construct the floating point number from its integer representation, either sequence or IEEE754. Also return the number in rational form. The integer may be either a signed integer produced by float-as-integer, or an IEEE754 integer. Acceptable float-type is either 'double-float or 'single-float." `````` liam committed Feb 18, 2008 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 `````` (let ((expbits (if (eq float-type 'double-float) 11 8)) (sigbits (if (eq float-type 'double-float) 52 23))) (let* ((pinteger (abs integer)) (bexponent (ldb (byte expbits sigbits) pinteger)) (significand (ldb (byte sigbits 0) pinteger)) (normalizedp (plusp bexponent)) ; normalized number (frac (if normalizedp 1 0))) (dotimes (pos sigbits) ;; compute the fractional part (incf frac (if (logbitp (- sigbits pos 1) significand) (expt 2 (- (+ (if normalizedp 1 0) pos))) 0))) (let ((result `````` lhealy committed Jul 24, 2008 121 122 123 124 125 126 `````` (* (signum integer) ; accept negative integer (if ; accept IEEE754 negative (ldb-test (ieee754-sign-bit float-type) pinteger) -1 1) `````` liam committed Feb 18, 2008 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 `````` (expt 2 (- bexponent (1- (expt 2 (1- expbits))))) frac))) (values (coerce result float-type) result))))) ;; (integer-as-float -2139095039 'single-float) ;; -3.4028235e38 ;; -340282346638528859811704183484516925440 ;; (integer-as-float -1 'single-float) ;; -1.4012985e-45 ;; -1/713623846352979940529142984724747568191373312 ;; (integer-as-float 0 'single-float) ;; 0.0 ;; 0 ;; (integer-as-float 1 'single-float) ;; 1.4012985e-45 ;; 1/713623846352979940529142984724747568191373312 ;; (integer-as-float 1065353215 'single-float) ;; 0.99999994 ;; 16777215/16777216 ;; (integer-as-float 1065353216 'single-float) ;; 1.0 ;; 1 ;; (integer-as-float 1065353217 'single-float) ;; 1.0000001 ;; 8388609/8388608 ;; (integer-as-float 2139095039 'single-float) ;; 3.4028235e38 ;; 340282346638528859811704183484516925440 (defun next-float (float &optional (increment 1)) (integer-as-float (+ increment (float-as-integer float)) (type-of float))) ``````