(in-package :maxima)

;; extended real interval, already in maxima as structure ri, printed as interval(..)

(defun $intervalp(r)(ri-p r))

(defmethod $intervalplus_2((r ri) (s ri))
  ($interval (+ (ri-lo r)(ri-lo s))  ;; fix for roundings
	     (+ (ri-hi r)(ri-hi s))))

;; just a few two-argument versions of plus and times

(defmethod $intervalplus_2((r ri) (s integer))  ;;[a,b]+c -> [a+c,b+c]
   ($interval (+ (ri-lo r)s)  ;; fix for roundings
	      (+ (ri-hi r) s)))

(defmethod $intervalplus_2((s integer)(r ri) )  ;; args reversed.
   ($interval (+ (ri-lo r)s)  ;; fix for roundings
	     (+ (ri-hi r) s)))

(defmethod $intervalplus_2((r t) (s t))  ;; any other rules?
  `((mplus simp) ,r ,s))

(defmethod $intervaltimes_2((r ri) (s integer))  
  
  (if (> s 0)
      ($interval 
       (* (ri-lo r)s)  ;; fix for roundings
       (* (ri-hi r) s))
   ($interval   ;; s is negative 
    (* (ri-hi r) s)
    (* (ri-lo r)s))))

  (defmethod $intervaltimes_2( (s integer) (r ri))  
    ($intervaltimes_2 r s))

(defmethod $intervalmexpt(  (r ri) (s integer))  
  (let ((a (expt (ri-lo r) s)) ;fix for roundings
	(b (expt (ri-hi r) s)))
    ($interval (min a b)(max a b))))

(defun $intervalsimp(r)
  (cond ((atom r) r)
	((member (caar r) '(mplus mtimes)) ;; add here other n-ary stuff

	 (let((intparts nil)
	      (otherparts nil)
	      (argsimp (mapcar '$intervalsimp (cdr r)))
	      (res (if (eq (caar r) 'mtimes) 1 0)) ;identity
	      (intop (get (caar r) '$interval_op)))

	   (if (null intop) ;; if no interval procedure for this op, do this
	       (return-from $intervalsimp
		 (simplifya  (cons (list (caar r))argsimp) nil)))
	      
	   (loop for i in argsimp do
		 ;; separate the pieces, putting all items that
		 ;; might combine with intervals ins the intparts stack
		 
		 (cond ((or (ri-p i)($numberp i)) (push i intparts))
		       (t (push i otherparts))))
	   
	   (if (or (null intparts)(null (cdr intparts)))
	       ;; maybe nothing to do at this level
	       ;;if only 0 or 1 intparts
	       (return-from $intervalsimp
		 (simplifya (cons (list (caar r))argsimp) nil)))
	   
	   (loop for i in intparts do
		 (setf res (funcall intop res i))) ;combine 2 at a time
	   (simplifya (cons (list(caar r)) (cons res otherparts)) nil)
	   ))
	(t (let 
	       ((argsimp (mapcar '$intervalsimp (cdr r)))
		(intop (get (caar r) '$interval_op)))
	     ;;(format t "~% intop=~s, argsimp=~s~%" intop argsimp)
	     (simplifya
	      (if intop (apply intop argsimp)
		;; if no interval procedure for this op, do this
		(cons (list (caar r))argsimp))
		
	      nil)))))


(setf (get 'mplus '$interval_op) '$intervalplus_2)
(setf (get 'mtimes '$interval_op) '$intervaltimes_2)
(setf (get 'mexpt '$interval_op) '$intervalmexpt)

;; need to add programs for every other operation
;; and every kind of argument and endpoints: float, rat, integer
;; and fix the rounding modes. 
;; expt, sin, cos, divide, max, comparisons

#| examples:
interval(1,2) -interval(1,2) 		;
intervalsimp(%)  /*should be interval(-1,1) */			;
x:interval(1,2) 			;
intervalsimp(x^2) 			;
x-x					;
intervalsimp(%)  /* should be 0 */ ;       

|#