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