1;; Maxima functions for finding the maximum or minimum
2;; Copyright (C) 2005, 2007 Barton Willis
3
4;; Barton Willis
5;; Department of Mathematics,
6;; University of Nebraska at Kearney
7;; Kearney NE 68847
8;; willisb@unk.edu
9
10;; This source code is licensed under the terms of the Lisp Lesser
11;; GNU Public License (LLGPL). The LLGPL consists of a preamble, published
12;; by Franz Inc. (http://opensource.franz.com/preamble.html), and the GNU
13;; Library General Public License (LGPL), version 2, or (at your option)
14;; any later version.  When the preamble conflicts with the LGPL,
15;; the preamble takes precedence.
16
17;; This library is distributed in the hope that it will be useful,
18;; but WITHOUT ANY WARRANTY; without even the implied warranty of
19;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20;; Library General Public License for details.
21
22;; You should have received a copy of the GNU Library General Public
23;; License along with this library; if not, write to the
24;; Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
25;; Boston, MA  02110-1301, USA.
26
27(in-package :maxima)
28(macsyma-module maxmin)
29
30(eval-when (:compile-toplevel :load-toplevel :execute)
31  ($put '$trylevel 1 '$maxmin)  ;; Default: only use basic simplification rules
32  ($put '$maxmin 1 '$version)   ;; Let's have version numbers 1,2,3,...
33  ;; Let's remove built-in symbols from list for user-defined properties.
34  (setq $props (remove '$trylevel $props))
35  (setq $props (remove '$maxmin $props)))
36
37;; Return true if there is pi in the CL list p and qi in the CL lisp q such that
38;; x is between pi and qi.  This means that either pi <= x <= qi or
39;; qi <= x <= pi. For example, 2x is between x and 3x.
40
41;; Strangely, sign((a-b)*(b-a)) --> pnz but sign(expand((a-b)*(b-a))) --> nz.
42;; This is the reason for the $expand.
43
44;; The betweenp simplification is done last; this has some interesting effects:
45;; max(x^2,x^4,x^6,x^2+1) (standard simplification) --> max(x^4,x^6,x^2+1)
46;; (betweenp) --> max(x^4,x^6,x^2+1). If the betweenp simplification were done
47;; 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).
48
49(defun betweenp (x p q)
50  (catch 'done
51      (dolist (pk p)
52	(dolist (qk q)
53	  (if (member (csign ($expand (mul (sub x pk) (sub qk x)))) '($pos $pz) :test #'eq) (throw 'done t))))
54      nil))
55
56;; Return true if y is the additive inverse of x.
57
58(defun add-inversep (x y)
59  (eq t (meqp x (neg y))))
60
61;; Define a simplim%function to handle a limit of $max.
62
63(defprop $max simplim$max simplim%function)
64
65(defun simplim$max (expr var val)
66  (cons '($max) (mapcar #'(lambda (e) (limit e var val 'think)) (cdr expr))))
67
68;; When get(trylevel,maxmin) is two or greater, max and min try additional
69;; O(n^2) and O(n^3) methods.
70
71;; Undone:  max(1-x,1+x) - max(x,-x) --> 1.
72
73(defprop $max simp-max operators)
74
75(defun simp-max (l tmp z)
76  (let ((acc nil) (sgn) (num-max nil) (issue-warning))
77    (setq l (margs (specrepcheck l)))
78    (dolist (li l)
79      (if (op-equalp li '$max) (setq acc (append acc (mapcar #'(lambda (s) (simplifya s z)) (margs li))))
80	(push (simplifya li z) acc)))
81
82    ;; First, delete duplicate members of l.
83
84    (setq l (sorted-remove-duplicates (sort acc '$orderlessp)))
85    (setq acc nil)
86
87    ;; Second, find the largest real number in l. Since (mnump '$%i) --> false, we don't
88    ;; have to worry that num-max is complex.
89
90    (dolist (li l)
91      (if (mnump li) (setq num-max (if (or (null num-max) (mgrp li num-max)) li num-max)) (push li acc)))
92    (setq l acc)
93    (setq acc (if (null num-max) num-max (list num-max)))
94
95    ;; Third, accumulate the maximum in the list acc. For each x in l, do:
96
97    ;; (a) if x is > or >= every member of acc, set acc to (x),
98    ;; (b) if x is < or <= to some member of acc, do nothing,
99    ;; (c) if neither 'a' or 'b', push x into acc,
100    ;; (d) if x cannot be compared to some member of acc, set issue-warning to true.
101
102    (dolist (x l)
103      (catch 'done
104	(dolist (ai acc)
105	  (setq sgn ($compare x ai))
106	  (cond ((member sgn '(">" ">=") :test #'equal)
107		 (setq acc (delete ai acc :test #'eq)))
108		((eq sgn '$notcomparable) (setq issue-warning t))
109		((member sgn '("<" "=" "<=") :test #'equal)
110		 (throw 'done t))))
111             (push x acc)))
112
113    ;; Fourth, when when trylevel is 2 or higher e and -e are members of acc, replace e by |e|.
114
115    (cond ((eq t (mgrp ($get '$trylevel '$maxmin) 1))
116           (let ((flag nil))
117             (setq sgn nil)
118             (dolist (ai acc)
119               (setq tmp (if (lenient-realp ai)
120                             (member-if #'(lambda (s) (add-inversep ai s)) sgn)
121                             nil))
122               (cond (tmp
123                      (setf (car tmp) (take '(mabs) ai))
124                      (setq flag t))
125                     (t (push ai sgn))))
126             (if flag
127                 ;; We have replaced -e and e with |e|. Call simp-max again.
128                 (return-from simp-max (simplify (cons '($max) sgn)))
129                 (setq acc sgn)))))
130
131    ;; Fifth, when trylevel is 3 or higher and issue-warning is false, try the
132    ;; betweenp simplification.
133
134    (cond ((and (not issue-warning) (eq t (mgrp ($get '$trylevel '$maxmin) 2)))
135	   (setq l nil)
136	   (setq sgn (cdr acc))
137	   (dolist (ai acc)
138	     (if (not (betweenp ai sgn sgn)) (push ai l))
139	     (setq sgn `(,@(cdr sgn) ,ai)))
140	   (setq acc l)))
141
142    ;; Finally, do a few clean ups:
143
144    (setq acc (if (not issue-warning) (delete '$minf acc) acc))
145    (cond ((null acc) '$minf)
146          ((and (not issue-warning) (member '$inf acc :test #'eq)) '$inf)
147          ((null (cdr acc)) (car acc))
148          (t  `(($max simp) ,@(sort acc '$orderlessp))))))
149
150(defun limitneg (x)
151  (cond ((eq x '$minf) '$inf)
152	((eq x '$inf) '$minf)
153	((member x '($und $ind $infinity) :test #'eq) x)
154	(t (neg x))))
155
156;; Define a simplim%function to handle a limit of $min.
157
158(defprop $min simplim$min simplim%function)
159
160(defun simplim$min (expr var val)
161  (cons '($min) (mapcar #'(lambda (e) (limit e var val 'think)) (cdr expr))))
162
163(defprop $min simp-min operators)
164
165(defun simp-min (l tmp z)
166  (declare (ignore tmp))
167  (let ((acc nil))
168    (setq l (margs (specrepcheck l)))
169    (dolist (li l)
170      (if (op-equalp li '$min) (setq acc (append acc (mapcar #'(lambda (s) (simplifya s z)) (margs li))))
171	(push (simplifya li z) acc)))
172    (setq l acc)
173    (setq l (mapcar #'limitneg acc))
174    (setq l (simplify `(($max) ,@l)))
175    (if (op-equalp l '$max)
176	`(($min simp) ,@(mapcar #'limitneg (margs l))) (limitneg l))))
177
178;; Several functions (derivdegree for example) use the maximin function. Here is
179;; a replacement that uses simp-min or simp-max.
180
181(defun maximin (l op) (simplify `((,op) ,@l)))
182
183(defmfun $lmax (e)
184  (simplify `(($max) ,@(require-list-or-set e '$lmax))))
185
186(defmfun $lmin (e)
187  (simplify `(($min) ,@(require-list-or-set e '$lmin))))
188
189;; Return the narrowest comparison operator op (<, <=, =, >, >=) such that
190;; a op b evaluates to true. Return 'unknown' when either there is no such
191;; operator or when  Maxima's sign function isn't powerful enough to determine
192;; such an operator; when Maxima is able to show that either argument is not
193;; real valued, return 'notcomparable.'
194
195;; The subtraction can be a problem--for example, compare(0.1, 1/10)
196;; evaluates to "=". But for flonum floats, I believe 0.1 > 1/10.
197;; If you want to convert flonum and big floats to exact rational
198;; numbers, use $rationalize.
199
200;; I think compare(asin(x), asin(x) + 1) should evaluate to < without
201;; being quizzed about the sign of x. Thus the call to lenient-extended-realp.
202
203(defmfun $compare (a b)
204  ;; Simplify expressions with infinities, indeterminates, or infinitesimals
205  (when (amongl '($ind $und $inf $minf $infinity $zeroa $zerob) a)
206    (setq a ($limit a)))
207  (when (amongl '($ind $und $inf $minf $infinity $zeroa $zerob) b)
208    (setq b ($limit b)))
209  (cond ((or (amongl '($infinity $ind $und) a)
210             (amongl '($infinity $ind $und) b))
211         ;; Expressions with $infinity, $ind, or $und are not comparable
212         '$notcomparable)
213        ((eq t (meqp a b)) "=")
214        ((or (not (lenient-extended-realp a))
215             (not (lenient-extended-realp b)))
216         '$notcomparable)
217	(t
218	 (let ((sgn (csign (specrepcheck (sub a b)))))
219	   (cond ((eq sgn '$neg) "<")
220		 ((eq sgn '$nz) "<=")
221		 ((eq sgn '$zero) "=")
222		 ((eq sgn '$pz) ">=")
223		 ((eq sgn '$pos) ">")
224		 ((eq sgn '$pn) "#")
225		 ((eq sgn '$pnz) '$unknown)
226		 (t '$unknown))))))
227
228;; When it's fairly likely that the real domain of e is nonempty, return true;
229;; otherwise, return false. Even if z has been declared complex, the real domain
230;; of z is nonempty; thus (lenient-extended-realp z) --> true.  When does this
231;; function lie?  One example is acos(abs(x) + 2). The real domain of this
232;; expression is empty, yet lenient-extended-realp returns true for this input.
233
234(defun lenient-extended-realp (e)
235  (and ($freeof '$infinity '$%i '$und '$ind '$false '$true t nil e) ;; what else?
236       (not (mbagp e))
237       (not ($featurep e '$nonscalar))
238       (not (mrelationp e))
239       (not ($member e $arrays))))
240
241(defun lenient-realp (e)
242  (and ($freeof '$inf '$minf e) (lenient-extended-realp e)))
243
244;; Convert all floats and big floats in e to an exact rational representation.
245
246(defmfun $rationalize (e)
247  (setq e (ratdisrep e))
248  (cond ((floatp e)
249	 (let ((significand) (expon) (sign))
250	   (multiple-value-setq (significand expon sign) (integer-decode-float e))
251	   (cl-rat-to-maxima (* sign significand (expt 2 expon)))))
252	(($bfloatp e) (cl-rat-to-maxima (* (cadr e)(expt 2 (- (caddr e) (third (car e)))))))
253	(($mapatom e) e)
254	(t (simplify (cons (list (mop e)) (mapcar #'$rationalize (margs e)))))))
255
256