floating point interval arithmetic bounds

[Nikodemus Siivola]

(defun foo (x)
  (assert (> x 0)))

(defun test ()
  (loop for f = most-positive-double-float then (/ f 2d0)
        while (> f 0d0)
        do (foo f)))

(test) -| assertion fails

PROPAGATE-FROM-REFS gets (DOUBLE-FLOAT (0.0D0)) as the type of F, which in turn leads (> F 0d0) being eliminated.

[Nikodemus Siivola (nikodemus)]

The real issue is BOUND-BINOP not realizing that (0d0) * 0.5 can round to zero. (And similarly for other ops and FP arguments.)

[Nikodemus Siivola (nikodemus)]

Fixed:

commit 26627f70af3ede3610e4092175965e0250456a4d
Author: Nikodemus Siivola <email address hidden>
Date: Tue Jun 7 14:14:21 2011 +0300

    more conservative bounds in FP interval arithmetic

     Make BOUND-BINOP return closed intervals when floating
     point rounding can cause an open bound to close.

     Consider:

      ;; Can return zero.
      (defun fii (x)
        (declare (type (single-float (0.0)) x))
        (/ x 2.0))

      ;; Can return 2.0.
      (defun fii (x y)
        (declare (type (single-float 2.0) x)
                 (type (single-float (0.0)) y))
        (+ x y))