;;; -*- Package: C; Log: C.Log -*-
;;;
;;; **********************************************************************

;;; This code was written as part of the CMU Common Lisp project at
;;; Carnegie Mellon University, and has been placed in the public domain.
;;;
(ext:file-comment

  "$Header: srctran.lisp $")

;;;

;;; **********************************************************************
;;;
;;;    This file contains macro-like source transformations which convert
;;; uses of certain functions into the canonical form desired within the
;;; compiler.  ### and other IR1 transforms and stuff.  Some code adapted from
;;; CLC, written by Wholey and Fahlman.
;;;
;;; Written by Rob MacLachlan
;;;

;;; Propagate-float-type extension by Raymond Toy.
;;;

(in-package "C")

(intl:textdomain "cmucl")

;;; Source transform for Not, Null  --  Internal
;;;
;;;    Convert into an IF so that IF optimizations will eliminate redundant
;;; negations.
;;;
(def-source-transform not (x) `(if ,x nil t))
(def-source-transform null (x) `(if ,x nil t))

;;; Source transform for Endp  --  Internal
;;;
;;;    Endp is just NULL with a List assertion.
;;;

(def-source-transform endp (x)
  `(null (the (values &optional list &rest t) ,x)))

;;; We turn Identity into Prog1 so that it is obvious that it just returns the
;;; first value of its argument.  Ditto for Values with one arg.
(def-source-transform identity (x) `(prog1 ,x))
(def-source-transform values (x) `(prog1 ,x))

;;; CONSTANTLY source transform  --  Internal
;;;
;;;    Bind the values and make a closure that returns them.
;;;
(def-source-transform constantly (value &rest values)
  (let ((temps (loop repeat (1+ (length values))
		     collect (gensym)))
	(dum (gensym)))
    `(let ,(loop for temp in temps and
	         value in (list* value values)
	         collect `(,temp ,value))
       #'(lambda (&rest ,dum)
	   (declare (ignore ,dum))
	   (values ,@temps)))))


;;; COMPLEMENT IR1 transform  --  Internal
;;;
;;;    If the function has a known number of arguments, then return a lambda
;;; with the appropriate fixed number of args.  If the destination is a
;;; FUNCALL, then do the &REST APPLY thing, and let MV optimization figure
;;; things out.
;;;

(deftransform complement ((fun) * * :node node :when :both)

  "open code"

  (multiple-value-bind (min max)
		       (function-type-nargs (continuation-type fun))
    (cond
     ((and min (eql min max))
      (let ((dums (loop repeat min collect (gensym))))
	`#'(lambda ,dums (not (funcall fun ,@dums)))))
     ((let* ((cont (node-cont node))
	     (dest (continuation-dest cont)))
	(and (combination-p dest)
	     (eq (combination-fun dest) cont)))
      '#'(lambda (&rest args)
	   (not (apply fun args))))
     (t

      (give-up (intl:gettext "Function doesn't have fixed argument count."))))))


;;;; List hackery:

;;;
;;; Translate CxxR into car/cdr combos.

