Skip to content
Snippets Groups Projects
Commit cc753512 authored by liam's avatar liam
Browse files

Update documentation, include floating-point and auxiliary material.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3296 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 679dc967
No related branches found
No related tags found
No related merge requests found
<!-- -*- mode: HTML; time-stamp-line-limit: -8; -*- -->
<?xml version="1.0" encoding="iso-8859-1"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<link rel="stylesheet" href="style.css" type="text/css"/>
<title>Comparison of floating point numbers in Common Lisp</title>
</head>
<body>
<div class="body">
<div class="header">
<h1>Comparison of floating point numbers in Common Lisp</h1>
</div>
<p>
Floating point numbers are a computer representation of the real
numbers. Unlike the real numbers, there are a finite
number of them. So there is a smallest and largest floating point
number, and all others have a predecessor and successor.
<p>
Because different compilers and platforms can reorder a calculation
and optimize in a way that is approximated differently and so do
not necessarily produce the same floating point number, it is
difficult to compare two floating point numbers and conclude that
they represent the same result. For the purposes of regression (or
unit) testing, we would like to do exactly this. Bruce Dawson
addressed this problem in "<a
href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">Comparing
Floating Point Numbers</a>." He makes the point that the best way
to do this correctly is to interpret each floating point number as
an integer. By taking advantage of the <a
href="http://en.wikipedia.org/wiki/IEEE_floating-point_standard">IEEE
754</a> standard for representation of floating point numbers, we
can construct a function that maps the floating point numbers to
the integers. The genius of the standard's inventor <a
href="http://http.cs.berkeley.edu/~wkahan/">W. Kahan</a> is that a
mapping derived from the standard, call it i(x), satisfies three
properties:</p>
<ul>
<li>If two floats a&lt;b, then i(a)&lt;i(b),
<li>if two floats are adjacent and a&lt;b, then i(b)=i(a)+1,</li>
<li>and finally i(0.0)=0.</li>
</ul>
<p>
Dawson provides some clever C constructs to read a floating point
number as an integer, and instruction on how to prevent the
compiler from complaining about your trickery in doing so. Common
Lisp instead provides us with
functions with which we can properly construct our own integers.
As a side benefit, we don't care what the <i>actual</i>
representation of the floating point number is; we will build our
own IEEE-like representation. We don't exactly want the full
IEEE754 word though; we leave off the most significant
bit, which is a sign bit, and instead make the sign of the integer agree with the
sign of the float.
<p>What we end up with is an <i>enumeration of the floats</i>.
That is, for every single precision float, there is one integer in
the range
[-2139095039, 2139095039], and vice versa, with the exception that
both positive and negative zero (allowed by the standard) map to
zero. Likewise, there is a one-to-one mapping of the
double precision floats to
[-9218868437227405311,9218868437227405311]. Floats with
special values (positive and negative infinity, and <a
href="http://en.wikipedia.org/wiki/NaN">NaN</a>) that are required
by the standard do not have integer values.
<p>
The following functions written in Common Lisp are given:
</p>
<ul>
<li><code>float-as-integer</code> which is the function i(x);</li>
<li><code>integer-as-float</code> which is the inverse function
(this isn't necessary but can be useful) and
also returns the rational form of the float;</li>
<li><code>decode-IEEE754</code> (used by other functions) that
returns five values:
significand, exponent, sign, bits in significand, bits in
exponent, all as integers;</li>
<li><code>format-IEEE754-bits</code> which prints out the
binary form of the IEEE word, separated into the three parts
(this isn't necessary but is nice for comparing with bit
expansions shown in references like the Wikipedia page).</li>
</ul>
<p>
Here are some interesting floats:
<pre>
(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 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
</pre>
<p>
A regression test would record not the floating point number, but
the integer produced by <code>float-as-integer</code>. Since
integers can be unambiguously formatted to and read from a text
file in a unique way, a subsequent recomputation would provide a
clear indication of how close the floats are. Of course, we must
decide how much error we're going to allow, because a correct
calculation may produce slightly different integers. As an added
bonus, these functions can be used to identify (in languages other
than Lisp) when a positive single float has been interpreted as a
double float.
<p>
The floating point numbers computed for the regression tests here
were done with IEEE 754 double precision compatibility. If these
tests are run on a platform that does not follow this standard, the
tests may fail.
<div class="footer">
<hr>
<address><a href="mailto:">Liam Healy</a></address>
<small>
Time-stamp: <2008-01-26 14:11:38EST floating-point.html>
</small>
</div>
</body>
</html>
;; Comparison of floating point numbers
;; Liam Healy 2008-01-22 19:00:17EST floating-point.lisp
;; Time-stamp: <2008-01-29 21:34:52EST floating-point.lisp>
;; $Id: $
(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)))
(defun float-as-integer (float)
"An integer corresponding to the float which satisfies three properties:
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
sign of the integer agrees with the sign of the float."
(multiple-value-bind (mant exp sign sigbits expbits)
(decode-ieee754 float)
(declare (ignore expbits))
(* sign (+ (ash exp sigbits) mant))))
;; (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)
"Construct the floating point number from its integer representation.
Also return the number in rational form."
(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
(* (signum integer) ; preserve the sign
(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)))
;; IEEE 754 Modes and masks
;; Liam Healy 2008-01-29 21:35:50EST ieee-modes.lisp
;; Time-stamp: <2008-02-17 18:36:37EST ieee-modes.lisp>
;; $Id: $
(in-package :gsl)
(cffi:defcenum ieee-types
"IEEE 754 types, from /usr/include/gsl/gsl_ieee_utils.h"
(:NAN 1) :INF :NORMAL :DENORMAL :ZERO)
(cffi:defcenum ieee-precisions
"IEEE 754 precisions, from /usr/include/gsl/gsl_ieee_utils.h"
(:SINGLE-PRECISION 1) :DOUBLE-PRECISION :EXTENDED-PRECISION)
(cffi:defcenum ieee-rounding
"IEEE 754 rounding, from /usr/include/gsl/gsl_ieee_utils.h"
(:TO-NEAREST 1) :DOWN :UP :TO-ZERO)
(cffi:defcenum ieee-mask
"IEEE 754 mask, from /usr/include/gsl/gsl_ieee_utils.h"
(:INVALID 1) (:DENORMALIZED 2) (:DIVISION-BY-ZERO 4)
(:OVERFLOW 8) (:UNDERFLOW 16) (:ALL 31) (:INEXACT 32))
#+ieee-floating-point
(defmfun set-floating-point-modes (precision rounding exception-mask)
"gsl_ieee_set_mode"
(((cffi:foreign-enum-value 'ieee-precisions precision) :int)
((cffi:foreign-enum-value 'ieee-rounding rounding) :int)
((cffi:foreign-enum-value 'ieee-mask exception-mask) :int))
:documentation
"Set the IEEE 754 precision, rounding mode, and exception mask.")
;; Definition of GSLL system ;; Definition of GSLL system
;; Liam Healy ;; Liam Healy
;; Time-stamp: <2008-02-17 18:41:29EST gsll.asd> ;; Time-stamp: <2008-02-18 11:27:00EST gsll.asd>
;; $Id: $ ;; $Id: $
(asdf:defsystem "gsll" (asdf:defsystem "gsll"
...@@ -23,6 +23,10 @@ ...@@ -23,6 +23,10 @@
;; http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html ;; http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html
(:file "lisp-unit") (:file "lisp-unit")
(:file "tests" :depends-on (init lisp-unit)))) (:file "tests" :depends-on (init lisp-unit))))
(:module floating-point
:depends-on (init)
:components
((:file "ieee-modes")))
(:file "mathematical" :depends-on (init)) (:file "mathematical" :depends-on (init))
;; complex numbers not necessary? Just make a struct. ;; complex numbers not necessary? Just make a struct.
(:module data (:module data
......
...@@ -74,22 +74,22 @@ debugger invoked on a GSL-ERROR in thread #<THREAD "initial thread" {10032258C1} ...@@ -74,22 +74,22 @@ debugger invoked on a GSL-ERROR in thread #<THREAD "initial thread" {10032258C1}
to find the equivalent GSLL function, for example to find the equivalent GSLL function, for example
<pre> <pre>
(gsl-lookup "gsl_sf_elljac_e") (gsl-lookup "gsl_sf_elljac_e")
((JACOBIAN-ELLIPTIC-FUNCTIONS . "gsl_sf_elljac_e")) JACOBIAN-ELLIPTIC-FUNCTIONS
T
</pre> </pre>
to find that the Lisp function name is <code>#'jacobian-elliptic-functions</code>. to find that the Lisp function name is <code>#'jacobian-elliptic-functions</code>.
</li> </li>
<li>Look at the documentation for that Lisp function, e.g. <li>Look at the documentation for that Lisp function, e.g.
<pre>(documentation #'jacobian-elliptic-functions 'function) <pre>(documentation #'jacobian-elliptic-functions 'function)
"The Jacobian elliptic functions @math{sn(u|m)}, "The Jacobian elliptic functions sn(u|m),
@math{cn(u|m)}, @math{dn(u|m)} computed by descending Landen cn(u|m), dn(u|m) computed by descending Landen transformations."
transformations."
</pre> </pre>
to get an explanation of the arguments etc. to get an explanation of the arguments etc.
</li> </li>
<li>Look at the end of the Lisp source file, or in an "example" file <li>Look at the end of the Lisp source file, or in an "example" file
in the same directory, for examples. Often, the examples are placed in the same directory, for examples. Often, the examples are placed
in a <code>lisp-unit:define-test</code> form, and therefore show the answer in a <code>make-tests</code> form (which will be commented out with
expected.</li> a #| |# block).</li>
</ul> </ul>
It is advisable to look at the examples first for calculations that It is advisable to look at the examples first for calculations that
require more complex setup (generally, the later chapters in the GSL require more complex setup (generally, the later chapters in the GSL
...@@ -339,7 +339,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p> ...@@ -339,7 +339,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p>
</table> </table>
</div> </div>
<h3>Regression Tests</h3> <h3>Regression tests</h3>
<div class="content"> <div class="content">
<p>Regression tests using <a <p>Regression tests using <a
href="http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html">lisp-unit</a> href="http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html">lisp-unit</a>
...@@ -352,12 +352,10 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p> ...@@ -352,12 +352,10 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p>
<pre> <pre>
(lisp-unit:run-tests) (lisp-unit:run-tests)
</pre> </pre>
and expect a few failures on amd64 (four), none? on i386. Often, the and expect a two failures on amd64, none? on i386.
failures aren't really failures, just floating point numbers that
aren't close enough to pass.
</p> </p>
<h3>Related projects from others</h3> <h3>Similar or related projects from others</h3>
<div class="content"> <div class="content">
A partial automatically-generated interface to GSL is provided in <a A partial automatically-generated interface to GSL is provided in <a
href="http://common-lisp.net/project/cl-gsl/">cl-gsl</a>. href="http://common-lisp.net/project/cl-gsl/">cl-gsl</a>.
...@@ -374,7 +372,7 @@ and expect a few failures on amd64 (four), none? on i386. Often, the ...@@ -374,7 +372,7 @@ and expect a few failures on amd64 (four), none? on i386. Often, the
<!-- Created: Feb 25 2005 --> <!-- Created: Feb 25 2005 -->
<!-- hhmts start --> <!-- hhmts start -->
<small> <small>
Time-stamp: <2008-02-10 23:25:02EST index.html> Time-stamp: <2008-02-18 00:25:04EST index.html>
</small> </small>
<!-- hhmts end --> <!-- hhmts end -->
</div> </div>
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment