CMUCL commit: src/compiler (float-tran.lisp)

Raymond Toy rtoy at common-lisp.net
Tue Aug 17 22:17:45 CEST 2010


    Date: Tuesday, August 17, 2010 @ 16:17:45
  Author: rtoy
    Path: /project/cmucl/cvsroot/src/compiler

Modified: float-tran.lisp

Add some smarts to the EXPT derive-type optimizer so it can handle the
case of a negative real to an integer power.  Previously, this case
wasn't handled.


-----------------+
 float-tran.lisp |  200 +++++++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 183 insertions(+), 17 deletions(-)


Index: src/compiler/float-tran.lisp
diff -u src/compiler/float-tran.lisp:1.139 src/compiler/float-tran.lisp:1.140
--- src/compiler/float-tran.lisp:1.139	Tue Apr 20 13:57:46 2010
+++ src/compiler/float-tran.lisp	Tue Aug 17 16:17:45 2010
@@ -5,7 +5,7 @@
 ;;; Carnegie Mellon University, and has been placed in the public domain.
 ;;;
 (ext:file-comment
-  "$Header: /project/cmucl/cvsroot/src/compiler/float-tran.lisp,v 1.139 2010-04-20 17:57:46 rtoy Exp $")
+  "$Header: /project/cmucl/cvsroot/src/compiler/float-tran.lisp,v 1.140 2010-08-17 20:17:45 rtoy Exp $")
 ;;;
 ;;; **********************************************************************
 ;;;
@@ -1033,7 +1033,7 @@
 	     (interval-expt-> x y+))))))
 
 ;;; Handle the case when x <= 1
