From ...
Path: archiver1.google.com!news2.google.com!news.maxwell.syr.edu!uio.no!nntp.uio.no!not-for-mail
From: Erik Naggum
Newsgroups: comp.lang.lisp
Subject: Re: Floating point arithmetic (epsilon)
Date: 20 Jan 2004 06:14:12 +0000
Organization: Naggum Software, Oslo, Norway
Lines: 89
Message-ID: <3283568052702575KL2065E@naggum.no>
References:
Mime-Version: 1.0
Content-Type: text/plain; charset=iso-8859-1
Content-Transfer-Encoding: 8bit
X-Trace: readme.uio.no 1074579253 16740 129.240.65.210 (20 Jan 2004 06:14:13 GMT)
X-Complaints-To: abuse@uio.no
NNTP-Posting-Date: Tue, 20 Jan 2004 06:14:13 +0000 (UTC)
Mail-Copies-To: never
User-Agent: Gnus/5.09 (Gnus v5.9.0) Emacs/21.3
Xref: archiver1.google.com comp.lang.lisp:10346
* 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.