1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
15 (load-macsyma-macros rzmac
)
17 (declare-top (special *a sum
*i
))
19 (defmvar $opproperties
21 ;; This list was obtained by using an existing version of
22 ;; maxima and printing out the value. It can also be
23 ;; obtained by examining the code below to see what is
24 ;; placed on the OPER variable.
25 '($linear $additive $multiplicative $outative $evenfun $oddfun
26 $commutative $symmetric $antisymmetric $nary $lassociative $rassociative
))
27 "List of the special operator properties recognized by the Maxima simplifier."
28 ;; Don't reset this. (This was originally a defparameter which
29 ;; wouldn't get reset.)
31 ;; We probably don't want the user to modify this except via
33 :properties
((assign 'neverset
)))
35 (loop for
(x y
) on
'(%cot %tan %csc %sin %sec %cos %coth %tanh %csch %sinh %sech %cosh
)
36 by
#'cddr do
(putprop x y
'recip
) (putprop y x
'recip
))
38 ;; polynomial predicates and other such things
40 (defun poly?
(exp var
)
41 (cond ((or (atom exp
) (free exp var
)))
42 ((member (caar exp
) '(mtimes mplus
) :test
#'eq
)
43 (do ((exp (cdr exp
) (cdr exp
)))
45 (and (null (poly?
(car exp
) var
)) (return nil
))))
46 ((and (eq (caar exp
) 'mexpt
)
47 (integerp (caddr exp
))
49 (poly?
(cadr exp
) var
))))
57 (defun smonogen (x var fl
) ; fl indicates whether to return *a *n
58 (cond ((free x var
) (and fl
(setq *n
0 *a x
)) t
)
59 ((atom x
) (and fl
(setq *n
(setq *a
1))) t
)
61 (eq (caar x
) 'mtimes
))
62 (do ((x (cdr x
) (cdr x
))
65 (and fl
(setq *n
(addn n nil
) *a
(muln a nil
))) t
)
67 (if (smonogen (car x
) var fl
)
68 (and fl
(setq a
(cons *a a
) n
(cons *n n
)))
72 (cond ((and (free (caddr x
) var
) (eq (cadr x
) var
))
73 (and fl
(setq *n
(caddr x
) *a
1)) t
)))))
78 (cond ((minusp %m
) (improper-arg-err %m
'$genfact
))
82 a
(if (= %m
1) (return ans
))
83 (setq n
(m- n i
) %m
(1- %m
) ans
(m* ans n
))
86 ;; From Richard Fateman's paper, "Comments on Factorial Programs",
87 ;; http://www.cs.berkeley.edu/~fateman/papers/factorial.pdf
89 ;; k(n,m) = n*(n-m)*(n-2*m)*...
93 ;; This is much faster (3-4 times) than the original factorial
101 (k (- n m
) (* 2 m
))))))
106 ;;; Factorial has mirror symmetry
108 (defprop mfactorial t commutes-with-conjugate
)
110 (defun simpfact (x y z
)
112 (setq y
(simpcheck (cadr x
) z
))
113 (cond ((and (mnump y
)
115 (zerop1 (sub (simplify (list '(%truncate
) y
)) y
)))
116 ;; Negative integer or a real representation of a negative integer.
117 (merror (intl:gettext
"factorial: factorial of negative integer ~:M not defined.") y
))
120 (and (not (integerp y
))
122 (or (and (complex-number-p y
'float-or-rational-p
)
124 (floatp ($realpart y
))
125 (floatp ($imagpart y
))))
126 (and (complex-number-p y
'bigfloat-or-number-p
)
128 ($bfloatp
($realpart y
))
129 ($bfloatp
($imagpart y
))))))
130 (and (not makef
) (ratnump y
) (equal (caddr y
) 2)))
131 ;; Numerically evaluate for real or complex argument in float or
132 ;; bigfloat precision using the Gamma function
133 (simplify (list '(%gamma
) (add 1 y
))))
135 ((and $factorial_expand
138 ;; factorial(n+m) and m integer. Expand.
140 (n (simplify (cons '(mplus) (cddr y
)))))
143 (simplify (list '($pochhammer
) (add n
1) m
))
144 (simplify (list '(mfactorial) n
))))
148 (mul (power -
1 m
) (simplify (list '(mfactorial) n
)))
149 ;; We factor to get the ordering (n-1)*(n-2)*...
151 (simplify (list '($pochhammer
) (mul -
1 n
) m
))))))))
152 ((or (not (fixnump y
)) (not (> y -
1)))
153 (eqtest (list '(mfactorial) y
) x
))
154 ((or (minusp $factlim
) (not (> y $factlim
)))
156 (t (eqtest (list '(mfactorial) y
) x
))))
158 (defun makegamma1 (e)
160 ((eq (caar e
) 'mfactorial
)
161 (list '(%gamma
) (list '(mplus) 1 (makegamma1 (cadr e
)))))
163 ;; Begin code copied from orthopoly/orthopoly-init.lisp
164 ;; Do pochhammer(x,n) ==> gamma(x+n)/gamma(x),
165 ;; but not if x is a negative integer or zero.
167 ((eq (caar e
) '$pochhammer
)
168 (let ((x (makegamma1 (nth 1 e
)))
169 (n (makegamma1 (nth 2 e
))))
170 (if (and (integerp x
) (<= x
0))
172 (div (take '(%gamma
) (add x n
)) (take '(%gamma
) x
)))))
174 ;; (gamma(x/z+1)*z^floor(y))/gamma(x/z-floor(y)+1)
176 ((eq (caar e
) '%genfact
)
177 (let ((x (makegamma1 (nth 1 e
)))
178 (y (makegamma1 (nth 2 e
)))
179 (z (makegamma1 (nth 3 e
))))
180 (setq y
(take '($floor
) y
))
183 (take '(%gamma
) (add (div x z
) 1))
185 (take '(%gamma
) (sub (add (div x z
) 1) y
)))))
186 ;; End code copied from orthopoly/orthopoly-init.lisp
190 ((eq (caar e
) '%double_factorial
)
191 (let ((x (makegamma1 (nth 1 e
))))
197 (sub 1 (simplify (list '(%cos
) (mul '$%pi x
))))))
199 (simplify (list '(%gamma
) (add 1 (div x
2)))))))
201 ((eq (caar e
) '%elliptic_kc
)
202 ;; Complete elliptic integral of the first kind
203 (cond ((alike1 (cadr e
) '((rat simp
) 1 2))
204 ;; K(1/2) = gamma(1/4)/4/sqrt(pi)
205 '((mtimes simp
) ((rat simp
) 1 4)
206 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
207 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 4)) 2)))
208 ((or (alike1 (cadr e
)
209 '((mtimes simp
) ((rat simp
) 1 4)
211 ((mexpt simp
) 3 ((rat simp
) 1 2)))))
213 '((mplus simp
) ((rat simp
) 1 2)
214 ((mtimes simp
) ((rat simp
) 1 4)
215 ((mexpt simp
) 3 ((rat simp
) 1 2)))))
221 ((mexpt simp
) 3 ((rat simp
) 1 2))))
224 '((mtimes simp
) ((rat simp
) 1 4)
225 ((mexpt simp
) 3 ((rat simp
) 1 4))
226 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
227 ((%gamma simp
) ((rat simp
) 1 6))
228 ((%gamma simp
) ((rat simp
) 1 3))))
229 ((or (alike1 (cadr e
)
231 '((mtimes simp
) ((rat simp
) 1 4)
234 ((mexpt simp
) 3 ((rat simp
) 1 2))))))
237 '((mplus simp
) ((rat simp
) 1 2)
238 ((mtimes simp
) ((rat simp
) -
1 4)
239 ((mexpt simp
) 3 ((rat simp
) 1 2)))))
245 ((mexpt simp
) 3 ((rat simp
) 1 2))))
248 '((mtimes simp
) ((rat simp
) 1 4)
249 ((mexpt simp
) 3 ((rat simp
) -
1 4))
250 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
251 ((%gamma simp
) ((rat simp
) 1 6))
252 ((%gamma simp
) ((rat simp
) 1 3))))
254 ;; (3-2*sqrt(2))/(3+2*sqrt(2))
259 ((mexpt simp
) 2 ((rat simp
) 1 2))))
263 ((mexpt simp
) 2 ((rat simp
) 1 2)))) -
1)))
268 ((mexpt simp
) 2 ((rat simp
) 1 2)))))
269 ;; (2*SQRT(2) - 3)/(2*SQRT(2) + 3)
274 ((mexpt simp
) 2 ((rat simp
) 1 2))))
278 ((mexpt simp
) 2 ((rat simp
) 1 2))))
280 '((mtimes simp
) ((rat simp
) 1 8)
281 ((mexpt simp
) 2 ((rat simp
) -
1 2))
282 ((mplus simp
) 1 ((mexpt simp
) 2 ((rat simp
) 1 2)))
283 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
284 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 4)) 2)))
288 ((eq (caar e
) '%elliptic_ec
)
289 ;; Complete elliptic integral of the second kind
290 (cond ((alike1 (cadr e
) '((rat simp
) 1 2))
291 ;; 2*E(1/2) - K(1/2) = 2*%pi^(3/2)*gamma(1/4)^(-2)
293 ((mtimes simp
) ((mexpt simp
) $%pi
((rat simp
) 3 2))
295 ((%gamma simp irreducible
) ((rat simp
) 1 4)) -
2))
296 ((mtimes simp
) ((rat simp
) 1 8)
297 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
298 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 4)) 2))))
299 ((or (alike1 (cadr e
)
300 '((mtimes simp
) ((rat simp
) 1 4)
303 ((mexpt simp
) 3 ((rat simp
) 1 2))))))
305 '((mplus simp
) ((rat simp
) 1 2)
306 ((mtimes simp
) ((rat simp
) -
1 4)
307 ((mexpt simp
) 3 ((rat simp
) 1 2))))))
310 ;; %pi/4/sqrt(3) = K*(E-(sqrt(3)+1)/2/sqrt(3)*K)
312 ((mtimes simp
) ((mexpt simp
) 3 ((rat simp
) -
1 4))
313 ((mexpt simp
) $%pi
((rat simp
) 3 2))
314 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 6)) -
1)
315 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 3)) -
1))
316 ((mtimes simp
) ((rat simp
) 1 8)
317 ((mexpt simp
) 3 ((rat simp
) -
3 4))
318 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
319 ((%gamma simp
) ((rat simp
) 1 6))
320 ((%gamma simp
) ((rat simp
) 1 3)))
321 ((mtimes simp
) ((rat simp
) 1 8)
322 ((mexpt simp
) 3 ((rat simp
) -
1 4))
323 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
324 ((%gamma simp
) ((rat simp
) 1 6))
325 ((%gamma simp
) ((rat simp
) 1 3)))))
326 ((or (alike1 (cadr e
)
327 '((mtimes simp
) ((rat simp
) 1 4)
329 ((mexpt simp
) 3 ((rat simp
) 1 2)))))
331 '((mplus simp
) ((rat simp
) 1 2)
332 ((mtimes simp
) ((rat simp
) 1 4)
333 ((mexpt simp
) 3 ((rat simp
) 1 2))))))
336 ;; %pi*sqrt(3)/4 = K1*(E1-(sqrt(3)-1)/2/sqrt(3)*K1)
338 ((mtimes simp
) 3 ((mexpt simp
) 3 ((rat simp
) -
3 4))
339 ((mexpt simp
) $%pi
((rat simp
) 3 2))
340 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 6)) -
1)
341 ((mexpt simp
) ((%gamma simp
) ((rat simp
) 1 3)) -
1))
342 ((mtimes simp
) ((rat simp
) 3 8)
343 ((mexpt simp
) 3 ((rat simp
) -
3 4))
344 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
345 ((%gamma simp
) ((rat simp
) 1 6))
346 ((%gamma simp
) ((rat simp
) 1 3)))
347 ((mtimes simp
) ((rat simp
) -
1 8)
348 ((mexpt simp
) 3 ((rat simp
) -
1 4))
349 ((mexpt simp
) $%pi
((rat simp
) -
1 2))
350 ((%gamma simp
) ((rat simp
) 1 6))
351 ((%gamma simp
) ((rat simp
) 1 3)))))
354 (t (recur-apply #'makegamma1 e
))))
356 (def-simplifier genfact
(x y z
)
358 (b (take '($floor
) y
))
360 (cond ((and (fixnump a
)
365 (or (<= c a
) (= b
0))
368 (merror (intl:gettext
"genfact: generalized factorial not defined for given arguments."))))
370 ;; Give up, we want to return a result with args that are
371 ;; different from the original. In particular, we want the
372 ;; floor of y if y was real number. Otherwise, we leave
375 (if (and (not (atom b
))
376 (eq (caar b
) '$floor
))
383 ;; These variables should be initialized where they belong.
385 ;; FIXME: maxtaydiff and *trunclist don't appear to used anywhere.
386 (defmvar $maxtaydiff
4)
387 (defvar *trunclist nil
)
389 (defmacro sum-arg
(sum)
392 (defmacro sum-index
(sum)
395 (defmacro sum-lower
(sum)
398 (defmacro sum-upper
(sum)
399 `(cadr (cdddr ,sum
)))
402 (arg-count-check 4 l
)
404 (dosum (car l
) (cadr l
) (meval (caddr l
)) (meval (cadddr l
)) t
:evaluate-summand t
))
408 (arg-count-check 3 l
)
410 ;;(or (= (length l) 3) (wna-err '$lsum))
413 (lis (meval (caddr l
)))
415 (or (symbolp ind
) (merror (intl:gettext
"lsum: second argument must be a variable; found ~M") ind
))
417 (loop for v in
(cdr lis
)
418 with lind
= (cons ind nil
)
421 (setq ans
(add* ans
(mbinding (lind w
) (meval form
)))))
423 (t `((%lsum
) ,form
,ind
,lis
)))))
425 (defun simpsum (x y z
)
426 (let (($ratsimpexpons t
))
427 (setq y
(simplifya (sum-arg x
) z
)))
428 (simpsum1 y
(sum-index x
) (simplifya (sum-lower x
) z
)
429 (simplifya (sum-upper x
) z
)))
431 ; This function was SIMPSUM1 until the sum/product code was revised Nov 2005.
432 ; The revised code punts back to this function since this code knows
433 ; some simplifications not handled by the revised code. -- Robert Dodier
435 (defun simpsum1-save (exp i lo hi
)
436 (cond ((not (symbolp i
)) (merror (intl:gettext
"sum: index must be a symbol; found ~M") i
))
437 ((equal lo hi
) (mbinding ((list i
) (list hi
)) (meval exp
)))
440 (getl '%sum
'($outative $linear
)))
441 (freesum exp lo hi
1))
442 ((null $simpsum
) (list (get '%sum
'msimpind
) exp i lo hi
))
443 ((and (or (eq lo
'$minf
)
444 (alike1 lo
'((mtimes simp
) -
1 $inf
)))
446 (let ((pos-part (simpsum2 exp i
0 '$inf
))
447 (neg-part (simpsum2 (maxima-substitute (m- i
) i exp
) i
1 '$inf
)))
449 ((or (eq neg-part
'$und
)
453 (if (eq neg-part
'$minf
) '$und
'$inf
))
454 ((eq pos-part
'$minf
)
455 (if (eq neg-part
'$inf
) '$und
'$minf
))
456 ((or (eq neg-part
'$inf
) (eq neg-part
'$minf
))
458 (t (m+ neg-part pos-part
)))))
460 (alike1 lo
'((mtimes simp
) -
1 '$inf
)))
461 (simpsum2 (maxima-substitute (m- i
) i exp
) i
(m- hi
) '$inf
))
462 (t (simpsum2 exp i lo hi
))))
464 ;; DOSUM, MEVALSUMARG, DO%SUM -- general principles
466 ;; - evaluate the summand/productand
467 ;; - substitute a gensym for the index variable and make assertions (via assume) about the gensym index
468 ;; - return 0/1 for empty sum/product. sumhack/prodhack are ignored
469 ;; - distribute sum/product over mbags when listarith = true
471 (defun dosum (expr ind low hi sump
&key
(evaluate-summand t
))
472 (setq low
(ratdisrep low
) hi
(ratdisrep hi
)) ;; UGH, GAG WITH ME A SPOON
473 (if (not (symbolp ind
))
474 (merror (intl:gettext
"~:M: index must be a symbol; found ~M") (if sump
'$sum
'$product
) ind
))
476 (prog (u *i lind l
*i
*hl
)
477 (setq lind
(cons ind nil
))
479 ((not (fixnump (setq *hl
(mfuncall '$floor
(m- hi low
)))))
480 (if evaluate-summand
(setq expr
(mevalsumarg expr ind low hi
)))
481 (return (cons (if sump
'(%sum
) '(%product
))
482 (list expr ind low hi
))))
484 (return (if sump
0 1))))
485 (setq *i low l
*i
(list *i
) u
(if sump
0 1))
488 (add u
(resimplify (let* ((foo (mbinding (lind l
*i
) (meval expr
)))
489 (bar (subst-if-not-freeof *i ind foo
)))
491 (mul u
(resimplify (let* ((foo (mbinding (lind l
*i
) (meval expr
)))
492 (bar (subst-if-not-freeof *i ind foo
)))
494 (when (zerop *hl
) (return u
))
496 (setq *i
(car (rplaca l
*i
(m+ *i
1))))
499 (defun subst-if-not-freeof (x y expr
)
501 ;; suppressing substitution here avoids substituting for
502 ;; local variables recognize by freeof, e.g., formal argument of lambda.
504 (let ($simp
) (maxima-substitute x y expr
))))
506 (defun mevalsumarg (expr ind low hi
)
507 (if (let (($prederror nil
))
508 (eq (mevalp `((mlessp) ,hi
,low
)) t
))
511 (let ((gensym-ind (gensym)))
512 (if (apparently-integer low
)
513 (meval `(($declare
) ,gensym-ind $integer
)))
514 (assume (list '(mgeqp) gensym-ind low
))
515 (if (not (eq hi
'$inf
))
516 (assume (list '(mgeqp) hi gensym-ind
)))
517 (let ((msump t
) (foo) (summand))
519 (if (and (not (atom expr
)) (get (caar expr
) 'mevalsumarg-macro
))
520 (funcall (get (caar expr
) 'mevalsumarg-macro
) expr
)
523 (setq summand
($substitute gensym-ind ind summand
)))
524 (setq foo
(mbinding ((list gensym-ind
) (list gensym-ind
))
525 (resimplify (meval summand
))))
526 ;; At this point we do not switch off simplification to preserve
527 ;; the achieved simplification of the summand (DK 02/2010).
529 (setq foo
($substitute ind gensym-ind foo
)))
530 (if (not (eq hi
'$inf
))
531 (forget (list '(mgeqp) hi gensym-ind
)))
532 (forget (list '(mgeqp) gensym-ind low
))
533 (if (apparently-integer low
)
534 (meval `(($remove
) ,gensym-ind $integer
)))
537 (defun apparently-integer (x)
538 (or ($integerp x
) ($featurep x
'$integer
)))
541 (if (not (= (length l
) 4)) (wna-err op
))
542 (let ((ind (cadr l
)))
543 (if (mquotep ind
) (setq ind
(cadr ind
)))
544 (if (not (symbolp ind
))
545 (merror (intl:gettext
"~:M: index must be a symbol; found ~M") op ind
))
546 (let ((low (caddr l
))
548 (list (mevalsumarg (car l
) ind low hi
)
549 ind
(meval (caddr l
)) (meval (cadddr l
))))))
551 (defun simpsum1 (e k lo hi
)
552 (with-new-context (context)
553 (let ((acc 0) (n) (sgn) ($prederror nil
) (i (gensym)) (ex))
554 (setq lo
($ratdisrep lo
))
555 (setq hi
($ratdisrep hi
))
557 (setq n
($limit
(add 1 (sub hi lo
))))
560 (if (not (eq t
(csign lo
))) (mfuncall '$assume
`((mgeqp) ,i
,lo
)))
561 (if (not (eq t
(csign hi
))) (mfuncall '$assume
`((mgeqp) ,hi
,i
)))
563 (setq ex
(subst i k e
))
564 (setq ex
(subst i k ex
))
567 (cond ((and (eq n
'$inf
) ($freeof i ex
))
568 (setq sgn
(csign ex
))
569 (cond ((eq sgn
'$pos
) '$inf
)
570 ((eq sgn
'$neg
) '$minf
)
572 (t `((%sum simp
) ,ex
,i
,lo
,hi
))))
574 ((and (mbagp e
) $listarith
)
576 `((,(caar e
)) ,@(mapcar #'(lambda (s) (mfuncall '$sum s k lo hi
)) (margs e
))) t
))
578 ((or (eq sgn
'$neg
) (eq sgn
'$zero
) (eq sgn
'$nz
)) 0)
582 (($freeof i ex
) (mult n ex
))
584 ((and (integerp n
) (eq sgn
'$pos
) $simpsum
)
586 (setq acc
(add acc
(resimplify (subst (add j lo
) i ex
))))))
589 (setq ex
(subst '%sum
'$sum ex
))
590 `((%sum simp
) ,(subst k i ex
) ,k
,lo
,hi
))))
592 (setq acc
(subst k i acc
))
594 ;; If expression is still a summation,
595 ;; punt to previous simplification code.
597 (if (and $simpsum
(op-equalp acc
'$sum
'%sum
))
598 (let* ((args (cdr acc
)) (e (first args
)) (i (second args
)) (i0 (third args
)) (i1 (fourth args
)))
599 (setq acc
(simpsum1-save e i i0 i1
))))
603 (defun simpprod1 (e k lo hi
)
604 (with-new-context (context)
605 (let ((acc 1) (n) (sgn) ($prederror nil
) (i (gensym)) (ex) (ex-mag) (realp))
607 (setq lo
($ratdisrep lo
))
608 (setq hi
($ratdisrep hi
))
609 (setq n
($limit
(add 1 (sub hi lo
))))
612 (if (not (eq t
(csign lo
))) (mfuncall '$assume
`((mgeqp) ,i
,lo
)))
613 (if (not (eq t
(csign hi
))) (mfuncall '$assume
`((mgeqp) ,hi
,i
)))
615 (setq ex
(subst i k e
))
616 (setq ex
(subst i k ex
))
622 ((and (eq n
'$inf
) ($freeof i ex
))
623 (setq ex-mag
(mfuncall '$cabs ex
))
624 (setq realp
(mfuncall '$imagpart ex
))
625 (setq realp
(mevalp `((mequal) 0 ,realp
)))
627 (cond ((eq t
(mevalp `((mlessp) ,ex-mag
1))) 0)
628 ((and (eq realp t
) (eq t
(mevalp `((mgreaterp) ,ex
1)))) '$inf
)
629 ((eq t
(mevalp `((mgreaterp) ,ex-mag
1))) '$infinity
)
630 ((eq t
(mevalp `((mequal) 1 ,ex-mag
))) '$und
)
631 (t `((%product
) ,e
,k
,lo
,hi
))))
633 ((or (eq sgn
'$neg
) (eq sgn
'$zero
) (eq sgn
'$nz
))
636 ((and (mbagp e
) $listarith
)
638 `((,(caar e
)) ,@(mapcar #'(lambda (s) (mfuncall '$product s k lo hi
)) (margs e
))) t
))
640 (($freeof i ex
) (power ex n
))
642 ((and (integerp n
) (eq sgn
'$pos
) $simpproduct
)
644 (setq acc
(mult acc
(resimplify (subst (add j lo
) i ex
))))))
647 (setq ex
(subst '%product
'$product ex
))
648 `((%product simp
) ,(subst k i ex
) ,k
,lo
,hi
))))
650 ;; Hmm, this is curious... don't call existing product simplifications
651 ;; if index range is infinite -- what's up with that??
653 (if (and $simpproduct
(op-equalp acc
'$product
'%product
) (not (like n
'$inf
)))
654 (let* ((args (cdr acc
)) (e (first args
)) (i (second args
)) (i0 (third args
)) (i1 (fourth args
)))
655 (setq acc
(simpprod1-save e i i0 i1
))))
657 (setq acc
(subst k i acc
))
658 (setq acc
(subst '%product
'$product acc
))
662 ; This function was SIMPPROD1 until the sum/product code was revised Nov 2005.
663 ; The revised code punts back to this function since this code knows
664 ; some simplifications not handled by the revised code. -- Robert Dodier
666 (defun simpprod1-save (exp i lo hi
)
668 (cond ((not (symbolp i
)) (merror (intl:gettext
"product: index must be a symbol; found ~M") i
))
670 (let ((valist (list i
)))
671 (mbinding (valist (list hi
))
673 ((eq ($sign
(setq u
(m- hi lo
))) '$neg
)
674 (cond ((eq ($sign
(add2 u
1)) '$zero
) 1)
675 (t (merror (intl:gettext
"product: lower bound ~M greater than upper bound ~M") lo hi
))))
677 (cond ((null (eq exp i
))
678 (power* exp
(list '(mplus) hi
1 (list '(mtimes) -
1 lo
))))
679 ((let ((lot (asksign lo
)))
680 (cond ((equal lot
'$zero
) 0)
682 (m// (list '(mfactorial) hi
)
683 (list '(mfactorial) (list '(mplus) lo -
1))))
684 ((m* (list '(mfactorial)
686 (cond ((member (asksign hi
) '($zero $positive
) :test
#'eq
)
690 (setq hi
(list '(mabs) hi
)))))
691 (list '(mfactorial) hi
))))))))
692 ((list '(%product simp
) exp i lo hi
)))))
695 ;; multiplication of sums
697 (defun gensumindex ()
699 (intern (format nil
"~S~D" $genindex
(incf $gensumnum
)))
700 (intern (format nil
"~S" $genindex
))))
702 (defun sumtimes (x y
)
705 ((or (safe-zerop x
) (safe-zerop y
)) 0)
706 ((or (atom x
) (not (eq (caar x
) '%sum
))) (sumultin x y
))
707 ((or (atom y
) (not (eq (caar y
) '%sum
))) (sumultin y x
))
709 (if (great (sum-arg x
) (sum-arg y
)) (setq u y v x
) (setq u x v y
))
710 (setq i
(let ((ind (gensumindex)))
711 (setq u
(subst ind
(sum-index u
) u
)) ind
))
712 (setq j
(let ((ind (gensumindex)))
713 (setq v
(subst ind
(sum-index v
) v
)) ind
))
714 (if (and $cauchysum
(eq (sum-upper u
) '$inf
)
715 (eq (sum-upper v
) '$inf
))
718 (sumtimes (maxima-substitute j i
(sum-arg u
))
719 (maxima-substitute (m- i j
) j
(sum-arg v
)))
720 j
(sum-lower u
) (m- i
(sum-lower v
)))
721 i
(m+ (sum-lower u
) (sum-lower v
)) '$inf
)
723 (list '(%sum
) (sumtimes (sum-arg u
) (sum-arg v
))
724 j
(sum-lower v
) (sum-upper v
))
725 i
(sum-lower u
) (sum-upper u
)))))))
727 (defun sumultin (x s
) ; Multiplies x into a sum adjusting indices.
728 (cond ((or (atom s
) (not (eq (caar s
) '%sum
))) (m* x s
))
729 ((free x
(sum-index s
))
730 (list* (car s
) (sumultin x
(sum-arg s
)) (cddr s
)))
731 (t (let ((ind (gensumindex)))
733 (sumultin x
(subst ind
(sum-index s
) (sum-arg s
)))
739 (defun sumpls (sum out
)
741 (if (null out
) (return (cons sum nil
)))
742 (setq out
(setq l
(cons nil out
)))
743 a
(if (null (cdr out
)) (return (cons sum
(cdr l
))))
744 (and (not (atom (cadr out
)))
746 (eq (caar (cadr out
)) '%sum
)
747 (alike1 (sum-arg (cadr out
)) (sum-arg sum
))
748 (alike1 (sum-index (cadr out
)) (sum-index sum
))
749 (cond ((onediff (sum-upper (cadr out
)) (sum-lower sum
))
750 (setq sum
(list (car sum
)
753 (sum-lower (cadr out
))
755 (rplacd out
(cddr out
))
757 ((onediff (sum-upper sum
) (sum-lower (cadr out
)))
758 (setq sum
(list (car sum
)
762 (sum-upper (cadr out
))))
763 (rplacd out
(cddr out
))
771 (defun freesum (e b a q
)
772 (m* e q
(m- (m+ a
1) b
)))
774 ;; linear operator stuff
776 (defun oper-apply (e z
)
777 (cond ((null opers-list
)
778 (let ((w (get (caar e
) 'operators
)))
779 (if w
(funcall w e
1 z
) (simpargs e z
))))
780 ((get (caar e
) (caar opers-list
))
781 (let ((opers-list (cdr opers-list
))
782 (fun (cdar opers-list
)))
784 (t (let ((opers-list (cdr opers-list
)))
787 ;; Define an operator simplification, the same as antisymmetric, commutative, linear, etc.
788 ;; Here OP = operator name, FN = function of 1 argument to carry out operator-specific simplification.
789 ;; 1. push operator name onto OPERS
790 ;; 2. update $OPPROPERTIES
791 ;; 3. push operator name and glue code onto *OPERS-LIST
792 ;; 4. declare operator name as a feature, so declare(..., <op>) is recognized
794 (defmfun $define_opproperty
(op fn
)
796 (merror "define_opproperty: first argument must be a symbol; found: ~M" op
))
797 (unless (or (symbolp fn
) (and (consp fn
) (eq (caar fn
) 'lambda
)))
798 (merror "define_opproperty: second argument must be a symbol or lambda expression; found: ~M" fn
))
800 (setq $opproperties
(cons '(mlist simp
) (reverse opers
)))
802 ((fn-glue (coerce (if (symbolp fn
)
804 (declare (ignorable z
))
805 (if (or (fboundp ',fn
) (mget ',fn
'mexpr
))
806 (let ((e1 (let ($simp
) (mfuncall ',fn e
))))
807 (if ($mapatom e1
) e1
(oper-apply e1 nil
)))
808 (list '(,fn
) (let ((*opers-list
(cdr *opers-list
))) (oper-apply e z
)))))
811 (let ((e1 (let ($simp
) (mfuncall ',fn e
))))
812 (if ($mapatom e1
) e1
(oper-apply e1 nil
)))))
814 (push `(,op .
,fn-glue
) *opers-list
))
815 (mfuncall '$declare op
'$feature
))
817 (defun linearize1 (e z
) ; z = t means args already simplified.
818 (linearize2 (cons (car e
) (mapcar #'(lambda (q) (simpcheck q z
)) (cdr e
)))
822 (cond ((eq op
'mplus
) 0)
823 ((eq op
'mtimes
) 1)))
825 (defun rem-const (e) ;removes constantp stuff
826 (do ((l (cdr e
) (cdr l
))
827 (a (list (opident (caar e
))))
828 (b (list (opident (caar e
)))))
830 (cons (simplifya (cons (list (caar e
)) a
) nil
)
831 (simplifya (cons (list (caar e
)) b
) nil
)))
832 (if ($constantp
(car l
))
833 (setq a
(cons (car l
) a
))
834 (setq b
(cons (car l
) b
)))))
836 (defun linearize2 (e times
)
837 (cond ((linearconst e
))
838 ((atom (cadr e
)) (oper-apply e t
))
839 ((eq (caar (cadr e
)) 'mplus
)
840 (addn (mapcar #'(lambda (q)
841 (linearize2 (list* (car e
) q
(cddr e
)) nil
))
844 ((and (eq (caar (cadr e
)) 'mtimes
) (null times
))
845 (let ((z (if (and (cddr e
)
847 ($subvarp
(caddr e
))))
848 (partition (cadr e
) (caddr e
) 1)
849 (rem-const (cadr e
))))
851 (setq w
(linearize2 (list* (car e
)
852 (simplifya (cdr z
) t
)
855 (linearize3 w e
(car z
))))
856 (t (oper-apply e t
))))
858 (defun linearconst (e)
859 (if (or (mnump (cadr e
))
862 (or (atom (caddr e
)) (member 'array
(cdar (caddr e
)) :test
#'eq
))
863 (free (cadr e
) (caddr e
))))
864 (if (or (zerop1 (cadr e
))
865 (and (member (caar e
) '(%sum %integrate
) :test
#'eq
)
867 (or (eq (cadddr e
) '$minf
)
868 (member (car (cddddr e
)) '($inf $infinity
) :test
#'eq
))
869 (eq ($asksign
(cadr e
)) '$zero
)))
871 (let ((w (oper-apply (list* (car e
) 1 (cddr e
)) t
)))
872 (linearize3 w e
(cadr e
))))))
874 (defun linearize3 (w e x
)
876 (if (and (member w
'($inf $minf $infinity
) :test
#'eq
) (safe-zerop x
))
877 (merror (intl:gettext
"LINEARIZE3: undefined form 0*inf: ~M") e
))
878 (setq w
(mul2 (simplifya x t
) w
))
879 (cond ((or (atom w
) (getl (caar w
) '($outative $linear
))) (setq w1
1))
880 ((eq (caar w
) 'mtimes
)
881 (setq w1
(cons '(mtimes) nil
))
882 (do ((w2 (cdr w
) (cdr w2
)))
883 ((null w2
) (setq w1
(nreverse w1
)))
884 (if (or (atom (car w2
))
885 (not (getl (caaar w2
) '($outative $linear
))))
886 (setq w1
(cons (car w2
) w1
)))))
888 (if (and (not (atom w1
)) (or (among '$inf w1
) (among '$minf w1
)))
892 (setq opers
(cons '$additive opers
)
893 *opers-list
(cons '($additive . additive
) *opers-list
))
895 (defun rem-opers-p (p)
896 (cond ((eq (caar opers-list
) p
)
897 (setq opers-list
(cdr p
)))
898 ((do ((l opers-list
(cdr l
)))
900 (if (eq (caar (cdr l
)) p
)
901 (return (rplacd l
(cddr l
))))))))
903 (defun additive (e z
)
904 (cond ((get (caar e
) '$outative
) ; Really a linearize!
905 (setq opers-list
(copy-list opers-list
))
906 (rem-opers-p '$outative
)
909 (addn (mapcar #'(lambda (q)
910 (let ((opers-list *opers-list
))
911 (oper-apply (list* (car e
) q
(cddr e
)) z
)))
914 (t (oper-apply e z
))))
916 (setq opers
(cons '$multiplicative opers
)
917 *opers-list
(cons '($multiplicative . multiplicative
) *opers-list
))
919 (defun multiplicative (e z
)
920 (cond ((mtimesp (cadr e
))
921 (muln (mapcar #'(lambda (q)
922 (let ((opers-list *opers-list
))
923 (oper-apply (list* (car e
) q
(cddr e
)) z
)))
926 (t (oper-apply e z
))))
928 (setq opers
(cons '$outative opers
)
929 *opers-list
(cons '($outative . outative
) *opers-list
))
931 (defun outative (e z
)
932 (setq e
(cons (car e
) (mapcar #'(lambda (q) (simpcheck q z
)) (cdr e
))))
933 (cond ((get (caar e
) '$additive
)
934 (setq opers-list
(copy-list opers-list
))
935 (rem-opers-p '$additive
)
939 (let ((u (if (and (cddr e
)
941 ($subvarp
(caddr e
))))
942 (partition (cadr e
) (caddr e
) 1)
943 (rem-const (cadr e
))))
945 (setq w
(oper-apply (list* (car e
)
946 (simplifya (cdr u
) t
)
949 (linearize3 w e
(car u
))))
950 (t (oper-apply e t
))))
952 (defprop %sum t $outative
)
953 (defprop %sum t opers
)
954 (defprop %integrate t $outative
)
955 (defprop %integrate t opers
)
956 (defprop %limit t $outative
)
957 (defprop %limit t opers
)
959 (setq opers
(cons '$evenfun opers
)
960 *opers-list
(cons '($evenfun . evenfun
) *opers-list
))
962 (setq opers
(cons '$oddfun opers
)
963 *opers-list
(cons '($oddfun . oddfun
) *opers-list
))
966 (if (or (null (cdr e
)) (cddr e
))
967 (merror (intl:gettext
"Function declared 'even' takes exactly one argument; found ~M") e
))
968 (let ((arg (simpcheck (cadr e
) z
)))
969 (oper-apply (list (car e
) (if (mminusp arg
) (neg arg
) arg
)) t
)))
972 (if (or (null (cdr e
)) (cddr e
))
973 (merror (intl:gettext
"Function declared 'odd' takes exactly one argument; found ~M") e
))
974 (let ((arg (simpcheck (cadr e
) z
)))
975 (if (mminusp arg
) (neg (oper-apply (list (car e
) (neg arg
)) t
))
976 (oper-apply (list (car e
) arg
) t
))))
978 (setq opers
(cons '$commutative opers
)
979 *opers-list
(cons '($commutative . commutative1
) *opers-list
))
981 (setq opers
(cons '$symmetric opers
)
982 *opers-list
(cons '($symmetric . commutative1
) *opers-list
))
984 (defun commutative1 (e z
)
985 (oper-apply (cons (car e
)
987 (sort (mapcar #'(lambda (q) (simpcheck q z
))
992 (setq opers
(cons '$antisymmetric opers
)
993 *opers-list
(cons '($antisymmetric . antisym
) *opers-list
))
996 (when (and $dotscrules
(mnctimesp e
))
998 (setq e
(simpnct e
1 nil
))))
999 (if ($atom e
) e
(antisym1 e z
)))
1001 (defun antisym1 (e z
)
1002 (let ((antisym-sign nil
)
1003 (l (mapcar #'(lambda (q) (simpcheck q z
)) (cdr e
))))
1004 (when (or (not (eq (caar e
) 'mnctimes
)) (freel l
'mnctimes
))
1005 (multiple-value-setq (l antisym-sign
) (bbsort1 l
)))
1006 (cond ((equal l
0) 0)
1009 (setq e
(oper-apply (cons (car e
) l
) t
)))
1014 (prog (sl sl1 antisym-sign
)
1015 (if (or (null l
) (null (cdr l
))) (return (values l antisym-sign
))
1016 (setq sl
(list nil
(car l
))))
1017 loop
(setq l
(cdr l
))
1018 (if (null l
) (return (values (nreverse (cdr sl
)) antisym-sign
)))
1020 loop1
(cond ((null (cdr sl1
)) (rplacd sl1
(cons (car l
) nil
)))
1021 ((alike1 (car l
) (cadr sl1
)) (return (values 0 nil
)))
1022 ((great (car l
) (cadr sl1
)) (rplacd sl1
(cons (car l
) (cdr sl1
))))
1023 (t (setq antisym-sign
(not antisym-sign
) sl1
(cdr sl1
)) (go loop1
)))
1026 (setq opers
(cons '$nary opers
)
1027 *opers-list
(cons '($nary . nary1
) *opers-list
))
1030 (oper-apply (nary2 e z
) z
))
1034 ((l (cdr e
) (cdr l
)) (ans) (some-change))
1038 (nary2 (cons (car e
) (nreverse ans
)) z
)
1042 ans
(if (and (not (atom (car l
))) (eq (caaar l
) (caar e
)))
1044 (setq some-change t
)
1045 (nconc (reverse (cdar l
)) ans
))
1046 (cons (car l
) ans
)))))
1048 (setq opers
(cons '$lassociative opers
)
1049 *opers-list
(cons '($lassociative . lassociative
) *opers-list
))
1051 (defun lassociative (e z
)
1053 ((ans0 (oper-apply (cons (car e
) (total-nary e
)) z
))
1054 (ans (if (consp ans0
) (cdr ans0
))))
1055 (cond ((or (null (cddr ans
)) (not (eq (caar ans0
) (caar e
)))) ans0
)
1056 ((do ((newans (list (car e
) (car ans
) (cadr ans
))
1057 (list (car e
) newans
(car ans
)))
1058 (ans (cddr ans
) (cdr ans
)))
1059 ((null ans
) newans
))))))
1061 (setq opers
(cons '$rassociative opers
)
1062 *opers-list
(cons '($rassociative . rassociative
) *opers-list
))
1064 (defun rassociative (e z
)
1066 ((ans0 (oper-apply (cons (car e
) (total-nary e
)) z
))
1067 (ans (if (consp ans0
) (cdr ans0
))))
1068 (cond ((or (null (cddr ans
)) (not (eq (caar ans0
) (caar e
)))) ans0
)
1069 (t (setq ans
(nreverse ans
))
1070 (do ((newans (list (car e
) (cadr ans
) (car ans
))
1071 (list (car e
) (car ans
) newans
))
1072 (ans (cddr ans
) (cdr ans
)))
1073 ((null ans
) newans
))))))
1075 (defun total-nary (e)
1076 (do ((l (cdr e
) (cdr l
)) (ans))
1077 ((null l
) (nreverse ans
))
1078 (setq ans
(if (and (not (atom (car l
))) (eq (caaar l
) (caar e
)))
1079 (nconc (reverse (total-nary (car l
))) ans
)
1080 (cons (car l
) ans
)))))
