From ...
Path: archiver1.google.com!news2.google.com!news.maxwell.syr.edu!uio.no!nntp.uio.no!not-for-mail
From: Erik Naggum <erik@naggum.no>
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: <q67smibr8b0.fsf@L75001820.us.ray.com>
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.