(in-package :maxima)

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

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

;; from Graham, On Lisp, macro hackery
(defun mkstr (&rest args)
  (with-output-to-string (s)(dolist (a args) (princ a s))))

(defun symb (&rest args) (values (intern (apply #'mkstr args))))


;; take 2 real intervals and grab their insides. Then
;; do something with them. sample usage...
;;(with-ri2 ri1 ri2 (lo1 hi1)(lo2 hi2) (ri (+ lo1 lo2)(+ hi1 h2)))

(defmacro with-ri2 (struct1 struct2 names1 names2 &body body)
  (let ((gs1 (gensym))
	(gs2 (gensym)))
    `(let ((,gs1 ,struct1)
	   (,gs2 ,struct2))
       (let ,(append 
	      (mapcar #'(lambda (f field)
			  `(,f (,(symb "RI-" field) ,gs1)))
		      names1
		      '(lo hi))
	      (mapcar #'(lambda (f field)
			  `(,f (,(symb "RI-" field) ,gs2)))
		      names2
		      '(lo hi)))
	 ,@body))))

(defmacro with-ri (struct1 names1  &body body)
  (let ((gs1 (gensym)))
    `(let ((,gs1 ,struct1))
       (let  ,(mapcar #'(lambda (f field)
			  `(,f (,(symb "RI-" field) ,gs1)))
		      names1
		      '(lo hi))
	 ,@body))))
(defun |ri-LO|(r)(ri-lo r)) ; for ascii CL if necessary
(defun |ri-HI|(r)(ri-hi 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((r ri) (s ri))  
  (with-ri2 r s (lo1 hi1)(lo2 hi2)
	    (let ((prods (sort (list (* lo1 lo2)(* lo1 hi2)(* hi1 lo2)(* hi1 hi2)) #'<)))
	      ($interval (first prods)(fourth prods)))))
  
  (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)))
    (if (and (< s 0)(<= (ri-lo r) 0 (ri-hi r)))
	;; we could mess with infinities and such but for now,
	;; just complain and signal an error
	(merror "Negative power of interval containing zero ~m" r))
    ($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)

(defun $interval_lo(r)(ri-lo r))
(defun $interval_hi(r)(ri-hi r))

(defun $installintervalprog(op prog)
  (setf (get op '$interval_op) #'(lambda(&rest r)(mapply prog r prog))))

#|
foo(x):=compute(interval_lo(x),interval_hi(x))$
installintervalprog(bar,foo)$
bar(interval(1,2)) 			;
intervalsimp(%) 			;
/* returns compute(1,2) */


|#

  
  






  

;; 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 */ ;       

|#