File: maxmin.lisp

package info (click to toggle)
maxima 5.27.0-3
  • links: PTS
  • area: main
  • in suites: wheezy
  • size: 120,648 kB
  • sloc: lisp: 322,503; fortran: 14,666; perl: 14,343; tcl: 11,031; sh: 4,146; makefile: 2,047; ansic: 471; awk: 24; sed: 10
file content (252 lines) | stat: -rw-r--r-- 9,612 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
;; 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 (member (csign ($expand (mul (sub x pk) (sub qk x)))) '($pos $pz) :test #'eq) (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))))

;; Define a simplim%function to handle a limit of $max.

(defprop $max simplim$max simplim%function)

(defun simplim$max (expr var val)
  (cons '($max) (mapcar #'(lambda (e) (limit e var val 'think)) (cdr expr))))

;; 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))
    (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
	(dolist (ai acc)
	  (setq sgn ($compare x ai))
	  (cond ((member sgn '(">" ">=") :test #'equal)
		 (setq acc (delete ai acc :test #'eq)))
		((eq sgn '$notcomparable) (setq issue-warning t))
		((member sgn '("<" "=" "<=") :test #'equal)
		 (throw 'done t))))
             (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))
           (let ((flag nil))
             (setq sgn nil)
             (dolist (ai acc)
               (setq tmp (if (lenient-realp ai)
                             (member-if #'(lambda (s) (add-inversep ai s)) sgn)
                             nil))
               (cond (tmp
                      (setf (car tmp) (take '(mabs) ai))
                      (setq flag t))
                     (t (push ai sgn))))
             (if flag
                 ;; We have replaced -e and e with |e|. Call simp-max again.
                 (return-from simp-max (simplify (cons '($max) 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 (if (not issue-warning) (delete '$minf acc) acc))
    (cond ((null acc) '$minf)
          ((and (not issue-warning) (member '$inf acc :test #'eq)) '$inf)
          ((null (cdr acc)) (car acc))
          (t  `(($max simp) ,@(sort acc '$orderlessp))))))

(defun limitneg (x)
  (cond ((eq x '$minf) '$inf)
	((eq x '$inf) '$minf)
	((member x '($und $ind $infinity) :test #'eq) x)
	(t (neg x))))

;; Define a simplim%function to handle a limit of $min.

(defprop $min simplim$min simplim%function)

(defun simplim$min (expr var val)
  (cons '($min) (mapcar #'(lambda (e) (limit e var val 'think)) (cdr expr))))

(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 flonum floats, I believe 0.1 > 1/10. 
;; If you want to convert flonum 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)
  ;; Simplify expressions with infinities, indeterminates, or infinitesimals
  (when (amongl '($ind $und $inf $minf $infinity $zeroa $zerob) a)
    (setq a ($limit a)))
  (when (amongl '($ind $und $inf $minf $infinity $zeroa $zerob) b)
    (setq b ($limit b)))
  (cond ((or (amongl '($infinity $ind $und) a)
             (amongl '($infinity $ind $und) b))
         ;; Expressions with $infinity, $ind, or $und are not comparable
         '$notcomparable)
        ((eq t (meqp a b)) "=")
        ((or (not (lenient-extended-realp a))
             (not (lenient-extended-realp 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)
  (setq e (ratdisrep e))
  (cond ((floatp e) 
	 (let ((significand) (expon) (sign))
	   (multiple-value-setq (significand expon sign) (integer-decode-float e))
	   (cl-rat-to-maxima (* sign significand (expt 2 expon)))))
	(($bfloatp e) (cl-rat-to-maxima (* (cadr e)(expt 2 (- (caddr e) (third (car e)))))))
	(($mapatom e) e)
	(t (simplify (cons (list (mop e)) (mapcar #'$rationalize (margs e)))))))