;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;

;;; A maxima file , 2015-2016, based on an older file ...
;;  A basis for interval arithmetic

;;  constructed by overloading generic arithmetic. 
;;  Richard Fateman, November, 2005
(in-package :maxima)

;;; mysteries. how many indefinite, undefined, empty, not-an-interval, must we have?
;;; can we have one type of object? what to return for empty + NaI ?
;;; caution -- to the extent that this works, we generally assume
;;; exact arithmetic is exact on the endpoints of intervals when the number
;;; types of the endpoints are themselves exact, i.e. integers, rationals.

;;;  We have paid less attention to doing the optimally right thing with 
;;; single/double and especially big-floats.  Roundings can be done
;;; in a slightly un-tight but machine-independent way by
;;; judicious use of bump up and bump down routines. Elsewhere we
;;; discuss computing elementary functions to extra-high precision
;;; in order to get tighter bounds.

;;; We do not provide proofs of validity of code for (for example) sin/cos/exp
;;; for limits on relative error that would be needed for tight bounds.

;;; There is now an IEEE standard for interval arithmetic, IEEE-1788
;;; and a subset of the standard is in process (11/2015). This file's program
;;; corresponds to a sub-sub set. It ignores, among other issues,
;;; computing with "decorations" and "flavors" of intervals.

;;; in some (all?) flavors,
;;; a decoration DAC, Defined And Continuous, is appropriate for
;;; a result f(y), if y is DAC and the function f is DAC over the range y.
;;; There are also bare (undecorated) intervals, which is what we 
;;; implemented here. The ri struct could easily be extended to include flavors.
;;; The computational cost is modest. It is a bother to write
;;; code to add these features for flavors and the very large number of
;;; function variations the standard requires. Many of them trivial.

;;; We should make sure that all the code works for the distinguished intervals
;;;  <-oo..oo> = whole real line, empty_interval, not_an_interval.
;;; Also, if we have machine floating-point numbers entering into the system
;;; as end-points, we must deal correctly with overflow-to-infinity.
;;; Underflow to zero is probably not hazardous.

;;; There is no conversion of a machine-float infinity to a bigfloat
;;; because there is no bigfloat infinity (exponents are unbounded!)


(eval-when (:load-toplevel :compile-toplevel) 
  (declare (proclaim (special $not_an_interval *op* *reducer* fpprec)))

  (aload (get '$printf 'autoload)))	; so $printf is autoloaded
;;  rint():=load("c:/lisp/maxima-interval.lisp")

;; define ri, real interval, as a structure, and define how to print it

(defstruct (ri (:constructor ri (lo hi))
	    (:print-function (lambda (p stream k)
			       (declare (ignore k))
			       (mfuncall '$printf stream "<~a..~a>"  
					 (pbg(ri-lo p))(pbg(ri-hi p))))))
    lo hi)

(setf (get 'ri 'formatter) #'(lambda(p) (mfuncall '$printf nil "<~a..~a>"  
					 (pbg(ri-lo p))(pbg(ri-hi p)))))

(defun $ri_p(x)(ri-p x))  ;; is x a real interval. ri-p defined via defstruct
(defun $ri_lo(x)(convert-badguy-to-maxima(ri-lo x))) ;; pick out endpoints
(defun $ri_hi(x)(convert-badguy-to-maxima(ri-hi x)))

			     

#+sbcl (defvar dfpi #.sb-ext:double-float-positive-infinity)
#+sbcl (defvar dfni #.sb-ext:double-float-negative-infinity)
#+allegro (defvar dfpi #.*infinity-double*)
#+allegro (defvar dfni #.*negative-infinity-double*)
#+allegro (defvar dfnan #.*nan-double*)


;; these two programs are inadequate because
;; (a) sbcl specific
;; (b) if x is an integer or rational or bigfloat, then float-infinity-p
;; produces an error.
(defun $infright(x)  ;;is x an interval with right endpoint of oo
  (and #+sbcl (sb-ext::float-infinity-p (ri-hi x))
	;; what else here
       (< 0 (ri-hi x))))

(defun $infleft(x)  ;;is x an interval with left endpoint of -oo
  (and  (maxima::float-infinity-p (ri-lo x))
	(<(ri-lo x) 0)))

