src/maxmin.lisp

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