Subject: Re: Floating point arithmetic (epsilon) From: Erik Naggum <email@example.com> 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.