From cc753512f287c1332a791c4ba6f5a50aa1b5d963 Mon Sep 17 00:00:00 2001
From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30>
Date: Mon, 18 Feb 2008 16:34:05 +0000
Subject: [PATCH] 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
---
 floating-point/floating-point.html | 132 +++++++++++++++++++++++++++++
 floating-point/floating-point.lisp | 131 ++++++++++++++++++++++++++++
 floating-point/ieee-modes.lisp     |  32 +++++++
 gsll.asd                           |   6 +-
 index.html                         |  22 +++--
 5 files changed, 310 insertions(+), 13 deletions(-)
 create mode 100644 floating-point/floating-point.html
 create mode 100644 floating-point/floating-point.lisp
 create mode 100644 floating-point/ieee-modes.lisp

diff --git a/floating-point/floating-point.html b/floating-point/floating-point.html
new file mode 100644
index 00000000..8bbef793
--- /dev/null
+++ b/floating-point/floating-point.html
@@ -0,0 +1,132 @@
+<!-- -*- 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>
diff --git a/floating-point/floating-point.lisp b/floating-point/floating-point.lisp
new file mode 100644
index 00000000..6bbd1a7d
--- /dev/null
+++ b/floating-point/floating-point.lisp
@@ -0,0 +1,131 @@
+;; 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)))
+
diff --git a/floating-point/ieee-modes.lisp b/floating-point/ieee-modes.lisp
new file mode 100644
index 00000000..f7c87852
--- /dev/null
+++ b/floating-point/ieee-modes.lisp
@@ -0,0 +1,32 @@
+;; 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.")
diff --git a/gsll.asd b/gsll.asd
index 29bb0b71..3c0b9731 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2008-02-17 18:41:29EST gsll.asd>
+;; Time-stamp: <2008-02-18 11:27:00EST gsll.asd>
 ;; $Id: $
 
 (asdf:defsystem "gsll"
@@ -23,6 +23,10 @@
 	     ;; http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html
 	     (:file "lisp-unit")
 	     (:file "tests" :depends-on (init lisp-unit))))
+   (:module floating-point
+	    :depends-on (init)
+	    :components
+	    ((:file "ieee-modes")))
    (:file "mathematical" :depends-on (init))
    ;; complex numbers not necessary?  Just make a struct.
    (:module data
diff --git a/index.html b/index.html
index be03c31d..b55304b1 100644
--- a/index.html
+++ b/index.html
@@ -74,22 +74,22 @@ debugger invoked on a GSL-ERROR in thread #<THREAD "initial thread" {10032258C1}
   to find the equivalent GSLL function, for example
 <pre>
 (gsl-lookup "gsl_sf_elljac_e")
-((JACOBIAN-ELLIPTIC-FUNCTIONS . "gsl_sf_elljac_e"))
+JACOBIAN-ELLIPTIC-FUNCTIONS
+T
 </pre>
 to find that the Lisp function name is <code>#'jacobian-elliptic-functions</code>.  
   </li>
  <li>Look at the documentation for that Lisp function, e.g.
    <pre>(documentation #'jacobian-elliptic-functions 'function)
-  "The Jacobian elliptic functions @math{sn(u|m)},
-  @math{cn(u|m)}, @math{dn(u|m)} computed by descending Landen
-  transformations."
+"The Jacobian elliptic functions sn(u|m),
+  cn(u|m), dn(u|m) computed by descending Landen transformations."
    </pre>
    to get an explanation of the arguments etc.
  </li>
  <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 a <code>lisp-unit:define-test</code> form, and therefore show the answer
- expected.</li>
+ in a <code>make-tests</code> form (which will be commented out with
+ a #| |# block).</li>
  </ul>
  It is advisable to look at the examples first for calculations that
  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>
     </table>
  </div>
 
-<h3>Regression Tests</h3>
+<h3>Regression tests</h3>
   <div class="content">
    <p>Regression tests using <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>
 <pre>
 (lisp-unit:run-tests)
 </pre>
-and expect a few failures on amd64 (four), none? on i386.  Often, the
-   failures aren't really failures, just floating point numbers that
-   aren't close enough to pass.
+and expect a two failures on amd64, none? on i386.
 </p>
 
- <h3>Related projects from others</h3>
+ <h3>Similar or related projects from others</h3>
   <div class="content">
   A partial automatically-generated interface to GSL is provided in <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
 <!-- Created: Feb 25 2005 -->
 <!-- hhmts start -->
     <small>
-       Time-stamp: <2008-02-10 23:25:02EST index.html>
+       Time-stamp: <2008-02-18 00:25:04EST index.html>
        </small>
 <!-- hhmts end -->
  </div>
-- 
GitLab