(defun source-transform-cxr (form)

  (if (or (byte-compiling) (/= (length form) 2))

      (values nil t)
      (let ((name (symbol-name (car form))))
	(do ((i (- (length name) 2) (1- i))
	     (res (cadr form)
		  `(,(ecase (char name i)
		       (#\A 'car)
		       (#\D 'cdr))
		    ,res)))
	    ((zerop i) res)))))

(do ((i 2 (1+ i))
     (b '(1 0) (cons i b)))
    ((= i 5))
  (dotimes (j (ash 1 i))
    (setf (info function source-transform
		(intern (format nil "C~{~:[A~;D~]~}R"
				(mapcar #'(lambda (x) (logbitp x j)) b))))
	  #'source-transform-cxr)))

;;;
;;; Turn First..Fourth and Rest into the obvious synonym, assuming whatever is
;;; right for them is right for us.  Fifth..Tenth turn into Nth, which can be
;;; expanded into a car/cdr later on if policy favors it.
(def-source-transform first (x) `(car ,x))
(def-source-transform rest (x) `(cdr ,x))
(def-source-transform second (x) `(cadr ,x))
(def-source-transform third (x) `(caddr ,x))
(def-source-transform fourth (x) `(cadddr ,x))
(def-source-transform fifth (x) `(nth 4 ,x))
(def-source-transform sixth (x) `(nth 5 ,x))
(def-source-transform seventh (x) `(nth 6 ,x))
(def-source-transform eighth (x) `(nth 7 ,x))
(def-source-transform ninth (x) `(nth 8 ,x))
(def-source-transform tenth (x) `(nth 9 ,x))


;;;
;;; Translate RPLACx to LET and SETF.
(def-source-transform rplaca (x y)
  (once-only ((n-x x))
    `(progn
       (setf (car ,n-x) ,y)
       ,n-x)))
;;;
(def-source-transform rplacd (x y)
  (once-only ((n-x x))
    `(progn
       (setf (cdr ,n-x) ,y)
       ,n-x)))


(def-source-transform nth (n l) `(car (nthcdr ,n ,l)))
  
(defvar *default-nthcdr-open-code-limit* 6)
(defvar *extreme-nthcdr-open-code-limit* 20)

(deftransform nthcdr ((n l) (unsigned-byte t) * :node node)

  "convert NTHCDR to CAxxR"

  (unless (constant-continuation-p n) (give-up))
  (let ((n (continuation-value n)))
    (when (> n
	     (if (policy node (= speed 3) (= space 0))
		 *extreme-nthcdr-open-code-limit*
		 *default-nthcdr-open-code-limit*))
      (give-up))

    (labels ((frob (n)
	       (if (zerop n)
		   'l
		   `(cdr ,(frob (1- n))))))
      (frob n))))


;;;; ARITHMETIC and NUMEROLOGY.

