Machine epsilon - find epsi 4.2 Machine Epsilon
  
  Find the floating point number epsi that has
  the the following properties:
  
  
  1.0+epsi is greater than 1.0 and
  Let m b e  any number less than epsi. Then 1.0+m is equal to 1.0.
  
  
  epsi is called machine epsilon. It is of great importance in understanding floating point numbers.


I have written the following program to try and find my machine epsilon value to a given search depth (10).  Do you have any ideas how I could write this program better?

(defpackage :find-epsi (:use cl))
(in-package :find-epsi)

(defun smaller-scale (&amp;OPTIONAL (epsi 1.0)) (if (&gt; (+ 1.0 epsi) 1.0) (smaller-scale (/ epsi 10)) epsi))
(defun bigger (epsi inc-unit) (if (&lt; (+ 1.0 epsi) 1.0) (bigger (+ epsi inc-unit) inc-unit) epsi))
(defun smaller (epsi dec-unit) (if (&gt; (+ 1.0 epsi) 1.0) (smaller (+ epsi dec-unit) dec-unit) epsi))

(defun find-epsi (&amp;OPTIONAL (search-depth 10) (epsi (smaller-scale)) (incdec-unit epsi))
  (if (= search-depth 0) epsi (find-epsi (1- search-depth) (bigger (smaller epsi incdec-unit) incdec-unit) incdec-unit)))

(format t "epsi: ~a ~%" (find-epsi))


It seems that it should be much simpler to find epsilon than I originally thought.  What do you think about the following program?

(defpackage :find-epsi (:use cl))
(in-package :find-epsi)

(defun find-epsi (&amp;OPTIONAL (epsi 1.0)) 
  (if (&gt; (+ 1.0 epsi) 1.0)  ; if the variable epsi is still significant
    (find-epsi (/ epsi 2)) ; halve it and try again
    epsi)) ; otherwise, we have found epsilon

(format t "epsi: ~a ~%" (find-epsi))
If we assume that a float is represented in memory as a (sign, mantissa, exponent) tuple, and assume a radix of 2, then we can find the machine epsilon exactly.  That is, if we can assume the machine stores floats using base-2 representations of the mantissa and exponent, then we know that:


The machine will store a value of 1 in floating point exactly - this would be stored as 1 for the mantissa, and 0 for the exponent, i.e. 1 * 2^0.
The machine will store all powers of two that it can represent using a single bit in the mantissa, and by varying the exponent.  E.g. 1/4 could be represented as 1 * (2 ^ -2).  Any representable power of two will be stored without losing information.
1 + epsi will be the smallest value greater than 1 that can be stored in the mantissa of the floating-point number.


EDIT

The second version looks much better than the first, but I believe there's an off-by-one error in the number of times you recurse in find-epsi.  I suggest that you create a test function, to see if your result is the machine epsilon:

(defun epsi-sig-p (epsi)
  (&gt; (+ 1.0 epsi) 1.0))


You'll probably find that (is-sig (find-epsi)) is #f...  This also suggests that you can refactor (under the DRY principle) find-epsi to use this test function in find-epsi:

(defun find-epsi (&amp;OPTIONAL (epsi 1.0))
  (if (epsi-sig-p epsi)
    (find-epsi (/ epsi 2))
    epsi))


but we didn't change the behavior to fix the calculation, yet.  For this, I'd suggest another routine, to check whether we should try the next possible epsi:

(defun next-epsi (epsi) (/ epsi 2))
(defun is-epsi-p (epsi)
  (and (epsi-sig-p epsi) (not (epsi-sig-p (next-epsi epsi)))))
(defun find-epsi (&amp;OPTIONAL (epsi 1.0))
  (if (is-epsi-p epsi)
    epsi
    (find-epsi (next-epsi epsi))))


is-epsi-p should return #t, now.