(defun positiveinfp(x)
  (or (and (floatp x)
	   (maxima::float-infinity-p x)
	    (> x 0.0))
      (eq x '$inf)
      ;; check for 1/0,  infinity[3]
      ))

(defun maxima::float-infinity-p(x)
  #+sbcl(sb-ext::float-infinity-p x)
   #+gcl(float-inf-p x)
    #+allegro (infinityp x) ;; what else?
     #+clisp nil)

(defun negativeinfp(x) 
  (or (and (floatp x) (maxima::float-infinity-p x)
	   (< x 0.0))
      (eq x '$minf)
      ;; check for -1/0,  infinity[3],  -1*inf, 
      (min-inf-p x)
      (min-inf-rat-p x) ))

(defmspec $ri(args)			; take zero or 1 or 2 args
  (if (null (cdr args))(make_empty_interval) ; zero args
    
    (let* ((arg1 (meval(cadr args)))
	   (arg2 (if (cddr args)	; a second arg is there. use it
		   (meval (caddr args))
		   arg1)))
    ;;  (format t "~% ri args= ~s , ~s" arg1 arg2)
      (cond ((not (and (mnump arg1)(mnump arg2))) ; both numbers
	     (merror  (intl:gettext "only numbers allowed in real intervals ~m")
		       (list '(mlist)arg1 arg2) )
	     (make_empty_interval))
					     
	    ((mfuncall 'mgrp arg1 arg2)	; require   a<=b
	     (merror  (intl:gettext "endpoints must be in order ~m") 
		       (list '(mlist)arg1 arg2))
	     (make_empty_interval))
	    (t  (ri arg1 arg2)))
  ;; do we want to convert args to floats, bigfloats, rationals? probably not
      )))


(defun badguyp(x) ;; this returns non-nil for single or double nans and infinities
  (and (floatp x)
  #+allegro (excl::exceptional-floating-point-number-p x)
  #+gcl  (or(float-inf-p x)(float-nan-p x))
  #+sbcl (or(sb-ext::float-infinity-p x)(sb-ext::float-nan-p x))
  ))

(defun sign-and-index-my-inf(x)
  (if (> x 0) (marrayref '$infinity (interval_counter) )
    (mul* -1  (marrayref '$infinity (interval_counter)  ))))  ; this makes it harder to check for negative infinity...

(defun index-my-nan(x)(declare (ignore x)) ; for now, anyway.
       (marrayref '$indefinite(interval_counter) ))

(defun $inf()(sign-and-index-my-inf 1))
(defun $minf()(sign-and-index-my-inf -1))
(defun $ind()(index-my-nan 0))

(defun convert-badguy-to-maxima(x)
  (if (not (badguyp x)) x
    #+allegro (cond((infinityp x) (sign-and-index-my-inf x))
		   ((nanp x)(index-my-nan x)))
    #+gcl(cond((float-inf-p x) (sign-and-index-my-inf x))
	      ((float-nan-p x)(index-my-nan x)))
    #+sbcl(cond((sb-ext::float-infinity-p x) (sign-and-index-my-inf x))
	       ((sb-ext::float-nan-p x)(index-my-nan x)))
    #-(or allegro gcl sbcl) x))  ;uh, just don't know what other implementations might do

#+allegro
(defun pbg(x) ;; print bad guy. This uses the mistake that infinities can be compared equal
  (if
      (badguyp x)
      (case x
	((#.excl::*infinity-double*   #.excl::*infinity-single*)
	 "oo")
	((#.excl::*negative-infinity-double* #.excl::*negative-infinity-single*) 
	 "-oo")
	((#.excl::*nan-double*  #.excl::*nan-single*) 
	 "NaN")
	(otherwise x)) ;; not really a bad guy
    x))

#+gcl 
(defun pbg(x) ;; print bad guy for gcl
  (cond
   ((float-inf-p x)(if (< x 0) "-oo" "oo"))

   ((float-nan-p x) "NaN")
   (t x)))

#+sbcl
(defun pbg(x) ;; print bad guy for sbcl.
  (cond((floatp x) ;sbcl barfs if x is not a float below
       (cond
	((sb-ext::float-infinity-p x)(if (< x 0) "-oo" "oo"))
	((sb-ext::float-nan-p x) "NaN")
	(t x)))
       (t x)))

#-(or gcl allegro sbcl) (defun pbg(x) x) ;; dunno about the other lisps

#+(or gcl allegro)
(defun make-inf-nan(n d)(declare(optimize (speed 3)(safety 0))(double-float n d))
       (the double-float (/ n d)))

#+sbcl
(defun make-inf-nan(n d) (sb-int:with-float-traps-masked(:overflow :divide-by-zero :invalid) (/ n d)))

(compile 'make-inf-nan)  ;; make sure the float divide is open coded

#+(or gcl allegro sbcl) (defparameter $infinity_float (make-inf-nan 1.0d0 0.0d0)) ; #<1.#inf00e+000> in gcl
#+(or gcl allegro sbcl) (defparameter $indefinite_float(make-inf-nan 0.0d0 0.0d0)) ; #<1.#ind00e+000> in gcl
;; note, this is something of a hack. There are many possible NaNs, since
;; the fraction part can be used to store data

#| There are a few distinguish interval-ish items of potential
use.  One is $empty_interval, which is (for example) the result of
computing the intersection of two disjoint intervals.  Another is
$not_an_interval, which is generally the result of invalid constructs
and propagates through all operations.  (NaI is like the floating NaN).
|#

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

(defmacro with-struct ((name . fields) struct &body body)
  (let ((gs (gensym)))
    `(let ((,gs ,struct))
      (let ,(mapcar #'(lambda (f)
			`(,f (,(symb name f) ,gs)))
		    fields)
	,@body))));;; e.g. (with-struct (ri- lo hi) r (f lo hi))
;;; based on...
;;;from Figure 18.3: Destructuring on structures. from On Lisp, P. Graham

;; 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 hi2)))
;(print 0)

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

(defmethod two-arg-+ ((r ri)(s ri))
  ;; adding 2 intervals, just add their parts.
  ;; to be more precise we should round down for lo, round up for hi.
  ;; this program does not include rounding directions
  ;; this program uses common lisp numbers like 1/2, not
  ;; Maxima numbers like ((rat) 1 2).
  
  (with-ri2 r s (lo1 hi1)(lo2 hi2) (ri (add* lo1 lo2)(add* hi1 hi2))))

;; maybe generalize above to (= (max r)(max s)(min r)(min s)) ??

(defmethod two-arg-+ (r (s ri)) ;adding num+interval
  (with-ri s (lo1 hi1) (ri (add* lo1 r)(add* hi1 r))))

(defmethod two-arg-+ ((s ri) r)
  (with-ri s (lo1 hi1) (ri (add* lo1 r)(add* hi1 r))))

;; backstop  
(defmethod two-arg-+ (s r) (add* s r))


(defmethod two-arg-* ((r ri)(s ri))
  ;; multiplying intervals, try all 4, taking min and max.
  ;; To be tighter, we should round down for lo, round up for hi.
  ;; Could be done with fewer ops if we check first, e.g. if intervals are 0<lo<hi.
  ;; Also, see comments for two-arg-+
  
  (with-ri2 r s (lo1 hi1)(lo2 hi2)
	    (let ((prods (sort 
			  (list (mul* lo1 lo2)(mul* lo1 hi2)
				(mul* hi1 lo2)(mul* hi1 hi2))
			  #'(lambda(a b)(not(mfuncall 'mgrp a b))))))
	      (ri (car prods)(fourth prods)))))

(defmethod two-arg-* (r (s ri))
  ;; multiplying num by interval
  (with-ri s (lo1 hi1)
	   (let ((prods (sort (list (mul* r lo1)(mul* r hi1))
			      #'(lambda(a b)(not(mfuncall 'mgrp a b))))))
	      (ri (car prods)(cadr prods)))))

(defmethod two-arg-* ((s ri) r) (two-arg-* r s))

(defmethod two-arg-* (s r) (list '(mtimes) s r))


;; ; rest of code from c:/lisp/generic/interval.lisp  left out.

#| The Simplified 1788 Standard requires the 
following 
Forward Elementary Functions (table 9.1)  on Reals

  basic operations
neg,add, sub, mul, div, recip sqr, sqrt,
fma(x,y,z)= (x*y)+z fused multiply-add

DONE: all but fma, in some form or other

  power functions
  pown, pow, exp exp2 exp10

  DONE: exp
  
  trig/hyperbolic functions

sin, cos tan asin, acos, atan, atan2, 
sinh, cosh, tanh, asinh, acosh, atanh

  DONE: sin, atan

  integer functions

  sign, ceil floor trunc 
  roundTiesToEven(x) = nearest integer to x rounding..
  , roundTiesToAway

  Absmax functions

  abs, min, max  

  From table 10.3

  equals  equal lower and upper bounds
  subset
  less    weakly, e.g. left(a)<=left(b), right(a)<=right(b)
  precedes  is less but might touch
  interior
  strictLess
  strictPrecedes  precedes and doesn't touch
  disjoint

  all relations except disjoint are transitive for nonempty

  whole bunch of other functions that sometimes
  make no sense. e.g. mid(x) is the midpoint but
  has "no value" if x is empty or unbounded.
  inf(x) = left, lower bound, is infinity if x is empty
  ........
 need set operations intersection, union, convex hull.
  Should arithmetic work on sets of intervals, e.g. {s1,s2}* {s3,s4} -> set of 4 items ?
  

|#

;; It would be preferable (tighter) to be more delicate with  rounding up/down.
;; This is crude over-expansion by 0.5 to 1.0 ulp per operation.
;; I hope it is not more than that. Machine-independent 
;; language-implementation, independent efficient access
;; to IEEE-754 rounding modes -- seems unavailable. Note that
;; if we could set rounding to "up", we can generally 
;; compute rounding "down" by suitable negations of values.

(defun partition_expression  ;; this runs in maxima.  No setting of res if nil
    (operator pred init combiner action res e)
  (let ((yes init) (no init))
    (cond ((and (not (atom e))(eq (caar e) operator))
           (map nil
               #'(lambda(r)   ;; maxima value false is lisp nil.
                   (if (funcall pred r)
		       (setf yes (funcall combiner r yes)) 
		     ;r fails predicate
                          (setf no (funcall combiner r no))))
               (cdr e))
           (if res (progn (set res (funcall action yes no)) t)
	     (funcall action yes no)))
          (t nil))))


(defun reduce2(op ll)(if (cdr ll)(reduce op ll) (car ll)))

(defun int_or_mnump (z)(or (typep z 'ri)(mnump z)))

(defun $doint(s)  ;; compute the expression s doing interval simplifications.
  (multiple-value-bind (ans condition)
      (ignore-errors   ;; catch arithmetic errors. easier than handler-case..
  (let ((opfun nil))
    (cond
     ((numberp s)(if (integerp s)s (ri s s))) ;; convert non-integer number n to [n..n]. Should we round up/down?
     ((atom s) s)			; just return atoms
     
     ((eq (caar s) 'rat)(/ (cadr s)(caddr s)))
     
     ((mnump s) s) ; bigfloats are not atoms but ARE still numbers. That's all?
     ;; handle commutative operators + and * by separating
     ;; intervals & numbers from  others.  need simplifya to reduce '((mtimes) interval) etc
     
     ;; here's a question: do we allow sets of intervals, e.g.
     ;; {ri(1,2),ri(3,4)}   * 5  --> {ri(5,10), ri(15,20)}
     ;; etc.  that is, in general a double loop {a,b}*{c,d} -> { a*b, a*c, b*c, b*d}, then
     ;; needing simplification by sorting, combining if overlapping, .
     ;; would need to be extended to, for example  sin({a,b})  etc. to be map(sin,{a,b}).
     ;; for addition, 
     ((memq (caar s)'(mplus mtimes))
      (let* ((*op* (caar s)) ;; set up for commutative operators + and *
	     (*reducer* (if (eq *op* 'mplus) #'two-arg-+ #'two-arg-*)))
	(declare (special *reducer* *op*))
	;(format t "~%partition ~s with op ~s" s *op*) 
	(partition_expression *op*
			      #'int_or_mnump  nil #'cons
			      #'int_combiner
			    
			      nil (cons (car s) (mapcar #'$doint (cdr s))))
      
	))
     ;; for others, do dispatch on operator if
     ;; we have a routine to compute interval version of the function
     ((setf opfun (get (caar s) 'interval_fun))
      (funcall opfun (cons (car s)(mapcar #'$doint (cdr s)))))  
     ;; e.g.  (setf (get '$foo 'interval_fun) #'foofun)
     ;; foofun is called on expression (($foo)  x)
   
     (t (mfuncall '$printf t "doint -- does not (yet) do operator  ~a. continuing.." (caar s))s))))
	
    (if condition (whole-real-line) ans))) ;; if condition non-nil then an error happened.

(defun whole-real-line ()(ri dfni dfpi))

(setf (get '$empty_interval 'interval_fun) #'(lambda(r) r)) ; unchanged
(setf (get '$not_an_interval 'interval_fun) #'(lambda(r) r)) ; unchanged? should be contagious.


(defun int_combiner (y n)
  (declare (special *op* *reducer* $not_an_interval))
  (cond ((null n)(reduce2 *reducer* y))
					; all intervs -> one
	;; check if any of these non-interval and non-number items are bad guys.
	((some #'(lambda(z)(or (empty_intervalp z)(not_an_intervalp z))
			) n)
	 $not_an_interval) ;; what to return for empty + NaI ?
	
	(t(simplifya 
	   (cons (list *op*)		; left over pieces
		 (cons (reduce2 *reducer* y)
		       n)) t))))



;; we should use the maxima version of +  for scalars

;;#-gcl (defmethod bu ((k double-float)) (if (> k 0)(* k #.(+ 1 double-float-epsilon))
;;				   (* k #.(- 1 double-float-epsilon))))

;;#-gcl (defmethod bu ((k single-float)) (if (> k 0)(* k #.(+ 1 single-float-epsilon))
;;					 (* k #.(- 1 single-float-epsilon))))
;; huh? gcl doesn't have class single-float, float, real
	

;; we should use the maxima version of +  for scalars

;;#-gcl (defmethod bu ((k double-float)) (if (> k 0)(* k #.(+ 1 double-float-epsilon))
;;				   (* k #.(- 1 double-float-epsilon))))

;;#-gcl (defmethod bu ((k single-float)) (if (> k 0)(* k #.(+ 1 single-float-epsilon))
;;					 (* k #.(- 1 single-float-epsilon))))
;; huh? gcl doesn't have class single-float, float, real
#+ignore 
(defmethod bu ((k t)) (if (> k 0)(* k (add* 1 double-float-epsilon))
			      (* k (add*  1 (mul* -1 double-float-epsilon)))))

(defun bu (k) 
  (declare (special fpprec))
   (cond ((bigfloatp k)	;;;
	(if (> (cadr k) 0)
	    (mul* k (marrayref '$bfbu fpprec))
	  (mul* k (marrayref  '$bfbd fpprec))))
	 
	 ((eql k 0.0d0) least-positive-double-float)

	 ((floatp k)
	  (if (> k 0)(* k #.(+ 1 double-float-epsilon))
	    (* k #.(- 1 double-float-epsilon))))
	 
	 (t k))) ; exact numbers don't bump up or down

(defun bd (k) 
  (declare (special fpprec))
   (cond ((bigfloatp k)	;;;
	(if (> (cadr k) 0)
	    (mul* k (marrayref '$bfbd fpprec))
	  (mul* k (marrayref  '$bfbu fpprec))))
	 ((eql k 0.0d0) least-negative-double-float)

	 ((floatp k)
	
	  (if (> k 0)(* k #.(- 1 double-float-epsilon))
	    (* k #.(+ 1 double-float-epsilon))))
	 
	 (t k)))

;;	patch from float.lisp
	
(defun intofp (l)
  (cond	;;((not (atom l)) ($bfloat l))
   
   ((floatp l) (floattofp l))
   ((and(numberp l)(not(integerp l)))  ; a rational cl number like 1/2
    (float-ratio (list (numerator l)(denominator l))))
   ((equal 0 l) '(0 0))
   ((eq l '$%pi) (fppi))
   ((eq l '$%e) (fpe))
   ((eq l '$%gamma) (fpgamma))
   (t (list (fpround l) (+ *m fpprec)))))

;(defmethod bu ((k (eql 0))) 0) ;exactly 0 stays zero same as rationals
;(defmethod bu ((k rational)) k)
;(defmethod bu ((k (eql 0.0))) least-positive-normalized-single-float)
;(defmethod bu ((k (eql 0.0d0)))least-positive-normalized-double-float)
;(defmethod bd ((k (eql 0.0))) least-negative-normalized-single-float)
;(defmethod bd ((k (eql 0.0d0))) least-negative-normalized-double-float)

;;#-gcl(defmethod bd ((k double-float)) (if (< k 0)(* k #.(+ 1 double-float-epsilon))
;;				   (* k #.(- 1 double-float-epsilon))))
;;#-gcl(defmethod bd ((k single-float)) (if (< k 0)(* k #.(+ 1 single-float-epsilon))
;;					(* k #.(- 1 single-float-epsilon))))
;#+gcl(defmethod bd ((k t)) (if (< k 0)(* k #.(+ 1 double-float-epsilon))
;				   (* k #.(- 1 double-float-epsilon))))
;(defmethod bd ((k rational)) k)

(defun widen-ri (a b) (ri (bd a)(bu b))) 

(defmspec $widen_ri(args) 
  ;; widen_ri(a)  is same as widen_ri(a,a), which bumps down and up
  ;; we can also widen an existing interval, e.g. widen_ri(ri(...));
  (let* ((arg1 (meval(cadr args)))
	 (arg2 (if (cddr args)		; a second arg is there. use it
		   (meval (caddr args))
		 arg1)))
    (if (ri-p arg1) (mfuncall '$ri (bd (ri-lo arg1))(bu (ri-hi arg1)))
    (mfuncall '$ri (bd arg1)(bu arg2)))))

(defun widen-ri(&rest args)
  (if (not (every #'mnump args))
      (merror "ri: only numbers allowed in real intervals ~a" args))
	
  (cond ((null args)(ri least-negative-normalized-double-float least-positive-normalized-double-float))
	((null (cdr args))		; one argument
	 (widen-ri (car args)(car args)))	; same lo and hi
	(t (widen-ri (car args)(cadr args)))))

(defmethod $left ((z ri)) (convert-badguy-to-maxima(ri-lo z)))
(defmethod $left((a t)) a)		;left(3) is 3.  left(ri(3,4)) is 3. left of -oo is -infinity[#]

(defmethod $right((z ri))(convert-badguy-to-maxima(ri-hi z))
)
(defmethod $right((a t)) a)  

;; if we extract an infinity as an endpoint, we should generate inf[new_index] or -inf[]
;; if the interval is empty or NaI, what do we return?
;; ???

#+ignore
(defun $interval_within (a b)   
  ;; a inside b, a,b are number or interval.
  ;; this works for common lisp numbers, won't work for infinities or Maxima numbers
  (<= ($left b) ($left a) ($right a)($right b)))

(defun $interval_within (a b)   
  (if (typep a 'ri)
  (with-ri a (lo hi)
	   (and ($interval_element_inside lo b)
		($interval_element_inside hi b)))
  ($interval_element_inside a b)
  ))


(defun $interval_element_inside(r a)
  (with-ri a (lo hi)
	   (cond ((numberp r)(<= lo r hi)) ; common lisp number - do comparison.
					; works even if lo=dfni, hi=dfpi
		 ((mnump r)		; looks like ((rat) 1 2) or bigfloat
		  
		  ;; this works .. eql for two float-infinities
		  ;; the other comparisons work for any 2 maxima numbers
		  (and (or (eql lo dfni)(mfuncall '$is `((mleqp) ,lo ,r)))
		       (or (eql hi dfpi)(mfuncall '$is `((mleqp) ,r ,hi)))))
		 (t (merror "interval_inside element: ~m not a real number" r)))))
	      
  
(defun empty_intervalp(q) ;; this means we are not using, e.g. [oo,oo] for empty.
  (and (consp q)(eq (caar q) '$empty_interval)))

(defun make_empty_interval() (marrayref '$empty_interval (interval_counter)))

(defun not_an_intervalp(q) ;; not an interval
  (and (consp q)(eq (caar q) '$not_an_interval)))

(defun make_not_an_interval() (marrayref '$not_an_interval (interval_counter)))

(defun $interval_intersection (a b)   
  ;; a intersect b, a,b are number or interval or empty
  (cond ((empty_intervalp a) a)
	((empty_intervalp b) b)
	(t
	 ;; lowest is higher of   left(a), left(b). etc
	 (let ((nlo (mfuncall '$max ($left b) ($left a)))
	       (nhi (mfuncall '$min ($right b) ($right a))))
    
	   (if (mfuncall 'mgrp nlo nhi)
	       (make_empty_interval)	     ;manufacture an empty interval 
	     (ri nlo nhi))))))

(defun $interval_union (a b)   
  ;; a union b, a,b are number or interval or empty
  (cond ((empty_intervalp a) b)
	((empty_intervalp b) a)
	(t
	 (with-ri2 a b 
		   (lo1 hi1)(lo2 hi2)
		       ;; we can't compare maxima bfloat or rat with float infinities ..
		       ;; maxima tries to convert inf. to other form, which fails.
		   ;; maxima also fails (oops) comparing 2^5000 to 2.0. bug
		   
	
		       
		   (if (every #'(lambda(z)(not(infinityp z))) (list lo1 hi1 lo2 hi2))
		       	   ;; next is easy case if all are rationals
		       ;; or floats that don't overflow

			   (let ((nlo (mfuncall '$min lo1 lo2)) ; lowest of the low
				 (nhi (mfuncall '$max hi1 hi2))) ; highest of the high
			     (ri nlo nhi))
		     ;otherwise
		     
		     (list '($interval_union) a b)))))) ;punt for now.


(defun $interval_intersection (a b)   
  ;; a intersect b, a,b are number or interval or empty
  (cond ((empty_intervalp a) b)
	((empty_intervalp b) a)
	(t
	 (with-ri2 a b 
		   (lo1 hi1)(lo2 hi2)
		       ;; we can't compare maxima bfloat or rat with float infinities ..
		       ;; maxima tries to convert inf. to other form, which fails.
		   ;; maxima also fails (oops) comparing 2^5000 to 2.0. bug
		   (if (every #'(lambda(z)(not(infinityp z))) (list lo1 hi1 lo2 hi2))
		       	   ;; next is easy case if all are rationals
		       ;; or floats that don't overflow

			   (let ((nlo (mfuncall '$max lo1 lo2)) ;max of lows
				 (nhi (mfuncall '$min hi1 hi2))) ;min of highs
			     (if (mfuncall 'mgrp nlo nhi) ; if emptyi..
				 (make_empty_interval) ;manufacture an empty interval 
			       (ri nlo nhi)))
		     
		     (list '($interval_intersection) a b)))))) ;punt, for now.


(defparameter interval_index 0)
(defun interval_counter()(incf interval_index)) ; make up a new number
    

;;   [x_1,x_2] ^n for odd integer n is [x_1^n,x_2^n]
;;   for even integer n..
;;  {[x_1, x_2]}^n = [x_1^n, x_2^n], if x_1 >= 0,
;;  {[x_1, x_2]}^n = [x_2^n, x_1^n], if x_2 < 0,
;;  {[x_1, x_2]}^n = [0, \max \{x_1^n, x_2^n \} ], otherwise.

;;  for negative integer n,  [ min(1/x_1,1/x_2),max(1/x_1, 1/x_2)] ^ (-n)

;;   for non-integer rational n .. ugh,   refuse?
;;   for interval n .. ugh

;; is sqrt([0,4])   [0,2]   or [-2,2] ??  Mathematica says [0,2]
;; is 1/[-1,1] an error or union ( [-oo,-1],[1,oo]);  Mathematica says the latter.
;; 




(defun mexpt-interval-top(s)
 (cond ((and (int_or_mnump (cadr s))(int_or_mnump (caddr s)))
	(mexpt-interval (cadr s)(caddr s)))
       ((and (eq (cadr s) '$%e)(int_or_mnump (caddr s)))
	(intexp (caddr s)))
       (t  s)))

(setf (get 'mexpt 'interval_fun) #'mexpt-interval-top)

(defmethod mexpt-interval((a ri) (n integer))
  (cond ((< n 0)  (mexpt-interval (invert-interval a) (- n)))
	((eql n 1) a)
	(t
	 (with-ri a
		  (lo hi)
		  (cond ((oddp n)  
			 (ri (expt lo n)(expt hi n)))
			((>= lo 0) (ri (expt lo n)(expt hi n)))	;; fix to make general compares
			((< hi 0)  (ri (expt hi n)(expt lo n)))
			(t (ri 0 (max (expt lo n)(expt hi n)))))))))

(defmethod mexpt-interval(a (n (eql 1))) a)
 
(defmethod mexpt-interval((a ri) (n (eql  1/2)))  ;; this is SQUARE ROOT
      (with-ri a
	(lo hi)
	(cond
	 ((>= lo 0) (ri ($left($enclose '$sqrt lo))
			($right($enclose '$sqrt hi))))
	;; lo<hi<0. empty?	    
	 ((< hi 0)  (make_empty_interval))
	 ;; lo < 0 < hi
	 (t (ri 0 ($right($enclose '$sqrt hi)))))))

;; info re maxima compares of bigfloats and lisp numbers.
;; compares can be (mfuncall 'mgrp a b)  (mfuncall 'mlsp a b)
;; lessthan or equal is (not (mfuncall 'mgrp a b))
;; greaterthan or equal is (not (mfuncall 'mlsp a b))

;; inta ^ intb =  exp(intb*log(inta))
#+ignore				; for now
(defmethod mexpt-interval((a ri) (b ri))
  (let*((loga (intlog a)))
    (if (typep loga 'ri)(intexp (two-arg-* loga b))
      (list '(mexpt) a b))))


(defmethod invert-interval ((a ri))
      (with-ri a
	(lo hi)
	
	;;
	;; To exclude 0 we  need either 0<lo<=hi or lo<=hi<0
	;; we assume that lo<=hi so we need only check  0<lo  or hi<0
	;; 0<lo == lo>0 == (mfuncall 'mgrp lo 0)
	;; hi<0 == 0>hi == (mfuncall 'mgrp 0 hi)

	(if (not (or (mfuncall 'mgrp lo 0)
		 (mfuncall 'mgrp 0 hi)))
	    ;; if we want to allow maxima SETS to be the result
	    ;; of division by an interval containing zero,
	    ;; we have to upgrade all the rest of the routines to
	    ;; add, multiply, etc SETS of intervals.  Here's a start, ignored..
	    
	    #+ignore
	  (list '($set) (ri (- $infinity_float) (mfuncall 'mexpt  lo -1))
		(ri (mfuncall 'mexpt hi -1) $infinity_float))
	  
	  (ri   (- $infinity_float) $infinity_float)
	
	(let* ((rlo (mfuncall 'mexpt  lo -1))
	       (rhi (mfuncall 'mexpt hi -1))
	       (nlo (mfuncall '$min rlo rhi))
	       (nhi (mfuncall '$max rlo rhi)))
	  (ri nlo nhi)))))

;(print 4)
  

(defmethod mexpt-interval((a ri) z)
  (if (ratnump z)(mexpt-interval a (/ (cadr z)(caddr z)))  ;; ((rat) 1 2) ->  1/2 common lisp
  (merror "I don't know yet how to compute interval to power ~a " z)))

(defmethod mexpt-interval((a ri) (z double-float))
        (with-ri a
	(lo hi)
	(if (< lo  0 hi) (merror "division by interval containing 0: ~m" a))
	
	(let* ((rlo (expt  lo z)) ;; these 4 expressions assume lisp number,
	       ;;not e.g. bigfloat
	       (rhi (expt hi z))
	       (nlo (min rlo rhi))
	       (nhi (max rlo rhi)))
	  (ri nlo nhi))))


#|See IEEE 1788 standard and Nathalie Revol implementation.
Loading this in to lisp will require  GMP, MPFR  (or equivalent)
and MPFI, Revol's code.  How hard is this? May be easier on
UNIX-like systems, not Windows.

|#


;;;................
;;; this patch should already be in newest wxmaxima
#+ignore
(defun wxxml-atom (x l r &aux tmp-x)
  (append l
          (list (cond ((numberp x) (wxxmlnumformat x))
                      ((and (symbolp x) (get x 'wxxmlword)))
                      ((and (symbolp x) (get x 'reversealias))
                       (wxxml-stripdollar (get x 'reversealias)))
		      ((stringp x)
		       (setq tmp-x (wxxml-fix-string x))
		       (if (and (boundp '$stringdisp) $stringdisp)
			   (setq tmp-x (format nil "\"~a\"" tmp-x)))
		       (concatenate 'string "<st>" tmp-x "</st>"))
		      ((arrayp x)
		       (format nil "<v>Lisp array [~{~a~^,~}]</v>"
			       (array-dimensions x)))
		      ((streamp x)
		       (format nil "<v>Stream [~A]</v>"
			       (stream-element-type x)))
		      ((member (type-of x) '(GRAPH DIGRAPH))
		       (format nil "<v>~a</v>" x))
				                      
		      ((typep x 'structure-object)
;; the patch		       
		       #+ignore   (format nil "<v>Structure [~A]</v>" (type-of x))
		       (format nil "<st>~a</st>" (wxxml-fix-string
						  (format nil "~s" x)))
		       ;; end of patch
		      )
		      
		      ((hash-table-p x)
		       (format nil "<v>HashTable</v>"))
                      (t (wxxml-stripdollar x))))
	  r))



;; here's an idea for sin(interval).
;;  we know that sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y).
;; so sin([lo,hi]) can be expressed as sin([0, hi-lo] + lo)
;; = sin([0,hi-lo])*cos(lo)+ cos([0,hi-lo])* sin(lo)
;; note that hi>=lo so hi-lo is >=0
;; which reduces the problem to computing sin/cos of an interval [0,v], v>=0
;; Is it easier to see what multiples of pi/2 are in [0,v]?
;; Is it useful that we know sin(0)=0, cos(0)=1 exactly?

;; here's old code, not using the idea above.

(defparameter piby2 (/ pi 2))
(defparameter twobypi (/ 2 pi))

(defun intsin (z)
  "real interval sin of a real interval of machine floats."
  ;; result endpoints are same precision as inputs, generally, except if
  ;; we note that extrema -1, 1 are reached, in which case they may be integers
  ;; extending this with bfloats would make good sense
  ;; extending to all other plausible  operations like cos, exp, log, ... too.
  
  ;; here, as for any interval function,  check for not_an_interval empty_interval
  (if (or (empty_intervalp z)(not_an_intervalp z)) z
  (with-ri z	   (low hi)
	   (cond 
	    #+allegro  
	    ((or (badguyp low)(badguyp hi)(> low hi))(ri -1 1))
	    ((eql low hi)(let ((v (sin low)))(widen-ri v v)))
	    ;; return a non-interval?? or maybe should widen interval a little?
	    (t(let (u v minval maxval 
		    ;;  (l (ceiling low piby2)) 
		    (l (mfuncall '$ceiling (mul* low twobypi)))
		    ;; (h (floor hi piby2))
		    (h (mfuncall '$floor (mul* hi twobypi)))
		       )
		(cond ((>= (- h l) 4) (return-from intsin (ri -1 1)))) ; h should be integer
		(setf u ($enclose '%sin low))
		(setf v ($enclose '%sin hi))
		(setf minval (mfuncall '$min ($left u)($left v)))
		(setf maxval (mfuncall '$max ($right u)($right v)))
		;(format t "~%l=~s, h=~s" l h )
		(do 
		    ((k  l (1+ k)))
		    ((> k h)(ri (bd minval)(bu maxval)))
		 ; (format t "~% intsin k=~s h=~s" k h)
		  (case (mod k 4)
		    (1 (setf maxval 1))
		    (3 (setf minval -1))))))))))


(setf (get '%sin 'interval_fun) #'(lambda(r)(intsin(cadr r))))

(defun intlog(r)
  (if (or (empty_intervalp r)(not_an_intervalp r)) r
  (with-ri r
     (low hi)
     (if (and (mfuncall 'mgrp low 0) (mfuncall 'mgrp hi 0))
	 (ri ($left($enclose '$log low)) ($right($enclose '$log hi)))
       (list '(%log) r)  ;; punt. we don't know real interval for log([-1,..]).))))
     ))))
     
(setf (get '%log 'interval_fun) #'(lambda(r)(intlog(cadr r))))


;; atan is monotonic . everywhere

(defun intatan(r)
  (if (or (empty_intervalp r)(not_an_intervalp r)) r
  (with-ri r  ;; works for any finite real bigfloat, float, integer interval.
     (low hi)
	 (ri ($left($enclose '$atan low)) ($right($enclose '$atan hi))))))

(setf (get '%atan 'interval_fun) #'(lambda(r)(intatan(cadr r))))

;; also monotonic asin,  (real, domain [-1,1}, clamped somehow) asinh
;; acos is monotonic with negative slop 

(defun intexp(r)
    (if (or (empty_intervalp r)(not_an_intervalp r)) r
   (with-ri r
     (low hi)
     (ri ($left($enclose '$exp low)) ($right($enclose '$exp hi))))))

;; intexp actually won't be called via property because exp(x) is parsed  into ((mexpt)$%e x)
;; (setf (get '%exp 'interval_fun) #'(lambda(r)(intexp(cadr r)))


;; here are tests: (intsin (ri 0 pi))
;;                 (intsin (ri pi (* 2 pi))
;;                 (intsin (ri (- pi) pi))
;;                 (intsin (ri -1.0 1.0))
;;                 (intsin (ri -1.0d0 1.0d0))

(defun intsinbf (z)
  "real interval sin of a real interval of big floats. Maybe."
  ;; incomplete.  in progress.
  
  (with-ri z
	   (low hi)
	   (cond 
	    ((equal low hi)
	     ($enclose '$sin low))	;so far so good
	    
	    ;; return a non-interval?? or maybe should widen interval a little?
	    (t(let ((u 0)(v 0)(minval 0)(maxval 0)
		    (l (mfuncall '$ceiling (/ low piby2))) ; need bf version, 2 arg
		    (h (floor hi piby2))) ;need bf version also
		(cond ((mfuncall 'mcond(list '(mgreaterp) h (mfuncall #'addk 4 l)))
		       (return-from intsinbf (ri -1 1))))
		(setf u (mfuncall '$sin low)) ;; should work for bigfloats
		(setf v (mfuncall 'sin hi))
		(setf minval (min u v))	;lower value. should round down. need bf
		(setf maxval (max u v))	;upper value. should round up. need bf
		(do 
		    ((k  l (1+ k)))
		    ((> k h)(ri (bd minval)(bu maxval))) ; need bf versions
		  (case (mod k 4)
		    (1 (setf maxval 1))
		    (3 (setf minval -1)))))))))


;;;bigfloat atan fpatan  takes cdr of a bigfloat
;;;e.g. (cons (car x)(fpatan(cdr x))) gets an atan(x), x a bigfloat.


;; a general framework for computing  widen-ri(bigfloat)

(defun widen-bfloat-fun(fun arg oldfpprec)
  ;; eg, %cos  bigfloat  precision,  returns  [a,b]  such that a<=cos(arg)<=b.
  
  (let* ((fpprec (+ 5 oldfpprec))
	 (preciseval (mfuncall fun arg)))
    (setf fpprec oldfpprec)
    (ri (bd preciseval) (bu preciseval))))

#|

maxima code

bfbu[i]:=block([?fpprec:i],bfloat(1+2^(1-i)))$ /*bump up factor x*bfbu[?fpprec] if x>0 */
bfbd[i]:=block([?fpprec:i],bfloat(1-2^(1-i)))$

enclose(fn,point):=
block([h:bfloat(point),
oldfpprec:?fpprec,
?fpprec:?fpprec+5,val],
 val: ?mfuncall (fn, h),
 ?fpprec:oldfpprec,
 if val>0 then ri(val*bfbd[?fpprec],val*bfbu[?fpprec])
   else ri(val*bfbu[?fpprec],val*bfbd[?fpprec]))


   /*e.g.  enclose (sin,1.0b0) returns an interval containing sin(1.0b0),
   where we assume the bigfloat sin () function routinely returns numbers that
   are correct within 32 ulp. */

   |#

(defun  $enclose (fn point)
  (let ((h ($bfloat point))
	(oldfpprec fpprec)
	(fpprec (+ 5 fpprec))
	(val 0))
    (setq val (mfuncall fn h)) ; compute fn(h) in  higher precision
    (setq fpprec oldfpprec)  ; reset precision
     (if
      (is-boole-check (mgrp val 0))
       (ri (mul* val
		 (marrayref '$bfbd fpprec))
	   (mul* val
		 (marrayref '$bfbu fpprec)))

       (ri (mul* val
		 (marrayref '$bfbu fpprec))
	   (mul* val
		 (marrayref '$bfbd fpprec))))))

;; This enclosure program makes sense if you believe that every double-float
;; numerical function produces a numerical result which is either the
;; closest correctly rounded answer or is wrong by less than one unit in last place (ulp).

;; This also means that if a numerical result is 0.0, it is essentially
;; exact, since the measure of any error (relative to zero) would be infinite,
;; and thus the numerical result would not fulfill our assumption.
;; another possibility would be to disbelieve that and have
;; bu(0.0d0) produce least-positive-double-float
;; bd(0.0d0) produce least-negative-double-float
;; which is what we do

(defun  $enclose_float (fn point)  ;; here we assume that point is a double-float
  (let* ((h ($float point))(val (mfuncall fn h)) )		; make sure it is float
    		; compute fn(h) in  higher precision
    (ri (bd val)(bu val))))

;; bfbu[i]:=block([?fpprec:i],bfloat(1+2^(1-i)))$ /*bump up factor x*bfbu[?fpprec] if x>0 */
;; bfbd[i]:=block([?fpprec:i],bfloat(1-2^(1-i)))$
;; ps.  if we really wanted to do arithmetic before simplification,
;; essential if we wanted 0*inf to come out as indefinite, we need
;; a function that receives its argument unsimplified.  That might
;; not be enough, since a*0 is simplified to 0, and it is too late
;; if you later set a to inf.   Anyway... 
(putprop '$noneval (lambda(r) (format t "~% arg to noneval is ~s" r)
			  (cadr r)) 'mfexpr*)

;; one might also wish to simplify oo[1]-oo[2] to indefinite[]
;; and other tricks.

;;  errorless max and min?
;; input could be common lisp numbers, maxima numbers.. bfloat and rat and decimal bfloat
;; input could be, I suppose, intervals, symbols and expressions.
;; problem we are trying to fix:  max(2.0, 2^5000) error: tries for float(2^5000),.

;; not yet written.  3.29/2016
#| ideas:
Assume there are no symbols, no complex numbers.
Or give error message in that case.

Consider that we can find min and max of any CL numbers, no sweat.  {A}
Consider that we can find min and max of bfloats, no sweat. {B}
Consider that we can find min and max of rationals, integers. no sweat. {C}
If we have numbers from 2 of these categories remaining, then
1. make a new list S of  ALL of them, converted to CL rational numbers
   if legit real intervals, choose left or right as appropriate
2. compare elements in S as CL numbers.
3. compute the index of the winner in list S but return number in its
type/form.

should work for:  float inf, float minf, signed zero? mnumps, 
possible OK symbols: $inf, $minf, (- $inf),  arggh.
fails for other symbolic stuff, not-an-interval, float not-a-number, indet, etc.



|#






(dsksetq $bfbu '$bfbu)  ;; this has to do with bigfloat bump up and down, 
(mdefprop $bfbu
    ((lambda) ((mlist) $i)
     ((mprog) ((mlist) ((msetq) fpprec $i))
      (($bfloat) ((mplus) 1 ((mexpt) 2 ((mplus) 1 ((mminus) $i)))))))
    aexpr) 

(add2lnc '$bfbu $props) 
(defvar aaaaa (gensym)) 
(remcompary '$bfbu) 
(mputprop '$bfbu aaaaa 'hashar) 
(*array aaaaa 't 7) 
(fillarray aaaaa '(4 0 1 nil nil nil nil)) 
(add2lnc '$bfbu $arrays) 

(dsksetq $bfbd '$bfbd) 
(mdefprop $bfbd
    ((lambda) ((mlist) $i)
     ((mprog) ((mlist) ((msetq) fpprec $i))
      (($bfloat)
       ((mplus) 1 ((mminus) ((mexpt) 2 ((mplus) 1 ((mminus) $i))))))))
    aexpr) 

(add2lnc '$bfbd $props) 
(setq aaaaa (gensym)) 
(remcompary '$bfbd) 
(mputprop '$bfbd aaaaa 'hashar) 
(*array aaaaa 't 7) 
(fillarray aaaaa '(4 0 1 nil nil nil nil)) 
(add2lnc '$bfbd $arrays) 

#| additional thoughts..
Since this file also provides inf() and ind() ... ($inf )  and ($ind) in lisp,
perhaps $doint should also combine such items when evident.

inf[i]-inf[j] is ind()
inf[i]/inf[j] is ind()  ; maybe positive?
inf[i]+inf[j]  COULD be replaced by a fresh  inf()
ind[i] +  any number is ind()
inf[i] + any number is inf()
inf[i]* any non-zero number x is  sign(x)*inf()
inf[i]*0 is ind().  This is, I fear, deeply problematical
1/inf[i] is zero, except inf[j]*(1/inf[i]) is not zero


|#

;; tests
;; Rump "polynomial" that evaluates the same in IEEE single or double precision,
;; about 1.1726   but
;; is wrong in both

;; rp(x,y):=1335/4*y^6+x^2*(11*x^2*y^2-y^6+(-121)*y^4-2)+11/2*y^8+x/(2*y)$

;; here's where to evaluate it..

;; float(rp(77617, 33096)); evaluates to -0.82739605994682. That's exact, reduced to a float
;; doint((fpprec:20,rp(widen_ri(77617.0b0,77617.0b0),widen_ri(33096.0b0, 33096.0b0))));
;; doint(rp(widen_ri(77617.0,77617.0),widen_ri(33096.0, 33096.0)));


;; for fpprec:10 step 5 thru 50 do print (rp(77617.0b0, 33096.0b0));
;;1.17260394b0
;;1.17260394005318b0
;;1.8014398509481985173b16
;;1.172603940053178631858835b0
;;1.17260394005317863185883490452b0
;;1.1726039400531786318588349045201837b0
;;-8.27396059946821368141165095479816291999b-1
;;-8.27396059946821368141165095479816291999033116b-1
;;-8.2739605994682136814116509547981629199903311578439b-1


;; samples
#|  load ("c:/lisp/maxima-interval.lisp");

doint( ri(1,2)+ 3.0b0);
doint(exp(ri(0,1)));
doint(sin(ri(1.0,1.1)));
doint(sin(ri(1.0b0,1.1b0)));
doint(1/2+ri(1/3,4)+x);
t1:ri(-1,1);
t2:ri(-1,1);
f(doint(t1*t1))+g(doint(t1*t2));
doint (sqrt(t1));

completesq(q,x):=block([c],local(c),q:ratexpand(q),
  c[i]:=coeff(q,x,i),
  c[2]*(x+ c[1]/c[2]/2)^2+c[0]-c[1]^2/4/c[2])$

  e1:x^2+x+1;
  e2:completesq(e1,x);  /*e2 evaluates tighter than e1, though they are "the same" */
  doint(subst(x=t1,e1));
doint(subst(x=t1,e2))			;

/*using tellsimp etc */

matchdeclare([ri1,ri2],ri_p)$ /* ri1 matches only a real interval */
tellsimp(addem(ri1,ri2),doint(ri1+ri2))$
yy:ri(3,4)$ addem(yy,yy) 		;

/*trickier */
simp:off$
tellsimp(min(ri1),left(ri1))$ /* min of one arg is ordinarily that arg */
simp:on$   
  
;; design issues
;; many items in IEEE 1788 just not implemented.

;; e.g. two division ops, 1/[-1,1] -->  [-oo,-1] union [1,oo]  or  [-oo,oo]

;; I think bu and bd are slightly too pessimistic, compared to hardware
;; rounding up/down for machine floats. Maybe we have this fixed for
;; bigfloats though.

;; we have generally tossed efficiency under the bus.

;; known bugs
doint(1/ri(0,1)) ; gives <-oo..oo>. maybe should be <1..oo> ?


|#

;; a very very simple pattern matcher

(defun expand-mat(r x)
  (cond  ((and (consp r)(eq (car r) '--)) t); wildcard match rest of list
	 ((consp r) `(and (consp ,x)
		      ,(expand-mat (car r) (list 'car x))
		      ,(expand-mat (cdr r) (list 'cdr x))))
	 
	 ((null r) `(null ,x))
	 ((eq r '*) t) ;; * means wildcard: match any one item.
	 (t `(eql ',r ,x))))

;; make-mat-fun constructs a function with given name that matches a pattern e.g. 
;; to match a minus infinity, try
;; (make-mat-fun 'min-inf-p '((mtimes simp) -1 (($infinity simp) *)))
;; where the index is unspecified.
;; Equivalent, since we don't care to match the "simp" flags  is
;; (make-mat-fun 'min-inf-p '((mtimes *) -1 (($infinity *) *)))

;; safer is (make-mat-fun 'min-inf-p '((mtimes --) -1 (($infinity --) *)))
;; which matches even '((mtimes) -1 (($infinity) 23)))
;;
(defun make-mat-fun (name template)
  (let ((g (gensym)))
    (compile name `(lambda(,g),(expand-mat template g)))))

(make-mat-fun '$min_inf_p '((mtimes --) -1 (($infinity --) *)))
(make-mat-fun '$min_inf_rat_p '((rat --) -1 0))

;; this will not work if certain ops with float-inf cause interrupts, e.g.
;; division  doint(ri(0,0)/ri(0,0))  gives a float invalid op in SBCL !



  

#+ignore(defun wxxml-stripdollar(r)
  (format nil "<fmn>~a</fmn>" (stripdollar r)))
  
#+ignore
(defmfun stripdollar (x)
  (cond ((not (atom x))
	 (cond ((and (eq (caar x) 'bigfloat) (not (minusp (cadr x)))) (implode (fpformat x)))
	       (t (merror (intl:gettext "STRIPDOLLAR: argument must be an atom; found: ~M") x))))
	((numberp x) x)
	((null x) 'false)
	((eq x t) 'true)
	((member (getcharn x 1) '(#\$ #\%))
	 (intern(maybe-invert-string-case (subseq (string x) 1)))) 
	(t (intern(maybe-invert-string-case(format nil"?~s" x))))))


;; from ellipt.lisp

#+ignore
(defun complex-number-p (u &optional (ntypep 'numberp))
  (or (and (numberp u)(complexp u))  ;; new line
  (labels ((a1 (x) (cadr x))
	   (a2 (x) (caddr x))
	   (a3+ (x) (cdddr x))
	   (N (x) (funcall ntypep x))	     ; N
	   (i (x) (and (eq x '$%i) (N 1)))   ; %i
	   (N+i (x) (and (null (a3+ x))	     ; mplus test is precondition
			 (N (a1 x))
			 (or (i (a2 x))
			     (and (mtimesp (a2 x)) (N*i (a2 x))))))
	   (N*i (x) (and (null (a3+ x))	     ; mtimes test is precondition
			 (N (a1 x))
			 (eq (a2 x) '$%i))))
    (declare (inline a1 a2 a3+ N i N+i N*i))
    (cond ((N u))			     ;2.3
	  ((atom u) (i u))		     ;%i
	  ((mplusp u) (N+i u))		     ;N+%i, N+N*%i
	  ((mtimesp u) (N*i u))		     ;N*%i
	  (t nil)))))

;(setf (gethash '%sin *float-ops*) #'sin)
;; from rpart.lisp
(defun $realpart (xx) (if (numberp xx)(realpart xx)
			  (car (trisplit xx))))

(defun $imagpart (xx)  (if (numberp xx)(imagpart xx)
			 (cdr (trisplit xx))))

;; from trigi.lisp
(defun flonum-eval (op z)
  (let ((flonum-op (gethash op *flonum-op*)))
    (if  (and (numberp z)(complexp z) (floatp (realpart z)))
	(funcall flonum-op z) ;; added
    (when (and flonum-op 

		   (complex-number-p z 'float-or-rational-p))
      (let ((x ($realpart z)) (y ($imagpart z)))
	(when (or $numer (floatp x) (floatp y))
	  (setq x ($float x))
	  (setq y ($float y))
          (bigfloat-alternative
           (complexify
            (funcall flonum-op (if (zerop y) x (complex x y))))
           (big-float-eval op z))))))))

;; from simp.lisp
(defmfun great (x y)
  (cond ((atom x)
	 (cond ((atom y)
		(cond ((numberp x)
		       (cond ((numberp y)
			      (setq y (- x y))
			      (cond ((zerop y) (floatp x))
				    ((and (realp y) (plusp y))) ;complex
				    ((zerop(realpart y)) (plusp (imagpart y)))
				    ((plusp(realpart y)) ) ;else nil     ;complex. end
				    ))))
		      ((constant x)
		       (cond ((constant y) (alphalessp y x)) (t (numberp y))))
		      ((mget x '$scalar)
		       (cond ((mget y '$scalar) (alphalessp y x))
		             (t (maxima-constantp y))))
		      ((mget x '$mainvar)
		       (cond ((mget y '$mainvar) (alphalessp y x)) (t t)))
		      (t (or (maxima-constantp y) (mget y '$scalar)
			     (and (not (mget y '$mainvar)) (alphalessp y x))))))
	       (t (not (ordfna y x)))))
	((atom y) (ordfna x y))
	((eq (caar x) 'rat)
	 (cond ((eq (caar y) 'rat)
		(> (* (caddr y) (cadr x)) (* (caddr x) (cadr y))))))
	((eq (caar y) 'rat))
	((member (caar x) '(mbox mlabox) :test #'eq) (great (cadr x) y))
	((member (caar y) '(mbox mlabox) :test #'eq) (great x (cadr y)))
	((or (member (caar x) '(mtimes mplus mexpt %del) :test #'eq)
	     (member (caar y) '(mtimes mplus mexpt %del) :test #'eq))
	 (ordfn x y))
	((and (eq (caar x) 'bigfloat) (eq (caar y) 'bigfloat)) (mgrp x y))
	((or (eq (caar x) 'mrat) (eq (caar y) 'mrat))
	 (error "GREAT: internal error: unexpected MRAT argument"))
	(t (do ((x1 (margs x) (cdr x1)) (y1 (margs y) (cdr y1))) (())
	     (cond ((null x1)
		    (return (cond (y1 nil)
				  ((not (alike1 (mop x) (mop y)))
				   (great (mop x) (mop y)))
				  ((member 'array (cdar x) :test #'eq) t))))
		   ((null y1) (return t))
		   ((not (alike1 (car x1) (car y1)))
		    (return (great (car x1) (car y1)))))))))

;; from float.lisp

(defmfun $bfloat (x)
  (let (y)
    (cond ((bigfloatp x)) ; possibly changes precision
	  ((numberp x)
	   (cond ((float-nan-p x) ($ind))
		 ((float-inf-p x) (sign-and-index-my-inf x))
		 ;; (t (decintofp x)) ;; if decimal bigfloats
		 (t (bcons(intofp x)))))
	  ((member x '($%e $%pi $%gamma) :test #'eq)
	   (bcons (intofp x)))
	  ((or (atom x) (member 'array (cdar x) :test #'eq))
	   (if (eq x '$%phi)
	       ($bfloat '((mtimes simp)
			  ((rat simp) 1 2)
			  ((mplus simp) 1 ((mexpt simp) 5 ((rat simp) 1 2)))))
	       x))
	  ((eq (caar x) 'mexpt)
	   (if (equal (cadr x) '$%e)
	       (*fpexp ($bfloat (caddr x)))
	       (exptbigfloat ($bfloat (cadr x)) (caddr x))))
	  ((eq (caar x) 'mncexpt)
	   (list '(mncexpt) ($bfloat (cadr x)) (caddr x)))
	  ((eq (caar x) 'rat)
	   (ratbigfloat (cdr x)))
	  ((setq y (safe-get (caar x) 'floatprog))
	   (funcall y (mapcar #'$bfloat (cdr x))))
	  ((or (trigp (caar x)) (arcp (caar x)) (eq (caar x) '$entier))
	   (setq y ($bfloat (cadr x)))
	   (if ($bfloatp y)
	       (cond ((eq (caar x) '$entier) ($entier y))
		     ((arcp (caar x))
		      (setq y ($bfloat (logarc (caar x) y)))
		      (if (free y '$%i)
			  y (let ($ratprint) (fparcsimp ($rectform y)))))
		     ((member (caar x) '(%cot %sec %csc) :test #'eq)
		      (invertbigfloat
		       ($bfloat (list (ncons (safe-get (caar x) 'recip)) y))))
		     (t ($bfloat (exponentialize (caar x) y))))
	       (subst0 (list (ncons (caar x)) y) x)))
	  (t (recur-apply #'$bfloat x)))))