1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module limit
)
16 ;;; **************************************************************
18 ;;; ** LIMIT PACKAGE **
20 ;;; **************************************************************
22 ;;; I believe a large portion of this file is described in Paul
23 ;;; Wang's thesis, "Evaluation of Definite Integrals by Symbolic
24 ;;; Manipulation," MIT/LCS/TR-92, Oct. 1971. This can be found at
25 ;;; https://web.archive.org/web/20191019131847/https://apps.dtic.mil/dtic/tr/fulltext/u2/732005.pdf
28 ;;; TOP LEVEL FUNCTION(S): $LIMIT $LDEFINT
30 (declare-top (special errorsw origval $lhospitallim low
*
31 *indicator half%pi nn
* dn
* numer denom exp var val varlist
32 *zexptsimp? $tlimswitch $logarc taylored logcombed
33 $exponentialize lhp? lhcount $ratfac genvar
34 loginprod? $limsubst $logabs a context limit-assumptions
35 limit-top limitp integer-info old-integer-info $keepfloat $logexpand
))
37 (defconstant +behavior-count
+ 4)
38 (defvar *behavior-count-now
*)
39 (defvar *getsignl-asksign-ok
* nil
)
41 (load-macsyma-macros rzmac
)
43 (defmvar infinities
'($inf $minf $infinity
)
44 "The types of infinities recognized by Maxima.
45 INFINITY is complex infinity")
47 (defmvar real-infinities
'($inf $minf
)
48 "The real infinities, `inf' is positive infinity, `minf' negative infinity")
50 (defmvar infinitesimals
'($zeroa $zerob
)
51 "The infinitesimals recognized by Maxima. ZEROA zero from above,
52 ZEROB zero from below")
54 (defmvar simplimplus-problems
()
55 "A list of all problems in the stack of recursive calls to simplimplus.")
57 (defmvar limit-answers
()
58 "An association list for storing limit answers.")
60 (defmvar limit-using-taylor
()
61 "Is the current limit computation using taylor expansion?")
63 (defmvar preserve-direction
() "Makes `limit' return Direction info.")
65 (unless (boundp 'integer-info
) (setq integer-info
()))
67 ;; This should be made to give more information about the error.
69 ;; (cond (errorsw (throw 'errorsw t))
70 ;; (t (merror "Discontinuity Encountered"))))
72 ;;(DEFUN PUTLIMVAL (E V)
73 ;; (let ((exp (cons '(%limit) (list e var val))))
74 ;; (cond ((not (assolike exp limit-answers))
75 ;; (setq limit-answers (cons (cons exp v) limit-answers))
79 (defun putlimval (e v
&aux exp
)
80 (setq exp
`((%limit
) ,e
,var
,val
))
81 (unless (assolike exp limit-answers
)
82 (push (cons exp v
) limit-answers
))
86 (let ((exp (cons '(%limit
) (list e var val
))))
87 (assolike exp limit-answers
)))
89 (defmacro limit-catch
(exp var val
)
91 (let ((ans (catch 'errorsw
92 (catch 'limit
(limit ,exp
,var
,val
'think
)))))
93 (if (or (null ans
) (eq ans t
))
97 (defmfun $limit
(&rest args
)
98 (let ((first-try (apply #'toplevel-$limit args
)))
99 (if (and (consp first-try
) (eq (caar first-try
) '%limit
))
100 (let ((*getsignl-asksign-ok
* t
))
101 (apply #'toplevel-$limit args
))
104 (defun toplevel-$limit
(&rest args
)
105 (let ((limit-assumptions ())
106 (old-integer-info ())
109 (declare (special limit-assumptions old-integer-info
110 $keepfloat limit-top
))
112 (setq old-integer-info integer-info
)
113 (setq integer-info
()))
116 (let ((exp1 ()) (lhcount $lhospitallim
) (*behavior-count-now
* 0)
117 (exp ()) (var ()) (val ()) (dr ())
118 (*indicator
()) (taylored ()) (origval ())
119 (logcombed ()) (lhp?
())
120 (varlist ()) (ans ()) (genvar ()) (loginprod?
())
121 (limit-answers ()) (limitp t
) (simplimplus-problems ())
122 (lenargs (length args
))
124 (declare (special lhcount
*behavior-count-now
* exp var val
*indicator
125 taylored origval logcombed lhp?
126 varlist genvar loginprod? limitp
))
128 (unless (or (= lenargs
3) (= lenargs
4) (= lenargs
1))
130 ;; Is it a LIST of Things?
131 (when (setq ans
(apply #'limit-list args
))
133 (setq exp1
(specrepcheck (first args
)))
134 (when (and (atom exp1
)
135 (member exp1
'(nil t
)))
136 ;; The expression is 'T or 'NIL. Return immediately.
139 (setq var
(setq genfoo
(gensym)) ; Use a gensym. Not foo.
142 (setq var
(second args
))
143 (when ($constantp var
)
144 (merror (intl:gettext
"limit: second argument must be a variable, not a constant; found: ~M") var
))
145 (unless (or ($subvarp var
) (atom var
))
146 (merror (intl:gettext
"limit: variable must be a symbol or subscripted symbol; found: ~M") var
))
147 (setq val
(infsimp (third args
)))
148 ;; infsimp converts -inf to minf. it also converts -infinity to
149 ;; infinity, although perhaps this should generate the error below.
150 (when (and (not (atom val
))
151 (some #'(lambda (x) (not (freeof x val
)))
153 (merror (intl:gettext
"limit: third argument must be a finite value or one of: inf, minf, infinity; found: ~M") val
))
154 (when (eq val
'$zeroa
) (setq dr
'$plus
))
155 (when (eq val
'$zerob
) (setq dr
'$minus
))))
157 (unless (member (fourth args
) '($plus $minus
) :test
#'eq
)
158 (merror (intl:gettext
"limit: direction must be either 'plus' or 'minus'; found: ~M") (fourth args
)))
159 (setq dr
(fourth args
))))
160 (if (and (atom var
) (not (among var val
)))
162 (let ((realvar var
)) ;; Var is funny so make it a gensym.
164 (setq exp
(maxima-substitute var realvar exp1
))
165 (putprop var realvar
'limitsub
)))
166 (unless (or $limsubst
(eq var genfoo
))
167 (when (limunknown exp
)
168 (return `((%limit
) ,@(cons exp1
(cdr args
))))))
169 (setq varlist
(ncons var
) genvar nil origval val
)
170 ;; Transform limits to minf to limits to inf by
171 ;; replacing var with -var everywhere.
172 (when (eq val
'$minf
)
175 exp
(subin (m* -
1 var
) exp
)))
177 ;; Hide noun form of %derivative, %integrate.
178 (setq exp
(hide exp
))
180 ;; Transform the limit value.
181 (unless (infinityp val
)
183 (let ((*atp
* t
) (realvar var
))
184 ;; *atp* prevents substitution from applying to vars
185 ;; bound by %sum, %product, %integrate, %limit
187 (putprop var t
'internal
)
188 (setq exp
(maxima-substitute (m+ val var
) realvar exp
))))
189 (setq val
(cond ((eq dr
'$plus
) '$zeroa
)
190 ((eq dr
'$minus
) '$zerob
)
194 ;; Make assumptions about limit var being very small or very large.
195 ;; Assumptions are forgotten upon exit.
196 (unless (= lenargs
1)
197 (limit-context var val dr
))
199 ;; Resimplify in light of new assumptions.
200 (setq exp
(resimplify
204 (limitsimp ($expand exp
1 0) var
))))))
206 (if (not (or (real-epsilonp val
) ;; if direction of limit not specified
208 (setq ans
(both-side exp var val
)) ;; compute from both sides
209 (let ((d (catch 'mabs
(mabs-subst exp var val
))))
210 (cond ;; otherwise try to remove absolute value
211 ((eq d
'$und
) (return '$und
))
214 (setq ans
(limit-catch exp var val
));; and find limit from one side
218 (setq ans
(catch 'taylor-catch
219 (let ((silent-taylor-flag t
))
220 (declare (special silent-taylor-flag
))
221 (gruntz1 exp var val
)))))
223 ;; try taylor series expansion if simple limit didn't work
224 (if (and (null ans
) ;; if no limit found and
225 $tlimswitch
;; user says ok to use taylor and
226 (not limit-using-taylor
));; not already doing taylor
227 (let ((limit-using-taylor t
))
228 (declare (special limit-using-taylor
))
229 (setq ans
(limit-catch exp var val
))))))
232 (return (clean-limit-exp ans
))
233 (return (cons '(%limit
) args
))))) ;; failure: return nounform
234 (restore-assumptions))))
236 (defun clean-limit-exp (exp)
237 (setq exp
(restorelim exp
))
238 (if preserve-direction exp
(ridofab exp
)))
240 (defun limit-list (exp1 &rest rest
)
242 `(,(car exp1
) ,@(mapcar #'(lambda (x) (apply #'toplevel-$limit
`(,x
,@rest
))) (cdr exp1
)))
245 (defun limit-context (var val direction
) ;Only works on entry!
247 (assume '((mgreaterp) lim-epsilon
0))
248 (assume '((mgreaterp) prin-inf
100000000))
249 (setq limit-assumptions
(make-limit-assumptions var val direction
))
254 (defun make-limit-assumptions (var val direction
)
255 (let ((new-assumptions))
256 (cond ((or (null var
) (null val
))
258 ((and (not (infinityp val
)) (null direction
))
261 `(,(assume `((mgreaterp) ,var
100000000)) ,@new-assumptions
))
263 `(,(assume `((mgreaterp) -
100000000 ,var
)) ,@new-assumptions
))
264 ((eq direction
'$plus
)
265 `(,(assume `((mgreaterp) ,var
0)) ,@new-assumptions
)) ;All limits around 0
266 ((eq direction
'$minus
)
267 `(,(assume `((mgreaterp) 0 ,var
)) ,@new-assumptions
))
271 (defun restore-assumptions ()
272 ;;;Hackery until assume and forget take reliable args. Nov. 9 1979.
274 (do ((assumption-list limit-assumptions
(cdr assumption-list
)))
275 ((null assumption-list
) t
)
276 (forget (car assumption-list
)))
277 (forget '((mgreaterp) lim-epsilon
0))
278 (forget '((mgreaterp) prin-inf
100000000))
279 (cond ((and (not (null integer-info
))
281 (do ((list integer-info
(cdr list
)))
283 (i-$remove
`(,(cadar list
) ,(caddar list
))))
284 (setq integer-info old-integer-info
))))
286 ;; The optional arg allows the caller to decide on the value of
287 ;; preserve-direction. Default is T, like it used to be.
288 (defun both-side (exp var val
&optional
(preserve t
))
289 (let* ((preserve-direction preserve
)
290 (la (toplevel-$limit exp var val
'$plus
)) lb
)
291 (when (eq la
'$und
) (return-from both-side
'$und
))
292 (setf lb
(toplevel-$limit exp var val
'$minus
))
293 (let ((ra (ridofab la
))
295 (cond ((eq t
(meqp ra rb
))
299 ; Maxima does not consider equal(ind,ind) to be true, but
300 ; if both one-sided limits are ind then we want to call
301 ; the two-sided limit ind (e.g., limit(sin(1/x),x,0)).
303 ((or (not (free la
'%limit
))
304 (not (free lb
'%limit
)))
306 ((and (infinityp la
) (infinityp lb
))
307 ; inf + minf => infinity
312 (defun limunknown (f)
313 (catch 'limunknown
(limunknown1 (specrepcheck f
))))
315 (defun limunknown1 (f)
316 (cond ((mapatom f
) nil
)
317 ((or (not (safe-get (caar f
) 'operators
))
318 (member (caar f
) '(%sum %product mncexpt
) :test
#'eq
)
319 ;;Special function code here i.e. for li[2](x).
320 (and (eq (caar f
) 'mqapply
)
321 (not (get (subfunname f
) 'specsimp
))))
322 (if (not (free f var
)) (throw 'limunknown t
)))
323 (t (mapc #'limunknown1
(cdr f
)) nil
)))
326 (if (involve e
'(%gamma
)) (setq e
($makefact e
)))
327 (cond ((involve e
'(mfactorial))
328 (setq e
(simplify ($minfactorial e
))))
332 ;; or nil if sign unknown or complex
334 (let ((z (ridofab z
)))
335 (if (not (free z var
)) (setq z
(toplevel-$limit z var val
)))
336 (let ((*complexsign
* t
))
337 (let ((sign (if *getsignl-asksign-ok
* ($asksign z
) ($sign z
))))
338 (cond ((eq sign
'$pos
) 1)
340 ((eq sign
'$zero
) 0))))))
342 (defun restorelim (exp)
343 (cond ((null exp
) nil
)
344 ((atom exp
) (or (and (symbolp exp
) (get exp
'limitsub
)) exp
))
345 ((and (consp (car exp
)) (eq (caar exp
) 'mrat
))
347 (cons (restorelim (cadr exp
))
348 (restorelim (cddr exp
)))))
349 (t (cons (car exp
) (mapcar #'restorelim
(cdr exp
))))))
352 (defun mabs-subst (exp var val
) ; RETURNS EXP WITH MABS REMOVED, OR THROWS.
353 (let ((d (involve exp
'(mabs)))
357 ((not (and (equal ($imagpart
(let ((v (limit-catch d var val
)))
358 ;; The above call might
359 ;; throw 'limit, so we
360 ;; need to catch it. If
362 ;; limit without ABS, we
363 ;; assume the limit is
364 ;; undefined. Is this
365 ;; right? Anyway, this
366 ;; fixes Bug 1548643.
371 (equal ($imagpart var
) 0)))
372 (cond ((eq arglim
'$infinity
)
373 ;; Check for $infinity as limit of argument.
376 (throw 'mabs
'retn
))))
377 (t (do ((ans d
(involve exp
'(mabs))) (a () ()))
379 (setq a
(mabs-subst ans var val
))
380 (setq d
(limit a var val t
))
384 (setq d
(behavior a var val
))
385 (if (zerop1 d
) (throw 'mabs
'retn
))))
388 (cond ((or (eq d
'$zeroa
) (eq d
'$inf
)
390 ;; fails on limit(abs(sin(x))/sin(x), x, inf)
391 (eq ($sign d
) '$pos
))
392 (setq exp
(maxima-substitute a
`((mabs) ,ans
) exp
)))
393 ((or (eq d
'$zerob
) (eq d
'$minf
)
394 (eq ($sign d
) '$neg
))
395 (setq exp
(maxima-substitute (m* -
1 a
) `((mabs) ,ans
) exp
)))
396 (t (throw 'mabs
'retn
))))
397 (t (throw 'mabs
'retn
))))))))))
399 ;; Called on an expression that might contain $INF, $MINF, $ZEROA, $ZEROB. Tries
400 ;; to simplify it to sort out things like inf^inf or inf+1.
402 (simpinf-ic exp
(count-general-inf exp
)))
404 (defun count-general-inf (expr)
405 (count-atoms-matching
406 (lambda (x) (or (infinityp x
) (real-epsilonp x
))) expr
))
408 (defun count-atoms-matching (predicate expr
)
409 "Count the number of atoms in the Maxima expression EXPR matching PREDICATE,
410 ignoring dummy variables and array indices."
412 ((atom expr
) (if (funcall predicate expr
) 1 0))
413 ;; Don't count atoms that occur as a limit of %integrate, %sum, %product,
415 ((member (caar expr
) dummy-variable-operators
)
416 (count-atoms-matching predicate
(cadr expr
)))
417 ;; Ignore array indices
418 ((member 'array
(car expr
)) 0)
420 for arg in
(cdr expr
)
421 summing
(count-atoms-matching predicate arg
)))))
423 (defun simpinf-ic (exp &optional infinity-count
)
425 ;; A very slow identity transformation...
428 ;; If there's only one infinity, we replace it by a variable and take the
429 ;; limit as that variable goes to infinity. Use $gensym in case we can't
430 ;; compute the answer and the limit leaks out.
431 (1 (let* ((val (or (inf-typep exp
) (epsilon-typep exp
)))
433 (expr (subst var val exp
))
434 (limit (toplevel-$limit expr var val
)))
436 ;; Now we look to see whether the computed limit is any simpler than
437 ;; what we shoved in (which we'll define as "doesn't contain EXPR as a
438 ;; subtree"). If so, return it.
439 ((not (subtree-p expr limit
:test
#'equal
))
442 ;; Otherwise, return the original form: apparently, we can't compute
443 ;; the limit we needed, and it's uglier than what we started with.
446 ;; If more than one infinity, we have to be a bit more careful.
448 (let* ((arguments (mapcar 'simpinf
(cdr exp
)))
449 (new-expression (cons (list (caar exp
)) arguments
))
452 ;; If any of the arguments are undefined, we are too.
453 ((among '$und arguments
) '$und
)
454 ;; If we ended up with something indeterminate, we punt and just return
456 ((amongl '(%limit $ind
) arguments
) exp
)
458 ;; Exponentiation & multiplication
459 ((mexptp exp
) (simpinf-expt (first arguments
) (second arguments
)))
460 ((mtimesp exp
) (simpinf-times arguments
))
462 ;; Down to at most one infinity? We do this after exponentiation to
463 ;; avoid zeroa^zeroa => 0^0, which will raise an error rather than just
464 ;; returning und. We do it after multiplication to avoid zeroa * inf =>
466 ((<= (setf infinities-left
(count-general-inf new-expression
)) 1)
467 (simpinf-ic new-expression infinities-left
))
470 ((mplusp exp
) (simpinf-plus arguments
))
473 (t new-expression
))))))
475 (defun simpinf-times (arguments)
476 (declare (special exp var val
))
477 ;; When we have a product, we need to spot that zeroa * zerob = zerob, zeroa *
478 ;; inf = und etc. Note that (SIMPINF '$ZEROA) => 0, so a nonzero atom is not
479 ;; an infinitesimal. Moreover, we can assume that each of ARGUMENTS is either
480 ;; a number, computed successfully by the recursive SIMPINF call, or maybe a
481 ;; %LIMIT noun-form (in which case, we aren't going to be able to tell the
484 ((member 0 arguments
)
486 ((find-if #'infinityp arguments
) '$und
)
487 ((every #'atom arguments
) 0)
490 ((member '$infinity arguments
)
491 (if (every #'atom arguments
)
495 (t (simplimit (cons '(mtimes) arguments
) var val
))))
497 (defun simpinf-expt (base exponent
)
498 ;; In the comments below, zero* represents one of 0, zeroa, zerob.
500 ;; TODO: In some cases we give up too early. E.g. inf^(2 + 1/inf) => inf^2
501 ;; (which should simplify to inf)
511 ((0 $zeroa $zerob
) '$und
)
512 (t (list '(mexpt) base exponent
))))
513 ;; minf^inf = infinity <== Or should it be und?
516 ;; minf^foo = minf^foo
521 ((0 $zeroa $zerob
) '$und
)
522 (t (list '(mexpt) base exponent
))))
526 ;; zero*^foo = zero*^foo
531 ((0 $zeroa $zerob
) '$und
)
532 (t (list '(mexpt) base exponent
))))
533 ;; a^b where a is pretty much anything except for a naked
534 ;; inf,minf,zeroa,zerob or 0.
537 ;; When a isn't crazy, try a^b = e^(b log(a))
538 ((not (amongl (append infinitesimals infinities
) base
))
539 (simpinf (m^
'$%e
(m* exponent
`((%log
) ,base
)))))
541 ;; No idea. Just return what we've found so far.
542 (t (list '(mexpt) base exponent
))))))
544 (defun simpinf-plus (arguments)
545 ;; We know that none of the arguments are infinitesimals, since SIMPINF never
546 ;; returns one of them. As such, we partition our arguments into infinities
547 ;; and everything else. The latter won't have any "hidden" infinities like
548 ;; lim(x,x,inf), since SIMPINF gave up on anything containing a %lim already.
549 (let ((bigs) (others))
550 (dolist (arg arguments
)
551 (cond ((infinityp arg
) (push arg bigs
))
552 (t (push arg others
))))
554 ;; inf + minf or the like
555 ((cdr (setf bigs
(delete-duplicates bigs
))) '$und
)
556 ;; inf + smaller + stuff
558 ;; I don't think this can happen, since SIMPINF goes back to the start if
559 ;; there are fewer than two infinities in the arguments, but let's be
561 (t (cons '(mplus) others
)))))
563 ;; Simplify expression with zeroa or zerob.
564 (defun simpab (small)
565 (cond ((null small
) ())
566 ((member small
'($zeroa $zerob $inf $minf $infinity
) :test
#'eq
) small
)
567 ((not (free small
'$ind
)) '$ind
) ;Not exactly right but not
568 ((not (free small
'$und
)) '$und
) ;causing trouble now.
569 ((mapatom small
) small
)
570 (t (let ((preserve-direction t
)
571 (new-small (subst (m^
'$inf -
1) '$zeroa
572 (subst (m^
'$minf -
1) '$zerob small
))))
573 (simpinf new-small
)))))
576 ;;;*I* INDICATES: T => USE LIMIT1,THINK, NIL => USE SIMPLIMIT.
577 (defun limit (exp var val
*i
*)
579 ((among '$und exp
) '$und
)
582 ((not (among var exp
))
583 (cond ((amongl '($inf $minf $infinity $ind
) exp
)
585 ((amongl '($zeroa $zerob
) exp
)
586 ;; Simplify expression with zeroa or zerob.
590 (t (putlimval exp
(cond ((and limit-using-taylor
593 (taylim exp var val
*i
*))
594 ((ratp exp var
) (ratlim exp
))
595 ((or (eq *i
* t
) (radicalp exp var
))
596 (limit1 exp var val
))
598 (cond ((or (mtimesp exp
) (mexptp exp
))
599 (limit1 exp var val
))
600 (t (simplimit exp var val
))))
601 (t (simplimit exp var val
)))))))
603 (defun limitsimp (exp var
)
604 (limitsimp-expt (sin-sq-cos-sq-sub exp
) var
))
605 ;;Hack for sin(x)^2+cos(x)^2.
607 ;; if var appears in base and power of expt,
608 ;; push var into power of of expt
609 (defun limitsimp-expt (exp var
)
610 (cond ((or (atom exp
)
612 (freeof var exp
)) exp
)
614 (not (freeof var
(cadr exp
)))
615 (not (freeof var
(caddr exp
))))
616 (m^
'$%e
(simplify `((%log
) ,exp
))))
617 (t (subst0 (cons (cons (caar exp
) ())
618 (mapcar #'(lambda (x)
619 (limitsimp-expt x var
))
623 (defun sin-sq-cos-sq-sub (exp) ;Hack ... Hack
624 (let ((arg (involve exp
'(%sin %cos
))))
627 (t (let ((new-exp ($substitute
(m+t
1 (m- (m^t
`((%sin simp
) ,arg
) 2)))
628 (m^t
`((%cos simp
) ,arg
) 2)
630 (m+t
1 (m- (m^t
`((%cos simp
) ,arg
) 2)))
631 (m^t
`((%sin simp
) ,arg
) 2)
633 (cond ((not (involve new-exp
'(%sin %cos
))) new-exp
)
636 (defun expand-trigs (x var
)
639 ((and (or (eq (caar x
) '%sin
)
641 (not (free (cadr x
) var
)))
643 ((member 'array
(car x
))
644 ;; Some kind of array reference. Return it.
646 (t (simplify (cons (ncons (caar x
))
647 (mapcar #'(lambda (x)
648 (expand-trigs x var
))
653 (cond ((not (involve e
654 '(%cot %csc %binomial
655 %sec %coth %sech %csch
656 %acot %acsc %asec %acoth
658 %jacobi_ns %jacobi_nc %jacobi_cs
659 %jacobi_ds %jacobi_dc
)))
661 (t ($ratsimp
(tansc1 e
)))))
663 (defun tansc1 (e &aux tem
)
665 ((and (setq e
(cons (car e
) (mapcar 'tansc1
(cdr e
)))) ()))
666 ((setq tem
(assoc (caar e
) '((%cot . %tan
) (%coth . %tanh
)
667 (%sec . %cos
) (%sech . %cosh
)
668 (%csc . %sin
) (%csch . %sinh
)) :test
#'eq
))
669 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
)) -
1.
)))
670 ((setq tem
(assoc (caar e
) '((%jacobi_nc . %jacobi_cn
)
671 (%jacobi_ns . %jacobi_sn
)
672 (%jacobi_cs . %jacobi_sc
)
673 (%jacobi_ds . %jacobi_sd
)
674 (%jacobi_dc . %jacobi_cd
)) :test
#'eq
))
675 ;; Converts Jacobi elliptic function to its reciprocal
677 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
) (third e
)) -
1.
)))
678 ((setq tem
(member (caar e
) '(%sinh %cosh %tanh
) :test
#'eq
))
679 (let (($exponentialize t
))
681 ((setq tem
(assoc (caar e
) '((%acsc . %asin
) (%asec . %acos
)
682 (%acot . %atan
) (%acsch . %asinh
)
683 (%asech . %acosh
) (%acoth . %atanh
)) :test
#'eq
))
684 (list (ncons (cdr tem
)) (m^t
(cadr e
) -
1.
)))
685 ((and (eq (caar e
) '%binomial
) (among var
(cdr e
)))
686 (m// `((mfactorial) ,(cadr e
))
687 (m* `((mfactorial) ,(m+t
(cadr e
) (m- (caddr e
))))
688 `((mfactorial) ,(caddr e
)))))
692 (cond ((not (involve ex
'(%sin %cos %tan %asin %acos %atan
693 %sinh %cosh %tanh %asinh %acosh %atanh
)))
699 ((eq (caar ex
) '%sinh
)
700 (m// (m+ (m^
'$%e
(cadr ex
)) (m- (m^
'$%e
(m- (cadr ex
)))))
702 ((eq (caar ex
) '%cosh
)
703 (m// (m+ (m^
'$%e
(cadr ex
)) (m^
'$%e
(m- (cadr ex
))))
705 ((and (member (caar ex
)
706 '(%sin %cos %tan %asin %acos %atan %sinh
707 %cosh %tanh %asinh %acosh %atanh
) :test
#'eq
)
710 (t (cons (car ex
) (mapcar #'hyperex0
(cdr ex
))))))
715 ;;Used by tlimit also.
716 (defun limit1 (exp var val
)
718 (let ((lhprogress? lhp?
)
721 (cond ((setq ans
(and (not (atom exp
)) (getlimval exp
)))
723 ((and (not (infinityp val
)) (setq ans
(simplimsubst val exp
)))
726 ;;;NUMDEN* => (values numerator denominator)
727 (multiple-value-bind (n dn
)
729 (cond ((not (among var dn
))
730 (return (simplimit (m// (simplimit n var val
) dn
) var val
)))
732 (return (simplimit (m* n
(simplimexpt dn -
1 (simplimit dn var val
) -
1)) var val
)))
734 (/#alike n
(car lhprogress?
))
735 (/#alike dn
(cdr lhprogress?
)))
736 (throw 'lhospital nil
)))
737 (return (limit2 n dn var val
))))))
742 (let ((deriv (sdiff (m// e f
) var
)))
744 ((=0 ($ratsimp deriv
)) t
)
747 (defun limit2 (n dn var val
)
748 (prog (n1 d1 lim-sign gcp sheur-ans
)
749 (setq n
(hyperex n
) dn
(hyperex dn
))
750 ;;;Change to uniform limit call.
751 (cond ((infinityp val
)
752 (setq d1
(limit dn var val nil
))
753 (setq n1
(limit n var val nil
)))
754 (t (cond ((setq n1
(simplimsubst val n
)) nil
)
755 (t (setq n1
(limit n var val nil
))))
756 (cond ((setq d1
(simplimsubst val dn
)) nil
)
757 (t (setq d1
(limit dn var val nil
))))))
758 (cond ((or (null n1
) (null d1
)) (return nil
))
759 (t (setq n1
(sratsimp n1
) d1
(sratsimp d1
))))
760 (cond ((or (involve n
'(mfactorial)) (involve dn
'(mfactorial)))
761 (let ((ans (limfact2 n dn var val
)))
762 (cond (ans (return ans
))))))
763 (cond ((and (zerop2 n1
) (zerop2 d1
))
764 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
765 (return (colexpt n dn gcp
)))
766 ((and (real-epsilonp val
)
768 (not (free dn
'%log
)))
769 (return (liminv (m// n dn
))))
770 ((setq n1
(try-lhospital-quit n dn nil
))
772 ((and (zerop2 n1
) (not (member d1
'($ind $und
) :test
#'eq
))) (return 0))
774 (setq n1
(ridofab n1
))
775 (return (simplimtimes `(,n1
,(simplimexpt dn -
1 d1 -
1))))))
776 (setq n1
(ridofab n1
))
777 (setq d1
(ridofab d1
))
778 (cond ((or (eq d1
'$und
)
779 (and (eq n1
'$und
) (not (real-infinityp d1
))))
782 ;; At this point we have n1/$ind. Look if n1 is one of the
783 ;; infinities or zero.
784 (cond ((and (infinityp n1
) (eq ($sign dn
) '$pos
))
786 ((and (infinityp n1
) (eq ($sign dn
) '$neg
))
787 (return (simpinf (m* -
1 n1
))))
788 ((and (not (eq n1
'$ind
))
789 (eq ($csign n1
) '$zero
))
792 ((eq n1
'$ind
) (return (cond ((infinityp d1
) 0)
795 ((and (real-infinityp d1
) (member n1
'($inf $und $minf
) :test
#'eq
))
796 (cond ((and (not (atom dn
)) (not (atom n
))
797 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
798 (return (colexpt n dn gcp
)))
800 (or (involve dn
'(mfactorial %gamma
))
801 (involve n
'(mfactorial %gamma
))))
802 (return (limfact n dn
))))))
803 ((eq n1 d1
) (setq lim-sign
1) (go cp
))
804 (t (setq lim-sign -
1) (go cp
))))
805 ((and (infinityp d1
) (infinityp n1
))
806 (setq lim-sign
(if (or (eq d1
'$minf
) (eq n1
'$minf
)) -
1 1))
808 (t (return (simplimtimes `(,n1
,(m^ d1 -
1))))))
809 cp
(setq n
($expand n
) dn
($expand dn
))
811 (let ((new-n (m+l
(maxi (cdr n
)))))
812 (cond ((not (alike1 new-n n
))
813 (return (limit (m// new-n dn
) var val
'think
))))
817 (let ((new-dn (m+l
(maxi (cdr dn
)))))
818 (cond ((not (alike1 new-dn dn
))
819 (return (limit (m// n new-dn
) var val
'think
))))
822 (setq sheur-ans
(sheur0 n1 d1
))
823 (cond ((or (member sheur-ans
'($inf $zeroa
) :test
#'eq
)
824 (free sheur-ans var
))
825 (return (simplimtimes `(,lim-sign
,sheur-ans
))))
826 ((and (alike1 sheur-ans dn
)
828 ((member (setq n1
(cond ((expfactorp n1 d1
) (expfactor n1 d1 var
))
830 '($inf $zeroa
) :test
#'eq
)
832 ((not (null (setq n1
(cond ((expfactorp n dn
) (expfactor n dn var
))
835 ((and (alike1 sheur-ans dn
) (not (mplusp n
))))
836 ((not (alike1 sheur-ans
(m// n dn
)))
837 (return (simplimit (m// ($expand
(m// n sheur-ans
))
838 ($expand
(m// dn sheur-ans
)))
841 (cond ((and (not (and (eq val
'$inf
) (expp n
) (expp dn
)))
842 (setq n1
(try-lhospital-quit n dn nil
))
847 ;; Test whether both n and dn have form
848 ;; product of poly^poly
849 (defun expfactorp (n dn
)
850 (do ((llist (append (cond ((mtimesp n
) (cdr n
))
852 (cond ((mtimesp dn
) (cdr dn
))
855 (exp? t
) ;IS EVERY ELEMENT SO FAR
856 (factor nil
)) ;A POLY^POLY?
860 (setq factor
(car llist
))
861 (setq exp?
(or (polyinx factor var
())
863 (polyinx (cadr factor
) var
())
864 (polyinx (caddr factor
) var
()))))))
866 (defun expfactor (n dn var
) ;Attempts to evaluate limit by grouping
867 (prog (highest-deg) ; terms with similar exponents.
868 (let ((new-exp (exppoly n
))) ;exppoly unrats expon
869 (setq n
(car new-exp
) ;and rtns deg of expons
870 highest-deg
(cdr new-exp
)))
871 (cond ((null n
) (return nil
))) ;nil means expon is not
872 (let ((new-exp (exppoly dn
))) ;a rat func.
873 (setq dn
(car new-exp
)
874 highest-deg
(max highest-deg
(cdr new-exp
))))
876 (= highest-deg
0)) ; prevent infinite recursion
880 (degree highest-deg
(1- degree
))
887 (limit (m// numerator denominator
)
891 (let ((newnumer-factor (get-newexp&factors
895 (setq numerator
(car newnumer-factor
)
896 numfactors
(cdr newnumer-factor
)))
897 (let ((newdenom-factor (get-newexp&factors
901 (setq denominator
(car newdenom-factor
)
902 denfactors
(cdr newdenom-factor
)))
903 (setq answer
(simplimit (list '(mexpt)
905 (m// numfactors denfactors
))
906 (cond ((> degree
0) var
)
910 (cond ((member answer
'($ind $und
) :test
#'equal
)
911 ;; cannot handle limit(exp(x*%i)*x, x, inf);
913 ((member answer
'($inf $minf
) :test
#'equal
)
914 ;; 0, zeroa, zerob are passed through to next iteration
915 (return (simplimtimes (list (m// numerator denominator
) answer
)))))))))
917 (defun exppoly (exp) ;RETURNS EXPRESSION WITH UNRATTED EXPONENTS
921 (exp (cond ((mtimesp exp
)
925 ((null exp
) (cons new-exp highest-deg
))
926 (setq factor
(car exp
))
928 (m* (cond ((or (not (mexptp factor
))
929 (not (ratp (caddr factor
) var
)))
933 (ratdegree (caddr factor
))))
934 (m^
(cadr factor
) (unrat (caddr factor
)))))
937 (defun unrat (exp) ;RETURNS UNRATTED EXPRESION
938 (multiple-value-bind (n d
)
940 (let ((tem ($divide n d
)))
942 (m// (caddr tem
) d
)))))
944 (defun get-newexp&factors
(exp degree var
) ;RETURNS (CONS NEWEXP FACTORS)
945 (do ((terms (cond ((mtimesp exp
)(cdr exp
)) ; SUCH THAT
946 (t (ncons exp
))) ; NEWEXP*FACTORS^(VAR^DEGREE)
947 (cdr terms
)) ; IS EQUAL TO EXP.
954 (setq factor
(car terms
))
955 (cond ((not (mexptp factor
))
957 (setq factors
(m* factor factors
)))
958 (t (setq newexp
(m* factor newexp
)))))
960 (= (ratdegree (caddr factor
))
962 (setq factors
(m* (m^
(cadr factor
)
963 (leading-coef (caddr factor
)))
965 newexp
(m* (m^
(cadr factor
)
967 (m* (leading-coef (caddr factor
))
970 (t (setq newexp
(m* factor newexp
))))))
972 (defun leading-coef (rat)
973 (ratlim (m// rat
(m^ var
(ratdegree rat
)))))
975 (defun ratdegree (rat)
976 (multiple-value-bind (n d
)
978 (- (deg n
) (deg d
))))
980 (defun limfact2 (n d var val
)
981 (let ((n1 (reflect0 n var val
))
982 (d1 (reflect0 d var val
)))
983 (cond ((and (alike1 n n1
)
986 (t (limit (m// n1 d1
) var val
'think
)))))
988 ;; takes expression and returns operator at front with all flags removed
989 ;; except array flag.
990 ;; array flag must match for alike1 to consider two things to be the same.
991 ;; ((MTIMES SIMP) ... ) => (MTIMES)
992 ;; ((PSI SIMP ARRAY) 0) => (PSI ARRAY)
993 (defun operator-with-array-flag (exp)
994 (cond ((member 'array
(car exp
) :test
#'eq
)
995 (list (caar exp
) 'array
))
996 (t (list (caar exp
)))))
998 (defun reflect0 (exp var val
)
999 (cond ((atom exp
) exp
)
1000 ((and (eq (caar exp
) 'mfactorial
)
1001 (let ((argval (limit (cadr exp
) var val
'think
)))
1002 (or (eq argval
'$minf
)
1003 (and (numberp argval
)
1005 (reflect (cadr exp
)))
1006 (t (cons (operator-with-array-flag exp
)
1009 (reflect0 term var val
)))
1012 (defun reflect (arg)
1015 (m^
(list (ncons 'mfactorial
)
1019 (m^
(list (ncons '%sin
)
1023 (defun limfact (n d
)
1025 (setq n
(stirling0 n
)
1027 (setq ans
(toplevel-$limit
(m// n d
) var
'$inf
))
1028 (cond ((and (atom ans
)
1029 (not (member ans
'(und ind
) :test
#'eq
))) ans
)
1030 ((eq (caar ans
) '%limit
) ())
1033 ;; substitute asymptotic approximations for gamma, factorial, and
1035 (defun stirling0 (e)
1037 ((and (setq e
(cons (car e
) (mapcar 'stirling0
(cdr e
))))
1039 ((and (eq (caar e
) '%gamma
)
1040 (eq (limit (cadr e
) var val
'think
) '$inf
))
1041 (stirling (cadr e
)))
1042 ((and (eq (caar e
) 'mfactorial
)
1043 (eq (limit (cadr e
) var val
'think
) '$inf
))
1044 (m* (cadr e
) (stirling (cadr e
))))
1045 ((and (eq (caar e
) 'mqapply
) ;; polylogarithm
1046 (eq (subfunname e
) '$li
)
1047 (integerp (car (subfunsubs e
))))
1048 (li-asymptotic-expansion (m- (car (subfunsubs e
)) 1)
1049 (car (subfunsubs e
))
1050 (car (subfunargs e
))))
1054 (maxima-substitute x
'$z
1056 ((mexpt simp
) 2 ((rat simp
) 1 2))
1057 ((mexpt simp
) $%pi
((rat simp
) 1 2))
1058 ((mexpt simp
) $z
((mplus simp
) ((rat simp
) -
1 2) $z
))
1059 ((mexpt simp
) $%e
((mtimes simp
) -
1 $z
)))))
1061 (defun no-err-sub (v e
&aux ans
)
1062 (let ((errorsw t
) (*zexptsimp? t
)
1064 ;; Don't print any error messages
1066 (declare (special errcatch
))
1067 ;; Should we just use IGNORE-ERRORS instead HANDLER-CASE here? I
1068 ;; (rtoy) am choosing the latter so that unexpected errors will
1069 ;; actually show up instead of being silently discarded.
1071 (setq ans
(catch 'errorsw
1073 (sratsimp (subin v e
)))))
1076 (cond ((null ans
) t
) ; Ratfun package returns NIL for failure.
1079 ;; substitute value v for var into expression e.
1080 ;; if result is defined and e is continuous, we have the limit.
1081 (defun simplimsubst (v e
)
1083 (cond ((involve e
'(mfactorial)) nil
)
1085 ;; functions that are defined at their discontinuities
1086 ((amongl '($atan2 $floor %round $ceiling %signum %integrate
1090 ;; substitute value into expression
1091 ((eq (setq ans
(no-err-sub (ridofab v
) e
)) t
)
1094 ((and (member v
'($zeroa $zerob
) :test
#'eq
) (=0 ($radcan ans
)))
1095 (setq ans
(behavior e var v
))
1096 (cond ((equal ans
1) '$zeroa
)
1097 ((equal ans -
1) '$zerob
)
1098 (t nil
))) ; behavior can't find direction
1101 ;;;returns (cons numerator denominator)
1103 (let ((e (factor (simplify e
)))
1109 (mapc #'forq
(cdr e
)))
1115 (setq numer
(car numer
)))
1117 (setq numer
(m*l numer
))))
1121 (setq denom
(car denom
)))
1123 (setq denom
(m*l denom
))))
1124 (values (factor numer
) (factor denom
))))
1126 ;;;FACTOR OR QUOTIENT
1127 ;;;Setq's the special vars numer and denom from numden*
1129 (cond ((and (mexptp e
)
1130 (not (freeof var e
))
1131 (null (pos-neg-p (caddr e
))))
1132 (push (m^
(cadr e
) (m* -
1.
(caddr e
))) denom
))
1133 (t (push e numer
))))
1135 ;;;Predicate to tell whether an expression is pos,zero or neg as var -> val.
1136 ;;;returns T if pos,zero. () if negative or don't know.
1137 (defun pos-neg-p (exp)
1138 (let ((ans (limit exp var val
'think
)))
1139 (cond ((and (not (member ans
'($und $ind $infinity
) :test
#'eq
))
1140 (equal ($imagpart ans
) 0))
1141 (let ((sign (getsignl ans
)))
1142 (cond ((or (equal sign
1)
1145 ((equal sign -
1) nil
))))
1148 (declare-top (unspecial n dn
))
1151 (cond ((radicalp e var
) nil
)
1152 ((member (caar e
) '(%log %sin %cos %tan %sinh %cosh %tanh mfactorial
1153 %asin %acos %atan %asinh %acosh %atanh
) :test
#'eq
) nil
)
1155 ((do ((e (cdr e
) (cdr e
)))
1157 (and (expp (car e
)) (return t
))))))
1161 (radicalp (cadr e
) var
)
1162 (among var
(caddr e
))
1163 (radicalp (caddr e
) var
)))
1166 (defun gcpower (a b
)
1167 ($gcd
(getexp a
) (getexp b
)))
1170 (cond ((and (mexptp exp
)
1171 (free (caddr exp
) var
)
1172 (eq (ask-integer (caddr exp
) '$integer
) '$yes
))
1174 ((mtimesp exp
) (getexplist (cdr exp
)))
1177 (defun getexplist (list)
1178 (cond ((null (cdr list
))
1179 (getexp (car list
)))
1180 (t ($gcd
(getexp (car list
))
1181 (getexplist (cdr list
))))))
1183 (defun limroot (exp power
)
1184 (cond ((or (atom exp
) (not (member (caar exp
) '(mtimes mexpt
) :test
#'eq
)))
1185 (limroot (list '(mexpt) exp
1) power
)) ;This is strange-JIM.
1186 ((mexptp exp
) (m^
(cadr exp
)
1187 (sratsimp (m* (caddr exp
) (m^ power -
1.
)))))
1188 (t (m*l
(mapcar #'(lambda (x)
1192 ;;NUMERATOR AND DENOMINATOR HAVE EXPONENTS WITH GCD OF GCP.
1193 ;;; Used to call simplimit but some of the transformations used here
1194 ;;; were not stable w.r.t. the simplifier, so try keeping exponent separate
1197 (defun colexpt (n dn gcp
)
1198 (let ((bas (m* (limroot n gcp
) (limroot dn
(m* -
1 gcp
))))
1201 (setq baslim
(limit bas var val
'think
))
1202 (setq expolim
(limit expo var val
'think
))
1203 (simplimexpt bas expo baslim expolim
)))
1205 ;;; This function will transform an expression such that either all logarithms
1206 ;;; contain arguments not becoming infinite or are of the form
1207 ;;; LOG(LOG( ... LOG(VAR))) This reduction takes place only over the operators
1208 ;;; MPLUS, MTIMES, MEXPT, and %LOG.
1210 (defun log-red-contract (facs)
1211 (do ((l facs
(cdr l
))
1215 (if log
(cons (cadr log
) (m*l consts
))
1217 (cond ((freeof var
(car l
)) (push (car l
) consts
))
1219 (if (null log
) (setq log
(car l
))
1223 (defun log-reduce (x)
1227 (do ((l (cdr x
) (cdr l
))
1232 ((null l
) (m+l
`(((%log
) ,(m*l strong-logs
))
1233 ((%log
) ,(m*l weak-logs
))
1235 (setq x
(log-reduce (car l
)))
1237 (if (infinityp (limit (cadr x
) var val
'think
))
1238 (push (cadr x
) strong-logs
)
1239 (push (cadr x
) weak-logs
)))
1240 ((and (mtimesp x
) (setq temp
(log-red-contract (cdr x
))))
1241 (if (infinityp (limit (car temp
) var val
'think
))
1242 (push (m^
(car temp
) (cdr temp
)) strong-logs
)
1243 (push (m^
(car temp
) (cdr temp
)) weak-logs
)))
1246 (do ((l (cdr x
) (cdr l
))
1249 (setq ans
($expand
(m* (log-reduce (car l
)) ans
)))))
1250 ((mexptp x
) (m^t
(log-reduce (cadr x
)) (caddr x
)))
1252 (cond ((not (infinityp (limit (cadr x
) var val
'think
))) x
)
1254 (cond ((eq (cadr x
) var
) x
)
1256 (let ((strongl (maxi (cdadr x
))))
1257 (m+ (log-reduce `((%log
) ,(car strongl
))) `((%log
) ,(m// (cadr x
) (car strongl
))))))
1259 (do ((l (cdadr x
) (cdr l
)) (ans 0)) ((null l
) ans
)
1260 (setq ans
(m+ (log-reduce (simplify `((%log
) ,(log-reduce (car l
))))) ans
))))
1262 (let ((red-log (simplify `((%log
) ,(log-reduce (cadr x
))))))
1263 (if (alike1 red-log x
) x
(log-reduce red-log
))))))))
1266 ;; this function is responsible for the following bug:
1267 ;; limit(x^2 + %i*x, x, inf) -> inf (should be infinity)
1269 (cond ((member val
'($inf $infinity
) :test
#'eq
)
1270 (setq e
(maxima-substitute (m^t
'x -
1) var e
)))
1272 (setq e
(maxima-substitute (m^t -
1 (m^t
'x -
1)) var e
)))
1274 (setq e
(maxima-substitute (m- 'x
) var e
)))
1276 (setq e
(maxima-substitute 'x var e
)))
1277 ((setq e
(maxima-substitute (m+t
'x val
) var e
))))
1278 (destructuring-let* ((e (let (($ratfac
()))
1279 ($rat
(sratsimp e
) 'x
)))
1287 (sratsimp (m// ($ratdisrep
`(,h
,(locoef n g
) .
1))
1288 ($ratdisrep
`(,h
,(locoef d g
) .
1))))))
1290 (cond ((not (member val
'($zerob $zeroa $inf $minf
) :test
#'eq
))
1292 ((not (equal ($imagpart e
) 0))
1294 ((null (setq e
(getsignl ($realpart e
))))
1296 ((equal e
1) '$zeroa
)
1297 ((equal e -
1) '$zerob
)
1300 ((not (member val
'($zerob $zeroa $infinity $inf $minf
) :test
#'eq
))
1302 ((eq val
'$infinity
) '$infinity
)
1303 ((not (equal ($imagpart e
) 0)) '$infinity
)
1304 ((null (setq e
(getsignl ($realpart e
))))
1307 ((equal e -
1) '$minf
)
1311 (if (or (atom n
) (not (eq (car n
) x
)))
1316 (if (or (atom n
) (not (eq (car n
) x
)))
1320 (defun behavior (exp var val
) ; returns either -1, 0, 1.
1321 (if (= *behavior-count-now
* +behavior-count
+)
1323 (let ((*behavior-count-now
* (1+ *behavior-count-now
*)) pair sign
)
1324 (cond ((real-infinityp val
)
1325 (setq val
(cond ((eq val
'$inf
) '$zeroa
)
1326 ((eq val
'$minf
) '$zerob
)))
1327 (setq exp
(sratsimp (subin (m^ var -
1) exp
)))))
1328 (cond ((eq val
'$infinity
) 0) ; Needs more hacking for complex.
1330 (prog2 (setq pair
(partition exp var
1))
1331 (not (mtimesp (cdr pair
)))))
1332 (setq sign
(getsignl (car pair
)))
1333 (if (not (fixnump sign
))
1335 (mul sign
(behavior (cdr pair
) var val
))))
1336 ((and (=0 (no-err-sub (ridofab val
) exp
))
1338 (free (caddr exp
) var
)
1339 (equal (getsignl (caddr exp
)) 1))
1340 (let ((bas (cadr exp
)) (expo (caddr exp
)))
1341 (behavior-expt bas expo
)))
1342 (t (behavior-by-diff exp var val
))))))
1344 (defun behavior-expt (bas expo
)
1345 (let ((behavior (behavior bas var val
)))
1346 (cond ((= behavior
1) 1)
1348 ((eq (ask-integer expo
'$integer
) '$yes
)
1349 (cond ((eq (ask-integer expo
'$even
) '$yes
) 1)
1352 (cond ((evenp (cadr expo
)) 1)
1353 ((oddp (caddr expo
)) behavior
)
1357 (defun behavior-by-diff (exp var val
)
1358 (cond ((not (or (eq val
'$zeroa
) (eq val
'$zerob
))) 0)
1359 (t (let ((old-val val
) (old-exp exp
))
1360 (setq val
(ridofab val
))
1362 (exp (sratsimp (sdiff exp var
)) (sratsimp (sdiff exp var
)))
1364 (ans ())) ; This do wins by a return.
1365 ((> ct
0) 0) ; This loop used to run up to 5 times,
1366 ;; but the size of some expressions would blow up.
1367 (setq ans
(no-err-sub val exp
)) ;Why not do an EVENFN and ODDFN
1370 (return (behavior-numden old-exp var old-val
)))
1371 ((=0 ans
) ()) ;Do it again.
1372 (t (setq ans
(getsignl ans
))
1373 (cond (n (return ans
))
1375 (return (if (eq old-val
'$zeroa
) 1 -
1)))
1377 (return (if (eq old-val
'$zeroa
) -
1 1)))
1378 (t (return 0))))))))))
1380 (defun behavior-numden (exp var val
)
1381 (let ((num ($num exp
)) (denom ($denom exp
)))
1382 (cond ((equal denom
1) 0) ;Could be hacked more from here.
1383 (t (let ((num-behav (behavior num var val
))
1384 (denom-behav (behavior denom var val
)))
1385 (cond ((or (= num-behav
0) (= denom-behav
0)) 0)
1386 ((= num-behav denom-behav
) 1)
1389 (defun try-lhospital (n d ind
)
1390 ;;Make one catch for the whole bunch of lhospital trials.
1391 (let ((ans (lhospital-catch n d ind
)))
1392 (cond ((null ans
) ())
1393 ((not (free-infp ans
)) (simpinf ans
))
1394 ((not (free-epsilonp ans
)) (simpab ans
))
1397 (defun try-lhospital-quit (n d ind
)
1398 (let ((ans (or (lhospital-catch n d ind
)
1399 (lhospital-catch (m^ d -
1) (m^ n -
1) ind
))))
1400 (cond ((null ans
) (throw 'limit t
))
1401 ((not (free-infp ans
)) (simpinf ans
))
1402 ((not (free-epsilonp ans
)) (simpab ans
))
1405 (defun lhospital-catch (n d ind
)
1406 (cond ((> 0 lhcount
)
1407 (setq lhcount $lhospitallim
)
1408 (throw 'lhospital nil
))
1409 ((equal lhcount $lhospitallim
)
1410 (let ((lhcount (m+ lhcount -
1)))
1411 (catch 'lhospital
(lhospital n d ind
))))
1412 (t (setq lhcount
(m+ lhcount -
1))
1413 (prog1 (lhospital n d ind
)
1414 (setq lhcount
(m+ lhcount
1))))))
1415 ;;If this succeeds then raise LHCOUNT.
1418 (defun lhospital (n d ind
)
1419 (declare (special val lhp?
))
1421 (setq n
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr n
)))))
1423 (setq d
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr d
)))))
1424 (multiple-value-bind (n d
)
1426 (let (const nconst dconst
)
1427 (setq lhp?
(and (null ind
) (cons n d
)))
1428 (multiple-value-setq (nconst n
) (var-or-const n
))
1429 (multiple-value-setq (dconst d
) (var-or-const d
))
1431 (setq n
(stirling0 n
)) ;; replace factorial and %gamma
1432 (setq d
(stirling0 d
)) ;; with approximations
1434 (setq n
(sdiff n var
) ;; take derivatives for l'hospital
1437 (if (or (not (free n
'%derivative
)) (not (free d
'%derivative
)))
1438 (throw 'lhospital
()))
1439 (setq n
(expand-trigs (tansc n
) var
))
1440 (setq d
(expand-trigs (tansc d
) var
))
1442 (multiple-value-setq (const n d
) (remove-singularities n d
))
1443 (setq const
(m* const
(m// nconst dconst
)))
1444 (simpinf (let ((ans (if ind
1445 (limit2 n d var val
)
1446 (limit-numden n d val
))))
1447 ;; When the limit function returns, it's possible that it will return NIL
1448 ;; (gave up without finding a limit). It's also possible that it will
1449 ;; return something containing UND. We treat that as a failure too.
1450 (when (and ans
(freeof '$und ans
))
1451 (m* const ans
)))))))
1453 ;; Try to compute the limit of a quotient NUM/DEN, trying to massage the input
1454 ;; into a convenient form for LIMIT on the way.
1455 (defun limit-numden (n d val
)
1457 ;; For general arguments, the best approach seems to be to use
1458 ;; sratsimp to simplify the quotient as much as we can, then
1459 ;; $multthru, which splits it up into a sum (presumably useful
1460 ;; because limit(a+b) = limit(a) + limit(b) if the limits exist, and
1461 ;; the right hand side might be easier to calculate)
1463 ($multthru
(sratsimp (m// n d
))))
1465 ;; If we've already got a sum in the numerator, it seems to be
1466 ;; better not to recombine it. Call LIMIT on the whole lot, though,
1467 ;; because terms with infinite limits might cancel to give a finite
1470 (m+l
(mapcar #'(lambda (x)
1471 (sratsimp (m// x d
)))
1474 (limit expr var val
'think
)))
1476 ;; Heuristics for picking the right way to express a LHOSPITAL problem.
1477 (defun lhop-numden (num denom
)
1478 (declare (special var
))
1479 (cond ((let ((log-num (involve num
'(%log
))))
1480 (cond ((null log-num
) ())
1481 ((lessthan (num-of-logs (factor (sratsimp (sdiff (m^ num -
1) var
))))
1482 (num-of-logs (factor (sratsimp (sdiff num var
)))))
1483 (psetq num
(m^ denom -
1) denom
(m^ num -
1))
1486 ((let ((log-denom (involve denom
'(%log
))))
1487 (cond ((null log-denom
) ())
1488 ((lessthan (num-of-logs (sratsimp (sdiff (m^ denom -
1) var
)))
1489 (num-of-logs (sratsimp (sdiff denom var
))))
1490 (psetq denom
(m^ num -
1) num
(m^ denom -
1))
1493 ((let ((exp-num (%einvolve num
)))
1495 (cond ((%e-right-placep exp-num
)
1497 (t (psetq num
(m^ denom -
1)
1498 denom
(m^ num -
1)) t
)))
1500 ((let ((exp-den (%einvolve denom
)))
1502 (cond ((%e-right-placep exp-den
)
1504 (t (psetq num
(m^ denom -
1)
1505 denom
(m^ num -
1)) t
)))
1507 ((let ((scnum (involve num
'(%sin
))))
1508 (cond (scnum (cond ((trig-right-placep '%sin scnum
) t
)
1509 (t (psetq num
(m^ denom -
1)
1510 denom
(m^ num -
1)) t
)))
1512 ((let ((scden (involve denom
'(%sin
))))
1513 (cond (scden (cond ((trig-right-placep '%sin scden
) t
)
1514 (t (psetq num
(m^ denom -
1)
1515 denom
(m^ num -
1)) t
)))
1517 ((let ((scnum (involve num
'(%asin %acos %atan
))))
1518 ;; If the numerator contains an inverse trig and the
1519 ;; denominator or reciprocal of denominator is polynomial,
1520 ;; leave everything as is. If the inverse trig is moved to
1521 ;; the denominator, things get messy, even if the numerator
1522 ;; becomes a polynomial. This is not perfect.
1523 (cond ((and scnum
(or (polyinx denom var
())
1524 (polyinx (m^ denom -
1) var
())))
1527 ((or (oscip num
) (oscip denom
)))
1529 (psetq num
(m^ denom -
1) denom
(m^ num -
1))))
1532 ;;i don't know what to do here for some cases, may have to be refined.
1533 (defun num-of-logs (exp)
1534 (cond ((mapatom exp
) 0)
1535 ((equal (caar exp
) '%log
)
1536 (m+ 1 (num-of-log-l (cdr exp
))))
1537 ((and (mexptp exp
) (mnump (caddr exp
)))
1538 (m* (simplify `((mabs) ,(caddr exp
)))
1539 (num-of-logs (cadr exp
))))
1540 (t (num-of-log-l (cdr exp
)))))
1542 (defun num-of-log-l (llist)
1543 (do ((temp llist
(cdr temp
)) (ans 0))
1545 (setq ans
(m+ ans
(num-of-logs (car temp
))))))
1547 (defun %e-right-placep
(%e-arg
)
1548 (let ((%e-arg-diff
(sdiff %e-arg var
)))
1550 ((free %e-arg-diff var
)) ;simple cases
1551 ((or (and (mexptp denom
)
1552 (equal (cadr denom
) -
1))
1553 (polyinx (m^ denom -
1) var
())) ())
1554 ((let ((%e-arg-diff-lim
(ridofab (limit %e-arg-diff var val
'think
)))
1555 (%e-arg-exp-lim
(ridofab (limit (m^
'$%e %e-arg
) var val
'think
))))
1558 (format t
"%e-arg-dif-lim = ~A~%" %e-arg-diff-lim
)
1559 (format t
"%e-arg-exp-lim = ~A~%" %e-arg-exp-lim
))
1560 (cond ((equal %e-arg-diff-lim %e-arg-exp-lim
)
1562 ((and (mnump %e-arg-diff-lim
) (mnump %e-arg-exp-lim
))
1564 ((and (mnump %e-arg-diff-lim
) (infinityp %e-arg-exp-lim
))
1565 ;; This is meant to make maxima handle bug 1469411
1566 ;; correctly. Undoubtedly, this needs work.
1570 (defun trig-right-placep (trig-type arg
)
1571 (let ((arglim (ridofab (limit arg var val
'think
)))
1572 (triglim (ridofab (limit `((,trig-type
) ,arg
) var val
'think
))))
1573 (cond ((and (equal arglim
0) (equal triglim
0)) t
)
1574 ((and (infinityp arglim
) (infinityp triglim
)) t
)
1577 ;;Takes a numerator and a denominator. If they tries all combinations of
1578 ;;products to try and make a simpler set of subproblems for LHOSPITAL.
1579 (defun remove-singularities (numer denom
)
1580 (cond ((or (null numer
) (null denom
)
1581 (atom numer
) (atom denom
)
1582 (not (mtimesp numer
)) ;Leave this here for a while.
1583 (not (mtimesp denom
)))
1584 (values 1 numer denom
))
1587 (multiple-value-bind (num-consts num-vars
)
1588 (var-or-const numer
)
1589 (multiple-value-bind (denom-consts denom-vars
)
1590 (var-or-const denom
)
1591 (if (not (mtimesp num-vars
))
1592 (setq num-vars
(list num-vars
))
1593 (setq num-vars
(cdr num-vars
)))
1594 (if (not (mtimesp denom-vars
))
1595 (setq denom-vars
(list denom-vars
))
1596 (setq denom-vars
(cdr denom-vars
)))
1597 (do ((nl num-vars
(cdr nl
))
1598 (num-list (copy-list num-vars
))
1599 (den-list denom-vars den-list-temp
)
1600 (den-list-temp (copy-list denom-vars
)))
1601 ((null nl
) (values (m* const
(m// num-consts denom-consts
))
1603 (m*l den-list-temp
)))
1604 (do ((dl den-list
(cdr dl
)))
1606 (if (or (%einvolve
(car nl
)) (%einvolve
(car nl
)))
1608 (let ((lim (catch 'limit
(simpinf (simpab (limit (m// (car nl
) (car dl
))
1609 var val
'think
))))))
1610 (cond ((or (eq lim t
)
1612 (equal (ridofab lim
) 0)
1614 (not (free lim
'$inf
))
1615 (not (free lim
'$minf
))
1616 (not (free lim
'$infinity
))
1617 (not (free lim
'$ind
))
1618 (not (free lim
'$und
)))
1621 (setq const
(m* lim const
))
1622 (setq num-list
(delete (car nl
) num-list
:count
1 :test
#'equal
))
1623 (setq den-list-temp
(delete (car dl
) den-list-temp
:count
1 :test
#'equal
))
1624 (return t
)))))))))))))
1626 ;; separate terms that contain var from constant terms
1627 ;; returns (const-terms . var-terms)
1628 (defun var-or-const (expr)
1629 (setq expr
($factor expr
))
1637 (do ((l (cdr expr
) (cdr l
))
1640 ((null l
) (values const varl
))
1641 (if (free (car l
) var
)
1642 (setq const
(m* (car l
) const
))
1643 (setq varl
(m* (car l
) varl
)))))
1647 ;; if term goes to non-zero constant, replace with constant
1648 (defun lhsimp (term var val
)
1649 (cond ((atom term
) term
)
1651 (let ((term-value (ridofab (limit term var val
'think
))))
1652 (cond ((not (member term-value
1653 '($inf $minf $und $ind $infinity
0)))
1657 (defun bylog (expo bas
)
1660 (try-lhospital-quit (simplify `((%log
) ,(tansc bas
)))
1665 (defun simplimexpt (bas expo bl el
)
1666 (cond ((or (eq bl
'$und
) (eq el
'$und
)) '$und
)
1668 (cond ((eq el
'$inf
) (if (eq bl
'$zeroa
) bl
0))
1669 ((eq el
'$minf
) (if (eq bl
'$zeroa
) '$inf
'$infinity
))
1670 ((eq el
'$ind
) '$ind
)
1671 ((eq el
'$infinity
) '$und
)
1672 ((zerop2 el
) (bylog expo bas
))
1673 (t (cond ((equal (getsignl el
) -
1)
1674 (cond ((eq bl
'$zeroa
) '$inf
)
1676 (cond ((even1 el
) '$inf
)
1677 ((eq (ask-integer el
'$integer
) '$yes
)
1678 (if (eq (ask-integer el
'$even
) '$yes
)
1680 '$minf
)))) ;Gotta be ODD.
1681 (t (setq bas
(behavior bas var val
))
1682 (cond ((equal bas
1) '$inf
)
1683 ((equal bas -
1) '$minf
)
1684 (t (throw 'limit t
))))))
1686 (member bl
'($zeroa $zerob
) :test
#'eq
))
1687 (cond ((even1 el
) '$zeroa
)
1688 ((and (eq bl
'$zerob
)
1690 (evenp (caddr el
))) 0)
1692 ((and (equal (getsignl el
) 1)
1693 (eq bl
'$zeroa
)) bl
)
1694 ((equal (getsignl el
) 0)
1698 (cond ((zerop2 el
) (bylog expo bas
))
1700 ((eq el
'$inf
) '$infinity
)
1701 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1702 ((equal (setq el
(getsignl el
)) 1) '$infinity
)
1705 (t (throw 'limit t
))))
1707 (cond ((eq el
'$inf
) '$inf
)
1708 ((equal el
'$minf
) 0)
1709 ((zerop2 el
) (bylog expo bas
))
1710 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1711 (t (cond ((eql 0 (getsignl el
)) 1)
1712 ((ratgreaterp 0 el
) '$zeroa
)
1713 ((ratgreaterp el
0) '$inf
)
1714 (t (throw 'limit t
))))))
1716 (cond ((zerop2 el
) (bylog expo bas
))
1717 ((eq el
'$inf
) '$und
)
1718 ((equal el
'$minf
) 0)
1719 ;;;Why not generalize this. We can ask about the number. -Jim 2/23/81
1720 ((mnump el
) (cond ((mnegp el
)
1723 (if (eq (ask-integer el
'$integer
) '$yes
)
1724 (if (eq (ask-integer el
'$even
) '$yes
)
1728 (t (cond ((even1 el
) '$inf
)
1729 ((eq (ask-integer el
'$integer
) '$yes
)
1730 (if (eq (ask-integer el
'$even
) '$yes
)
1734 (loginprod?
(throw 'lip?
'lip
!))
1736 ((equal (simplify (ratdisrep (ridofab bl
))) 1)
1737 (if (infinityp el
) (bylog expo bas
) 1))
1738 ((and (equal (ridofab bl
) -
1)
1739 (infinityp el
)) '$ind
) ;LB
1740 ((eq bl
'$ind
) (cond ((or (zerop2 el
) (infinityp el
)) '$und
)
1741 ((not (equal (getsignl el
) -
1)) '$ind
)
1743 ((eq el
'$inf
) (cond ((abeq1 bl
)
1744 (if (equal (getsignl bl
) 1) 1 '$ind
))
1746 (if (equal (getsignl bl
) 1) '$zeroa
0))
1747 ((equal (getsignl (m1- bl
)) 1) '$inf
)
1748 ((equal (getsignl (m1- `((mabs) ,bl
))) 1) '$infinity
)
1749 (t (throw 'limit t
))))
1750 ((eq el
'$minf
) (cond ((abeq1 bl
)
1751 (if (equal (getsignl bl
) 1) 1 '$ind
))
1753 (if (equal (getsignl bl
) 1) '$zeroa
0))
1754 ((ratgreaterp 0 bl
) '$infinity
)
1757 (if (equal val
'$infinity
)
1758 '$und
;Not enough info to do anything.
1759 (destructuring-bind (real-el . imag-el
)
1761 (setq real-el
(limit real-el var origval nil
))
1762 (cond ((eq real-el
'$minf
)
1764 ((and (eq real-el
'$inf
)
1765 (not (equal (ridofab (limit imag-el var origval nil
)) 0)))
1767 ((eq real-el
'$infinity
)
1768 (throw 'limit t
)) ;; don't really know real component
1772 ((eq el
'$ind
) '$ind
)
1777 (cond ((numberp x
) (and (integerp x
) (evenp x
)))
1778 ((and (mnump x
) (evenp (cadr x
))))))
1780 ;; is absolute value less than one?
1784 (and (ratgreaterp 1. bl
) (ratgreaterp bl -
1.
)))
1785 (t (equal (getsignl (m1- `((mabs) ,bl
))) -
1.
))))
1787 ;; is absolute value equal to one?
1791 (or (equal 1. bl
) (equal bl -
1.
)))
1792 (t (equal (getsignl (m1- `((mabs) ,bl
))) 0))))
1794 (defun simplimit (exp var val
&aux op
)
1797 ((or (atom exp
) (mnump exp
)) exp
)
1798 ((and (not (infinityp val
))
1799 (not (amongl '(%sin %cos %atanh %cosh %sinh %tanh mfactorial %log
)
1801 (not (inf-typep exp
))
1802 (simplimsubst val exp
)))
1803 ((eq (caar exp
) '%limit
) (throw 'limit t
))
1804 ((mplusp exp
) (simplimplus exp
))
1805 ((mtimesp exp
) (simplimtimes (cdr exp
)))
1806 ((mexptp exp
) (simplimexpt (cadr exp
) (caddr exp
)
1807 (limit (cadr exp
) var val
'think
)
1808 (limit (caddr exp
) var val
'think
)))
1809 ((mlogp exp
) (simplimln exp var val
))
1810 ((member (caar exp
) '(%sin %cos
) :test
#'eq
)
1811 (simplimsc exp
(caar exp
) (limit (cadr exp
) var val
'think
)))
1812 ((eq (caar exp
) '%tan
) (simplim%tan
(cadr exp
)))
1813 ((eq (caar exp
) '%atan
) (simplim%atan
(limit (cadr exp
) var val
'think
)))
1814 ((eq (caar exp
) '$atan2
) (simplim%atan2 exp
))
1815 ((member (caar exp
) '(%sinh %cosh
) :test
#'eq
)
1816 (simplimsch (caar exp
) (limit (cadr exp
) var val
'think
)))
1817 ((eq (caar exp
) 'mfactorial
)
1818 (simplimfact exp var val
))
1819 ((member (caar exp
) '(%erf %tanh
) :test
#'eq
)
1820 (simplim%erf-%tanh
(caar exp
) (cadr exp
)))
1821 ((member (caar exp
) '(%acos %asin
) :test
#'eq
)
1822 (simplim%asin-%acos
(caar exp
) (limit (cadr exp
) var val
'think
)))
1823 ((eq (caar exp
) '%atanh
)
1824 (simplim%atanh
(limit (cadr exp
) var val
'think
) val
))
1825 ((eq (caar exp
) '%acosh
)
1826 (simplim%acosh
(limit (cadr exp
) var val
'think
)))
1827 ((eq (caar exp
) '%asinh
)
1828 (simplim%asinh
(limit (cadr exp
) var val
'think
)))
1829 ((eq (caar exp
) '%inverse_jacobi_ns
)
1830 (simplim%inverse_jacobi_ns
(limit (cadr exp
) var val
'think
) (third exp
)))
1831 ((eq (caar exp
) '%inverse_jacobi_nc
)
1832 (simplim%inverse_jacobi_nc
(limit (cadr exp
) var val
'think
) (third exp
)))
1833 ((eq (caar exp
) '%inverse_jacobi_sc
)
1834 (simplim%inverse_jacobi_sc
(limit (cadr exp
) var val
'think
) (third exp
)))
1835 ((eq (caar exp
) '%inverse_jacobi_cs
)
1836 (simplim%inverse_jacobi_cs
(limit (cadr exp
) var val
'think
) (third exp
)))
1837 ((eq (caar exp
) '%inverse_jacobi_dc
)
1838 (simplim%inverse_jacobi_dc
(limit (cadr exp
) var val
'think
) (third exp
)))
1839 ((eq (caar exp
) '%inverse_jacobi_ds
)
1840 (simplim%inverse_jacobi_ds
(limit (cadr exp
) var val
'think
) (third exp
)))
1841 ((and (eq (caar exp
) 'mqapply
)
1842 (eq (subfunname exp
) '$li
))
1843 (simplim$li
(subfunsubs exp
) (subfunargs exp
) val
))
1844 ((and (eq (caar exp
) 'mqapply
)
1845 (eq (subfunname exp
) '$psi
))
1846 (simplim$psi
(subfunsubs exp
) (subfunargs exp
) val
))
1847 ((and (eq (caar exp
) var
)
1848 (member 'array
(car exp
) :test
#'eq
)
1849 (every #'(lambda (sub-exp)
1852 exp
) ;LIMIT(B[I],B,INF); -> B[I]
1853 ((setq op
(safe-get (mop exp
) 'simplim%function
))
1854 ;; Lookup a simplim%function from the property list
1855 (funcall op exp var val
))
1857 (simplify (cons (operator-with-array-flag exp
)
1858 (mapcar #'(lambda (a)
1859 (limit a var val
'think
))
1861 (throw 'limit t
)))))
1864 (setq e
(resimplify (subst (m// 1 var
) var e
)))
1865 (let ((new-val (cond ((eq val
'$zeroa
) '$inf
)
1866 ((eq val
'$zerob
) '$minf
))))
1867 (if new-val
(let ((preserve-direction t
))
1868 (toplevel-$limit e var new-val
)) (throw 'limit t
))))
1870 (defun simplimtimes (exp)
1871 ;; The following test
1872 ;; handles (-1)^x * 2^x => (-2)^x => $infinity
1873 ;; wants to avoid (-1)^x * 2^x => $ind * $inf => $und
1875 (and (expfactorp (cons '(mtimes) exp
) 1)
1876 (expfactor (cons '(mtimes) exp
) 1 var
))))
1877 (when try
(return-from simplimtimes try
)))
1879 (let ((prod 1) (num 1) (denom 1)
1880 (zf nil
) (ind-flag nil
) (inf-type nil
)
1881 (constant-zero nil
) (constant-infty nil
))
1883 (let* ((loginprod?
(involve term
'(%log
)))
1884 (y (catch 'lip?
(limit term var val
'think
))))
1886 ;; limit failed due to log in product
1888 (return-from simplimtimes
(liminv (cons '(mtimes simp
) exp
))))
1890 ;; If the limit is infinitesimal or zero
1892 (setf num
(m* num term
)
1893 constant-zero
(or constant-zero
(not (among var term
))))
1896 (unless zf
(setf zf
1)))
1898 (setf zf
(* -
1 (or zf
1))))))
1900 ;; If the limit is not some form of infinity or
1901 ;; undefined/indeterminate.
1902 ((not (member y
'($inf $minf $infinity $ind $und
) :test
#'eq
))
1903 (setq prod
(m* prod y
)))
1905 ((eq y
'$und
) (return-from simplimtimes
'$und
))
1906 ((eq y
'$ind
) (setq ind-flag t
))
1908 ;; Some form of infinity
1910 (setf denom
(m* denom term
)
1911 constant-infty
(or constant-infty
(not (among var term
))))
1912 (unless (eq inf-type
'$infinity
)
1914 ((eq y
'$infinity
) (setq inf-type
'$infinity
))
1915 ((null inf-type
) (setf inf-type y
))
1916 ;; minf * minf or inf * inf
1917 ((eq y inf-type
) (setf inf-type
'$inf
))
1919 (t (setf inf-type
'$minf
))))))))
1922 ;; If there are zeros and infinities among the terms that are free of
1923 ;; VAR, then we have an expression like "inf * zeroa * f(x)" or
1924 ;; similar. That gives an undefined result. Note that we don't
1925 ;; necessarily have something undefined if only the zeros have a term
1926 ;; free of VAR. For example "zeroa * exp(-1/x) * 1/x" as x -> 0. And
1927 ;; similarly for the infinities.
1928 ((and constant-zero constant-infty
) '$und
)
1930 ;; If num=denom=1, we didn't find any explicit infinities or zeros, so we
1931 ;; either return the simplified product or ind
1932 ((and (eql num
1) (eql denom
1))
1933 (if ind-flag
'$ind prod
))
1934 ;; If denom=1 (and so num != 1), we have some form of zero
1938 (let ((sign (getsignl prod
)))
1939 (if (or (not sign
) (eq sign
'complex
))
1945 ;; If num=1 (and so denom != 1), we have some form of infinity
1947 (let ((sign ($csign prod
)))
1950 ((eq sign
'$pos
) inf-type
)
1951 ((eq sign
'$neg
) (case inf-type
1955 ((member sign
'($complex $imaginary
)) '$infinity
)
1956 ; sign is '$zero, $pnz, $pz, etc
1957 (t (throw 'limit t
)))))
1958 ;; Both zeros and infinities
1960 ;; All bets off if there are some infinities or some zeros, but it
1961 ;; needn't be undefined (see above)
1962 (when (or constant-zero constant-infty
) (throw 'limit t
))
1964 (let ((ans (limit2 num
(m^ denom -
1) var val
)))
1966 (simplimtimes (list prod ans
))
1967 (throw 'limit t
)))))))
1969 ;;;PUT CODE HERE TO ELIMINATE FAKE SINGULARITIES??
1971 (defun simplimplus (exp)
1972 (cond ((memalike exp simplimplus-problems
)
1975 (progn (push exp simplimplus-problems
)
1976 (let ((ans (catch 'limit
(simplimplus1 exp
))))
1977 (cond ((or (eq ans
())
1979 (among '%limit ans
))
1980 (let ((new-exp (sratsimp exp
)))
1981 (cond ((not (alike1 exp new-exp
))
1983 (limit new-exp var val
'think
))))
1984 (cond ((or (eq ans
())
1989 (pop simplimplus-problems
)))))
1991 (defun simplimplus1 (exp)
1992 (prog (sum y infl infinityl minfl indl
)
1994 (do ((exp (cdr exp
) (cdr exp
)) (f))
1995 ((or y
(null exp
)) nil
)
1996 (setq f
(limit (car exp
) var val
'think
))
1999 ((eq f
'$und
) (setq y t
))
2000 ((not (member f
'($inf $minf $infinity $ind
) :test
#'eq
))
2001 (setq sum
(m+ sum f
)))
2002 ((eq f
'$ind
) (push (car exp
) indl
))
2003 (infinityl (throw 'limit t
))
2004 ;;;Don't know what to do with an '$infinity and an $inf or $minf
2005 ((eq f
'$inf
) (push (car exp
) infl
))
2006 ((eq f
'$minf
) (push (car exp
) minfl
))
2008 (cond ((or infl minfl
)
2010 (t (push (car exp
) infinityl
))))))
2011 (cond ((not (or infl minfl indl infinityl
))
2012 (return (cond ((atom sum
) sum
)
2013 ((or (not (free sum
'$zeroa
))
2014 (not (free sum
'$zerob
)))
2017 (t (cond ((null infinityl
)
2018 (cond (infl (cond ((null minfl
) (return '$inf
))
2020 (minfl (return '$minf
))
2021 ((> (length indl
) 1)
2022 ;; At this point we have a sum of '$ind. We factor
2023 ;; the sum and try again. This way we get the limit
2024 ;; of expressions like (a-b)*ind, where (a-b)--> 0.
2025 (cond ((not (alike1 (setq y
($factorsum exp
)) exp
))
2026 (return (limit y var val
'think
)))
2029 (t (return '$ind
))))
2030 (t (setq infl
(append infl infinityl
))))))
2032 oon
(setq y
(m+l
(append minfl infl
)))
2033 (cond ((alike1 exp
(setq y
(sratsimp (log-reduce (hyperex y
)))))
2034 (cond ((not (infinityp val
))
2035 (setq infl
(cnv infl val
)) ;THIS IS HORRIBLE!!!!
2036 (setq minfl
(cnv minfl val
))))
2038 (cond ((every #'(lambda (j) (radicalp j var
))
2039 (append infl minfl
))
2040 (setq y
(rheur infl minfl
)))
2041 (t (setq y
(sheur infl minfl
))))))
2042 (t (setq y
(limit y var val
'think
))))
2043 (cond ((or (eq y
())
2044 (eq y t
)) (return ()))
2045 ((infinityp y
) (return y
))
2046 (t (return (m+ sum y
))))))
2048 ;; Limit n/d, using heuristics on the order of growth.
2051 (cond ((and (free n var
)
2054 (t (setq n
(cpa n d nil
))
2056 (cond ((oscip orig-n
) '$und
)
2058 ((equal n -
1) '$zeroa
)
2059 ((equal n
0) (m// orig-n d
)))))))
2062 ;;;L1 is a list of INF's and L2 is a list of MINF's. Added together
2063 ;;;it is indeterminate.
2064 (defun sheur (l1 l2
)
2065 (let ((term (sheur1 l1 l2
)))
2066 (cond ((equal term
'$inf
) '$inf
)
2067 ((equal term
'$minf
) '$minf
)
2068 (t (let ((new-num (m+l
(mapcar #'(lambda (num-term)
2069 (m// num-term
(car l1
)))
2071 (cond ((limit2 new-num
(m// 1 (car l1
)) var val
))))))))
2074 (cond ((atom exp
) nil
)
2075 (t (setq exp
(nformat exp
))
2076 (cond ((and (eq (caar exp
) 'mquotient
)
2077 (among var
(caddr exp
)))
2080 (defun zerop2 (z) (=0 (ridofab z
)))
2082 (defun raise (a) (m+ a
'$zeroa
))
2084 (defun lower (a) (m+ a
'$zerob
))
2086 (defun sincoshk (exp1 l sc
)
2087 (cond ((equal l
1) (lower l
))
2088 ((equal l -
1) (raise l
))
2090 ((member val
'($zeroa $zerob
) :test
#'eq
) (spangside exp1 l
))
2093 (defun spangside (e l
)
2094 (setq e
(behavior e var val
))
2095 (cond ((equal e
1) (raise l
))
2096 ((equal e -
1) (lower l
))
2099 ;; get rid of zeroa and zerob
2101 (if (among '$zeroa e
) (setq e
(maxima-substitute 0 '$zeroa e
)))
2102 (if (among '$zerob e
) (setq e
(maxima-substitute 0 '$zerob e
)))
2106 ;; returns true if exp is a polynomial raised to a numeric power
2107 (defun simplerd (exp)
2109 (mnump (caddr exp
)) ;; exponent must be a number - no variables
2110 (polyp (cadr exp
))))
2112 (defun branch1 (exp val
)
2113 (cond ((polyp exp
) nil
)
2114 ((simplerd exp
) (zerop2 (subin val
(cadr exp
))))
2116 (loop for v on
(cdr exp
)
2117 when
(branch1 (car v
) val
)
2120 (defun branch (exp val
)
2121 (cond ((polyp exp
) nil
)
2122 ((or (simplerd exp
) (mtimesp exp
))
2125 (every #'(lambda (j) (branch j val
)) (the list
(cdr exp
))))))
2127 (defun ser0 (e n d val
)
2128 (cond ((and (branch n val
) (branch d val
))
2132 ;;;NN* gets set by POFX, called by SER1, to get a list of exponents.
2133 (setq nn
* (ratmin nn
*))
2134 (setq n
(sratsimp (m^ n
(m^ var nn
*))))
2135 (setq d
(sratsimp (m^ d
(m^ var nn
*))))
2136 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2139 (t (try-lhospital-quit n d nil
))))
2141 (defun rheur (l1 l2
)
2143 (setq m1
(mapcar (function asymredu
) l1
))
2144 (setq m2
(mapcar (function asymredu
) l2
))
2145 (setq ans
(m+l
(append m1 m2
)))
2146 (cond ((rptrouble (m+l
(append l1 l2
)))
2147 (return (limit (simplify (rdsget (m+l
(append l1 l2
))))
2151 ((mplusp ans
) (return (sheur m1 m2
)))
2152 (t (return (limit ans var val t
))))))
2154 (defun rptrouble (rp)
2155 (not (equal (rddeg rp nil
) (rddeg (asymredu rp
) nil
))))
2157 (defun radicalp (exp var
)
2158 (cond ((polyinx exp var
()))
2159 ((mexptp exp
) (cond ((equal (caddr exp
) -
1.
)
2160 (radicalp (cadr exp
) var
))
2162 ((member (caar exp
) '(mplus mtimes
) :test
#'eq
)
2163 (every #'(lambda (j) (radicalp j var
))
2166 (defun involve (e nn
*)
2167 (declare (special var
))
2168 (cond ((atom e
) nil
)
2170 ((and (member (caar e
) nn
* :test
#'eq
) (among var
(cdr e
))) (cadr e
))
2171 (t (some #'(lambda (j) (involve j nn
*)) (cdr e
)))))
2173 (defun notinvolve (exp nn
*)
2174 (cond ((atom exp
) t
)
2176 ((member (caar exp
) nn
* :test
#'eq
) (not (among var
(cdr exp
))))
2177 ((every #'(lambda (j) (notinvolve j nn
*))
2180 (defun sheur1 (l1 l2
)
2182 (setq l1
(m+l
(maxi l1
)))
2183 (setq l2
(m+l
(maxi l2
)))
2184 (setq ans
(cpa l1 l2 t
))
2185 (return (cond ((=0 ans
) (m+ l1 l2
))
2186 ((equal ans
1.
) '$inf
)
2189 (defun zero-lim (cpa-list)
2190 (do ((l cpa-list
(cdr l
)))
2192 (and (eq (caar l
) 'gen
)
2193 (zerop2 (limit (cadar l
) var val
'think
))
2196 ;; Compare order of growth for R1 and R2. The result is 0, -1, +1
2197 ;; depending on the relative order of growth. 0 is returned if R1 and
2198 ;; R2 have the same growth; -1 if R1 grows much more slowly than R2;
2199 ;; +1 if R1 grows much more quickly than R2.
2200 (defun cpa (r1 r2 flag
)
2203 (cond ((alike1 t1 t2
) 0.
)
2205 (cond ((free t2 var
) 0.
)
2206 (t (let ((lim-ans (limit1 t2 var val
)))
2207 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2210 (let ((lim-ans (limit1 t1 var val
)))
2211 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2214 ;; Make T1 and T2 be a list of terms that are multiplied
2216 (cond ((mtimesp t1
) (setq t1
(cdr t1
)))
2217 (t (setq t1
(list t1
))))
2218 (cond ((mtimesp t2
) (setq t2
(cdr t2
)))
2219 (t (setq t2
(list t2
))))
2220 ;; Find the strengths of each term of T1 and T2
2221 (setq t1
(mapcar (function istrength
) t1
))
2222 (setq t2
(mapcar (function istrength
) t2
))
2223 ;; Compute the max of the strengths of the terms.
2224 (let ((ans (ismax t1
))
2226 (cond ((or (null ans
) (null d
)
2227 (eq (car ans
) 'gen
) (eq (car d
) 'gen
)) 0.
))
2228 (if (eq (car ans
) 'var
) (setq ans
(add-up-deg t1
)))
2229 (if (eq (car d
) 'var
) (setq d
(add-up-deg t2
)))
2230 ;; Can't just just compare dominating terms if there are
2231 ;; indeterm-inates present; e.g. X-X^2*LOG(1+1/X). So
2233 (cond ((or (zero-lim t1
)
2235 (cpa-indeterm ans d t1 t2 flag
))
2236 ((isgreaterp ans d
) 1.
)
2237 ((isgreaterp d ans
) -
1.
)
2240 (defun cpa-indeterm (ans d t1 t2 flag
)
2241 (cond ((not (eq (car ans
) 'var
))
2242 (setq ans
(gather ans t1
) d
(gather d t2
))))
2243 (let ((*indicator
(and (eq (car ans
) 'exp
)
2246 (setq test
(cpa1 ans d
))
2247 (cond ((and (zerop1 test
)
2248 (or (equal ($radcan
(m// (cadr ans
) (cadr d
))) 1.
)
2249 (and (polyp (cadr ans
))
2251 (equal (limit (m// (cadr ans
) (cadr d
)) var val
'think
)
2253 (let ((new-term1 (m// t1
(cadr ans
)))
2254 (new-term2 (m// t2
(cadr d
))))
2255 (cpa new-term1 new-term2 flag
)))
2258 (defun add-up-deg (strengthl)
2259 (do ((stl strengthl
(cdr stl
))
2262 ((null stl
) (list 'var
(m*l poxl
) (m+l degl
)))
2263 (cond ((eq (caar stl
) 'var
)
2264 (push (cadar stl
) poxl
)
2265 (push (caddar stl
) degl
)))))
2269 (cond ((eq (car p1
) 'gen
) (return 0.
)))
2270 (setq flag
(car p1
))
2275 (setq s1
(istrength p1
))
2276 (setq s2
(istrength p2
))
2279 ((isgreaterp s1 s2
) 1.
)
2280 ((isgreaterp s2 s1
) -
1.
)
2282 (setq *indicator nil
)
2284 ((and (poly? p1 var
) (poly? p2 var
))
2285 (setq p1
(m- p1 p2
))
2286 (cond ((zerop1 p1
) 0.
)
2287 (t (getsignl (hot-coef p1
)))))
2291 (list (m*t -
1 p2
))))
2292 (cond ((zerop2 s1
) 0.
)
2293 ((ratgreaterp s1
0.
) 1.
)
2297 (setq p1
(caddr p1
))
2298 (setq p2
(caddr p2
))
2299 (cond ((and (poly? p1 var
) (poly? p2 var
))
2300 (setq p1
(m- p1 p2
))
2301 (return (cond ((or (zerop1 p1
)
2302 (not (among var p1
)))
2304 (t (getsignl (hot-coef p1
))))))
2305 ((and (radicalp p1 var
) (radicalp p2 var
))
2308 (list (m*t -
1 p2
))))
2309 (return (cond ((eq s1
'$inf
) 1.
)
2310 ((eq s1
'$minf
) -
1.
)
2312 (cond ((ratgreaterp s1
0.
) 1.
)
2313 ((ratgreaterp 0. s1
) -
1.
)
2316 (t (return (cpa p1 p2 t
)))))
2318 (setq p1
(try-lhospital (asymredu p1
) (asymredu p2
) nil
))
2319 (return (cond ((zerop2 p1
) -
1.
)
2320 ((real-infinityp p1
) 1.
)
2323 ;;;EXPRESSIONS TO ISGREATERP ARE OF THE FOLLOWING FORMS
2324 ;;; ("VAR" POLY DEG)
2326 ;;; ("LOG" LOG(EXP))
2327 ;;; ("FACT" <A FACTORIAL EXPRESSION>)
2328 ;;; ("GEN" <ANY OTHER TYPE OF EXPRESSION>)
2330 (defun isgreaterp (a b
)
2333 (cond ((or (eq ta
'gen
)
2335 ((and (eq ta tb
) (eq ta
'var
))
2336 (ratgreaterp (caddr a
) (caddr b
)))
2337 ((and (eq ta tb
) (eq ta
'exp
))
2338 ;; Both are exponential order of infinity. Check the
2339 ;; exponents to determine which exponent is bigger.
2340 (eq (limit (m- `((%log
) ,(second a
)) `((%log
) ,(second b
)))
2343 ((member ta
(cdr (member tb
'(num log var exp fact gen
) :test
#'eq
)) :test
#'eq
)))))
2346 ;; Preprocess the list of products. Separate the terms that
2347 ;; exponentials and those that don't. Actually multiply the
2348 ;; exponential terms together to form a single term. Pass this and
2349 ;; the rest to ismax-core to find the max.
2350 (let (exp-terms non-exp-terms
)
2352 (if (eq 'exp
(car term
))
2353 (push term exp-terms
)
2354 (push term non-exp-terms
)))
2355 ;; Multiply the exp-terms together
2358 ;;(format t "exp-terms = ~A~%" exp-terms)
2359 (dolist (term exp-terms
)
2360 (setf product
(simplify (mul product
(second term
)))))
2361 ;;(format t "product = ~A~%" product)
2362 (setf product
`(exp ,($logcontract product
)))
2363 ;;(format t "product = ~A~%" product)
2364 (ismax-core (cons product non-exp-terms
)))
2367 (defun ismax-core (l)
2370 ((= (length l
) 1) (car l
)) ;If there is only 1 thing give it back.
2371 ((every #'(lambda (x)
2372 (not (eq (car x
) 'gen
))) l
)
2374 (do ((l1 (cdr l
) (cdr l1
))
2378 (cond ((isgreaterp temp-ans
(car l1
))
2379 (setq ans temp-ans
))
2380 ((isgreaterp (car l1
) temp-ans
)
2381 (setq temp-ans
(car l1
))
2382 (setq ans temp-ans
))
2383 (t (setq ans
())))))
2386 ;RETURNS LIST OF HIGH TERMS
2388 (cond ((atom all
) nil
)
2389 (t (do ((l (cdr all
) (cdr l
))
2391 (total 1) ; running total constant factor of hi-term
2392 (hi-terms (ncons (car all
)))
2394 ((null l
) (if (zerop2 total
) ; if high-order terms cancel each other
2395 all
; keep everything
2396 hi-terms
)) ; otherwise return list of high terms
2397 (setq compare
(limit (m// (car l
) hi-term
) var val
'think
))
2399 ((or (infinityp compare
)
2400 (and (eq compare
'$und
)
2401 (zerop2 (limit (m// hi-term
(car l
)) var val
'think
))))
2402 (setq total
1) ; have found new high term
2403 (setq hi-terms
(ncons (setq hi-term
(car l
)))))
2404 ((zerop2 compare
) nil
)
2405 ;; COMPARE IS IND, FINITE-VALUED, or und in both directions
2406 (t ; add to list of high terms
2407 (setq total
(m+ total compare
))
2408 (setq hi-terms
(append hi-terms
(ncons (car l
))))))))))
2412 (cond ((atom l
) (return nil
)))
2413 l1
(setq ans
(car l
))
2415 (cond ((null l
) (return ans
))
2416 ((ratgreaterp ans
(car l
)) (go l2
))
2421 (cond ((atom l
) (return nil
)))
2422 l1
(setq ans
(car l
))
2424 (cond ((null l
) (return ans
))
2425 ((ratgreaterp (car l
) ans
) (go l2
))
2433 ((or (mnump e
) (not (among var e
))) nil
)
2434 ((and (mexptp e
) (eq (cadr e
) var
))
2435 (push (caddr e
) nn
*))
2437 (t (mapc #'pofx
(cdr e
)))))
2440 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2441 (t (setq e
(subin (m+ var val
) e
))))
2443 (cond ((pofx e
) e
)))
2445 (defun gather (ind l
)
2447 (setq ind
(car ind
))
2448 loop
(cond ((null l
)
2449 (return (list ind
(m*l ans
))))
2450 ((equal (caar l
) ind
)
2451 (push (cadar l
) ans
)))
2455 ; returns rough class-of-growth of term
2456 (defun istrength (term)
2457 (cond ((mnump term
) (list 'num term
))
2458 ((atom term
) (cond ((eq term var
)
2460 (t (list 'num term
))))
2461 ((not (among var term
)) (list 'num term
))
2463 (let ((temp (ismax (mapcar #'istrength
(cdr term
)))))
2464 (cond ((not (null temp
)) temp
)
2467 (let ((temp (mapcar #'istrength
(cdr term
)))
2469 (setq temp1
(ismax temp
))
2470 (cond ((null temp1
) `(gen ,term
))
2471 ((eq (car temp1
) 'log
) `(log ,temp
))
2472 ((eq (car temp1
) 'var
) (add-up-deg temp
))
2475 (real-infinityp (limit term var val t
)))
2476 (let ((logterm (logred term
)))
2477 (cond ((and (among var
(caddr term
))
2478 (member (car (istrength logterm
))
2479 '(var exp fact
) :test
#'eq
)
2480 (real-infinityp (limit logterm var val t
)))
2481 (list 'exp
(m^
'$%e logterm
)))
2482 ((not (among var
(caddr term
)))
2483 (let ((temp (istrength (cadr term
))))
2484 (cond ((not (alike1 temp term
))
2485 (rplaca (cdr temp
) term
)
2486 (and (eq (car temp
) 'var
)
2488 (m* (caddr temp
) (caddr term
))))
2492 ((and (eq (caar term
) '%log
)
2493 (real-infinityp (limit term var val t
)))
2494 (let ((stren (istrength (cadr term
))))
2495 (cond ((member (car stren
) '(log var
) :test
#'eq
)
2497 ((and (eq (car stren
) 'exp
)
2498 (eq (caar (second stren
)) 'mexpt
))
2499 (istrength (logred (second stren
))))
2501 ((eq (caar term
) 'mfactorial
)
2503 ((let ((temp (hyperex term
)))
2504 (and (not (alike1 term temp
))
2506 (t (list 'gen term
))))
2508 ;; log reduce - returns log of s1
2510 (or (and (eq (cadr s1
) '$%e
) (caddr s1
))
2511 (m* (caddr s1
) `((%log
) ,(cadr s1
)))))
2513 (defun asymredu (rd)
2514 (cond ((atom rd
) rd
)
2516 ((not (among var rd
)) rd
)
2517 ((polyinx rd var t
))
2519 (cond ((eq (cadr rd
) var
) rd
)
2521 (factor ($expand
(m^
(polyinx (cadr rd
) var t
)
2525 (t (simplify (cons (list (caar rd
))
2526 (mapcar #'asymredu
(cdr rd
)))))))
2529 (let ((dn** ()) (nn** ()))
2530 (cond ((atom rd
) rd
)
2532 ((not (among var rd
)) rd
)
2536 (cond ((eq (cadr rd
) var
) rd
)
2537 (t (setq dn
** (caddr rd
))
2538 (setq nn
** (factor (cadr rd
)))
2539 (cond ((mtimesp nn
**)
2540 (m*l
(mapcar #'(lambda (j) (m^ j dn
**))
2543 (t (simplify (cons (ncons (caar rd
))
2544 (mapcar #'rdfact
(cdr rd
))))))))
2546 (defun cnv (expl val
)
2547 (mapcar #'(lambda (e)
2548 (maxima-substitute (cond ((eq val
'$zerob
)
2549 (m* -
1 (m^ var -
1)))
2554 (t (m^
(m+ var
(m* -
1 val
)) -
1.
)))
2559 (defun pwtaylor (exp var l terms
)
2560 (prog (coef ans c mc
)
2561 (cond ((=0 terms
) (return nil
)) ((=0 l
) (setq mc t
)))
2564 loop
(setq c
(1+ c
))
2565 (cond ((or (> c
10.
) (equal c terms
))
2567 (t (setq exp
(sdiff exp var
))))
2568 tag1
(setq coef
($radcan
(subin l exp
)))
2569 (cond ((=0 coef
) (setq terms
(1+ terms
)) (go loop
)))
2576 (m^
`((mfactorial) ,c
) -
1)
2577 (m^
(if mc var
(m+t
(m*t -
1 l
) var
)) c
)))))
2582 ((simplerd e
) (rdtay e
))
2583 (t (cons (list (caar e
))
2584 (mapcar #'rdsget
(cdr e
))))))
2587 (cond (limit-using-taylor ($ratdisrep
($taylor rd var val
1.
)))
2591 (prog (varlist p c e d $ratfac
)
2592 (setq varlist
(ncons var
))
2593 (setq p
(ratnumerator (cdr (ratrep* (cadr rd
)))))
2594 (cond ((< (length p
) 3.
) (return rd
)))
2597 (setq c
(m^ var
(m* d e
)))
2598 (setq d
($ratsimp
(varinvert (m* (pdis p
) (m^ var
(m- d
)))
2600 (setq d
(pwtaylor (m^ d e
) var
0.
3.
))
2601 (return (m* c
(varinvert d var
)))))
2603 (defun varinvert (e var
) (subin (m^t var -
1.
) e
))
2606 (prog ((varlist (list var
)))
2607 (return (let (($ratfac nil
))
2609 (pdegr (cadr (ratrep* p
)))))))
2611 (defun rat-no-ratfac (e)
2612 (let (($ratfac nil
))
2617 (defun rddeg (rd low
*)
2618 (cond ((or (mnump rd
)
2619 (not (among var rd
)))
2624 (m* (deg (cadr rd
)) (caddr rd
)))
2626 (addn (mapcar #'(lambda (j)
2630 (setq rd
(andmapcar #'(lambda (j) (rddeg j low
*))
2632 (cond (low* (ratmin rd
))
2636 (cond ((or (atom pf
) (not (eq (caadr (ratf var
)) (car pf
))))
2638 (low* (cadr (reverse pf
)))
2640 ;;There is some confusion here. We need to be aware of Branch cuts etc....
2641 ;;when doing this section of code. It is not very carefully done so there
2642 ;;are bugs still lurking. Another misfortune is that LIMIT or its inferiors
2643 ;;sometimes decides to change the limit VAL in midstream. This must be corrected
2644 ;;since LIMIT's interaction with the data base environment must be maintained.
2645 ;;I'm not sure that this code can ever be called with VAL other than $INF but
2646 ;;there is a hook in the first important cond clause to cathc them anyway.
2649 (let ((num-power (rddeg n nil
))
2650 (den-power (rddeg d nil
))
2651 (coef ()) (coef-sign ()) (power ()))
2652 (setq coef
(m// ($ratcoef
($expand n
) var num-power
)
2653 ($ratcoef
($expand d
) var den-power
)))
2654 (setq coef-sign
(getsignl coef
))
2655 (setq power
(m// num-power den-power
))
2656 (cond ((eq (ask-integer power
'$integer
) '$integer
)
2657 (cond ((eq (ask-integer power
'$even
) '$even
) '$even
)
2658 (t '$odd
)))) ;Can be extended from here.
2659 (cond ((or (eq val
'$minf
)
2662 (equal val
0)) ()) ;Can be extended to cover some these.
2663 ((ratgreaterp den-power num-power
)
2664 (cond ((equal coef-sign
1.
) '$zeroa
)
2665 ((equal coef-sign -
1) '$zerob
)
2666 ((equal coef-sign
0) 0)
2668 ((ratgreaterp num-power den-power
)
2669 (cond ((equal coef-sign
1.
) '$inf
)
2670 ((equal coef-sign -
1) '$minf
)
2671 ((equal coef-sign
0) nil
) ; should never be zero
2672 ((null coef-sign
) '$infinity
)))
2675 (defun radlim (e n d
)
2677 (cond ((eq val
'$infinity
) (throw 'limit nil
))
2679 (setq nl
(m* var -
1))
2680 (setq n
(subin nl n
))
2681 (setq d
(subin nl d
))
2682 (setq val
'$inf
))) ;This is the Culprit. Doesn't tell the DATABASE.
2683 (cond ((eq val
'$inf
)
2684 (setq nl
(asymredu n
))
2685 (setq dl
(asymredu d
))
2687 ((or (rptrouble n
) (rptrouble d
))
2688 (return (limit (m* (rdsget n
) (m^
(rdsget d
) -
1.
)) var val t
)))
2689 (t (return (asy nl dl
))))))
2690 (setq nl
(limit n var val t
))
2691 (setq dl
(limit d var val t
))
2692 (cond ((and (zerop2 nl
) (zerop2 dl
))
2693 (cond ((or (polyp n
) (polyp d
))
2694 (return (try-lhospital-quit n d t
)))
2695 (t (return (ser0 e n d val
)))))
2696 (t (return ($radcan
(ratrad (m// n d
) n d nl dl
)))))))
2698 (defun ratrad (e n d nl dl
)
2700 (cond ((equal nl
0) (return 0))
2703 (cond ((equal dl
0) (setq d1
'$infinity
)) ;No direction Info.
2706 ((equal (setq d
(behavior d var val
)) 1)
2708 ((equal d -
1) (setq d1
'$minf
))
2709 (t (throw 'limit nil
))))
2712 (cond ((equal (setq n
(behavior n var val
)) 1)
2714 ((equal n -
1) (setq n1
'$zerob
))
2716 (t (return ($radcan
(ridofab (subin val e
))))))
2717 (return (simplimtimes (list n1 d1
)))))
2719 ;;; Limit of the Logarithm function
2721 (defun simplimln (expr var val
)
2722 ;; We need to be careful with log because of the branch cut on the
2723 ;; negative real axis. So we look at the imagpart of the argument. If
2724 ;; it's not identically zero, we compute the limit of the real and
2725 ;; imaginary parts and combine them. Otherwise, we can use the
2726 ;; original method for real limits.
2727 (let ((arglim (limit (cadr expr
) var val
'think
)))
2728 (cond ((eq arglim
'$inf
) '$inf
)
2729 ((member arglim
'($minf $infinity
) :test
#'eq
)
2731 ((member arglim
'($ind $und
) :test
#'eq
) '$und
)
2732 ((eq arglim
'$zeroa
) '$minf
)
2733 ((eq arglim
'$zerob
) '$infinity
)
2734 ((equalp arglim
0) '$infinity
)
2736 (let ((dir (behavior (cadr expr
) var val
)))
2737 (cond ((equal dir
1) '$zeroa
)
2738 ((equal dir -
1) '$zerob
)
2740 ((equalp ($imagpart
(cadr expr
)) 0)
2741 ;; argument is real.
2742 (simplify `((%log
) ,arglim
)))
2743 (t ;; argument is complex.
2744 (destructuring-bind (rp . ip
)
2746 (if (eq (setq rp
(limit rp var val
'think
)) '$minf
)
2747 ;; Realpart is minf, do not return minf+%i*ip but infinity.
2749 ;; Return a complex limit value.
2750 (add rp
(mul '$%i
(limit ip var val
'think
)))))))))
2752 ;;; Limit of the Factorial function
2754 (defun simplimfact (expr var val
)
2755 (let* ((arglim (limit (cadr expr
) var val
'think
)) ; Limit of the argument.
2757 (cond ((eq arglim
'$inf
) '$inf
)
2758 ((member arglim
'($minf $infinity $und $ind
) :test
#'eq
) '$und
)
2759 ((and (or (maxima-integerp arglim
)
2760 (setq arg2
(integer-representation-p arglim
)))
2761 (eq ($sign arg2
) '$neg
))
2762 ;; A negative integer or float or bigfloat representation.
2763 (let ((dir (limit (add (cadr expr
) (mul -
1 arg2
)) var val
'think
))
2764 (even (mevenp arg2
)))
2765 (cond ((or (and even
2775 (t (throw 'limit nil
)))))
2777 ;; Call simplifier to get value at the limit of the argument.
2778 (simplify (list '(mfactorial) arglim
))))))
2780 (defun simplim%erf-%tanh
(fn arg
)
2781 (let ((arglim (limit arg var val
'think
))
2784 (cond ((eq arglim
'$inf
) 1)
2785 ((eq arglim
'$minf
) -
1)
2786 ((eq arglim
'$infinity
)
2787 (destructuring-bind (rpart . ipart
)
2789 (setq rlim
(limit rpart var origval
'think
))
2790 (cond ((eq fn
'%tanh
)
2791 (cond ((equal rlim
'$inf
) 1)
2792 ((equal rlim
'$minf
) -
1)))
2794 (setq ans
(limit (m* rpart
(m^t ipart -
1)) var origval
'think
))
2795 (setq ans
($asksign
(m+ `((mabs) ,ans
) -
1)))
2796 (cond ((or (eq ans
'$pos
) (eq ans
'$zero
))
2797 (cond ((eq rlim
'$inf
) 1)
2798 ((eq rlim
'$minf
) -
1)
2801 ((eq arglim
'$und
) '$und
)
2802 ((member arglim
'($zeroa $zerob $ind
) :test
#'eq
) arglim
)
2803 ;;;Ignore tanh(%pi/2*%I) and multiples of the argument.
2805 ;; erf (or tanh) of a known value is just erf(arglim).
2806 (simplify (list (ncons fn
) arglim
))))))
2808 (defun simplim%atan
(exp1)
2809 (cond ((zerop2 exp1
) exp1
)
2810 ((member exp1
'($und $infinity
) :test
#'eq
)
2812 ((eq exp1
'$inf
) half%pi
)
2815 (t `((%atan
) ,exp1
))))
2817 ;; Most instances of atan2 are simplified to expressions in atan
2818 ;; by simpatan2 before we get to this point. This routine handles
2819 ;; tricky cases such as limit(atan2((x^2-2), x^3-2*x), x, sqrt(2), minus).
2820 ;; Taylor and Gruntz cannot handle the discontinuity at atan(0, -1)
2821 (defun simplim%atan2
(exp)
2822 (let* ((exp1 (cadr exp
))
2824 (lim1 (limit (cadr exp
) var val
'think
))
2825 (lim2 (limit (caddr exp
) var val
'think
))
2826 (sign2 ($csign lim2
)))
2827 (cond ((and (zerop2 lim1
) ;; atan2( 0+, + )
2829 lim1
) ;; result is zeroa or zerob
2830 ((and (eq lim1
'$zeroa
)
2833 ((and (eq lim1
'$zerob
) ;; atan2( 0-, - )
2836 ((and (eq lim1
'$zeroa
) ;; atan2( 0+, 0 )
2838 (simplim%atan
(limit (m// exp1 exp2
) var val
'think
)))
2839 ((and (eq lim1
'$zerob
) ;; atan2( 0-, 0 )
2841 (m+ (porm (eq lim2
'$zeroa
) '$%pi
)
2842 (simplim%atan
(limit (m// exp1 exp2
) var val
'think
))))
2843 ((member lim1
'($und $infinity
) :test
#'eq
)
2845 ((eq lim1
'$inf
) half%pi
)
2848 (t (take '($atan2
) lim1 lim2
)))))
2850 (defun simplimsch (sch arg
)
2851 (cond ((real-infinityp arg
)
2852 (cond ((eq sch
'%sinh
) arg
) (t '$inf
)))
2853 ((eq arg
'$infinity
) '$infinity
)
2854 ((eq arg
'$ind
) '$ind
)
2855 ((eq arg
'$und
) '$und
)
2856 (t (let (($exponentialize t
))
2857 (resimplify (list (ncons sch
) (ridofab arg
)))))))
2859 ;; simple limit of sin and cos
2860 (defun simplimsc (exp fn arg
)
2861 (cond ((member arg
'($inf $minf $ind
) :test
#'eq
) '$ind
)
2862 ((member arg
'($und $infinity
) :test
#'eq
)
2864 ((member arg
'($zeroa $zerob
) :test
#'eq
)
2865 (cond ((eq fn
'%sin
) arg
)
2866 (t (m+ 1 '$zerob
))))
2868 (simplify (list (ncons fn
) (ridofab arg
)))
2871 (defun simplim%tan
(arg)
2872 (let ((arg1 (ridofab (limit arg var val
'think
))))
2874 ((member arg1
'($inf $minf $infinity $ind $und
) :test
#'eq
) '$und
)
2876 (let ((c (trigred (pip arg1
))))
2877 (cond ((not (equal ($imagpart arg1
) 0)) '$infinity
)
2878 ((and (eq (caar c
) 'rat
)
2881 (setq arg1
(behavior arg var val
))
2882 (cond ((= arg1
1) '$inf
)
2883 ((= arg1 -
1) '$minf
)
2885 ((and (eq (caar c
) 'rat
)
2888 (setq arg1
(behavior arg var val
))
2889 (cond ((= arg1
1) '$minf
)
2892 (t (throw 'limit
())))))
2894 (setq arg1
(behavior arg var val
))
2895 (cond ((equal arg1
1) '$zeroa
)
2896 ((equal arg1 -
1) '$zerob
)
2898 (t (simp-%tan
(list '(%tan
) arg1
) 1 nil
)))))
2900 (defun simplim%asinh
(arg)
2901 (cond ((member arg
'($inf $minf $zeroa $zerob $ind $und
) :test
#'eq
)
2903 ((eq arg
'$infinity
) '$und
)
2904 (t (simplify (list '(%asinh
) (ridofab arg
))))))
2906 (defun simplim%acosh
(arg)
2907 (cond ((equal (ridofab arg
) 1) '$zeroa
)
2908 ((eq arg
'$inf
) arg
)
2909 ((eq arg
'$minf
) '$infinity
)
2910 ((member arg
'($und $ind $infinity
) :test
#'eq
) '$und
)
2911 (t (simplify (list '(%acosh
) (ridofab arg
))))))
2913 (defun simplim%atanh
(arg dir
)
2914 ;; Compute limit(atanh(x),x,arg). If ARG is +/-1, we need to take
2915 ;; into account which direction we're approaching ARG.
2916 (cond ((zerop2 arg
) arg
)
2917 ((member arg
'($ind $und $infinity $minf $inf
) :test
#'eq
)
2919 ((equal (setq arg
(ridofab arg
)) 1.
)
2920 ;; The limit at 1 should be complex infinity because atanh(x)
2921 ;; is complex for x > 1, but inf if we're approaching 1 from
2923 (if (eq dir
'$zerob
)
2927 ;; Same as above, except for the limit is at -1.
2928 (if (eq dir
'$zeroa
)
2931 (t (simplify (list '(%atanh
) arg
)))))
2933 (defun simplim%asin-%acos
(fn arg
)
2934 (cond ((member arg
'($und $ind $inf $minf $infinity
) :test
#'eq
)
2936 ((and (eq fn
'%asin
)
2937 (member arg
'($zeroa $zerob
) :test
#'eq
))
2939 (t (simplify (list (ncons fn
) (ridofab arg
))))))
2941 (defun simplim$li
(order arg val
)
2942 (if (and (not (equal (length order
) 1))
2943 (not (equal (length arg
) 1)))
2945 (setq order
(car order
)
2947 (if (not (equal order
2))
2949 (destructuring-bind (rpart . ipart
)
2951 (cond ((not (equal ipart
0))
2954 (setq rpart
(limit rpart var val
'think
))
2955 (cond ((eq rpart
'$zeroa
) '$zeroa
)
2956 ((eq rpart
'$zerob
) '$zerob
)
2957 ((eq rpart
'$minf
) '$minf
)
2958 ((eq rpart
'$inf
) '$infinity
)
2959 (t (simplify (subfunmake '$li
(list order
)
2960 (list rpart
))))))))))
2962 (defun simplim$psi
(order arg val
)
2963 (if (and (not (equal (length order
) 1))
2964 (not (equal (length arg
) 1)))
2966 (setq order
(car order
)
2968 (cond ((equal order
0)
2969 (destructuring-bind (rpart . ipart
)
2971 (cond ((not (equal ipart
0)) (throw 'limit
()))
2972 (t (setq rpart
(limit rpart var val
'think
))
2973 (cond ((eq rpart
'$zeroa
) '$minf
)
2974 ((eq rpart
'$zerob
) '$inf
)
2975 ((eq rpart
'$inf
) '$inf
)
2976 ((eq rpart
'$minf
) '$und
)
2977 ((equal (getsignl rpart
) -
1) (throw 'limit
()))
2978 (t (simplify (subfunmake '$psi
(list order
)
2979 (list rpart
)))))))))
2980 ((and (integerp order
) (> order
0)
2981 (equal (limit arg var val
'think
) '$inf
))
2982 (cond ((mevenp order
) '$zerob
)
2983 ((moddp order
) '$zeroa
)
2984 (t (throw 'limit
()))))
2985 (t (throw 'limit
()))))
2987 (defun simplim%inverse_jacobi_ns
(arg m
)
2988 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
2990 `((%inverse_jacobi_ns
) ,arg
,m
)))
2992 (defun simplim%inverse_jacobi_nc
(arg m
)
2993 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
2994 `((%elliptic_kc
) ,m
)
2995 `((%inverse_jacobi_nc
) ,arg
,m
)))
2997 (defun simplim%inverse_jacobi_sc
(arg m
)
2998 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
2999 `((%elliptic_kc
) ,m
)
3000 `((%inverse_jacobi_sc
) ,arg
,m
)))
3002 (defun simplim%inverse_jacobi_dc
(arg m
)
3003 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3004 `((%elliptic_kc
) ,m
)
3005 `((%inverse_jacobi_dc
) ,arg
,m
)))
3007 (defun simplim%inverse_jacobi_cs
(arg m
)
3008 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3010 `((%inverse_jacobi_cs
) ,arg
,m
)))
3012 (defun simplim%inverse_jacobi_ds
(arg m
)
3013 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3015 `((%inverse_jacobi_ds
) ,arg
,m
)))
3017 (setf (get '%signum
'simplim%function
) 'simplim%signum
)
3019 (defun simplim%signum
(e x pt
)
3020 (let* ((e (limit (cadr e
) x pt
'think
)) (sgn (mnqp e
0)))
3021 (cond ((eq t sgn
) (take '(%signum
) e
)) ;; limit of argument of signum is not zero
3022 ((eq nil sgn
) '$und
) ;; limit of argument of signum is zero (noncontinuous)
3023 (t (throw 'limit nil
))))) ;; don't know
3025 ;; more functions for limit to handle
3027 (defun lfibtophi (e)
3028 (cond ((not (involve e
'($fib
))) e
)
3029 ((eq (caar e
) '$fib
)
3030 (let ((lnorecurse t
))
3031 ($fibtophi
(list '($fib
) (lfibtophi (cadr e
))) lnorecurse
)))
3033 (mapcar #'lfibtophi
(cdr e
))))))
3035 ;;; FOLLOWING CODE MAKES $LDEFINT WORK
3037 (defmfun $ldefint
(exp var ll ul
&aux $logabs ans a1 a2
)
3038 (setq $logabs t ans
(sinint exp var
)
3039 a1
(toplevel-$limit ans var ul
'$minus
)
3040 a2
(toplevel-$limit ans var ll
'$plus
))
3041 (and (member a1
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3042 (setq a1
(nounlimit ans var ul
)))
3043 (and (member a2
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3044 (setq a2
(nounlimit ans var ll
)))
3045 ($expand
(m- a1 a2
)))
3047 (defun nounlimit (exp var val
)
3048 (setq exp
(restorelim exp
))
3049 (nconc (list '(%limit
) exp var
(ridofab val
))
3050 (cond ((eq val
'$zeroa
) '($plus
))
3051 ((eq val
'$zerob
) '($minus
)))))
3053 ;; replace noun form of %derivative and indefinite %integrate with gensym.
3054 ;; prevents substitution x -> x+1 for limit('diff(x+2,x), x, 1)
3056 ;; however, this doesn't work for limit('diff(x+2,x)/x, x, inf)
3057 ;; because the rest of the limit code thinks the gensym is const wrt x.
3059 (cond ((atom exp
) exp
)
3060 ((or (eq '%derivative
(caar exp
))
3061 (and (eq '%integrate
(caar exp
)) ; indefinite integral
3062 (null (cdddr exp
))))
3063 (hidelim exp
(caar exp
)))
3064 (t (cons (car exp
) (mapcar 'hide
(cdr exp
))))))
3066 (defun hidelim (exp func
)
3067 (setq func
(gensym))
3075 (nounlimit e var val
)
3078 ;;;Used by Defint also.
3080 (or (involve e
'(%sin %cos %tan
))
3081 (among '$%i
(%einvolve e
))))
3083 (defun %einvolve
(e)
3084 (when (among '$%e e
) (%einvolve01 e
)))
3086 (defun %einvolve01
(e)
3087 (cond ((atom e
) nil
)
3091 (among var
(caddr e
)))
3093 (t (some #'%einvolve
(cdr e
)))))
3095 (declare-top (unspecial *indicator nn
* dn
* exp var val origval taylored
3096 $tlimswitch logcombed lhp? lhcount $ratfac
))
3102 ;; "On Computing Limits in a Symbolic Manipulation System"
3103 ;; PhD Dissertation ETH Zurich 1996
3105 ;; The algorithm identifies the most rapidly varying (MRV) subexpression,
3106 ;; replaces it with a new variable w, rewrites the expression in terms
3107 ;; of the new variable, and then repeats.
3109 ;; The algorithm doesn't handle oscillating functions, so it can't do things like
3110 ;; limit(sin(x)/x, x, inf).
3112 ;; To handle limits involving functions like gamma(x) and erf(x), the
3113 ;; gruntz algorithm requires them to be written in terms of asymptotic
3114 ;; expansions, which maxima cannot currently do.
3116 ;; The algorithm assumes that everything is real, so it can't
3117 ;; currently handle limit((-2)^x, x, inf).
3119 ;; This is one of the methods used by maxima's $limit.
3120 ;; It is also directly available to the user as $gruntz.
3123 ;; most rapidly varying subexpression of expression exp with respect to limit variable var.
3124 ;; returns a list of subexpressions which are in the same MRV equivalence class.
3125 (defun mrv (exp var
)
3126 (cond ((freeof var exp
)
3131 (mrv-max (mrv (cadr exp
) var
)
3132 (mrv (m*l
(cddr exp
)) var
)
3135 (mrv-max (mrv (cadr exp
) var
)
3136 (mrv (m+l
(cddr exp
)) var
)
3139 (cond ((freeof var
(caddr exp
))
3140 (mrv (cadr exp
) var
))
3141 ((member (limitinf (logred exp
) var
) '($inf $minf
) :test
#'eq
)
3142 (mrv-max (list exp
) (mrv (caddr exp
) var
) var
))
3143 (t (mrv-max (mrv (cadr exp
) var
) (mrv (caddr exp
) var
) var
))))
3145 (mrv (cadr exp
) var
))
3146 ((equal (length (cdr exp
)) 1)
3147 (mrv (cadr exp
) var
))
3148 ((equal (length (cdr exp
)) 2)
3149 (mrv-max (mrv (cadr exp
) var
)
3150 (mrv (caddr exp
) var
)
3152 (t (tay-error "mrv not implemented" exp
))))
3154 ;; takes two lists of expresions, f and g, and limit variable var.
3155 ;; members in each list are assumed to be in same MRV equivalence
3156 ;; class. returns MRV set of the union of the inputs - either f or g
3157 ;; or the union of f and g.
3158 (defun mrv-max (f g var
)
3165 (return (union f g
))))
3166 (let ((c (mrv-compare (car f
) (car g
) var
)))
3172 (return (union f g
)))
3173 (t (merror "MRV-MAX: expected '>' '<' or '='; found: ~M" c
))))))
3175 (defun mrv-compare (a b var
)
3176 (let ((c (limitinf (m// `((%log
) ,a
) `((%log
) ,b
)) var
)))
3179 ((member c
'($inf $minf
) :test
#'eq
)
3183 ;; rewrite expression exp by replacing members of MRV set omega with
3184 ;; expressions in terms of new variable wsym. return cons pair of new
3185 ;; version of exp and the log of the new variable wsym.
3186 (defun mrv-rewrite (exp omega var wsym
)
3187 (setq omega
(sort omega
(lambda (x y
) (> (length (mrv x var
))
3188 (length (mrv y var
))))))
3189 (let* ((g (car (last omega
)))
3191 (sig (equal (mrv-sign logg var
) 1))
3192 (w (if sig
(m// 1 wsym
) wsym
))
3193 (logw (if sig
(m* -
1 logg
) logg
)))
3194 (mapcar (lambda (x y
)
3195 ;;(mtell "y:~M x:~M exp:~M~%" y x exp)
3196 (setq exp
(syntactic-substitute y x exp
)))
3198 (mapcar (lambda (f) ;; rewrite each element of omega
3199 (let* ((logf (logred f
))
3200 (c (mrv-leadterm (m// logf logg
) var nil
)))
3201 (cond ((not (equal (cadr c
) 0))
3202 (merror "MRV-REWRITE: expected leading term to be constant in ~M" c
)))
3203 ;;(mtell "logg: ~M logf: ~M~%" logg logf)
3206 (m* (car c
) logg
))))))
3210 ;;; if log w(x) = h(x), rewrite all subexpressions of the form
3211 ;;; log(f(x)) as log(w^-c f(x)) + c h(x) with c the unique constant
3212 ;;; such that w^-c f(x) is strictly less rapidly varying than w.
3213 (defun mrv-rewrite-logs (exp wsym logw
)
3214 (cond ((atom exp
) exp
)
3216 (not (freeof wsym exp
)))
3217 (let* ((f (cadr exp
))
3218 (c ($lopow
(calculate-series f wsym
)
3221 (m* (m^ wsym
(m- c
))
3222 (mrv-rewrite-logs f wsym logw
)))
3227 (mrv-rewrite-logs e wsym logw
))
3230 ;; returns list of two elements: coeff and exponent of leading term of exp,
3231 ;; after rewriting exp in term of its MRV set omega.
3232 (defun mrv-leadterm (exp var omega
)
3233 (prog ((new-omega ()))
3234 (cond ((freeof var exp
)
3235 (return (list exp
0))))
3236 (dolist (term omega
)
3237 (cond ((subexp exp term
)
3238 (push term new-omega
))))
3239 (setq omega new-omega
)
3241 (setq omega
(mrv exp var
))))
3242 (cond ((member var omega
:test
#'eq
)
3243 (let* ((omega-up (mrv-moveup omega var
))
3244 (e-up (car (mrv-moveup (list exp
) var
)))
3245 (mrv-leadterm-up (mrv-leadterm e-up var omega-up
)))
3246 (return (mrv-movedown mrv-leadterm-up var
)))))
3247 (destructuring-let* ((wsym (gensym "w"))
3250 ((f . logw
) (mrv-rewrite exp omega var wsym
))
3251 (series (calculate-series (mrv-rewrite-logs f wsym logw
)
3253 (setq series
(maxima-substitute logw
`((%log
) ,wsym
) series
))
3254 (setq lo
($lopow series wsym
))
3255 (when (or (not ($constantp lo
))
3256 (not (free series
'%derivative
)))
3257 ;; (mtell "series: ~M lo: ~M~%" series lo)
3258 (tay-error "error in series expansion" f
))
3259 (setq coef
($coeff series wsym lo
))
3260 ;;(mtell "exp: ~M f: ~M~%" exp f)
3261 ;;(mtell "series: ~M~%coeff: ~M~%pow: ~M~%" series coef lo)
3262 (return (list coef lo
)))))
3264 (defun mrv-moveup (l var
)
3265 (mapcar (lambda (exp)
3266 (simplify-log-of-exp
3267 (syntactic-substitute `((mexpt) $%e
,var
) var exp
)))
3270 (defun mrv-movedown (l var
)
3271 (mapcar (lambda (exp) (syntactic-substitute `((%log simp
) ,var
) var exp
))
3274 ;; test whether sub is a subexpression of exp
3275 (defun subexp (exp sub
)
3276 (not (equal (maxima-substitute 'dummy
3281 ;; Generate $lhospitallim terms of taylor expansion.
3282 ;; Ideally we would use a lazy series representation that generates
3283 ;; more terms as higher order terms cancel.
3284 (defun calculate-series (exp var
)
3285 (assume `((mgreaterp) ,var
0))
3286 (putprop var t
'internal
);; keep var from appearing in questions to user
3287 (let ((series ($taylor exp var
0 $lhospitallim
)))
3288 (forget `((mgreaterp) ,var
0))
3291 (defun mrv-sign (exp var
)
3292 (cond ((freeof var exp
)
3293 (let ((sign ($sign
($radcan exp
))))
3294 (cond ((eq sign
'$zero
)
3300 (t (tay-error " cannot determine mrv-sign" exp
)))))
3304 (* (mrv-sign (cadr exp
) var
)
3305 (mrv-sign (m*l
(cddr exp
)) var
)))
3307 (equal (mrv-sign (cadr exp
) var
) 1))
3310 (cond ((equal (mrv-sign (cadr exp
) var
) -
1)
3311 (tay-error " complex expression in gruntz limit" exp
)))
3312 (mrv-sign (m+ -
1 (cadr exp
)) var
))
3314 (mrv-sign (limitinf exp var
) var
))
3315 (t (tay-error " cannot determine mrv-sign" exp
))))
3317 ;; gruntz algorithm for limit of exp as var goes to positive infinity
3318 (defun limitinf (exp var
)
3319 (prog (($exptsubst nil
))
3320 (cond ((freeof var exp
)
3322 (destructuring-let* ((c0-e0 (mrv-leadterm exp var nil
))
3325 (sig (mrv-sign e0 var
)))
3326 (cond ((equal sig
1)
3329 (cond ((equal (mrv-sign c0 var
) 1)
3331 ((equal (mrv-sign c0 var
) -
1)
3334 (return (limitinf c0 var
)))))))
3336 ;; user-level function equivalent to $limit.
3337 ;; direction must be specified if limit point is not infinite
3338 ;; The arguments are checked and a failure of taylor is catched.
3340 (defmfun $gruntz
(expr var val
&rest rest
)
3342 (when (> (length rest
) 1)
3344 (intl:gettext
"gruntz: too many arguments; expected just 3 or 4")))
3345 (setq dir
(car rest
))
3346 (when (and (not (member val
'($inf $minf $zeroa $zerob
)))
3347 (not (member dir
'($plus $minus
))))
3349 (intl:gettext
"gruntz: direction must be 'plus' or 'minus'")))
3351 (catch 'taylor-catch
3352 (let ((silent-taylor-flag t
))
3353 (declare (special silent-taylor-flag
))
3354 (gruntz1 expr var val dir
))))
3355 (if (or (null ans
) (eq ans t
))
3357 `(($gruntz simp
) ,expr
,var
, val
,dir
)
3358 `(($gruntz simp
) ,expr
,var
,val
))
3361 ;; This function is for internal use in $limit.
3362 (defun gruntz1 (exp var val
&rest rest
)
3363 (cond ((> (length rest
) 1)
3364 (merror (intl:gettext
"gruntz: too many arguments; expected just 3 or 4"))))
3365 (let (($logexpand t
) ; gruntz needs $logexpand T
3366 (newvar (gensym "w"))
3368 (cond ((eq val
'$inf
)
3371 (setq exp
(maxima-substitute (m* -
1 newvar
) var exp
)))
3373 (setq exp
(maxima-substitute (m// 1 newvar
) var exp
)))
3375 (setq exp
(maxima-substitute (m// -
1 newvar
) var exp
)))
3377 (setq exp
(maxima-substitute (m+ val
(m// 1 newvar
)) var exp
)))
3379 (setq exp
(maxima-substitute (m+ val
(m// -
1 newvar
)) var exp
)))
3380 (t (merror (intl:gettext
"gruntz: direction must be 'plus' or 'minus'; found: ~M") dir
)))
3381 (limitinf exp newvar
)))
3383 ;; substitute y for x in exp
3384 ;; similar to maxima-substitute but does not simplify result
3385 (defun syntactic-substitute (y x exp
)
3386 (cond ((alike1 x exp
) y
)
3389 (mapcar (lambda (exp)
3390 (syntactic-substitute y x exp
))
3393 ;; log(exp(subexpr)) -> subexpr
3394 ;; without simplifying entire exp
3395 (defun simplify-log-of-exp (exp)
3396 (cond ((atom exp
) exp
)
3399 (eq '$%e
(cadadr exp
)))
3402 (mapcar #'simplify-log-of-exp
