Testing a floating point error agnostic = in Common Lisp

user_12421304 0
Viewed 117

For doing some testing of a system where lots of floating arithmetic is involved, I have defined a deviation range of floating arithmetic error, so if the difference between two floats lies within the deviation range they are considered as mathematically equal:

;;; Floating Error Agnostic =
;;; F2 is = to F1 within the deviation range +-DEV

(defparameter *flerag-devi* .0005
  "Allowed deviation range, within which two floats 
should be considered as mathematically equal.")

(defun flerag= (f1 f2 &optional (devi *flerag-devi*))
  ""
  (<= (abs (- f1 f2)) devi))

Now when I test the function by comparing a float with a random fraction added to it, the responses are (at least as for my testings of the function) positive, e.g.:

(loop repeat 100000000
      with f = 1.0
      always (flerag= f (+ f (random *flerag-devi*)))) ;T

This also is the case when I subtract a random fraction from the original float, e.g.:


(loop repeat 100000000
      with f = 19.0
      always (flerag= f (- f (random *flerag-devi*)))) ;T

(loop repeat 100000000
      with f = 3
      always (flerag= f (- f (random *flerag-devi*)))) ;T

However when I set the original float to be 1.0 (or 1):

(loop repeat 100000000
      with f = 1.0
      always (flerag= f (- f (random *flerag-devi*)))) ;NIL

it always returns NIL, which seems to me more strange as it returns NIL immediately after evaluation (in other cases it took some seconds to compute the 100000000 times). Until now this has happened only with f = 1.0 and no other numbers, regardless of addition or subtraction of the fraction to them. Can any one reproduce this behavior on his/her Lisp too? Any explanation and help would be appreciated.

Update The case with f = 1.0 also works as with the other ones when I use a double-float 1, i.e. 1d0 or by doing:

(setf *read-default-float-format* 'double-float)
common-lisp
1 Answers

Welcome to floating point arithmetic: a world where horrible traps await even the wary and the unwary are immediately eaten by grues. See this which you can get a PDF copy of here (this copy might be questionably legal, but people like Oracle have freely available versions as well) and it's kind of stupid that this paper is not freely available to everyone.

What is going on here is that, assuming single floats, at some point (random 0.0005) produces a number r (for instance 4.9999706E-4) that is close enough to 0.0005 that (- 1.0 r) is 0.9995. But (- 1.0 0.9995) is greater than 0.0005: in particular for single floats (and perhaps for doubles too, but definitely for singles) there exists single floats r such that

(and (< r 0.0005)
     (> (- 1.0 (- 1.0 r)) 0.0005)))

It is probably possible to systematically enumerate such single floats, but I got too lazy.

The definition of what it means to be 'close enough' also depends on what the semantics of what you are measuring is. If you're comparing diameters of objects you may not want to have a case arise where both x is 'close enough' to y and x >> y: 1.5E-18 is within 0.0005 of 1.5E-8, but if you were measuring the radius of a proton you'd be out by ten orders of magnitude, which is a not insignificant error.

A definition like the below (this is amended from a previous one which was buggy) might be suitable for that.

(defun close-enough-p (f1 f2 &optional (epsilon 0.0005))
  (declare (type real f1 f2 epsilon))
  (let ((delta (* (abs f1) epsilon)))
    (<= (- f1 delta)
        f2
        (+ f1 delta))))

But if you're measuring temperatures in centigrade, then you don't want things to get particularly fussy close to zero (unless you care a lot about ice forming, in which case maybe you do).

However I'm not a floating point expert: I just know it's all a nightmare.

As an unrelated note: your loop syntax isn't legal. with needs to be before for or other iteration constructs.

Admin1
Related