Testing a floating point error agnostic = in Common Lisp

Question

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)

Bail

(1) you say random "fraction" but this is bit confusing (2) regarding execution times "always" will terminate the loop early on the first NIL, whereas for result T all tests must pass – Piane

1) With random fraction I just meant the (random flerag-devi) 2) I am aware of that, I am just curious why the early termination happens with f = 1.0 (and under `*read-default-float-format* = single) and not with other reals. – Bail

Answer 1

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.

Achromatic

Thanks for your comment! I could figure out some of those horrible rs with some integers in a simple loop. Unless I am totally mistaken, the fastest and most secure way for doing my testing tasks is to convert floats as strings to a certain decimal point, something like (string= (format nil "~D" 0.465188032387971243567) (format nil "~D" 0.46567714) :end1 5 :end2 5) to check to the third digit after the point. – Bail

(close-enough-p 0.0 0.000000000000000001) ;NIL obviously since (<= 0.0 1e-20 0.0) is false for any non-normalized float (multiplying by zero as for f1). – Bail

@Student: my close-enough-p was aimed at cases where magnitude matters: a is not 'close enough' to be if a >> b. I am trying to think of a worse way of comparing floats than by printing them, but I can't: if you want to test the floating point implementation then you want to decode the floats and look at the bits. If you want to know that x = y +/- delta y, then write code to check that and live with the fact that there will be odd boundary cases, which is fine: floats are not useful for exact numbers. – Achromatic

Actually I am comparing pixels, for which tiny fractional shifts over a long computation can lead to possibly considerable changes in the graphic. Since pixels don't have a physical substance, I think doing the comparison with a scaled delta is not a good idea in my case (since single pixels shouldn't lose importance depending on other criteria). Thank you anyway for the discussion! – Bail

@Bail then I think your equality test is fine, but you can't expect it always to give the mathematically right answer very near the edge of the interval. – Achromatic