open intervals and type derivation
authorChristophe Rhodes <csr21@cantab.net>
Sat, 12 May 2012 12:40:24 +0000 (13:40 +0100)
committerChristophe Rhodes <csr21@cantab.net>
Sat, 12 May 2012 12:40:24 +0000 (13:40 +0100)
When dealing with open intervals, preserving the openness of a bound in
a result depends on the operation being strictly monotonic, rather than
merely monotonic.  However, many operations which at first sight seem
strictly monotonic are in fact not, including:
- squaring (expt least-positive-double-float 2) = (expt 0d0 2)
- coercion (coerce (1+ double-float-epsilon) 'single-float)
  = (coerce 1d0 'single-float)

Modify various coercion and interval functions to include conservatism
in these situations.  (Fixes lp#997528)

NEWS
src/compiler/float-tran.lisp
src/compiler/srctran.lisp
tests/compiler.impure.lisp

diff --git a/NEWS b/NEWS
index 2c3b5f6..51b81fa 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -61,6 +61,9 @@ changes relative to sbcl-1.0.56:
     allowing faster release of memory back to the OS. (lp#991293)
   * bug fix: WITH-DEADLINE (:SECONDS NIL :OVERRIDE T) now drops any
     existing deadline for the dynamic scope of its body.
+  * bug fix: compiler-internal interval arithmetic needed to be more
+    conservative about open intervals when operated on by monotonic but not
+    strictly-monotonic functions.  (lp#975528)
   * documentation:
     ** improved docstrings: REPLACE (lp#965592)
 
index b5cd536..1877cef 100644 (file)
                                        (if (< x ,most-negative)
                                            ,most-negative
                                            (coerce x ',type)))
-                                     (numeric-type-low num)))
+                                     (numeric-type-low num)
+                                     nil))
                      (hi (bound-func (lambda (x)
                                        (if (< ,most-positive x )
                                            ,most-positive
                                            (coerce x ',type)))
-                                     (numeric-type-high num))))
+                                     (numeric-type-high num)
+                                     nil)))
                 (specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
 
             (defoptimizer (,fun derive-type) ((num))
                  ;; Process the intersection.
                  (let* ((low (interval-low intersection))
                         (high (interval-high intersection))
-                        (res-lo (or (bound-func fun (if increasingp low high))
+                        (res-lo (or (bound-func fun (if increasingp low high) nil)
                                     default-low))
-                        (res-hi (or (bound-func fun (if increasingp high low))
+                        (res-hi (or (bound-func fun (if increasingp high low) nil)
                                     default-high))
                         (format (case (numeric-type-class arg)
                                   ((integer rational) 'single-float)
                    (int-hi (if hi
                                (ceiling (type-bound-number hi))
                                '*))
-                   (f-lo (or (bound-func #'float lo)
+                   (f-lo (or (bound-func #'float lo nil)
                              '*))
-                   (f-hi (or (bound-func #'float hi)
+                   (f-hi (or (bound-func #'float hi nil)
                              '*)))
               (specifier-type `(or (rational ,int-lo ,int-hi)
                                 (single-float ,f-lo, f-hi)))))
                    (int-hi (if hi
                                (ceiling (type-bound-number hi))
                                '*))
-                   (f-lo (or (bound-func #'float lo)
+                   (f-lo (or (bound-func #'float lo nil)
                              '*))
-                   (f-hi (or (bound-func #'float hi)
+                   (f-hi (or (bound-func #'float hi nil)
                              '*)))
               (specifier-type `(or (rational ,int-lo ,int-hi)
                                 (single-float ,f-lo, f-hi)))))
               ;; exactly the same way as the functions themselves do
               ;; it.
               (if (csubtypep arg domain)
-                  (let ((res-lo (bound-func fun (numeric-type-low arg)))
-                        (res-hi (bound-func fun (numeric-type-high arg))))
+                  (let ((res-lo (bound-func fun (numeric-type-low arg) nil))
+                        (res-hi (bound-func fun (numeric-type-high arg) nil)))
                     (unless increasingp
                       (rotatef res-lo res-hi))
                     (make-numeric-type
index ad7abbc..c2f1aec 100644 (file)
 (defun set-bound (x open-p)
   (if (and x open-p) (list x) x))
 
-;;; Apply the function F to a bound X. If X is an open bound, then
-;;; the result will be open. IF X is NIL, the result is NIL.
-(defun bound-func (f x)
+;;; Apply the function F to a bound X. If X is an open bound and the
+;;; function is declared strictly monotonic, then the result will be
+;;; open. IF X is NIL, the result is NIL.
+(defun bound-func (f x strict)
   (declare (type function f))
   (and x
        (handler-case
              (if (and (floatp y)
                       (float-infinity-p y))
                  nil
-                 (set-bound y (consp x)))))
+                 (set-bound y (and strict (consp x))))))
          ;; Some numerical operations will signal SIMPLE-TYPE-ERROR, e.g.
          ;; in the course of converting a bignum to a float.  Default to
          ;; NIL in that case.
                                   `(and (not (fp-zero-p ,xb))
                                         (not (fp-zero-p ,yb))))))))))))
 
+(defun coercion-loses-precision-p (val type)
+  (typecase val
+    (single-float)
+    (double-float (subtypep type 'single-float))
+    (rational (subtypep type 'float))
+    (t (bug "Unexpected arguments to bounds coercion: ~S ~S" val type))))
+
 (defun coerce-for-bound (val type)
   (if (consp val)
-      (list (coerce-for-bound (car val) type))
+      (let ((xbound (coerce-for-bound (car val) type)))
+        (if (coercion-loses-precision-p (car val) type)
+            xbound
+            (list xbound)))
       (cond
         ((subtypep type 'double-float)
          (if (<= most-negative-double-float val most-positive-double-float)
 (defun coerce-and-truncate-floats (val type)
   (when val
     (if (consp val)
-        (list (coerce-and-truncate-floats (car val) type))
+        (let ((xbound (coerce-for-bound (car val) type)))
+          (if (coercion-loses-precision-p (car val) type)
+              xbound
+              (list xbound)))
         (cond
           ((subtypep type 'double-float)
            (if (<= most-negative-double-float val most-positive-double-float)
 ;;; the negative of an interval
 (defun interval-neg (x)
   (declare (type interval x))
-  (make-interval :low (bound-func #'- (interval-high x))
-                 :high (bound-func #'- (interval-low x))))
+  (make-interval :low (bound-func #'- (interval-high x) t)
+                 :high (bound-func #'- (interval-low x) t)))
 
 ;;; Add two intervals.
 (defun interval-add (x y)
 
 ;;; Apply the function F to the interval X. If X = [a, b], then the
 ;;; result is [f(a), f(b)]. It is up to the user to make sure the
-;;; result makes sense. It will if F is monotonic increasing (or
-;;; non-decreasing).
-(defun interval-func (f x)
+;;; result makes sense. It will if F is monotonic increasing (or, if
+;;; the interval is closed, non-decreasing).
+;;;
+;;; (Actually most uses of INTERVAL-FUNC are coercions to float types,
+;;; which are not monotonic increasing, so default to calling
+;;; BOUND-FUNC with a non-strict argument).
+(defun interval-func (f x &optional increasing)
   (declare (type function f)
            (type interval x))
-  (let ((lo (bound-func f (interval-low x)))
-        (hi (bound-func f (interval-high x))))
+  (let ((lo (bound-func f (interval-low x) increasing))
+        (hi (bound-func f (interval-high x) increasing)))
     (make-interval :low lo :high hi)))
 
 ;;; Return T if X < Y. That is every number in the interval X is
 ;;; Compute the square of an interval.
 (defun interval-sqr (x)
   (declare (type interval x))
-  (interval-func (lambda (x) (* x x))
-                 (interval-abs x)))
+  (interval-func (lambda (x) (* x x)) (interval-abs x)))
 \f
 ;;;; numeric DERIVE-TYPE methods
 
index b5d6512..4361794 100644 (file)
             (assert (= 7 (funcall fun 15 3))))))
     (assert (string= "" trace-output))))
 
+(test-util:with-test (:name :bug-997528)
+  (let ((fun (compile nil '(lambda (x)
+                            (declare (optimize (speed 0) (space 0))
+                             (type (integer -228645653448155482 -228645653447928749) x))
+                            (floor 1.0 (the (integer -228645653448151677 -228645653448150900) x))))))
+    (multiple-value-bind (quo rem)
+        (funcall fun -228645653448151381)
+      (assert (= quo -1))
+      (assert (= rem (float -228645653448151381))))))
+
 ;;; success