-(defun interval-expt-< (x y)
+(defun interval-expt-< (x y &optional integer-power-p)
   (case (c::interval-range-info x 0d0)
     ('+
      ;; The case of 0 <= x <= 1 is easy
@@ -1067,34 +1067,200 @@
 	;; Split the interval in half
 	(destructuring-bind (y- y+)
 	    (c::interval-split 0 y t)
-	  (list (interval-expt-< x y-)
-		(interval-expt-< x y+))))))
+	  (list (interval-expt-< x y- integer-power-p)
+		(interval-expt-< x y+ integer-power-p))))))
     ('-
-     ;; The case where x <= 0.  Y MUST be an INTEGER for this to
-     ;; work!  The calling function must insure this! For now we'll
-     ;; just return the appropriate unbounded float type.
-     (list (c::make-interval :low nil :high nil)))
+     ;; The case where x <= 0.
+     (cond (integer-power-p
+	    ;; Y is an integer, so we can do something useful.  But we
+	    ;; need to split into x < -1 and -1 <= x <= 0, first
+	    #+(or)
+	    (progn
+	      (format t "integer-power-p = ~A~%" integer-power-p)
+	      (format t "x = ~A~%" x)
+	      (format t "range-info y (~A) = ~A~%" y (interval-range-info y)))
+	    (flet ((handle-positive-power-0 (x y)
+		     ;; -1 <= X <= 0 and Y is positive.  We need to
+		     ;; consider if Y contains an odd integer or not.
+		     ;;
+		     (let* ((y-lo (bound-value (interval-low y)))
+			    (min-odd (if (oddp y-lo)
+					 y-lo
+					 (let ((y-odd (1+ y-lo)))
+					   (if (interval-contains-p y-odd y)
+					       y-odd
+					       nil))))
+			    (min-even (if (evenp y-lo)
+					  y-lo
+					  (let ((y-even (1+ y-lo)))
+					    (if (interval-contains-p y-even y)
+						y-even
+						nil)))))
+		       ;; At least one of min-odd and min-even must be non-NIL!
+		       (assert (or min-odd min-even))
+		       (cond ((and min-odd min-even)
+			      ;; The Y interval contains both even and odd
+			      ;; integers.  Then the lower bound is (least
+			      ;; x)^(least positive odd), because this
+			      ;; creates the most negative value.  The upper
+			      ;; is (most x)^(least positive even), because
+			      ;; this is the most positive number.
+			      ;;
+			      (let ((lo (safe-expt (bound-value (interval-low x))
+						   min-odd))
+				    (hi (safe-expt (bound-value (interval-high x))
+						   min-even)))
+				(list (make-interval :low lo :high hi))))
+			     (min-odd
+			      ;; Y consists of just one odd integer.
+			      (assert (oddp min-odd))
+			      (let ((lo (safe-expt (bound-value (interval-low x))
+						   min-odd))
+				    (hi (safe-expt (bound-value (interval-high x))
+						   min-odd)))
+				(list (make-interval :low lo :high hi))))
+			     (min-even
+			      ;; Y consists of just one even integer.
+			      (assert (evenp min-even))
+			      (let ((lo (safe-expt (bound-value (interval-high x))
+						   min-even))
+				    (hi (safe-expt (bound-value (interval-low x))
+						   min-even)))
+				(list (make-interval :low lo :high hi)))))))
+		   (handle-positive-power-1 (x y)
+		     ;; X <= -1, Y is a positive integer.
+		     (let* ((y-hi (bound-value (interval-high y)))
+			    (max-odd (if (oddp y-hi)
+					 y-hi
+					 (let ((y-odd (1- y-hi)))
+					   (if (interval-contains-p y-odd y)
+					       y-odd
+					       nil))))
+			    (max-even (if (evenp y-hi)
+					  y-hi
+					  (let ((y-even (1- y-hi)))
+					    (if (interval-contains-p y-even y)
+						y-even
+						nil)))))
+		       ;; At least one of max-odd and max-even must be non-NIL!
+		       (assert (or max-odd max-even))
+		       (cond ((and max-odd max-even)
+			      ;; The Y interval contains both even and odd
+			      ;; integers.  Then the lower bound is (least
+			      ;; x)^(most positive odd), because this
+			      ;; creates the most negative value.  The upper
+			      ;; is (least x)^(most positive even), because
+			      ;; this is the most positive number.
+			      ;;
+			      (let ((lo (safe-expt (bound-value (interval-low x))
+						   max-odd))
+				    (hi (safe-expt (bound-value (interval-low x))
+						   max-even)))
+				(list (make-interval :low lo :high hi))))
+			     (max-odd
+			      ;; Y consists of just one odd integer.
+			      (assert (oddp max-odd))
+			      (let ((lo (safe-expt (bound-value (interval-low x))
+						   max-odd))
+				    (hi (safe-expt (bound-value (interval-high x))
+						   max-odd)))
+				(list (make-interval :low lo :high hi))))
+			     (max-even
+			      ;; Y consists of just one even integer.
+			      (assert (evenp max-even))
+			      (let ((lo (safe-expt (bound-value (interval-high x))
+						   max-even))
+				    (hi (safe-expt (bound-value (interval-low x))
+						   max-even)))
+				(list (make-interval :low lo :high hi))))))))
+	      (case (interval-range-info x -1)
+		('+
+		 ;; -1 <= x <= 0
+		 #+(or)
+		 (format t "x range +~%")
+		 (case (interval-range-info y 0)
+		   ('+
+		    (handle-positive-power-0 x y))
+		   ('-
+		    ;; Y is negative.  We should do something better
+		    ;; than this because there's an extra rounding which
+		    ;; we shouldn't do.
+		    #+(or)
+		    (format t "Handle y neg~%")
+		    (let ((unit (make-interval :low 1 :high 1))
+			  (result (handle-positive-power-0 x (interval-neg y))))
+		      #+(or)
+		      (format t "result = ~A~%" result)
+		      (mapcar #'(lambda (r)
+				  (interval-div unit r))
+			      result)))
+		   (t
+		    ;; Split the interval and try again.
+		    (multiple-value-bind (y- y+)
+			(values (make-interval :low (interval-low y)
+					       :high -1)
+				(make-interval :low 1
+					       :high (interval-high y)))
+		      (append (list (make-interval :low 1 :high 1))
+			      (interval-expt-< x y- integer-power-p)
+			      (interval-expt-< x y+ integer-power-p))))))
+		('-
+		 ;; x < -1
+		 (case (c::interval-range-info y)
+		   ('+
+		    ;; Y is positive.  We need to consider if Y contains an
+		    ;; odd integer or not.
+		    ;;
+		    (handle-positive-power-1 x y))
+		   ('-
+		    ;; Y is negative.  Do this in a better way
+		    (let ((unit (make-interval :low 1 :high 1))
+			  (result (handle-positive-power-1 x (interval-neg y))))
+		      (mapcar #'(lambda (r)
+				  (interval-div unit r))
+			      result)))
+		   (t
+		    ;; Split the interval and try again.
+		    #+(or)
+		    (format t "split y ~A~%" y)
+		    (multiple-value-bind (y- y+)
+			(values (make-interval :low (interval-low y) :high -1)
+				(make-interval :low 1 :high (interval-high y)))
+		      (append (list (make-interval :low 1 :high 1))
+			      (interval-expt-< x y- integer-power-p)
+			      (interval-expt-< x y+ integer-power-p))))))
+		(t
+		 #+(or)
+		 (format t "splitting x ~A~%" x)
+		 (destructuring-bind (neg pos)
+		     (interval-split -1 x t t)
+		   (append (interval-expt-< neg y integer-power-p)
+			   (interval-expt-< pos y integer-power-p)))))))
+	    (t
+	     ;; Y is not an integer.  Just give up and return an
+	     ;; unbounded interval.
+	     (list (c::make-interval :low nil :high nil)))))
     (t
      (destructuring-bind (neg pos)
 	 (interval-split 0 x t t)
-       (list (interval-expt-< neg y)
-	     (interval-expt-< pos y))))))
+       (append (interval-expt-< neg y integer-power-p)
+	       (interval-expt-< pos y integer-power-p))))))
 
 ;;; Compute bounds for (expt x y)
 