(def-source-transform plusp (x) `(> ,x 0))
(def-source-transform minusp (x) `(< ,x 0))
(def-source-transform zerop (x) `(= ,x 0))

(def-source-transform 1+ (x) `(+ ,x 1))
(def-source-transform 1- (x) `(- ,x 1))

(def-source-transform oddp (x) `(not (zerop (logand ,x 1))))
(def-source-transform evenp (x) `(zerop (logand ,x 1)))

;;; Note that all the integer division functions are available for inline
;;; expansion.

(macrolet ((frob (fun)
	     `(def-source-transform ,fun (x &optional (y nil y-p))
		(declare (ignore y))
		(if y-p
		    (values nil t)
		    `(,',fun ,x 1)))))
  (frob truncate)

  (frob round)
  (frob floor)
  (frob ceiling))

;; Some of these source transforms are not needed when modular
;; arithmetic is available.  When modular arithmetic is available, the
;; various backends need to define them there.
#-modular-arith
(progn

(def-source-transform lognand (x y) `(lognot (logand ,x ,y)))
(def-source-transform lognor (x y) `(lognot (logior ,x ,y)))
(def-source-transform logandc1 (x y) `(logand (lognot ,x) ,y))
(def-source-transform logandc2 (x y) `(logand ,x (lognot ,y)))
(def-source-transform logorc1 (x y) `(logior (lognot ,x) ,y))
(def-source-transform logorc2 (x y) `(logior ,x (lognot ,y)))

)

(def-source-transform logtest (x y) `(not (zerop (logand ,x ,y))))
(def-source-transform byte (size position) `(cons ,size ,position))
(def-source-transform byte-size (spec) `(car ,spec))
(def-source-transform byte-position (spec) `(cdr ,spec))
(def-source-transform ldb-test (bytespec integer)

  `(not (zerop (mask-field ,bytespec ,integer))))

;;; With the ratio and complex accessors, we pick off the "identity" case, and
;;; use a primitive to handle the cell access case.
;;;

(def-source-transform numerator (num)

  (once-only ((n-num `(the (values rational &rest t) ,num)))

    `(if (ratiop ,n-num)

	 (%numerator ,n-num)

	 ,n-num)))
;;;
(def-source-transform denominator (num)

  (once-only ((n-num `(the (values rational &rest t) ,num)))

    `(if (ratiop ,n-num)

	 (%denominator ,n-num)

	 1)))

(deftransform logbitp ((index integer)

		       (integer (or (signed-byte #.vm:word-bits)
				    (unsigned-byte #.vm:word-bits)))

		       (member nil t))
  `(if (>= index #.vm:word-bits)
       (minusp integer)

       (not (zerop (logand integer (ash 1 index))))))


;;;; Interval arithmetic for computing bounds
;;;; (toy@rtp.ericsson.se)
;;;;

;;;; This is a set of routines for operating on intervals.  It implements a
;;;; simple interval arithmetic package.  Although CMUCL has an interval type
;;;; in numeric-type, we choose to use our own for two reasons:

;;;;
;;;;   1.  This package is simpler than numeric-type
;;;;

;;;;   2.  It makes debugging much easier because you can just strip out these
;;;;   routines and test them independently of CMUCL.  (A big win!)

;;;;

;;;; One disadvantage is a probable increase in consing because we have to
;;;; create these new interval structures even though numeric-type has
;;;; everything we want to know.  Reason 2 wins for now.

;;; The basic interval type.  It can handle open and closed intervals.  A
;;; bound is open if it is a list containing a number, just like Lisp says.
;;; NIL means unbounded.
;;;

(defstruct (interval
	     (:constructor %make-interval))

  low high)

(defun make-interval (&key low high)
  (labels ((normalize-bound (val)
	     (cond ((and (floatp val)
			 (float-infinity-p val))
		    ;; Handle infinities
		    nil)
		   ((or (numberp val)
			(eq val nil))
		    ;; Handle any closed bounds
		    val)
		   ((listp val)

		    ;; We have an open bound.  Normalize the numeric bound.
		    ;; If the normalized bound is still a number (not nil),
		    ;; keep the bound open.  Otherwise, the bound is really
		    ;; unbounded, so drop the openness.

		    (let ((new-val (normalize-bound (first val))))
		      (when new-val
			;; Bound exists, so keep it open still
			(list new-val))))
		   (t

		    (error (intl:gettext "Unknown bound type in make-interval!"))))))

    (%make-interval :low (normalize-bound low)
		    :high (normalize-bound high))))

(declaim (inline bound-value set-bound))

;;; Extract the numeric value of a bound.  Return NIL, if X is NIL.

;;;

(defun bound-value (x)
  (if (consp x) (car x) x))

;;; Given a number X, create a form suitable as a bound for an interval.
;;; Make the bound open if OPEN-P is T.  NIL remains NIL.
;;;

(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)
  (and x

       (with-float-traps-masked (:underflow :overflow :inexact :divide-by-zero)

	 ;; With these traps masked, we might get things like infinity or
	 ;; negative infinity returned.  Check for this and return NIL to
	 ;; indicate unbounded.

	 ;;
	 ;; We also ignore any errors that funcall might cause and
	 ;; return NIL instead to indicate infinity.
	 (let ((y (ignore-errors (funcall f (bound-value x)))))

	   (if (and (floatp y)
		    (float-infinity-p y))
	       nil

	       (set-bound y (consp x)))))))

;;; Apply a binary operator OP to two bounds X and Y.  The result is NIL if
;;; either is NIL.  Otherwise bound is computed and the result is open if
;;; either X or Y is open.
;;;

(defmacro bound-binop (op x y)
  `(and ,x ,y

       (with-float-traps-masked (:underflow :overflow :inexact :divide-by-zero)

	 (set-bound (,op (bound-value ,x)
			 (bound-value ,y))
	            (or (consp ,x) (consp ,y))))))

;;; NUMERIC-TYPE->INTERVAL
;;;
;;; Convert a numeric-type object to an interval object.

;;;

(defun numeric-type->interval (x)
  (declare (type numeric-type x))
  (make-interval :low (numeric-type-low x)
		 :high (numeric-type-high x)))

(defun copy-interval-limit (limit)
  (if (numberp limit)
      limit
      (copy-list limit)))

(defun copy-interval (x)
  (declare (type interval x))
  (make-interval :low (copy-interval-limit (interval-low x))
		 :high (copy-interval-limit (interval-high x))))

;;; INTERVAL-SPLIT
;;;

;;; Given a point P contained in the interval X, split X into two interval at
;;; the point P.  If CLOSE-LOWER is T, then the left interval contains P.  If
;;; CLOSE-UPPER is T, the right interval contains P. You can specify both to
;;; be T or NIL.

;;;

(defun interval-split (p x &optional close-lower close-upper)
  (declare (type number p)
	   (type interval x))

  ;; Need to be careful if the lower limit is -0.0 and the split point is 0.
  (let ((low (interval-low x)))
    (cond ((and (zerop p)
		(floatp (bound-value low))
		(member (bound-value low) '(-0f0 -0d0)))
	   (list (make-interval :low (copy-interval-limit low)
				:high (float -0d0 (bound-value low)))
		 (make-interval :low (if close-upper (list p) p)
				:high (copy-interval-limit (interval-high x)))))
	  (t
	   (list (make-interval :low (copy-interval-limit (interval-low x))
				:high (if close-lower p (list p)))
		 (make-interval :low (if close-upper (list p) p)
				:high (copy-interval-limit (interval-high x))))))))

;;; INTERVAL-CLOSURE
;;;

;;; Return the closure of the interval.  That is, convert open bounds to
;;; closed bounds.

;;;

(defun interval-closure (x)
  (declare (type interval x))
  (make-interval :low (bound-value (interval-low x))
		 :high (bound-value (interval-high x))))

;;; INTERVAL-RANGE-INFO
;;;

;;; For an interval X, if X >= POINT, return '+.  If X <= POINT, return

;;; '-. Otherwise return NIL.

;;;
(defun interval-range-info (x &optional (point 0))

  (declare (type interval x))

  (labels ((signed->= (x y)

	     ;; If one of the args is a float, we need to do a float
	     ;; comparison to get the correct value when testing for a
	     ;; signed-zero.  That is, we want (>= -0.0 0) to be false.
	     (if (and (zerop x) (zerop y)
		      (or (floatp x) (floatp y)))
		 (>= (float-sign (float x)) (float-sign (float y)))

		 (>= x y))))
    (let ((lo (interval-low x))
	  (hi (interval-high x)))

      ;; FIXME!  We get confused if X is the interval -0d0 to 0d0.
      ;; Special case that.  What else could we be missing?
      (cond ((and (zerop point)
		  (numberp (bound-value lo))
		  (numberp (bound-value hi))
		  (floatp (bound-value lo))
		  (zerop (bound-value lo))
		  (= (bound-value lo) (bound-value hi)))
	     ;; At this point lo = hi = +/- 0.0.

	     (cond ((or (eql (bound-value lo) (bound-value hi))
			(integerp (bound-value hi)))
		    ;; Both bounds are the same kind of signed 0.  Or
		    ;; the high bound is an exact 0.  The sign of the
		    ;; zero tells us the sign of the interval.

		    (if (= (float-sign (bound-value lo)) -1)
			'-
			'+))
		   (t
		    ;; They have different signs
		    nil)))
	    ((and lo (signed->= (bound-value lo) point))

	     '+)
	    ((and hi (signed->= point (bound-value hi)))
	     '-)
	    (t
	     nil)))))

;;; INTERVAL-BOUNDED-P
;;;

;;; Test to see if the interval X is bounded.  HOW determines the test, and
;;; should be either ABOVE, BELOW, or BOTH.

;;;

(defun interval-bounded-p (x how)
  (declare (type interval x))
  (ecase how
    ('above
     (interval-high x))
    ('below
     (interval-low x))
    ('both
     (and (interval-low x) (interval-high x)))))

;;; Return the sign of the number, taking into account the sign of
;;; signed-zeros.  An integer 0 has a positive sign here.
(declaim (inline number-sign))
(defun number-sign (x)
  (declare (real x))

  (if (floatp x)
      (float-sign x)
      (if (minusp x) -1.0 1.0)))

;;; Signed zero comparison functions.  Use these functions if we need
;;; to distinguish between signed zeroes.  Thus -0.0 < 0.0, which not
;;; true with normal Lisp comparison functions.

(defun signed-zero-= (x y)
  (declare (real x y))
  (and (= x y)

       (= (number-sign x)
	  (number-sign y))))


(macrolet ((frob (name op1 op2)
	     `(defun ,name (x y)
		(declare (real x y))
		;; Comparison (op1) is true, or the numbers are EQUAL
		;; so we need to compare the signs of the numbers
		;; appropriately.
		(or (,op1 x y)
		    (and (= x y)
			 ;; Convert the numbers to long-floats so we
			 ;; don't get problems converting to
			 ;; shorter floats from longer. 
			 (,op2 (number-sign x)
			       (number-sign y)))))))
  (frob signed-zero-< < <)
  (frob signed-zero-> > >)
  (frob signed-zero-<= < <=)
  (frob signed-zero->= > >=))

;;; INTERVAL-CONTAINS-P
;;;

;;; See if the interval X contains the number P, taking into account that the
;;; interval might not be closed.

;;;

(defun interval-contains-p (p x)
  (declare (type number p)
	   (type interval x))

  ;; Does the interval X contain the number P?  This would be a lot easier if
  ;; all intervals were closed!

  (let ((lo (interval-low x))
	(hi (interval-high x)))
    (cond ((and lo hi)
	   ;; The interval is bounded

	   (if (and (signed-zero-<= (bound-value lo) p)
		    (signed-zero-<= p (bound-value hi)))

	       ;; P is definitely in the closure of the interval.
	       ;; We just need to check the end points now.

	       (cond ((signed-zero-= p (bound-value lo))

		      (numberp lo))

		     ((signed-zero-= p (bound-value hi))

		      (numberp hi))
		     (t t))
	       nil))
	  (hi
	   ;; Interval with upper bound

	   (if (signed-zero-< p (bound-value hi))

	       t

	       (and (numberp hi) (signed-zero-= p hi))))

	  (lo
	   ;; Interval with lower bound

	   (if (signed-zero-> p (bound-value lo))

	       t

	       (and (numberp lo) (signed-zero-= p lo))))

	  (t
	   ;; Interval with no bounds
	   t))))

;;; INTERVAL-INTERSECT-P
;;;
;;; Determine if two intervals X and Y intersect.  Return T if so.  If

;;; CLOSED-INTERVALS-P is T, the treat the intervals as if they were closed.
;;; Otherwise the intervals are treated as they are.

;;;

;;; Thus if X = [0, 1) and Y = (1, 2), then they do not intersect because no
;;; element in X is in Y.  However, if CLOSED-INTERVALS-P is T, then they do
;;; intersect because we use the closure of X = [0, 1] and Y = [1, 2] to
;;; determine intersection.

;;;

(defun interval-intersect-p (x y &optional closed-intervals-p)
  (declare (type interval x y))

  (multiple-value-bind (intersect diff)
      (interval-intersection/difference (if closed-intervals-p
					    (interval-closure x)
					    x)
					(if closed-intervals-p
					    (interval-closure y)
					    y))
    (declare (ignore diff))
    intersect))

;;; Are the two intervals adjacent?  That is, is there a number between the
;;; two intervals that is not an element of either interval?  If so, they are
;;; not adjacent.  For example [0, 1) and [1, 2] are adjacent but [0, 1) and
;;; (1, 2] are not because 1 lies between both intervals.
;;;

(defun interval-adjacent-p (x y)
  (declare (type interval x y))
  (flet ((adjacent (lo hi)
	   ;; Check to see if lo and hi are adjacent.  If either is
	   ;; nil, they can't be adjacent.
	   (when (and lo hi (= (bound-value lo) (bound-value hi)))

	     ;; The bounds are equal.  They are adjacent if one of them is
	     ;; closed (a number).  If both are open (consp), then there is a
	     ;; number that lies between them.

	     (or (numberp lo) (numberp hi)))))
    (or (adjacent (interval-low y) (interval-high x))
	(adjacent (interval-low x) (interval-high y)))))

;;; INTERVAL-INTERSECTION/DIFFERENCE
;;;
;;; Compute the intersection and difference between two intervals.
;;; Two values are returned: the intersection and the difference.
;;;

;;; Let the two intervals be X and Y, and let I and D be the two values
;;; returned by this function.  Then I = X intersect Y.  If I is NIL (the
;;; empty set), then D is X union Y, represented as the list of X and Y.  If I
;;; is not the empty set, then D is (X union Y) - I, which is a list of two
;;; intervals.

;;;

;;; For example, let X = [1,5] and Y = [-1,3).  Then I = [1,3) and D = [-1,1)
;;; union [3,5], which is returned as a list of two intervals.

;;;
(defun interval-intersection/difference (x y)
  (declare (type interval x y))

  (let ((x-lo (interval-low x))

	(x-hi (interval-high x))
	(y-lo (interval-low y))
	(y-hi (interval-high y)))
    (labels

	((test-lower-bound (p int)
	   ;; Test if the low bound P is in the interval INT.

	   (if p

	       (if (interval-contains-p (bound-value p)
					(interval-closure int))
		   (let ((lo (interval-low int))
			 (hi (interval-high int)))
		     ;; Check for endpoints
		     (cond ((and lo (= (bound-value p) (bound-value lo)))
			    (not (and (numberp p) (consp lo))))
			   ((and hi (= (bound-value p) (bound-value hi)))
			    (and (numberp p) (numberp hi)))
			   (t t))))

	       (not (interval-bounded-p int 'below))))
	 (test-upper-bound (p int)

	   ;; Test if the upper bound P is in the interval INT.

	   (if p

	       (if (interval-contains-p (bound-value p)
					(interval-closure int))
		   (let ((lo (interval-low int))
			 (hi (interval-high int)))
		     ;; Check for endpoints
		     (cond ((and lo (= (bound-value p) (bound-value lo)))
			    (and (numberp p) (numberp lo)))
			   ((and hi (= (bound-value p) (bound-value hi)))
			    (not (and (numberp p) (consp hi))))
			   (t t))))
	       (not (interval-bounded-p int 'above))))
	 (opposite-bound (p)

	   ;; If P is an open bound, make it closed.  If P is a closed bound,
	   ;; make it open.

	   (if (listp p)
	       (first p)
	       (list p))))

      (let ((x-lo-in-y (test-lower-bound x-lo y))
	    (x-hi-in-y (test-upper-bound x-hi y))
	    (y-lo-in-x (test-lower-bound y-lo x))
	    (y-hi-in-x (test-upper-bound y-hi x)))
	(cond ((or x-lo-in-y x-hi-in-y y-lo-in-x y-hi-in-x)

	       ;; Intervals intersect.  Let's compute the intersection and the
	       ;; difference.

	       (multiple-value-bind (lo left-lo left-hi)
		   (cond (x-lo-in-y
			  (values x-lo y-lo (opposite-bound x-lo)))
			 (y-lo-in-x
			  (values y-lo x-lo (opposite-bound y-lo))))
		 (multiple-value-bind (hi right-lo right-hi)
		     (cond (x-hi-in-y
			    (values x-hi (opposite-bound x-hi) y-hi))
			   (y-hi-in-x
			    (values y-hi (opposite-bound y-hi) x-hi)))
		   (values (make-interval :low lo :high hi)
			   (list (make-interval :low left-lo :high left-hi)
				 (make-interval :low right-lo :high right-hi))))))
	      (t
	       (values nil (list x y))))))))

;;; INTERVAL-MERGE-PAIR
;;;

;;; If intervals X and Y intersect, return a new interval that is the union of
;;; the two.  If they do not intersect, return NIL.

;;;

(defun interval-merge-pair (x y)
  (declare (type interval x y))
  ;; If x and y intersect or are adjacent, create the union.
  ;; Otherwise return nil
  (when (or (interval-intersect-p x y)

	    (interval-adjacent-p x y))

    (flet ((select-bound (x1 x2 min-op max-op)
	     (let ((x1-val (bound-value x1))
		   (x2-val (bound-value x2)))
	       (cond ((and x1 x2)
		      ;; Both bounds are finite.  Select the right one.
		      (cond ((funcall min-op x1-val x2-val)
			     ;; x1 definitely better
			     x1)
			    ((funcall max-op x1-val x2-val)
			     ;; x2 definitely better
			     x2)
			    (t

			     ;; Bounds are equal.  Select either value and
			     ;; make it open only if both were open.

			     (set-bound x1-val (and (consp x1) (consp x2))))))
		     (t

		      ;; At least one bound is not finite.  The non-finite
		      ;; bound always wins.

		      nil)))))

      (let* ((x-lo (copy-interval-limit (interval-low x)))
	     (x-hi (copy-interval-limit (interval-high x)))
	     (y-lo (copy-interval-limit (interval-low y)))
	     (y-hi (copy-interval-limit (interval-high y))))

	(make-interval :low (select-bound x-lo y-lo
					  #'signed-zero-< #'signed-zero->)
		       :high (select-bound x-hi y-hi
					   #'signed-zero-> #'signed-zero-<))))))

;; Wrap a handler-case around BODY so that any errors in the body
;; will return a doubly-infinite interval.
;;
;; This is intended to catch things like (* 0f0 n) where n is an
;; integer that won't fit in a single-float.
;;
;; This is a bit heavy-handed since ANY error gets converted to an
;; unbounded interval.  Perhaps some more fine-grained control would
;; be appropriate?

(defmacro with-unbounded-interval-on-error (() &body body)
  `(handler-case
       (progn ,@body)
     (error ()
       (make-interval :low nil :high nil))))
     

;;; Basic arithmetic operations on intervals.  We probably should do true
;;; interval arithmetic here, but it's complicated because we have float and
;;; integer types and bounds can be open or closed.

;;; INTERVAL-NEG
;;;
;;; 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))))

		       
;;; INTERVAL-ADD
;;;
;;; Add two intervals

;;;

(defun interval-add (x y)
  (declare (type interval x y))

  (with-unbounded-interval-on-error ()
    (make-interval :low (bound-binop + (interval-low x) (interval-low y))
		   :high (bound-binop + (interval-high x) (interval-high y)))))

;;; INTERVAL-SUB
;;;
;;; Subtract two intervals

;;;

(defun interval-sub (x y)
  (declare (type interval x y))

  (with-unbounded-interval-on-error ()
    (make-interval :low (bound-binop - (interval-low x) (interval-high y))
		   :high (bound-binop - (interval-high x) (interval-low y)))))

;;; INTERVAL-MUL
;;;
;;; Multiply two intervals

;;;

(defun interval-mul (x y)
  (declare (type interval x y))
  (flet ((bound-mul (x y)
	   (cond ((or (null x) (null y))
		  ;; Multiply by infinity is infinity
		  nil)
		 ((or (and (numberp x) (zerop x))
		      (and (numberp y) (zerop y)))

		  ;; Multiply by closed zero is special.  The result is always
		  ;; a closed bound.  But don't replace this with zero; we
		  ;; want the multiplication to produce the correct signed
		  ;; zero, if needed.

		  (* (bound-value x) (bound-value y)))

		 ((or (and (floatp x) (float-infinity-p x))
		      (and (floatp y) (float-infinity-p y)))
		  ;; Infinity times anything is infinity
		  nil)
		 (t
		  ;; General multiply.  The result is open if either is open.
		  (bound-binop * x y)))))
    (let ((x-range (interval-range-info x))
	  (y-range (interval-range-info y)))
      (cond ((null x-range)
	     ;; Split x into two and multiply each separately
	     (destructuring-bind (x- x+)
		 (interval-split 0 x t t)
	       (interval-merge-pair (interval-mul x- y)
				    (interval-mul x+ y))))
	    ((null y-range)
	     ;; Split y into two and multiply each separately
	     (destructuring-bind (y- y+)
		 (interval-split 0 y t t)
	       (interval-merge-pair (interval-mul x y-)
				    (interval-mul x y+))))
	    ((eq x-range '-)
	     (interval-neg (interval-mul (interval-neg x) y)))
	    ((eq y-range '-)
	     (interval-neg (interval-mul x (interval-neg y))))
	    ((and (eq x-range '+) (eq y-range '+))
	     ;; If we are here, X and Y are both positive