;; Maxima functions for finding the maximum or minimum
;; Copyright (C) 2005, 2007 Barton Willis

;; Barton Willis
;; Department of Mathematics, 
;; University of Nebraska at Kearney
;; Kearney NE 68847
;; willisb@unk.edu

;; This source code is licensed under the terms of the Lisp Lesser 
;; GNU Public License (LLGPL). The LLGPL consists of a preamble, published
;; by Franz Inc. (http://opensource.franz.com/preamble.html), and the GNU 
;; Library General Public License (LGPL), version 2, or (at your option)
;; any later version.  When the preamble conflicts with the LGPL, 
;; the preamble takes precedence. 

;; This library is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;; Library General Public License for details.

;; You should have received a copy of the GNU Library General Public
;; License along with this library; if not, write to the
;; Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
;; Boston, MA  02110-1301, USA.

(in-package :maxima)
(macsyma-module maxmin)

(eval-when (:compile-toplevel :load-toplevel :execute)
  ($put '$trylevel 1 '$maxmin)  ;; Default: only use basic simplification rules
  ($put '$maxmin 1 '$version))  ;; Let's have version numbers 1,2,3,...

;; Return true if there is pi in the CL list p and qi in the CL lisp q such that
;; x is between pi and qi.  This means that either pi <= x <= qi or
;; qi <= x <= pi. For example, 2x is between x and 3x.

;; Strangely, sign((a-b)*(b-a)) --> pnz but sign(expand((a-b)*(b-a))) --> nz.
;; This is the reason for the $expand.

;; The betweenp simplification is done last; this has some interesting effects:
;; max(x^2,x^4,x^6,x^2+1) (standard simplification) --> max(x^4,x^6,x^2+1) 
;; (betweenp) --> max(x^4,x^6,x^2+1). If the betweenp simplification were done 
;; first, we'd have max(x^2,x^4,x^6,x^2+1) --> max(x^2,x^6,x^2+1) --> max(x^6,x^2+1).

(defun betweenp (x p q)
  (catch 'done
      (dolist (pk p)
	(dolist (qk q)
	  (if (memq (csign ($expand (mul (sub x pk) (sub qk x)))) '($pos $pz)) (throw 'done t))))
      nil))
	  	       
;; Return true if y is the additive inverse of x. 

(defun add-inversep (x y)
  (eq t (meqp x (neg y))))
	    
;; When get(trylevel,maxmin) is two or greater, max and min try additional 
;; O(n^2) and O(n^3) methods.
 
;; Undone:  max(1-x,1+x) - max(x,-x) --> 1.

(defprop $max simp-max operators)

(defun simp-max (l tmp z)
  (let ((acc nil) (sgn) (num-max nil) (issue-warning) (new-max))
    (setq l (margs (specrepcheck l)))
    (dolist (li l)
      (if (op-equalp li '$max) (setq acc (append acc (mapcar #'(lambda (s) (simplifya s z)) (margs li))))
	(push (simplifya li z) acc)))
	  
    ;; First, delete duplicate members of l.
    
    (setq l (sorted-remove-duplicates (sort acc '$orderlessp)))
    (setq acc nil)
        	   		 
    ;; Second, find the largest real number in l. Since (mnump '$%i) --> false, we don't 
    ;; have to worry that num-max is complex. 
	  	    
    (dolist (li l)
      (if (mnump li) (setq num-max (if (or (null num-max) (mgrp li num-max)) li num-max)) (push li acc)))
    (setq l acc)
    (setq acc (if (null num-max) num-max (list num-max)))
    
    ;; Third, accumulate the maximum in the list acc. For each x in l, do:
    
    ;; (a) if x is > or >= every member of acc, set acc to (x),
    ;; (b) if x is < or <= to some member of acc, do nothing,
    ;; (c) if neither 'a' or 'b', push x into acc,
    ;; (d) if x cannot be compared to some member of acc, set issue-warning to true.
    
    (dolist (x l)
      (catch 'done
	(setq new-max t)
	(dolist (ai acc)
	  (setq sgn ($compare x ai))
	  (setq new-max (and new-max (memq sgn '(&> &>=))))
	  (if (eq sgn '$notcomparable) (setq issue-warning t))
	  (if (memq sgn '(&< &= &<=)) (throw 'done t)))
	(if new-max (setq acc (list x)) (push x acc))))
    
    ;; Fourth, when when trylevel is 2 or higher e and -e are members of acc, replace e by |e|.
    
    (cond ((eq t (mgrp ($get '$trylevel '$maxmin) 1))
	   (setq sgn nil)
	   (dolist (ai acc)
	     (setq tmp (if (lenient-realp ai) (member-if #'(lambda (s) (add-inversep ai s)) sgn) nil))
	     (if tmp (setf (car tmp) (take '(mabs) ai)) (push ai sgn)))
	   (setq acc sgn)))
 
    ;; Fifth, when trylevel is 3 or higher and issue-warning is false, try the
    ;; betweenp simplification.

    (cond ((and (not issue-warning) (eq t (mgrp ($get '$trylevel '$maxmin) 2)))
	   (setq l nil)
	   (setq sgn (cdr acc))
	   (dolist (ai acc)
	     (if (not (betweenp ai sgn sgn)) (push ai l))
	     (setq sgn `(,@(cdr sgn) ,ai)))
	   (setq acc l)))

    ;; Finally, do a few clean ups:
	   
    (setq acc (delete '$minf acc))
    (cond ((null acc) '$minf)
	  ((memq '$inf acc) '$inf)
	  ((null (cdr acc)) (car acc))
	  (t  `(($max simp) ,@(sort acc '$orderlessp))))))

(defun limitneg (x)
  (cond ((eq x '$minf) '$inf)
	((eq x '$inf) '$minf)
	((memq x '($und $ind $infinity)) x)
	(t (neg x))))

(defprop $min simp-min operators)
	  
(defun simp-min (l tmp z)
  (declare (ignore tmp))
  (let ((acc nil))
    (setq l (margs (specrepcheck l)))
    (dolist (li l)
      (if (op-equalp li '$min) (setq acc (append acc (mapcar #'(lambda (s) (simplifya s z)) (margs li))))
	(push (simplifya li z) acc)))
    (setq l acc)
    (setq l (mapcar #'limitneg acc))
    (setq l (simplify `(($max) ,@l)))
    (if (op-equalp l '$max)
	`(($min simp) ,@(mapcar #'limitneg (margs l))) (limitneg l))))

;; Several functions (derivdegree for example) use the maximin function. Here is 
;; a replacement that uses simp-min or simp-max.

(defun maximin (l op) (simplify `((,op) ,@l)))
 
(defmfun $lmax (e)
  (simplify `(($max) ,@(require-list-or-set e "$lmax")))) 

(defmfun $lmin (e)
  (simplify `(($min) ,@(require-list-or-set e "$lmin"))))

;; Return the narrowest comparison operator op (<, <=, =, >, >=) such that
;; a op b evaluates to true. Return 'unknown' when either there is no such 
;; operator or when  Maxima's sign function isn't powerful enough to determine
;; such an operator; when Maxima is able to show that either argument is not 
;; real valued, return 'notcomparable.'

;; The subtraction can be a problem--for example, compare(0.1, 1/10)
;; evaluates to "=". But for double floats, I believe 0.1d0 > 1/10. 
;; If you want to convert double and big floats to exact rational
;; numbers, use $rationalize.

;; I think compare(asin(x), asin(x) + 1) should evaluate to < without
;; being quizzed about the sign of x. Thus the call to lenient-extended-realp.

(defun $compare (a b)
  (cond ((or (not (lenient-extended-realp a)) (not (lenient-extended-realp b)))
	 (if (eq t (meqp a b)) '&= '$notcomparable))
	(t
	 (let ((sgn (csign (specrepcheck (sub a b)))))
	   (cond ((eq sgn '$neg) '&<)
		 ((eq sgn '$nz) '&<=)
		 ((eq sgn '$zero) '&=)
		 ((eq sgn '$pz) '&>=)
		 ((eq sgn '$pos) '&>)
		 ((eq sgn '$pn) '&#)
		 ((eq sgn '$pnz) '$unknown)
		 (t '$unknown))))))

;; When it's fairly likely that the real domain of e is nonempty, return true; 
;; otherwise, return false. Even if z has been declared complex, the real domain
;; of z is nonempty; thus (lenient-extended-realp z) --> true.  When does this
;; function lie?  One example is acos(abs(x) + 2). The real domain of this 
;; expression is empty, yet lenient-extended-realp returns true for this input.

(defun lenient-extended-realp (e)
  (and ($freeof '$infinity '$%i '$und '$ind '$false '$true t nil e) ;; what else?
       (not (mbagp e))
       (not ($featurep e '$nonscalarp))
       (not (mrelationp e))
       (not ($member e $arrays))))

(defun lenient-realp (e)
  (and ($freeof '$inf '$minf e) (lenient-extended-realp e)))

;; Convert all floats and big floats in e to an exact rational representation. 

(defun $rationalize (e)
  (cond ((floatp e) (cl-rat-to-maxima (rationalize e)))
	(($bfloatp e) (cl-rat-to-maxima (* (cadr e)(expt 2 (- (caddr e) (third (car e)))))))
	(($mapatom e) e)
	(t (mfuncall '$fullmap '$rationalize e))))