Machine epsilon - find epsi

Question

4.2 Machine Epsilon

Find the floating point number epsi that has the the following properties:

1. 1.0+epsi is greater than 1.0 and
2. 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 (&OPTIONAL (epsi 1.0)) (if (> (+ 1.0 epsi) 1.0) (smaller-scale (/ epsi 10)) epsi))
(defun bigger (epsi inc-unit) (if (< (+ 1.0 epsi) 1.0) (bigger (+ epsi inc-unit) inc-unit) epsi))
(defun smaller (epsi dec-unit) (if (> (+ 1.0 epsi) 1.0) (smaller (+ epsi dec-unit) dec-unit) epsi))

(defun find-epsi (&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 (&OPTIONAL (epsi 1.0)) 
  (if (> (+ 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))

Answer 1

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)
  (> (+ 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 (&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 (&OPTIONAL (epsi 1.0))
  (if (is-epsi-p epsi)
    epsi
    (find-epsi (next-epsi epsi))))

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