floating-point.lisp 5.68 KB
Newer Older
1
;; Comparison of floating point numbers using sequence integers
2
;; Liam Healy 2008-01-22 19:00:17EST floating-point.lisp
3
;; Time-stamp: <2008-04-04 11:35:20EDT floating-point.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5

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.

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

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: 
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
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."
74
75
76
  (multiple-value-bind (mant exp sign sigbits expbits)
      (decode-ieee754 float)
    (declare (ignore expbits))
77
    (when (eq ieee754 t) (setf ieee754 (type-of float)))
78
79
80
    (if ieee754
	(* (dpb 1 (ieee754-sign-bit ieee754) (+ (ash exp sigbits) mant)))
	(* sign (+ (ash exp sigbits) mant)))))
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)
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."
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
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)
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)))