```Subject: Re: Floating point arithmetic (epsilon)
From: Erik Naggum <erik@naggum.no>
Date: 20 Jan 2004 06:14:12 +0000
Newsgroups: comp.lang.lisp
Message-ID: <3283568052702575KL2065E@naggum.no>

* Steven E. Harris
| So 2.1597653E-7, which we'll call "delta," is greater than
| single-float(-negative)-epsilon, suggesting that delta is large
| enough to do useful arithmetic with on single-floats.

No, this is very wrong.  1.0±epsilon is the smallest representable
number different from 1.0.  To calculate the smallest representable
number different from any other value, you have to scale epsilon.

This should be instructive:

(rational (- 1.0 single-float-negative-epsilon))
=> 16777215/16777216

(rational (+ 1.0 single-float-epsilon))
=> 8388609/8388608

To find the epsilon of any given floating-point number, you have to
find the number that would be a change to the least significant bit
away from the number you look at.  This will be a scaled +epsilon when
the exponent does not change, and a scaled -epsilon when it does.

The function DECODE-FLOAT returns a number of type (float 1/b (1.0)),
where b is the floating point representation base, which you can
always find with (/ (decode-float 1.0)), which means that in order to
scale the epsilon correctly, you need to compute the scaling factor:

(defun scaled-epsilon (float)
(let ((scale (- (nth-value 1 (decode-float float))
(nth-value 1 (decode-float 1.0)))))
(scale-float single-float-epsilon scale)))

To show how this works, consider the following dialogue:

(scaled-epsilon 8.5)
=> 4.768372e-7

(integer-decode-float (- 8.5 *))
=> 8912895
=> -20
=> 1

(integer-decode-float (+ 8.5 **))
=> 8912897
=> -20
=> 1

(integer-decode-float 8.5)
=> 8912896
=> -20
=> 1

This takes care of +epsilon, but -epsilon is a different matter, since
it applies only when the value in question is exactly representable
with a single bit and the operation yields a smaller value.  Let me
therefore consider the full range:

(defun scaled-epsilon (float &optional (operation '+))
"Return the smallest number that would return a value different from
FLOAT if OPERATION were applied to FLOAT and this number.  OPERATION
should be either + or -, and defauls to +."
(multiple-value-bind (significand exponent)
(decode-float float)
(multiple-value-bind (1.0-significand 1.0-exponent)
(decode-float (float 1.0 float))
(if (and (eq operation '-)
(= significand 1.0-significand))
(scale-float (typecase float
(short-float short-float-negative-epsilon)
(single-float single-float-negative-epsilon)
(double-float double-float-negative-epsilon)
(long-float long-float-negative-epsilon))
(- exponent 1.0-exponent))
(scale-float (typecase float
(short-float short-float-epsilon)
(single-float single-float-epsilon)
(double-float double-float-epsilon)
(long-float long-float-epsilon))
(- exponent 1.0-exponent))))))

| I must be misunderstanding how to interpret these epsilon constants.

I hope this helps.

--
Erik Naggum | Oslo, Norway

Act from reason, and failure makes you rethink and study harder.
Act from faith, and failure makes you blame someone and push harder.

```