-(defun interval-expt (x y)
+(defun interval-expt (x y &optional integer-power-p)
   (case (interval-range-info x 1)
     ('+
      ;; X >= 1
 	 (interval-expt-> x y))
     ('-
      ;; X <= 1
-     (interval-expt-< x y))
+     (interval-expt-< x y integer-power-p))
     (t
      (destructuring-bind (left right)
 	 (interval-split 1 x t t)
-       (list (interval-expt left y)
-	     (interval-expt right y))))))
+       (append (interval-expt left y integer-power-p)
+	       (interval-expt right y integer-power-p))))))
 
 (defun fixup-interval-expt (bnd x-int y-int x-type y-type)
   (declare (ignore x-int))
@@ -1228,12 +1394,12 @@
 	   ;; A number to some power is a number.
 	   (specifier-type 'number)))))
   
-(defun merged-interval-expt (x y)
+(defun merged-interval-expt (x y &optional integer-power-p)
   (let* ((x-int (numeric-type->interval x))
 	 (y-int (numeric-type->interval y)))
     (mapcar #'(lambda (type)
 		(fixup-interval-expt type x-int y-int x y))
-	    (flatten-list (interval-expt x-int y-int)))))
+	    (flatten-list (interval-expt x-int y-int integer-power-p)))))
 
 (defun expt-derive-type-aux (x y same-arg)
   (declare (ignore same-arg))
@@ -1243,7 +1409,7 @@
 	 (numeric-contagion x y))
 	((csubtypep y (specifier-type 'integer))
 	 ;; A real raised to an integer power is well-defined
-	 (merged-interval-expt x y))
+	 (merged-interval-expt x y t))
 	(t
 	 ;; A real raised to a non-integral power is complicated....
 	 (cond ((or (csubtypep x (specifier-type '(rational 0)))



More information about the cmucl-commit mailing list