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
))
103 ;; The function toplevel-$limit sets $numer, $%enumer, %emode, and
104 ;; $%e_to_numlog to false. When any of these option variables are true,
105 ;; we resimplify the limit in the current context.
107 (when (or $numer $%enumer $%emode $%e_to_numlog
)
108 (setq first-try
(resimplify first-try
)))
111 (defun toplevel-$limit
(&rest args
)
112 (let ((limit-assumptions ())
113 (old-integer-info ())
120 (declare (special limit-assumptions old-integer-info
121 $keepfloat limit-top
))
123 (setq old-integer-info integer-info
)
124 (setq integer-info
()))
127 (let ((exp1 ()) (lhcount $lhospitallim
) (*behavior-count-now
* 0)
128 (exp ()) (var ()) (val ()) (dr ())
129 (*indicator
()) (taylored ()) (origval ())
130 (logcombed ()) (lhp?
())
131 (varlist ()) (ans ()) (genvar ()) (loginprod?
())
132 (limit-answers ()) (limitp t
) (simplimplus-problems ())
133 (lenargs (length args
))
135 (declare (special lhcount
*behavior-count-now
* exp var val
*indicator
136 taylored origval logcombed lhp?
137 varlist genvar loginprod? limitp
))
139 (unless (or (= lenargs
3) (= lenargs
4) (= lenargs
1))
141 ;; Is it a LIST of Things?
142 (when (setq ans
(apply #'limit-list args
))
144 (setq exp1
(specrepcheck (first args
)))
145 (when (and (atom exp1
)
146 (member exp1
'(nil t
)))
147 ;; The expression is 'T or 'NIL. Return immediately.
150 (setq var
(setq genfoo
(gensym)) ; Use a gensym. Not foo.
153 (setq var
(second args
))
154 (when ($constantp var
)
155 (merror (intl:gettext
"limit: second argument must be a variable, not a constant; found: ~M") var
))
156 (unless (or ($subvarp var
) (atom var
))
157 (merror (intl:gettext
"limit: variable must be a symbol or subscripted symbol; found: ~M") var
))
158 (setq val
(infsimp (third args
)))
159 ;; infsimp converts -inf to minf. it also converts -infinity to
160 ;; infinity, although perhaps this should generate the error below.
161 (when (and (not (atom val
))
162 (some #'(lambda (x) (not (freeof x val
)))
164 (merror (intl:gettext
"limit: third argument must be a finite value or one of: inf, minf, infinity; found: ~M") val
))
165 (when (eq val
'$zeroa
) (setq dr
'$plus
))
166 (when (eq val
'$zerob
) (setq dr
'$minus
))))
168 (unless (member (fourth args
) '($plus $minus
) :test
#'eq
)
169 (merror (intl:gettext
"limit: direction must be either 'plus' or 'minus'; found: ~M") (fourth args
)))
170 (setq dr
(fourth args
))))
171 (if (and (atom var
) (not (among var val
)))
173 (let ((realvar var
)) ;; Var is funny so make it a gensym.
175 (setq exp
(maxima-substitute var realvar exp1
))
176 (putprop var realvar
'limitsub
)))
177 (unless (or $limsubst
(eq var genfoo
))
178 (when (limunknown exp
)
179 (return `((%limit
) ,@(cons exp1
(cdr args
))))))
180 (setq varlist
(ncons var
) genvar nil origval val
)
181 ;; Transform limits to minf to limits to inf by
182 ;; replacing var with -var everywhere.
183 (when (eq val
'$minf
)
186 exp
(subin (m* -
1 var
) exp
)))
188 ;; Transform the limit value.
189 (unless (infinityp val
)
191 (let ((*atp
* t
) (realvar var
))
192 ;; *atp* prevents substitution from applying to vars
193 ;; bound by %sum, %product, %integrate, %limit
195 (putprop var t
'internal
)
196 (setq exp
(derivative-subst exp val var realvar
))
197 (setq exp
(maxima-substitute (m+ val var
) realvar exp
))))
198 (setq val
(cond ((eq dr
'$plus
) '$zeroa
)
199 ((eq dr
'$minus
) '$zerob
)
203 ;; Make assumptions about limit var being very small or very large.
204 ;; Assumptions are forgotten upon exit.
205 (unless (= lenargs
1)
206 (limit-context var val dr
))
208 ;; Resimplify in light of new assumptions.
209 (setq exp
(resimplify
213 (limitsimp ($expand exp
1 0) var
))))))
215 (if (not (or (real-epsilonp val
) ;; if direction of limit not specified
217 (setq ans
(both-side exp var val
)) ;; compute from both sides
218 (let ((d (catch 'mabs
(mabs-subst exp var val
))))
219 (cond ;; otherwise try to remove absolute value
220 ((eq d
'$und
) (return '$und
))
223 (setq ans
(limit-catch exp var val
));; and find limit from one side
227 (setq ans
(catch 'taylor-catch
228 (let ((silent-taylor-flag t
))
229 (declare (special silent-taylor-flag
))
230 (gruntz1 exp var val
)))))
232 ;; try taylor series expansion if simple limit didn't work
233 (if (and (null ans
) ;; if no limit found and
234 $tlimswitch
;; user says ok to use taylor and
235 (not limit-using-taylor
));; not already doing taylor
236 (let ((limit-using-taylor t
))
237 (declare (special limit-using-taylor
))
238 (setq ans
(limit-catch exp var val
))))))
241 (return (clean-limit-exp ans
))
242 (return (cons '(%limit
) args
))))) ;; failure: return nounform
243 (restore-assumptions))))
245 (defun clean-limit-exp (exp)
246 (setq exp
(restorelim exp
))
247 (if preserve-direction exp
(ridofab exp
)))
249 ;; Users who want limit to map over equality (mequal) will need to do that
251 (defun limit-list (exp1 &rest rest
)
252 (if (and (mbagp exp1
) (not (mequalp exp1
)))
253 `(,(car exp1
) ,@(mapcar #'(lambda (x) (apply #'toplevel-$limit
`(,x
,@rest
))) (cdr exp1
)))
256 (defun limit-context (var val direction
) ;Only works on entry!
258 (assume '((mgreaterp) lim-epsilon
0))
259 (assume '((mgreaterp) prin-inf
100000000))
260 (setq limit-assumptions
(make-limit-assumptions var val direction
))
265 (defun make-limit-assumptions (var val direction
)
266 (let ((new-assumptions))
267 (cond ((or (null var
) (null val
))
269 ((and (not (infinityp val
)) (null direction
))
272 `(,(assume `((mgreaterp) ,var
100000000)) ,@new-assumptions
))
274 `(,(assume `((mgreaterp) -
100000000 ,var
)) ,@new-assumptions
))
275 ((eq direction
'$plus
)
276 `(,(assume `((mgreaterp) ,var
0)) ,@new-assumptions
)) ;All limits around 0
277 ((eq direction
'$minus
)
278 `(,(assume `((mgreaterp) 0 ,var
)) ,@new-assumptions
))
282 (defun restore-assumptions ()
283 ;;;Hackery until assume and forget take reliable args. Nov. 9 1979.
285 (do ((assumption-list limit-assumptions
(cdr assumption-list
)))
286 ((null assumption-list
) t
)
287 (forget (car assumption-list
)))
288 (forget '((mgreaterp) lim-epsilon
0))
289 (forget '((mgreaterp) prin-inf
100000000))
290 (cond ((and (not (null integer-info
))
292 (do ((list integer-info
(cdr list
)))
294 (i-$remove
`(,(cadar list
) ,(caddar list
))))
295 (setq integer-info old-integer-info
))))
297 ;; The optional arg allows the caller to decide on the value of
298 ;; preserve-direction. Default is nil, since we immediately ridofab.
299 (defun both-side (exp var val
&optional
(preserve nil
))
300 (let* ((preserve-direction preserve
)
301 (la (toplevel-$limit exp var val
'$plus
)) lb
)
302 ; Immediately propagate an und without trying the
304 (when (eq la
'$und
) (return-from both-side
'$und
))
305 (setf lb
(toplevel-$limit exp var val
'$minus
))
306 ; Immediately propagate an und
307 (when (eq lb
'$und
) (return-from both-side
'$und
))
308 (let ((ra (ridofab la
))
310 (cond ((eq t
(meqp ra rb
))
314 ; Maxima does not consider equal(ind,ind) to be true, but
315 ; if both one-sided limits are ind then we want to call
316 ; the two-sided limit ind (e.g., limit(sin(1/x),x,0)).
318 ((or (not (free la
'%limit
))
319 (not (free lb
'%limit
)))
322 (let ((infa (infinityp la
))
323 (infb (infinityp lb
)))
324 (cond ((and infa infb
)
325 ; inf + minf => infinity
332 (defun limunknown (f)
333 (catch 'limunknown
(limunknown1 (specrepcheck f
))))
335 (defun limunknown1 (f)
336 (cond ((mapatom f
) nil
)
337 ((or (not (safe-get (caar f
) 'operators
))
338 (member (caar f
) '(%sum %product mncexpt
) :test
#'eq
)
339 ;;Special function code here i.e. for li[2](x).
340 (and (eq (caar f
) 'mqapply
)
341 (not (get (subfunname f
) 'specsimp
))))
342 (if (not (free f var
)) (throw 'limunknown t
)))
343 (t (mapc #'limunknown1
(cdr f
)) nil
)))
346 (if (involve e
'(%gamma
)) (setq e
($makefact e
)))
347 (cond ((involve e
'(mfactorial))
348 (setq e
(simplify ($minfactorial e
))))
352 ;; or nil if sign unknown or complex
354 (let ((z (ridofab z
)))
355 (if (not (free z var
)) (setq z
(toplevel-$limit z var val
)))
356 (let ((*complexsign
* t
))
357 (let ((sign (if *getsignl-asksign-ok
* ($asksign z
) ($sign z
))))
358 (cond ((eq sign
'$pos
) 1)
360 ((eq sign
'$zero
) 0))))))
362 (defun restorelim (exp)
363 (cond ((null exp
) nil
)
364 ((atom exp
) (or (and (symbolp exp
) (get exp
'limitsub
)) exp
))
365 ((and (consp (car exp
)) (eq (caar exp
) 'mrat
))
367 (cons (restorelim (cadr exp
))
368 (restorelim (cddr exp
)))))
369 (t (cons (car exp
) (mapcar #'restorelim
(cdr exp
))))))
372 (defun mabs-subst (exp var val
) ; RETURNS EXP WITH MABS REMOVED, OR THROWS.
373 (let ((d (involve exp
'(mabs)))
377 ((not (and (equal ($imagpart
(let ((v (limit-catch d var val
)))
378 ;; The above call might
379 ;; throw 'limit, so we
380 ;; need to catch it. If
382 ;; limit without ABS, we
383 ;; assume the limit is
384 ;; undefined. Is this
385 ;; right? Anyway, this
386 ;; fixes Bug 1548643.
391 (equal ($imagpart var
) 0)))
392 (cond ((eq arglim
'$infinity
)
393 ;; Check for $infinity as limit of argument.
396 (throw 'mabs
'retn
))))
397 (t (do ((ans d
(involve exp
'(mabs))) (a () ()))
399 (setq a
(mabs-subst ans var val
))
400 (setq d
(limit a var val t
))
404 (setq d
(behavior a var val
))
405 (if (zerop1 d
) (throw 'mabs
'retn
))))
408 (cond ((or (eq d
'$zeroa
) (eq d
'$inf
)
410 ;; fails on limit(abs(sin(x))/sin(x), x, inf)
411 (eq ($sign d
) '$pos
))
412 (setq exp
(maxima-substitute a
`((mabs) ,ans
) exp
)))
413 ((or (eq d
'$zerob
) (eq d
'$minf
)
414 (eq ($sign d
) '$neg
))
415 (setq exp
(maxima-substitute (m* -
1 a
) `((mabs) ,ans
) exp
)))
416 (t (throw 'mabs
'retn
))))
417 (t (throw 'mabs
'retn
))))))))))
419 ;; Called on an expression that might contain $INF, $MINF, $ZEROA, $ZEROB. Tries
420 ;; to simplify it to sort out things like inf^inf or inf+1.
422 (simpinf-ic exp
(count-general-inf exp
)))
424 (defun count-general-inf (expr)
425 (count-atoms-matching
426 (lambda (x) (or (infinityp x
) (real-epsilonp x
))) expr
))
428 (defun count-atoms-matching (predicate expr
)
429 "Count the number of atoms in the Maxima expression EXPR matching PREDICATE,
430 ignoring dummy variables and array indices."
432 ((atom expr
) (if (funcall predicate expr
) 1 0))
433 ;; Don't count atoms that occur as a limit of %integrate, %sum, %product,
435 ((member (caar expr
) dummy-variable-operators
)
436 (count-atoms-matching predicate
(cadr expr
)))
437 ;; Ignore array indices
438 ((member 'array
(car expr
)) 0)
440 for arg in
(cdr expr
)
441 summing
(count-atoms-matching predicate arg
)))))
443 (defun simpinf-ic (exp &optional infinity-count
)
445 ;; A very slow identity transformation...
448 ;; If there's only one infinity, we replace it by a variable and take the
449 ;; limit as that variable goes to infinity. Use $gensym in case we can't
450 ;; compute the answer and the limit leaks out.
451 (1 (let* ((val (or (inf-typep exp
) (epsilon-typep exp
)))
453 (expr (subst var val exp
))
454 (limit (toplevel-$limit expr var val
)))
456 ;; Now we look to see whether the computed limit is any simpler than
457 ;; what we shoved in (which we'll define as "doesn't contain EXPR as a
458 ;; subtree"). If so, return it.
459 ((not (subtree-p expr limit
:test
#'equal
))
462 ;; Otherwise, return the original form: apparently, we can't compute
463 ;; the limit we needed, and it's uglier than what we started with.
466 ;; If more than one infinity, we have to be a bit more careful.
468 (let* ((arguments (mapcar 'simpinf
(cdr exp
)))
469 (new-expression (cons (list (caar exp
)) arguments
))
472 ;; If any of the arguments are undefined, we are too.
473 ((among '$und arguments
) '$und
)
474 ;; If we ended up with something indeterminate, we punt and just return
476 ((amongl '(%limit $ind
) arguments
) exp
)
478 ;; Exponentiation & multiplication
479 ((mexptp exp
) (simpinf-expt (first arguments
) (second arguments
)))
480 ((mtimesp exp
) (simpinf-times arguments
))
482 ;; Down to at most one infinity? We do this after exponentiation to
483 ;; avoid zeroa^zeroa => 0^0, which will raise an error rather than just
484 ;; returning und. We do it after multiplication to avoid zeroa * inf =>
486 ((<= (setf infinities-left
(count-general-inf new-expression
)) 1)
487 (simpinf-ic new-expression infinities-left
))
490 ((mplusp exp
) (simpinf-plus arguments
))
493 (t new-expression
))))))
495 (defun simpinf-times (arguments)
496 (declare (special exp var val
))
497 ;; When we have a product, we need to spot that zeroa * zerob = zerob, zeroa *
498 ;; inf = und etc. Note that (SIMPINF '$ZEROA) => 0, so a nonzero atom is not
499 ;; an infinitesimal. Moreover, we can assume that each of ARGUMENTS is either
500 ;; a number, computed successfully by the recursive SIMPINF call, or maybe a
501 ;; %LIMIT noun-form (in which case, we aren't going to be able to tell the
504 ((member 0 arguments
)
506 ((find-if #'infinityp arguments
) '$und
)
507 ((every #'atom arguments
) 0)
510 ((member '$infinity arguments
)
511 (if (every #'atom arguments
)
515 (t (simplimit (cons '(mtimes) arguments
) var val
))))
517 (defun simpinf-expt (base exponent
)
518 ;; In the comments below, zero* represents one of 0, zeroa, zerob.
520 ;; TODO: In some cases we give up too early. E.g. inf^(2 + 1/inf) => inf^2
521 ;; (which should simplify to inf)
531 ((0 $zeroa $zerob
) '$und
)
532 (t (list '(mexpt) base exponent
))))
533 ;; minf^inf = infinity <== Or should it be und?
536 ;; minf^foo = minf^foo
541 ((0 $zeroa $zerob
) '$und
)
542 (t (list '(mexpt) base exponent
))))
546 ;; zero*^foo = zero*^foo
551 ((0 $zeroa $zerob
) '$und
)
552 (t (list '(mexpt) base exponent
))))
553 ;; a^b where a is pretty much anything except for a naked
554 ;; inf,minf,zeroa,zerob or 0.
557 ;; When a isn't crazy, try a^b = e^(b log(a))
558 ((not (amongl (append infinitesimals infinities
) base
))
559 (simpinf (m^
'$%e
(m* exponent
`((%log
) ,base
)))))
561 ;; No idea. Just return what we've found so far.
562 (t (list '(mexpt) base exponent
))))))
564 (defun simpinf-plus (arguments)
565 ;; We know that none of the arguments are infinitesimals, since SIMPINF never
566 ;; returns one of them. As such, we partition our arguments into infinities
567 ;; and everything else. The latter won't have any "hidden" infinities like
568 ;; lim(x,x,inf), since SIMPINF gave up on anything containing a %lim already.
569 (let ((bigs) (others))
570 (dolist (arg arguments
)
571 (cond ((infinityp arg
) (push arg bigs
))
572 (t (push arg others
))))
574 ;; inf + minf or the like
575 ((cdr (setf bigs
(delete-duplicates bigs
))) '$und
)
576 ;; inf + smaller + stuff
578 ;; I don't think this can happen, since SIMPINF goes back to the start if
579 ;; there are fewer than two infinities in the arguments, but let's be
581 (t (cons '(mplus) others
)))))
583 ;; Simplify expression with zeroa or zerob.
584 (defun simpab (small)
585 (cond ((null small
) ())
586 ((member small
'($zeroa $zerob $inf $minf $infinity
) :test
#'eq
) small
)
587 ((not (free small
'$ind
)) '$ind
) ;Not exactly right but not
588 ((not (free small
'$und
)) '$und
) ;causing trouble now.
589 ((mapatom small
) small
)
590 (t (let ((preserve-direction t
)
591 (new-small (subst (m^
'$inf -
1) '$zeroa
592 (subst (m^
'$minf -
1) '$zerob small
))))
593 (simpinf new-small
)))))
596 ;;;*I* INDICATES: T => USE LIMIT1,THINK, NIL => USE SIMPLIMIT.
597 (defun limit (exp var val
*i
*)
599 ((among '$und exp
) '$und
)
602 ((not (among var exp
))
603 (cond ((amongl '($inf $minf $infinity $ind
) exp
)
605 ((amongl '($zeroa $zerob
) exp
)
606 ;; Simplify expression with zeroa or zerob.
610 (t (putlimval exp
(cond ((and limit-using-taylor
613 (taylim exp var val
*i
*))
614 ((ratp exp var
) (ratlim exp
))
615 ((or (eq *i
* t
) (radicalp exp var
))
616 (limit1 exp var val
))
618 (cond ((or (mtimesp exp
) (mexptp exp
))
619 (limit1 exp var val
))
620 (t (simplimit exp var val
))))
621 (t (simplimit exp var val
)))))))
623 (defun limitsimp (exp var
)
624 (limitsimp-expt (sin-sq-cos-sq-sub exp
) var
))
625 ;;Hack for sin(x)^2+cos(x)^2.
627 ;; if var appears in base and power of expt,
628 ;; push var into power of of expt
629 (defun limitsimp-expt (exp var
)
630 (cond ((or (atom exp
)
632 (freeof var exp
)) exp
)
634 (not (freeof var
(cadr exp
)))
635 (not (freeof var
(caddr exp
))))
636 (m^
'$%e
(simplify `((%log
) ,exp
))))
637 (t (subst0 (cons (cons (caar exp
) ())
638 (mapcar #'(lambda (x)
639 (limitsimp-expt x var
))
643 (defun sin-sq-cos-sq-sub (exp) ;Hack ... Hack
644 (let ((arg (involve exp
'(%sin %cos
))))
647 (t (let ((new-exp ($substitute
(m+t
1 (m- (m^t
`((%sin simp
) ,arg
) 2)))
648 (m^t
`((%cos simp
) ,arg
) 2)
650 (m+t
1 (m- (m^t
`((%cos simp
) ,arg
) 2)))
651 (m^t
`((%sin simp
) ,arg
) 2)
653 (cond ((not (involve new-exp
'(%sin %cos
))) new-exp
)
656 (defun expand-trigs (x var
)
659 ((and (or (eq (caar x
) '%sin
)
661 (not (free (cadr x
) var
)))
663 ((member 'array
(car x
))
664 ;; Some kind of array reference. Return it.
666 (t (simplify (cons (ncons (caar x
))
667 (mapcar #'(lambda (x)
668 (expand-trigs x var
))
673 (cond ((not (involve e
674 '(%cot %csc %binomial
675 %sec %coth %sech %csch
676 %acot %acsc %asec %acoth
678 %jacobi_ns %jacobi_nc %jacobi_cs
679 %jacobi_ds %jacobi_dc
)))
681 (t ($ratsimp
(tansc1 e
)))))
683 (defun tansc1 (e &aux tem
)
685 ((and (setq e
(cons (car e
) (mapcar 'tansc1
(cdr e
)))) ()))
686 ((setq tem
(assoc (caar e
) '((%cot . %tan
) (%coth . %tanh
)
687 (%sec . %cos
) (%sech . %cosh
)
688 (%csc . %sin
) (%csch . %sinh
)) :test
#'eq
))
689 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
)) -
1.
)))
690 ((setq tem
(assoc (caar e
) '((%jacobi_nc . %jacobi_cn
)
691 (%jacobi_ns . %jacobi_sn
)
692 (%jacobi_cs . %jacobi_sc
)
693 (%jacobi_ds . %jacobi_sd
)
694 (%jacobi_dc . %jacobi_cd
)) :test
#'eq
))
695 ;; Converts Jacobi elliptic function to its reciprocal
697 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
) (third e
)) -
1.
)))
698 ((setq tem
(member (caar e
) '(%sinh %cosh %tanh
) :test
#'eq
))
699 (let (($exponentialize t
))
701 ((setq tem
(assoc (caar e
) '((%acsc . %asin
) (%asec . %acos
)
702 (%acot . %atan
) (%acsch . %asinh
)
703 (%asech . %acosh
) (%acoth . %atanh
)) :test
#'eq
))
704 (list (ncons (cdr tem
)) (m^t
(cadr e
) -
1.
)))
705 ((and (eq (caar e
) '%binomial
) (among var
(cdr e
)))
706 (m// `((mfactorial) ,(cadr e
))
707 (m* `((mfactorial) ,(m+t
(cadr e
) (m- (caddr e
))))
708 `((mfactorial) ,(caddr e
)))))
712 (cond ((not (involve ex
'(%sin %cos %tan %asin %acos %atan
713 %sinh %cosh %tanh %asinh %acosh %atanh
)))
719 ((eq (caar ex
) '%sinh
)
720 (m// (m+ (m^
'$%e
(cadr ex
)) (m- (m^
'$%e
(m- (cadr ex
)))))
722 ((eq (caar ex
) '%cosh
)
723 (m// (m+ (m^
'$%e
(cadr ex
)) (m^
'$%e
(m- (cadr ex
))))
725 ((and (member (caar ex
)
726 '(%sin %cos %tan %asin %acos %atan %sinh
727 %cosh %tanh %asinh %acosh %atanh
) :test
#'eq
)
730 (t (cons (car ex
) (mapcar #'hyperex0
(cdr ex
))))))
735 ;;Used by tlimit also.
736 (defun limit1 (exp var val
)
738 (let ((lhprogress? lhp?
)
741 (cond ((setq ans
(and (not (atom exp
)) (getlimval exp
)))
743 ((and (not (infinityp val
)) (setq ans
(simplimsubst val exp
)))
746 ;;;NUMDEN* => (values numerator denominator)
747 (multiple-value-bind (n dn
)
749 (cond ((not (among var dn
))
750 (return (simplimit (m// (simplimit n var val
) dn
) var val
)))
752 (return (simplimit (m* n
(simplimexpt dn -
1 (simplimit dn var val
) -
1)) var val
)))
754 (/#alike n
(car lhprogress?
))
755 (/#alike dn
(cdr lhprogress?
)))
756 (throw 'lhospital nil
)))
757 (return (limit2 n dn var val
))))))
762 (let ((deriv (sdiff (m// e f
) var
)))
764 ((=0 ($ratsimp deriv
)) t
)
767 ;; The standard /#alike function is possibly somewhat inefficient. Here is
768 ;; a possible replacement:
770 ;; If e/f is free of var (a special), return true. This code
771 ;; assumes that f is not zero. First, we test if e & f
772 ;; are alike--this check is relatively fast; second, we check
773 ;; if e/f is free of var.
775 ;;(defun /#alike (e f)
776 ;; (or (alike1 e f) (freeof var (sratsimp (div e f)))))
778 (defun limit2 (n dn var val
)
779 (prog (n1 d1 lim-sign gcp sheur-ans
)
780 (setq n
(hyperex n
) dn
(hyperex dn
))
781 ;;;Change to uniform limit call.
782 (cond ((infinityp val
)
783 (setq d1
(limit dn var val nil
))
784 (setq n1
(limit n var val nil
)))
785 (t (cond ((setq n1
(simplimsubst val n
)) nil
)
786 (t (setq n1
(limit n var val nil
))))
787 (cond ((setq d1
(simplimsubst val dn
)) nil
)
788 (t (setq d1
(limit dn var val nil
))))))
789 (cond ((or (null n1
) (null d1
)) (return nil
))
790 (t (setq n1
(sratsimp n1
) d1
(sratsimp d1
))))
791 (cond ((or (involve n
'(mfactorial)) (involve dn
'(mfactorial)))
792 (let ((ans (limfact2 n dn var val
)))
793 (cond (ans (return ans
))))))
794 (cond ((and (zerop2 n1
) (zerop2 d1
))
795 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
796 (return (colexpt n dn gcp
)))
797 ((and (real-epsilonp val
)
799 (not (free dn
'%log
)))
800 (return (liminv (m// n dn
))))
801 ((setq n1
(try-lhospital-quit n dn nil
))
803 ((and (zerop2 n1
) (not (member d1
'($ind $und
) :test
#'eq
))) (return 0))
805 (setq n1
(ridofab n1
))
806 (return (simplimtimes `(,n1
,(simplimexpt dn -
1 d1 -
1))))))
807 (setq n1
(ridofab n1
))
808 (setq d1
(ridofab d1
))
809 (cond ((or (eq d1
'$und
)
810 (and (eq n1
'$und
) (not (real-infinityp d1
))))
813 ;; At this point we have n1/$ind. Look if n1 is one of the
814 ;; infinities or zero.
815 (cond ((and (infinityp n1
) (eq ($sign dn
) '$pos
))
817 ((and (infinityp n1
) (eq ($sign dn
) '$neg
))
818 (return (simpinf (m* -
1 n1
))))
820 (or (eq ($sign dn
) '$pos
)
821 (eq ($sign dn
) '$neg
)))
824 ((eq n1
'$ind
) (return (cond ((infinityp d1
) 0)
827 ((and (real-infinityp d1
) (member n1
'($inf $und $minf
) :test
#'eq
))
828 (cond ((and (not (atom dn
)) (not (atom n
))
829 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
830 (return (colexpt n dn gcp
)))
832 (or (involve dn
'(mfactorial %gamma
))
833 (involve n
'(mfactorial %gamma
))))
834 (return (limfact n dn
))))))
835 ((eq n1 d1
) (setq lim-sign
1) (go cp
))
836 (t (setq lim-sign -
1) (go cp
))))
837 ((and (infinityp d1
) (infinityp n1
))
838 (setq lim-sign
(if (or (eq d1
'$minf
) (eq n1
'$minf
)) -
1 1))
840 (t (return (simplimtimes `(,n1
,(m^ d1 -
1))))))
841 cp
(setq n
($expand n
) dn
($expand dn
))
843 (let ((new-n (m+l
(maxi (cdr n
)))))
844 (cond ((not (alike1 new-n n
))
845 (return (limit (m// new-n dn
) var val
'think
))))
849 (let ((new-dn (m+l
(maxi (cdr dn
)))))
850 (cond ((not (alike1 new-dn dn
))
851 (return (limit (m// n new-dn
) var val
'think
))))
854 (setq sheur-ans
(sheur0 n1 d1
))
855 (cond ((or (member sheur-ans
'($inf $zeroa
) :test
#'eq
)
856 (free sheur-ans var
))
857 (return (simplimtimes `(,lim-sign
,sheur-ans
))))
858 ((and (alike1 sheur-ans dn
)
860 ((member (setq n1
(cond ((expfactorp n1 d1
) (expfactor n1 d1 var
))
862 '($inf $zeroa
) :test
#'eq
)
864 ((not (null (setq n1
(cond ((expfactorp n dn
) (expfactor n dn var
))
867 ((and (alike1 sheur-ans dn
) (not (mplusp n
))))
868 ((not (alike1 sheur-ans
(m// n dn
)))
869 (return (simplimit (m// ($expand
(m// n sheur-ans
))
870 ($expand
(m// dn sheur-ans
)))
873 (cond ((and (not (and (eq val
'$inf
) (expp n
) (expp dn
)))
874 (setq n1
(try-lhospital-quit n dn nil
))
879 ;; Test whether both n and dn have form
880 ;; product of poly^poly
881 (defun expfactorp (n dn
)
882 (do ((llist (append (cond ((mtimesp n
) (cdr n
))
884 (cond ((mtimesp dn
) (cdr dn
))
887 (exp? t
) ;IS EVERY ELEMENT SO FAR
888 (factor nil
)) ;A POLY^POLY?
892 (setq factor
(car llist
))
893 (setq exp?
(or (polyinx factor var
())
895 (polyinx (cadr factor
) var
())
896 (polyinx (caddr factor
) var
()))))))
898 (defun expfactor (n dn var
) ;Attempts to evaluate limit by grouping
899 (prog (highest-deg) ; terms with similar exponents.
900 (let ((new-exp (exppoly n
))) ;exppoly unrats expon
901 (setq n
(car new-exp
) ;and rtns deg of expons
902 highest-deg
(cdr new-exp
)))
903 (cond ((null n
) (return nil
))) ;nil means expon is not
904 (let ((new-exp (exppoly dn
))) ;a rat func.
905 (setq dn
(car new-exp
)
906 highest-deg
(max highest-deg
(cdr new-exp
))))
908 (= highest-deg
0)) ; prevent infinite recursion
912 (degree highest-deg
(1- degree
))
919 (limit (m// numerator denominator
)
923 (let ((newnumer-factor (get-newexp&factors
927 (setq numerator
(car newnumer-factor
)
928 numfactors
(cdr newnumer-factor
)))
929 (let ((newdenom-factor (get-newexp&factors
933 (setq denominator
(car newdenom-factor
)
934 denfactors
(cdr newdenom-factor
)))
935 (setq answer
(simplimit (list '(mexpt)
937 (m// numfactors denfactors
))
938 (cond ((> degree
0) var
)
942 (cond ((member answer
'($ind $und
) :test
#'equal
)
943 ;; cannot handle limit(exp(x*%i)*x, x, inf);
945 ((member answer
'($inf $minf
) :test
#'equal
)
946 ;; 0, zeroa, zerob are passed through to next iteration
947 (return (simplimtimes (list (m// numerator denominator
) answer
)))))))))
949 (defun exppoly (exp) ;RETURNS EXPRESSION WITH UNRATTED EXPONENTS
953 (exp (cond ((mtimesp exp
)
957 ((null exp
) (cons new-exp highest-deg
))
958 (setq factor
(car exp
))
960 (m* (cond ((or (not (mexptp factor
))
961 (not (ratp (caddr factor
) var
)))
965 (ratdegree (caddr factor
))))
966 (m^
(cadr factor
) (unrat (caddr factor
)))))
969 (defun unrat (exp) ;RETURNS UNRATTED EXPRESION
970 (multiple-value-bind (n d
)
972 (let ((tem ($divide n d
)))
974 (m// (caddr tem
) d
)))))
976 (defun get-newexp&factors
(exp degree var
) ;RETURNS (CONS NEWEXP FACTORS)
977 (do ((terms (cond ((mtimesp exp
)(cdr exp
)) ; SUCH THAT
978 (t (ncons exp
))) ; NEWEXP*FACTORS^(VAR^DEGREE)
979 (cdr terms
)) ; IS EQUAL TO EXP.
986 (setq factor
(car terms
))
987 (cond ((not (mexptp factor
))
989 (setq factors
(m* factor factors
)))
990 (t (setq newexp
(m* factor newexp
)))))
992 (= (ratdegree (caddr factor
))
994 (setq factors
(m* (m^
(cadr factor
)
995 (leading-coef (caddr factor
)))
997 newexp
(m* (m^
(cadr factor
)
999 (m* (leading-coef (caddr factor
))
1002 (t (setq newexp
(m* factor newexp
))))))
1004 (defun leading-coef (rat)
1005 (ratlim (m// rat
(m^ var
(ratdegree rat
)))))
1007 (defun ratdegree (rat)
1008 (multiple-value-bind (n d
)
1010 (- (deg n
) (deg d
))))
1012 (defun limfact2 (n d var val
)
1013 (let ((n1 (reflect0 n var val
))
1014 (d1 (reflect0 d var val
)))
1015 (cond ((and (alike1 n n1
)
1018 (t (limit (m// n1 d1
) var val
'think
)))))
1020 ;; takes expression and returns operator at front with all flags removed
1021 ;; except array flag.
1022 ;; array flag must match for alike1 to consider two things to be the same.
1023 ;; ((MTIMES SIMP) ... ) => (MTIMES)
1024 ;; ((PSI SIMP ARRAY) 0) => (PSI ARRAY)
1025 (defun operator-with-array-flag (exp)
1026 (cond ((member 'array
(car exp
) :test
#'eq
)
1027 (list (caar exp
) 'array
))
1028 (t (list (caar exp
)))))
1030 (defun reflect0 (exp var val
)
1031 (cond ((atom exp
) exp
)
1032 ((and (eq (caar exp
) 'mfactorial
)
1033 (let ((argval (limit (cadr exp
) var val
'think
)))
1034 (or (eq argval
'$minf
)
1035 (and (numberp argval
)
1037 (reflect (cadr exp
)))
1038 (t (cons (operator-with-array-flag exp
)
1041 (reflect0 term var val
)))
1044 (defun reflect (arg)
1047 (m^
(list (ncons 'mfactorial
)
1051 (m^
(list (ncons '%sin
)
1055 (defun limfact (n d
)
1057 (setq n
(stirling0 n
)
1059 (setq ans
(toplevel-$limit
(m// n d
) var
'$inf
))
1060 (cond ((and (atom ans
)
1061 (not (member ans
'(und ind
) :test
#'eq
))) ans
)
1062 ((eq (caar ans
) '%limit
) ())
1065 ;; substitute asymptotic approximations for gamma, factorial, and
1067 (defun stirling0 (e)
1069 ((and (setq e
(cons (car e
) (mapcar 'stirling0
(cdr e
))))
1071 ((and (eq (caar e
) '%gamma
)
1072 (eq (limit (cadr e
) var val
'think
) '$inf
))
1073 (stirling (cadr e
)))
1074 ((eq (caar e
) 'mfactorial
)
1075 (let ((n (limit (cadr e
) var val
'think
)))
1077 (m* (cadr e
) (stirling (cadr e
))))
1078 ((and (integerp n
) (< n
0))
1080 (div (power -
1 n
) (mul (ftake 'mfactorial n
) (add var n
))))
1082 ((and (eq (caar e
) 'mqapply
) ;; polylogarithm
1083 (eq (subfunname e
) '$li
)
1084 (integerp (car (subfunsubs e
))))
1085 (li-asymptotic-expansion (m- (car (subfunsubs e
)) 1)
1086 (car (subfunsubs e
))
1087 (car (subfunargs e
))))
1091 (maxima-substitute x
'$z
1093 ((mexpt simp
) 2 ((rat simp
) 1 2))
1094 ((mexpt simp
) $%pi
((rat simp
) 1 2))
1095 ((mexpt simp
) $z
((mplus simp
) ((rat simp
) -
1 2) $z
))
1096 ((mexpt simp
) $%e
((mtimes simp
) -
1 $z
)))))
1098 (defun no-err-sub (v e
&aux ans
)
1099 (let ((errorsw t
) (*zexptsimp? t
)
1101 ;; Don't print any error messages
1103 (declare (special errcatch
))
1104 ;; Should we just use IGNORE-ERRORS instead HANDLER-CASE here? I
1105 ;; (rtoy) am choosing the latter so that unexpected errors will
1106 ;; actually show up instead of being silently discarded.
1108 (setq ans
(catch 'errorsw
1110 (sratsimp (subin v e
)))))
1113 (cond ((null ans
) t
) ; Ratfun package returns NIL for failure.
1116 (defun extended-real-p (e)
1117 (member e
(list '$minf
'$zerob
'$zeroa
'$ind
'$und
'$inf
'$infinity
)))
1119 ;; Evaluate the limit of e as var (special) approaches v using direct substitution.
1120 ;; This function is the last of the chain of methods tried by limit. As such it
1121 ;; shouldn't call higher level functions such as simplimit or limit--doing so
1122 ;; would risk creating an infinite loop.
1124 ;; This function special cases sums, products, and powers. It declines to use
1125 ;; direct subsitution on any expression whose main operator has a simplim%function
1126 ;; function. Generally, limits of functions that have a simplim%function should be
1127 ;; handled by those specialized functions, not by simplimsubst. Having said this,
1128 ;; it's OK for simplimsubst to special case an operator and do direct substitution
1129 ;; when its OK. The cases of log, cosine, and sine happen often enough that these
1130 ;; functions are special cased.
1132 ;; Possibly simplimsubst should decline to do direct substituion on functions
1133 ;; that have a limit function, for example inverse_jacobi_ns.
1135 ;; For all other kinds of expressions, this function assumes that if substiution
1136 ;; doesn't result in an error (division by zero, for example), the function
1137 ;; is continuous at the limit point. That's dodgy. This dodginess underscores the
1138 ;; need to define a simplim%function functions that are not continous on their
1141 ;; This function sometimes recieves unsimplifed expressions. Maybe they should be
1142 ;; simplified, maybe not.
1144 ;; Locally setting $numer and $%enumer to nil keeps some limits, for example
1145 ;; limit(sin(x)/x,x,0) from returning 1.0 when numer is true. (Barton Willis)
1147 (defun simplimsubst (v e
)
1148 (let ((ans nil
) ($numer nil
) ($%enumer nil
) (ee))
1150 ;; When e is a number, return it.
1152 ;; When e is a mapatom, substitute v for var and return.
1153 (($mapatom e
) ($substitute v var e
))
1154 ;; Special case mexpt expressions. Decline direct subsitution for
1156 ((and (mexptp e
) (not (extended-real-p v
)))
1157 (let ((x (simplimsubst v
(second e
)))
1158 (n (simplimsubst v
(third e
))))
1159 ;; Decline direct substitution (DS) for 0^negative. Also decline
1160 ;; DS when x is on the negative real axis and n isn't an integer.
1161 ;; Additionally, we require that DS is OK for both x & n.
1162 (if (and x n
(not (and (zerop2 x
) (eq t
(mgqp 0 n
)))) ; not 0^negative
1163 (or (off-negative-real-axisp x
) (integerp n
)))
1164 (ftake 'mexpt x n
) nil
)))
1165 ;; Special case product and sum expressions. Again, we decline direct
1166 ;; subsitution for extended reals.
1167 ((and (or (mplusp e
) (mtimesp e
)) (not (extended-real-p v
)))
1168 (setq ee
(mapcar #'(lambda(q) (simplimsubst v q
)) (cdr e
)))
1169 (if (some #'(lambda (q) (eq q nil
)) ee
) nil
1170 (simplifya (cons (list (caar e
)) ee
) t
)))
1171 ;; Decline direct substitution for sums, products, and powers for
1172 ;; the extended real case.
1173 ((and (or (mexptp e
) (mplusp e
) (mtimesp e
)) (extended-real-p v
))
1175 ;; Special case log expressions--possibly the log case happens
1176 ;; often enough to make this worthwile.
1177 ((and (consp e
) (consp (car e
)) (eq '%log
(caar e
)))
1178 (let ((w (simplimsubst v
(cadr e
))))
1179 (if (and w
(off-negative-real-axisp w
)) (ftake '%log w
) nil
)))
1180 ;; Special case %cos and %sin expressions--we could special case others.
1181 ;; The lenient-realp check declines direct substitution for complex arguments--
1182 ;; this check could be relaxed.
1183 ((and (consp e
) (consp (car e
)) (or (member (caar e
) (list '%cos
'%sin
))))
1184 (let ((op (caar e
)))
1185 (setq e
(simplimsubst v
(second e
)))
1186 (if (and e
(lenient-realp e
) (not (extended-real-p e
))) (ftake op e
) nil
)))
1187 ;; Don't use direct substitution on expressions whose main operator has
1188 ;; a simplim%function.
1189 ((and (consp e
) (consp (car e
)) (get (caar e
) 'simplim%function
)) nil
)
1190 ;; The function no-err-sub returns true when there is an error
1191 ((not (eq t
(setq ans
(no-err-sub (ridofab v
) e
))))
1192 ;; Previously the condition (zerop2 ans) was (=0 ($radcan ans)). In
1193 ;; December 2021, the testsuite + the share testsuite only gives ans = 0,
1194 ;; making the radcan unneeded.
1195 (cond ((and (member v
'($zeroa $zerob
) :test
#'eq
) (zerop2 ans
))
1196 (setq ans
(behavior e var v
))
1197 (cond ((eql ans
1) '$zeroa
)
1198 ((eql ans -
1) '$zerob
)
1199 (t nil
))) ;behavior can't find direction
1201 ;; direct substitution fails, so return nil.
1204 ;;;returns (cons numerator denominator)
1206 (let ((e (factor (simplify e
)))
1212 (mapc #'forq
(cdr e
)))
1218 (setq numer
(car numer
)))
1220 (setq numer
(m*l numer
))))
1224 (setq denom
(car denom
)))
1226 (setq denom
(m*l denom
))))
1227 (values (factor numer
) (factor denom
))))
1229 ;;;FACTOR OR QUOTIENT
1230 ;;;Setq's the special vars numer and denom from numden*
1232 (cond ((and (mexptp e
)
1233 (not (freeof var e
))
1234 (null (pos-neg-p (caddr e
))))
1235 (push (m^
(cadr e
) (m* -
1.
(caddr e
))) denom
))
1236 (t (push e numer
))))
1238 ;;;Predicate to tell whether an expression is pos,zero or neg as var -> val.
1239 ;;;returns T if pos,zero. () if negative or don't know.
1240 (defun pos-neg-p (exp)
1241 (let ((ans (limit exp var val
'think
)))
1242 (cond ((and (not (member ans
'($und $ind $infinity
) :test
#'eq
))
1243 (equal ($imagpart ans
) 0))
1244 (let ((sign (getsignl ans
)))
1245 (cond ((or (equal sign
1)
1248 ((equal sign -
1) nil
))))
1251 (declare-top (unspecial n dn
))
1254 (cond ((radicalp e var
) nil
)
1255 ((member (caar e
) '(%log %sin %cos %tan %sinh %cosh %tanh mfactorial
1256 %asin %acos %atan %asinh %acosh %atanh
) :test
#'eq
) nil
)
1258 ((do ((e (cdr e
) (cdr e
)))
1260 (and (expp (car e
)) (return t
))))))
1264 (radicalp (cadr e
) var
)
1265 (among var
(caddr e
))
1266 (radicalp (caddr e
) var
)))
1269 (defun gcpower (a b
)
1270 ($gcd
(getexp a
) (getexp b
)))
1273 (cond ((and (mexptp exp
)
1274 (free (caddr exp
) var
)
1275 (eq (ask-integer (caddr exp
) '$integer
) '$yes
))
1277 ((mtimesp exp
) (getexplist (cdr exp
)))
1280 (defun getexplist (list)
1281 (cond ((null (cdr list
))
1282 (getexp (car list
)))
1283 (t ($gcd
(getexp (car list
))
1284 (getexplist (cdr list
))))))
1286 (defun limroot (exp power
)
1287 (cond ((or (atom exp
) (not (member (caar exp
) '(mtimes mexpt
) :test
#'eq
)))
1288 (limroot (list '(mexpt) exp
1) power
)) ;This is strange-JIM.
1289 ((mexptp exp
) (m^
(cadr exp
)
1290 (sratsimp (m* (caddr exp
) (m^ power -
1.
)))))
1291 (t (m*l
(mapcar #'(lambda (x)
1295 ;;NUMERATOR AND DENOMINATOR HAVE EXPONENTS WITH GCD OF GCP.
1296 ;;; Used to call simplimit but some of the transformations used here
1297 ;;; were not stable w.r.t. the simplifier, so try keeping exponent separate
1300 (defun colexpt (n dn gcp
)
1301 (let ((bas (m* (limroot n gcp
) (limroot dn
(m* -
1 gcp
))))
1304 (setq baslim
(limit bas var val
'think
))
1305 (setq expolim
(limit expo var val
'think
))
1306 (simplimexpt bas expo baslim expolim
)))
1308 ;; this function is responsible for the following bug:
1309 ;; limit(x^2 + %i*x, x, inf) -> inf (should be infinity)
1311 (cond ((member val
'($inf $infinity
) :test
#'eq
)
1312 (setq e
(maxima-substitute (m^t
'x -
1) var e
)))
1314 (setq e
(maxima-substitute (m^t -
1 (m^t
'x -
1)) var e
)))
1316 (setq e
(maxima-substitute (m- 'x
) var e
)))
1318 (setq e
(maxima-substitute 'x var e
)))
1319 ((setq e
(maxima-substitute (m+t
'x val
) var e
))))
1320 (destructuring-let* ((e (let (($ratfac
()))
1321 ($rat
(sratsimp e
) 'x
)))
1329 (sratsimp (m// ($ratdisrep
`(,h
,(locoef n g
) .
1))
1330 ($ratdisrep
`(,h
,(locoef d g
) .
1))))))
1332 (cond ((not (member val
'($zerob $zeroa $inf $minf
) :test
#'eq
))
1334 ((not (equal ($imagpart e
) 0))
1336 ((null (setq e
(getsignl ($realpart e
))))
1338 ((equal e
1) '$zeroa
)
1339 ((equal e -
1) '$zerob
)
1342 ((not (member val
'($zerob $zeroa $infinity $inf $minf
) :test
#'eq
))
1344 ((eq val
'$infinity
) '$infinity
)
1345 ((not (equal ($imagpart e
) 0)) '$infinity
)
1346 ((null (setq e
(getsignl ($realpart e
))))
1349 ((equal e -
1) '$minf
)
1353 (if (or (atom n
) (not (eq (car n
) x
)))
1358 (if (or (atom n
) (not (eq (car n
) x
)))
1362 (defun behavior (exp var val
) ; returns either -1, 0, 1.
1363 (if (= *behavior-count-now
* +behavior-count
+)
1365 (let ((*behavior-count-now
* (1+ *behavior-count-now
*)) pair sign
)
1366 (cond ((real-infinityp val
)
1367 (setq val
(cond ((eq val
'$inf
) '$zeroa
)
1368 ((eq val
'$minf
) '$zerob
)))
1369 (setq exp
(sratsimp (subin (m^ var -
1) exp
)))))
1370 (cond ((eq val
'$infinity
) 0) ; Needs more hacking for complex.
1372 (prog2 (setq pair
(partition exp var
1))
1373 (not (mtimesp (cdr pair
)))))
1374 (setq sign
(getsignl (car pair
)))
1375 (if (not (fixnump sign
))
1377 (mul sign
(behavior (cdr pair
) var val
))))
1378 ((and (=0 (no-err-sub (ridofab val
) exp
))
1380 (free (caddr exp
) var
)
1381 (equal (getsignl (caddr exp
)) 1))
1382 (let ((bas (cadr exp
)) (expo (caddr exp
)))
1383 (behavior-expt bas expo
)))
1384 (t (behavior-by-diff exp var val
))))))
1386 (defun behavior-expt (bas expo
)
1387 (let ((behavior (behavior bas var val
)))
1388 (cond ((= behavior
1) 1)
1390 ((eq (ask-integer expo
'$integer
) '$yes
)
1391 (cond ((eq (ask-integer expo
'$even
) '$yes
) 1)
1394 (cond ((evenp (cadr expo
)) 1)
1395 ((oddp (caddr expo
)) behavior
)
1399 (defun behavior-by-diff (exp var val
)
1400 (cond ((not (or (eq val
'$zeroa
) (eq val
'$zerob
))) 0)
1401 (t (let ((old-val val
) (old-exp exp
))
1402 (setq val
(ridofab val
))
1404 (exp (sratsimp (sdiff exp var
)) (sratsimp (sdiff exp var
)))
1406 (ans ())) ; This do wins by a return.
1407 ((> ct
1) 0) ; This loop used to run up to 5 times,
1408 ;; but the size of some expressions would blow up.
1409 (setq ans
(no-err-sub val exp
)) ;Why not do an EVENFN and ODDFN
1412 (return (behavior-numden old-exp var old-val
)))
1413 ((=0 ans
) ()) ;Do it again.
1414 (t (setq ans
(getsignl ans
))
1415 (cond (n (return ans
))
1417 (return (if (eq old-val
'$zeroa
) 1 -
1)))
1419 (return (if (eq old-val
'$zeroa
) -
1 1)))
1420 (t (return 0))))))))))
1422 (defun behavior-numden (exp var val
)
1423 (let ((num ($num exp
)) (denom ($denom exp
)))
1424 (cond ((equal denom
1) 0) ;Could be hacked more from here.
1425 (t (let ((num-behav (behavior num var val
))
1426 (denom-behav (behavior denom var val
)))
1427 (cond ((or (= num-behav
0) (= denom-behav
0)) 0)
1428 ((= num-behav denom-behav
) 1)
1431 (defun try-lhospital (n d ind
)
1432 ;;Make one catch for the whole bunch of lhospital trials.
1433 (let ((ans (lhospital-catch n d ind
)))
1434 (cond ((null ans
) ())
1435 ((not (free-infp ans
)) (simpinf ans
))
1436 ((not (free-epsilonp ans
)) (simpab ans
))
1439 (defun try-lhospital-quit (n d ind
)
1440 (let ((ans (or (lhospital-catch n d ind
)
1441 (lhospital-catch (m^ d -
1) (m^ n -
1) ind
))))
1442 (cond ((null ans
) (throw 'limit t
))
1443 ((not (free-infp ans
)) (simpinf ans
))
1444 ((not (free-epsilonp ans
)) (simpab ans
))
1447 (defun lhospital-catch (n d ind
)
1448 (cond ((> 0 lhcount
)
1449 (setq lhcount $lhospitallim
)
1450 (throw 'lhospital nil
))
1451 ((equal lhcount $lhospitallim
)
1452 (let ((lhcount (m+ lhcount -
1)))
1453 (catch 'lhospital
(lhospital n d ind
))))
1454 (t (setq lhcount
(m+ lhcount -
1))
1455 (prog1 (lhospital n d ind
)
1456 (setq lhcount
(m+ lhcount
1))))))
1457 ;;If this succeeds then raise LHCOUNT.
1460 (defun lhospital (n d ind
)
1461 (declare (special val lhp?
))
1463 (setq n
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr n
)))))
1465 (setq d
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr d
)))))
1466 (multiple-value-bind (n d
)
1468 (let (const nconst dconst
)
1469 (setq lhp?
(and (null ind
) (cons n d
)))
1470 (multiple-value-setq (nconst n
) (var-or-const n
))
1471 (multiple-value-setq (dconst d
) (var-or-const d
))
1473 (setq n
(stirling0 n
)) ;; replace factorial and %gamma
1474 (setq d
(stirling0 d
)) ;; with approximations
1476 (setq n
(sdiff n var
) ;; take derivatives for l'hospital
1479 (if (or (not (free n
'%derivative
)) (not (free d
'%derivative
)))
1480 (throw 'lhospital
()))
1481 (setq n
(expand-trigs (tansc n
) var
))
1482 (setq d
(expand-trigs (tansc d
) var
))
1484 (multiple-value-setq (const n d
) (remove-singularities n d
))
1485 (setq const
(m* const
(m// nconst dconst
)))
1486 (simpinf (let ((ans (if ind
1487 (limit2 n d var val
)
1488 (limit-numden n d val
))))
1489 ;; When the limit function returns, it's possible that it will return NIL
1490 ;; (gave up without finding a limit). It's also possible that it will
1491 ;; return something containing UND. We treat that as a failure too.
1492 (when (and ans
(freeof '$und ans
))
1493 (m* const ans
)))))))
1495 ;; Try to compute the limit of a quotient NUM/DEN, trying to massage the input
1496 ;; into a convenient form for LIMIT on the way.
1497 (defun limit-numden (n d val
)
1499 ;; For general arguments, the best approach seems to be to use
1500 ;; sratsimp to simplify the quotient as much as we can, then
1501 ;; $multthru, which splits it up into a sum (presumably useful
1502 ;; because limit(a+b) = limit(a) + limit(b) if the limits exist, and
1503 ;; the right hand side might be easier to calculate)
1505 ($multthru
(sratsimp (m// n d
))))
1507 ;; If we've already got a sum in the numerator, it seems to be
1508 ;; better not to recombine it. Call LIMIT on the whole lot, though,
1509 ;; because terms with infinite limits might cancel to give a finite
1512 (m+l
(mapcar #'(lambda (x)
1513 (sratsimp (m// x d
)))
1516 (limit expr var val
'think
)))
1518 ;; Heuristics for picking the right way to express a LHOSPITAL problem.
1519 (defun lhop-numden (num denom
)
1520 (declare (special var
))
1521 (cond ((let ((log-num (involve num
'(%log
))))
1522 (cond ((null log-num
) ())
1523 ((lessthan (num-of-logs (factor (sratsimp (sdiff (m^ num -
1) var
))))
1524 (num-of-logs (factor (sratsimp (sdiff num var
)))))
1525 (psetq num
(m^ denom -
1) denom
(m^ num -
1))
1528 ((let ((log-denom (involve denom
'(%log
))))
1529 (cond ((null log-denom
) ())
1530 ((lessthan (num-of-logs (sratsimp (sdiff (m^ denom -
1) var
)))
1531 (num-of-logs (sratsimp (sdiff denom var
))))
1532 (psetq denom
(m^ num -
1) num
(m^ denom -
1))
1535 ((let ((exp-num (%einvolve num
)))
1537 (cond ((%e-right-placep exp-num
)
1539 (t (psetq num
(m^ denom -
1)
1540 denom
(m^ num -
1)) t
)))
1542 ((let ((exp-den (%einvolve denom
)))
1544 (cond ((%e-right-placep exp-den
)
1546 (t (psetq num
(m^ denom -
1)
1547 denom
(m^ num -
1)) t
)))
1549 ((let ((scnum (involve num
'(%sin
))))
1550 (cond (scnum (cond ((trig-right-placep '%sin scnum
) t
)
1551 (t (psetq num
(m^ denom -
1)
1552 denom
(m^ num -
1)) t
)))
1554 ((let ((scden (involve denom
'(%sin
))))
1555 (cond (scden (cond ((trig-right-placep '%sin scden
) t
)
1556 (t (psetq num
(m^ denom -
1)
1557 denom
(m^ num -
1)) t
)))
1559 ((let ((scnum (involve num
'(%asin %acos %atan
))))
1560 ;; If the numerator contains an inverse trig and the
1561 ;; denominator or reciprocal of denominator is polynomial,
1562 ;; leave everything as is. If the inverse trig is moved to
1563 ;; the denominator, things get messy, even if the numerator
1564 ;; becomes a polynomial. This is not perfect.
1565 (cond ((and scnum
(or (polyinx denom var
())
1566 (polyinx (m^ denom -
1) var
())))
1569 ((or (oscip num
) (oscip denom
)))
1571 (psetq num
(m^ denom -
1) denom
(m^ num -
1))))
1574 ;;i don't know what to do here for some cases, may have to be refined.
1575 (defun num-of-logs (exp)
1576 (cond ((mapatom exp
) 0)
1577 ((equal (caar exp
) '%log
)
1578 (m+ 1 (num-of-log-l (cdr exp
))))
1579 ((and (mexptp exp
) (mnump (caddr exp
)))
1580 (m* (simplify `((mabs) ,(caddr exp
)))
1581 (num-of-logs (cadr exp
))))
1582 (t (num-of-log-l (cdr exp
)))))
1584 (defun num-of-log-l (llist)
1585 (do ((temp llist
(cdr temp
)) (ans 0))
1587 (setq ans
(m+ ans
(num-of-logs (car temp
))))))
1589 (defun %e-right-placep
(%e-arg
)
1590 (let ((%e-arg-diff
(sdiff %e-arg var
)))
1592 ((free %e-arg-diff var
)) ;simple cases
1593 ((or (and (mexptp denom
)
1594 (equal (cadr denom
) -
1))
1595 (polyinx (m^ denom -
1) var
())) ())
1596 ((let ((%e-arg-diff-lim
(ridofab (limit %e-arg-diff var val
'think
)))
1597 (%e-arg-exp-lim
(ridofab (limit (m^
'$%e %e-arg
) var val
'think
))))
1600 (format t
"%e-arg-dif-lim = ~A~%" %e-arg-diff-lim
)
1601 (format t
"%e-arg-exp-lim = ~A~%" %e-arg-exp-lim
))
1602 (cond ((equal %e-arg-diff-lim %e-arg-exp-lim
)
1604 ((and (mnump %e-arg-diff-lim
) (mnump %e-arg-exp-lim
))
1606 ((and (mnump %e-arg-diff-lim
) (infinityp %e-arg-exp-lim
))
1607 ;; This is meant to make maxima handle bug 1469411
1608 ;; correctly. Undoubtedly, this needs work.
1612 (defun trig-right-placep (trig-type arg
)
1613 (let ((arglim (ridofab (limit arg var val
'think
)))
1614 (triglim (ridofab (limit `((,trig-type
) ,arg
) var val
'think
))))
1615 (cond ((and (equal arglim
0) (equal triglim
0)) t
)
1616 ((and (infinityp arglim
) (infinityp triglim
)) t
)
1619 ;;Takes a numerator and a denominator. If they tries all combinations of
1620 ;;products to try and make a simpler set of subproblems for LHOSPITAL.
1621 (defun remove-singularities (numer denom
)
1622 (cond ((or (null numer
) (null denom
)
1623 (atom numer
) (atom denom
)
1624 (not (mtimesp numer
)) ;Leave this here for a while.
1625 (not (mtimesp denom
)))
1626 (values 1 numer denom
))
1629 (multiple-value-bind (num-consts num-vars
)
1630 (var-or-const numer
)
1631 (multiple-value-bind (denom-consts denom-vars
)
1632 (var-or-const denom
)
1633 (if (not (mtimesp num-vars
))
1634 (setq num-vars
(list num-vars
))
1635 (setq num-vars
(cdr num-vars
)))
1636 (if (not (mtimesp denom-vars
))
1637 (setq denom-vars
(list denom-vars
))
1638 (setq denom-vars
(cdr denom-vars
)))
1639 (do ((nl num-vars
(cdr nl
))
1640 (num-list (copy-list num-vars
))
1641 (den-list denom-vars den-list-temp
)
1642 (den-list-temp (copy-list denom-vars
)))
1643 ((null nl
) (values (m* const
(m// num-consts denom-consts
))
1645 (m*l den-list-temp
)))
1646 (do ((dl den-list
(cdr dl
)))
1648 (if (or (%einvolve
(car nl
)) (%einvolve
(car nl
)))
1650 (let ((lim (catch 'limit
(simpinf (simpab (limit (m// (car nl
) (car dl
))
1651 var val
'think
))))))
1652 (cond ((or (eq lim t
)
1654 (equal (ridofab lim
) 0)
1656 (not (free lim
'$inf
))
1657 (not (free lim
'$minf
))
1658 (not (free lim
'$infinity
))
1659 (not (free lim
'$ind
))
1660 (not (free lim
'$und
)))
1663 (setq const
(m* lim const
))
1664 (setq num-list
(delete (car nl
) num-list
:count
1 :test
#'equal
))
1665 (setq den-list-temp
(delete (car dl
) den-list-temp
:count
1 :test
#'equal
))
1666 (return t
)))))))))))))
1668 ;; separate terms that contain var from constant terms
1669 ;; returns (const-terms . var-terms)
1670 (defun var-or-const (expr)
1671 (setq expr
($factor expr
))
1679 (do ((l (cdr expr
) (cdr l
))
1682 ((null l
) (values const varl
))
1683 (if (free (car l
) var
)
1684 (setq const
(m* (car l
) const
))
1685 (setq varl
(m* (car l
) varl
)))))
1689 ;; if term goes to non-zero constant, replace with constant
1690 (defun lhsimp (term var val
)
1691 (cond ((atom term
) term
)
1693 (let ((term-value (ridofab (limit term var val
'think
))))
1694 (cond ((not (member term-value
1695 '($inf $minf $und $ind $infinity
0)))
1699 (defun bylog (expo bas
)
1702 (try-lhospital-quit (simplify `((%log
) ,(tansc bas
)))
1707 (defun simplimexpt (bas expo bl el
)
1708 (cond ((or (eq bl
'$und
) (eq el
'$und
)) '$und
)
1710 (cond ((eq el
'$inf
) (if (eq bl
'$zeroa
) bl
0))
1711 ((eq el
'$minf
) (if (eq bl
'$zeroa
) '$inf
'$infinity
))
1712 ((eq el
'$ind
) '$ind
)
1713 ((eq el
'$infinity
) '$und
)
1714 ((zerop2 el
) (bylog expo bas
))
1715 (t (cond ((equal (getsignl el
) -
1)
1716 (cond ((eq bl
'$zeroa
) '$inf
)
1718 (cond ((even1 el
) '$inf
)
1719 ((eq (ask-integer el
'$integer
) '$yes
)
1720 (if (eq (ask-integer el
'$even
) '$yes
)
1722 '$minf
)))) ;Gotta be ODD.
1723 (t (setq bas
(behavior bas var val
))
1724 (cond ((equal bas
1) '$inf
)
1725 ((equal bas -
1) '$minf
)
1726 (t (throw 'limit t
))))))
1728 (member bl
'($zeroa $zerob
) :test
#'eq
))
1729 (cond ((even1 el
) '$zeroa
)
1730 ((and (eq bl
'$zerob
)
1732 (evenp (caddr el
))) 0)
1734 ((and (equal (getsignl el
) 1)
1735 (eq bl
'$zeroa
)) bl
)
1736 ((equal (getsignl el
) 0)
1740 (cond ((zerop2 el
) (bylog expo bas
))
1742 ((eq el
'$inf
) '$infinity
)
1743 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1744 ((equal (setq el
(getsignl el
)) 1) '$infinity
)
1747 (t (throw 'limit t
))))
1749 (cond ((eq el
'$inf
) '$inf
)
1750 ((equal el
'$minf
) 0)
1751 ((zerop2 el
) (bylog expo bas
))
1752 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1753 (t (cond ((eql 0 (getsignl el
)) 1)
1754 ((ratgreaterp 0 el
) '$zeroa
)
1755 ((ratgreaterp el
0) '$inf
)
1756 (t (throw 'limit t
))))))
1758 (cond ((zerop2 el
) (bylog expo bas
))
1759 ((eq el
'$inf
) '$und
)
1760 ((equal el
'$minf
) 0)
1761 ;;;Why not generalize this. We can ask about the number. -Jim 2/23/81
1762 ((mnump el
) (cond ((mnegp el
)
1765 (if (eq (ask-integer el
'$integer
) '$yes
)
1766 (if (eq (ask-integer el
'$even
) '$yes
)
1770 (t (cond ((even1 el
) '$inf
)
1771 ((eq (ask-integer el
'$integer
) '$yes
)
1772 (if (eq (ask-integer el
'$even
) '$yes
)
1776 (loginprod?
(throw 'lip?
'lip
!))
1778 ((equal (simplify (ratdisrep (ridofab bl
))) 1)
1779 (if (infinityp el
) (bylog expo bas
) 1))
1780 ((and (equal (ridofab bl
) -
1)
1781 (infinityp el
)) '$ind
) ;LB
1782 ((eq bl
'$ind
) (cond ((or (zerop2 el
) (infinityp el
)) '$und
)
1783 ((not (equal (getsignl el
) -
1)) '$ind
)
1785 ((eq el
'$inf
) (cond ((abeq1 bl
)
1786 (if (equal (getsignl bl
) 1) 1 '$ind
))
1788 (if (equal (getsignl bl
) 1) '$zeroa
0))
1789 ((equal (getsignl (m1- bl
)) 1) '$inf
)
1790 ((equal (getsignl (m1- `((mabs) ,bl
))) 1) '$infinity
)
1791 (t (throw 'limit t
))))
1792 ((eq el
'$minf
) (cond ((abeq1 bl
)
1793 (if (equal (getsignl bl
) 1) 1 '$ind
))
1795 (if (equal (getsignl bl
) 1) '$zeroa
0))
1796 ((ratgreaterp 0 bl
) '$infinity
)
1799 (if (equal val
'$infinity
)
1800 '$und
;Not enough info to do anything.
1801 (destructuring-bind (real-el . imag-el
)
1803 (setq real-el
(limit real-el var origval nil
))
1804 (cond ((eq real-el
'$minf
)
1806 ((and (eq real-el
'$inf
)
1807 (not (equal (ridofab (limit imag-el var origval nil
)) 0)))
1809 ((eq real-el
'$infinity
)
1810 (throw 'limit t
)) ;; don't really know real component
1814 ((eq el
'$ind
) '$ind
)
1819 (cond ((numberp x
) (and (integerp x
) (evenp x
)))
1820 ((and (mnump x
) (evenp (cadr x
))))))
1822 ;; is absolute value less than one?
1826 (and (ratgreaterp 1. bl
) (ratgreaterp bl -
1.
)))
1827 (t (equal (getsignl (m1- `((mabs) ,bl
))) -
1.
))))
1829 ;; is absolute value equal to one?
1833 (or (equal 1. bl
) (equal bl -
1.
)))
1834 (t (equal (getsignl (m1- `((mabs) ,bl
))) 0))))
1836 (defun simplimit (exp var val
&aux op
)
1839 ((or (atom exp
) (mnump exp
)) exp
)
1840 ;; Lookup and dispatch a simplim%function from the property list
1841 ((setq op
(safe-get (mop exp
) 'simplim%function
))
1842 (funcall op exp var val
))
1844 ;; And do the same for subscripted:
1845 ((and (consp exp
) (consp (car exp
)) (eq (caar exp
) 'mqapply
)
1846 (setq op
(safe-get (subfunname exp
) 'simplim%function
)))
1847 (funcall op exp var val
))
1849 ;; Without the call to rootscontract,
1850 ;; limit(((-4)*x^2-10*x+24)/((4*x+8)^(1/3)+2),x,-4)
1851 ;; returns zero (should be 66). This bug is due to the fact that
1852 ;; 4^(1/3)*2^(1/3)-2 does not simplify to zero. Although the call
1853 ;; to rootscontract takes care of this case, almost surely there are
1854 ;; many other limit problems that need more than rootscontract.
1855 ((mplusp exp
) (let (($rootsconmode nil
)) ($rootscontract
(simplimplus exp
))))
1856 ((mtimesp exp
) (simplimtimes (cdr exp
)))
1857 ((mexptp exp
) (simplimexpt (cadr exp
) (caddr exp
)
1858 (limit (cadr exp
) var val
'think
)
1859 (limit (caddr exp
) var val
'think
)))
1860 ((member (caar exp
) '(%sin %cos
) :test
#'eq
)
1861 (simplimsc exp
(caar exp
) (limit (cadr exp
) var val
'think
)))
1862 ((eq (caar exp
) '%tan
) (simplim%tan
(cadr exp
)))
1863 ((member (caar exp
) '(%sinh %cosh
) :test
#'eq
)
1864 (simplimsch (caar exp
) (limit (cadr exp
) var val
'think
)))
1865 ((member (caar exp
) '(%erf %tanh
) :test
#'eq
)
1866 (simplim%erf-%tanh
(caar exp
) (cadr exp
)))
1867 ((eq (caar exp
) '%atanh
)
1868 (simplim%atanh
(limit (cadr exp
) var val
'think
) val
))
1869 ((eq (caar exp
) '%acosh
)
1870 (simplim%acosh
(limit (cadr exp
) var val
'think
)))
1871 ((eq (caar exp
) '%asinh
)
1872 (simplim%asinh
(limit (cadr exp
) var val
'think
)))
1873 ((eq (caar exp
) '%inverse_jacobi_ns
)
1874 (simplim%inverse_jacobi_ns
(limit (cadr exp
) var val
'think
) (third exp
)))
1875 ((eq (caar exp
) '%inverse_jacobi_nc
)
1876 (simplim%inverse_jacobi_nc
(limit (cadr exp
) var val
'think
) (third exp
)))
1877 ((eq (caar exp
) '%inverse_jacobi_sc
)
1878 (simplim%inverse_jacobi_sc
(limit (cadr exp
) var val
'think
) (third exp
)))
1879 ((eq (caar exp
) '%inverse_jacobi_cs
)
1880 (simplim%inverse_jacobi_cs
(limit (cadr exp
) var val
'think
) (third exp
)))
1881 ((eq (caar exp
) '%inverse_jacobi_dc
)
1882 (simplim%inverse_jacobi_dc
(limit (cadr exp
) var val
'think
) (third exp
)))
1883 ((eq (caar exp
) '%inverse_jacobi_ds
)
1884 (simplim%inverse_jacobi_ds
(limit (cadr exp
) var val
'think
) (third exp
)))
1885 ((and (eq (caar exp
) 'mqapply
)
1886 (eq (subfunname exp
) '$psi
))
1887 (simplim$psi
(subfunsubs exp
) (subfunargs exp
) val
))
1888 ((and (eq (caar exp
) var
)
1889 (member 'array
(car exp
) :test
#'eq
)
1890 (every #'(lambda (sub-exp)
1893 exp
) ;LIMIT(B[I],B,INF); -> B[I]
1894 ;; When limsubst is true, limit(f(n+1)/f(n),n,inf) = 1. The user documentation
1895 ;; warns against setting limsubst to true.
1897 (simplify (cons (operator-with-array-flag exp
)
1898 (mapcar #'(lambda (a)
1899 (limit a var val
'think
))
1901 (throw 'limit t
)))))
1904 (setq e
(resimplify (subst (m// 1 var
) var e
)))
1905 (let ((new-val (cond ((eq val
'$zeroa
) '$inf
)
1906 ((eq val
'$zerob
) '$minf
))))
1907 (if new-val
(let ((preserve-direction t
))
1908 (toplevel-$limit e var new-val
)) (throw 'limit t
))))
1910 (defun simplimtimes (exp)
1911 ;; The following test
1912 ;; handles (-1)^x * 2^x => (-2)^x => $infinity
1913 ;; wants to avoid (-1)^x * 2^x => $ind * $inf => $und
1915 (and (expfactorp (cons '(mtimes) exp
) 1)
1916 (expfactor (cons '(mtimes) exp
) 1 var
))))
1917 (when try
(return-from simplimtimes try
)))
1919 (let ((prod 1) (num 1) (denom 1)
1920 (zf nil
) (ind-flag nil
) (inf-type nil
)
1921 (constant-zero nil
) (constant-infty nil
))
1923 (let* ((loginprod?
(involve term
'(%log
)))
1924 (y (catch 'lip?
(limit term var val
'think
))))
1926 ;; limit failed due to log in product
1928 (return-from simplimtimes
(liminv (cons '(mtimes simp
) exp
))))
1930 ;; If the limit is infinitesimal or zero
1932 (setf num
(m* num term
)
1933 constant-zero
(or constant-zero
(not (among var term
))))
1936 (unless zf
(setf zf
1)))
1938 (setf zf
(* -
1 (or zf
1))))))
1940 ;; If the limit is not some form of infinity or
1941 ;; undefined/indeterminate.
1942 ((not (member y
'($inf $minf $infinity $ind $und
) :test
#'eq
))
1943 (setq prod
(m* prod y
)))
1945 ((eq y
'$und
) (return-from simplimtimes
'$und
))
1946 ((eq y
'$ind
) (setq ind-flag t
))
1948 ;; Some form of infinity
1950 (setf denom
(m* denom term
)
1951 constant-infty
(or constant-infty
(not (among var term
))))
1952 (unless (eq inf-type
'$infinity
)
1954 ((eq y
'$infinity
) (setq inf-type
'$infinity
))
1955 ((null inf-type
) (setf inf-type y
))
1956 ;; minf * minf or inf * inf
1957 ((eq y inf-type
) (setf inf-type
'$inf
))
1959 (t (setf inf-type
'$minf
))))))))
1962 ;; If there are zeros and infinities among the terms that are free of
1963 ;; VAR, then we have an expression like "inf * zeroa * f(x)" or
1964 ;; similar. That gives an undefined result. Note that we don't
1965 ;; necessarily have something undefined if only the zeros have a term
1966 ;; free of VAR. For example "zeroa * exp(-1/x) * 1/x" as x -> 0. And
1967 ;; similarly for the infinities.
1968 ((and constant-zero constant-infty
) '$und
)
1970 ;; If num=denom=1, we didn't find any explicit infinities or zeros, so we
1971 ;; either return the simplified product or ind
1972 ((and (eql num
1) (eql denom
1))
1973 (if ind-flag
'$ind prod
))
1974 ;; If denom=1 (and so num != 1), we have some form of zero
1978 (let ((sign (getsignl prod
)))
1979 (if (or (not sign
) (eq sign
'complex
))
1985 ;; If num=1 (and so denom != 1), we have some form of infinity
1987 (let ((sign ($csign prod
)))
1990 ((eq sign
'$pos
) inf-type
)
1991 ((eq sign
'$neg
) (case inf-type
1995 ((member sign
'($complex $imaginary
)) '$infinity
)
1996 ; sign is '$zero, $pnz, $pz, etc
1997 (t (throw 'limit t
)))))
1998 ;; Both zeros and infinities
2000 ;; All bets off if there are some infinities or some zeros, but it
2001 ;; needn't be undefined (see above)
2002 (when (or constant-zero constant-infty
) (throw 'limit t
))
2004 (let ((ans (limit2 num
(m^ denom -
1) var val
)))
2006 (simplimtimes (list prod ans
))
2007 (throw 'limit t
)))))))
2009 ;;;PUT CODE HERE TO ELIMINATE FAKE SINGULARITIES??
2011 (defun simplimplus (exp)
2012 (cond ((memalike exp simplimplus-problems
)
2015 (progn (push exp simplimplus-problems
)
2016 (let ((ans (catch 'limit
(simplimplus1 exp
))))
2017 (cond ((or (eq ans
())
2019 (among '%limit ans
))
2020 (let ((new-exp (sratsimp exp
)))
2021 (cond ((not (alike1 exp new-exp
))
2023 (limit new-exp var val
'think
))))
2024 (cond ((or (eq ans
())
2029 (pop simplimplus-problems
)))))
2031 (defun simplimplus1 (exp)
2032 (prog (sum y infl infinityl minfl indl undl
)
2034 (do ((exp (cdr exp
) (cdr exp
)) (f))
2035 ((or y
(null exp
)) nil
)
2036 (setq f
(limit (car exp
) var val
'think
))
2039 ((eq f
'$und
) (push (car exp
) undl
))
2040 ((not (member f
'($inf $minf $infinity $ind
) :test
#'eq
))
2041 (setq sum
(m+ sum f
)))
2042 ((eq f
'$ind
) (push (car exp
) indl
))
2043 (infinityl (throw 'limit t
))
2044 ;;;Don't know what to do with an '$infinity and an $inf or $minf
2045 ((eq f
'$inf
) (push (car exp
) infl
))
2046 ((eq f
'$minf
) (push (car exp
) minfl
))
2048 (cond ((or infl minfl
)
2050 (t (push (car exp
) infinityl
))))))
2052 (if (or infl minfl indl infinityl
)
2053 (setq infinityl
(append undl infinityl
)); x^2 + x*sin(x)
2054 (return '$und
))) ; 1 + x*sin(x)
2055 ((not (or infl minfl indl infinityl
))
2056 (return (cond ((atom sum
) sum
)
2057 ((or (not (free sum
'$zeroa
))
2058 (not (free sum
'$zerob
)))
2061 (t (cond ((null infinityl
)
2062 (cond (infl (cond ((null minfl
) (return '$inf
))
2064 (minfl (return '$minf
))
2065 ((> (length indl
) 1)
2066 ;; At this point we have a sum of '$ind. We factor
2067 ;; the sum and try again. This way we get the limit
2068 ;; of expressions like (a-b)*ind, where (a-b)--> 0.
2069 (cond ((not (alike1 (setq y
($factorsum exp
)) exp
))
2070 (return (limit y var val
'think
)))
2073 (t (return '$ind
))))
2074 (t (setq infl
(append infl infinityl
))))))
2076 oon
(setq y
(m+l
(append minfl infl
)))
2077 (cond ((alike1 exp
(setq y
(sratsimp (hyperex y
))))
2078 (cond ((not (infinityp val
))
2079 (setq infl
(cnv infl val
)) ;THIS IS HORRIBLE!!!!
2080 (setq minfl
(cnv minfl val
))))
2082 (cond ((every #'(lambda (j) (radicalp j var
))
2083 (append infl minfl
))
2084 (setq y
(rheur infl minfl
)))
2085 (t (setq y
(sheur infl minfl
))))))
2086 (t (setq y
(limit y var val
'think
))))
2087 (cond ((or (eq y
())
2088 (eq y t
)) (return ()))
2089 ((infinityp y
) (return y
))
2090 (t (return (m+ sum y
))))))
2092 ;; Limit n/d, using heuristics on the order of growth.
2095 (cond ((and (free n var
)
2098 (t (setq n
(cpa n d nil
))
2100 (cond ((oscip orig-n
) '$und
)
2102 ((equal n -
1) '$zeroa
)
2103 ((equal n
0) (m// orig-n d
)))))))
2106 ;;;L1 is a list of INF's and L2 is a list of MINF's. Added together
2107 ;;;it is indeterminate.
2108 (defun sheur (l1 l2
)
2109 (let ((term (sheur1 l1 l2
)))
2110 (cond ((equal term
'$inf
) '$inf
)
2111 ((equal term
'$minf
) '$minf
)
2112 (t (let ((new-num (m+l
(mapcar #'(lambda (num-term)
2113 (m// num-term
(car l1
)))
2115 (cond ((limit2 new-num
(m// 1 (car l1
)) var val
))))))))
2118 (cond ((atom exp
) nil
)
2119 (t (setq exp
(nformat exp
))
2120 (cond ((and (eq (caar exp
) 'mquotient
)
2121 (among var
(caddr exp
)))
2124 (defun zerop2 (z) (=0 (ridofab z
)))
2126 (defun raise (a) (m+ a
'$zeroa
))
2128 (defun lower (a) (m+ a
'$zerob
))
2130 (defun sincoshk (exp1 l sc
)
2131 (cond ((equal l
1) (lower l
))
2132 ((equal l -
1) (raise l
))
2134 ((member val
'($zeroa $zerob
) :test
#'eq
) (spangside exp1 l
))
2137 (defun spangside (e l
)
2138 (setq e
(behavior e var val
))
2139 (cond ((equal e
1) (raise l
))
2140 ((equal e -
1) (lower l
))
2143 ;; get rid of zeroa and zerob
2145 (if (among '$zeroa e
) (setq e
(maxima-substitute 0 '$zeroa e
)))
2146 (if (among '$zerob e
) (setq e
(maxima-substitute 0 '$zerob e
)))
2150 ;; returns true if exp is a polynomial raised to a numeric power
2151 (defun simplerd (exp)
2153 (mnump (caddr exp
)) ;; exponent must be a number - no variables
2154 (polyp (cadr exp
))))
2156 (defun branch1 (exp val
)
2157 (cond ((polyp exp
) nil
)
2158 ((simplerd exp
) (zerop2 (subin val
(cadr exp
))))
2160 (loop for v on
(cdr exp
)
2161 when
(branch1 (car v
) val
)
2164 (defun branch (exp val
)
2165 (cond ((polyp exp
) nil
)
2166 ((or (simplerd exp
) (mtimesp exp
))
2169 (every #'(lambda (j) (branch j val
)) (the list
(cdr exp
))))))
2171 (defun ser0 (e n d val
)
2172 (cond ((and (branch n val
) (branch d val
))
2176 ;;;NN* gets set by POFX, called by SER1, to get a list of exponents.
2177 (setq nn
* (ratmin nn
*))
2178 (setq n
(sratsimp (m^ n
(m^ var nn
*))))
2179 (setq d
(sratsimp (m^ d
(m^ var nn
*))))
2180 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2183 (t (try-lhospital-quit n d nil
))))
2185 (defun rheur (l1 l2
)
2187 (setq m1
(mapcar (function asymredu
) l1
))
2188 (setq m2
(mapcar (function asymredu
) l2
))
2189 (setq ans
(m+l
(append m1 m2
)))
2190 (cond ((rptrouble (m+l
(append l1 l2
)))
2191 (return (limit (simplify (rdsget (m+l
(append l1 l2
))))
2195 ((mplusp ans
) (return (sheur m1 m2
)))
2196 (t (return (limit ans var val t
))))))
2198 (defun rptrouble (rp)
2199 (not (equal (rddeg rp nil
) (rddeg (asymredu rp
) nil
))))
2201 (defun radicalp (exp var
)
2202 (cond ((polyinx exp var
()))
2203 ((mexptp exp
) (cond ((equal (caddr exp
) -
1.
)
2204 (radicalp (cadr exp
) var
))
2206 ((member (caar exp
) '(mplus mtimes
) :test
#'eq
)
2207 (every #'(lambda (j) (radicalp j var
))
2210 (defun involve (e nn
*)
2211 (declare (special var
))
2212 (cond ((atom e
) nil
)
2214 ((and (member (caar e
) nn
* :test
#'eq
) (among var
(cdr e
))) (cadr e
))
2215 (t (some #'(lambda (j) (involve j nn
*)) (cdr e
)))))
2217 (defun notinvolve (exp nn
*)
2218 (cond ((atom exp
) t
)
2220 ((member (caar exp
) nn
* :test
#'eq
) (not (among var
(cdr exp
))))
2221 ((every #'(lambda (j) (notinvolve j nn
*))
2224 (defun sheur1 (l1 l2
)
2226 (setq l1
(mapcar #'stirling0 l1
))
2227 (setq l2
(mapcar #'stirling0 l2
))
2228 (setq l1
(m+l
(maxi l1
)))
2229 (setq l2
(m+l
(maxi l2
)))
2230 (setq ans
(cpa l1 l2 t
))
2231 (return (cond ((=0 ans
) (m+ l1 l2
))
2232 ((equal ans
1.
) '$inf
)
2235 (defun zero-lim (cpa-list)
2236 (do ((l cpa-list
(cdr l
)))
2238 (and (eq (caar l
) 'gen
)
2239 (zerop2 (limit (cadar l
) var val
'think
))
2242 ;; Compare order of growth for R1 and R2. The result is 0, -1, +1
2243 ;; depending on the relative order of growth. 0 is returned if R1 and
2244 ;; R2 have the same growth; -1 if R1 grows much more slowly than R2;
2245 ;; +1 if R1 grows much more quickly than R2.
2246 (defun cpa (r1 r2 flag
)
2249 (cond ((alike1 t1 t2
) 0.
)
2251 (cond ((free t2 var
) 0.
)
2252 (t (let ((lim-ans (limit1 t2 var val
)))
2253 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2256 (let ((lim-ans (limit1 t1 var val
)))
2257 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2260 ;; Make T1 and T2 be a list of terms that are multiplied
2262 (cond ((mtimesp t1
) (setq t1
(cdr t1
)))
2263 (t (setq t1
(list t1
))))
2264 (cond ((mtimesp t2
) (setq t2
(cdr t2
)))
2265 (t (setq t2
(list t2
))))
2266 ;; Find the strengths of each term of T1 and T2
2267 (setq t1
(mapcar (function istrength
) t1
))
2268 (setq t2
(mapcar (function istrength
) t2
))
2269 ;; Compute the max of the strengths of the terms.
2270 (let ((ans (ismax t1
))
2272 (cond ((or (null ans
) (null d
))
2273 ;;(eq (car ans) 'gen) (eq (car d) 'gen))
2274 ;; ismax couldn't find highest term; give up
2277 (if (eq (car ans
) 'var
) (setq ans
(add-up-deg t1
)))
2278 (if (eq (car d
) 'var
) (setq d
(add-up-deg t2
)))
2279 ;; Can't just just compare dominating terms if there are
2280 ;; indeterm-inates present; e.g. X-X^2*LOG(1+1/X). So
2282 (cond ((or (zero-lim t1
)
2284 (cpa-indeterm ans d t1 t2 flag
))
2285 ((isgreaterp ans d
) 1.
)
2286 ((isgreaterp d ans
) -
1.
)
2289 (defun cpa-indeterm (ans d t1 t2 flag
)
2290 (cond ((not (eq (car ans
) 'var
))
2291 (setq ans
(gather ans t1
) d
(gather d t2
))))
2292 (let ((*indicator
(and (eq (car ans
) 'exp
)
2295 (setq test
(cpa1 ans d
))
2296 (cond ((and (zerop1 test
)
2297 (or (equal ($radcan
(m// (cadr ans
) (cadr d
))) 1.
)
2298 (and (polyp (cadr ans
))
2300 (equal (limit (m// (cadr ans
) (cadr d
)) var val
'think
)
2302 (let ((new-term1 (m// t1
(cadr ans
)))
2303 (new-term2 (m// t2
(cadr d
))))
2304 (cpa new-term1 new-term2 flag
)))
2307 (defun add-up-deg (strengthl)
2308 (do ((stl strengthl
(cdr stl
))
2311 ((null stl
) (list 'var
(m*l poxl
) (m+l degl
)))
2312 (cond ((eq (caar stl
) 'var
)
2313 (push (cadar stl
) poxl
)
2314 (push (caddar stl
) degl
)))))
2318 (cond ((eq (car p1
) 'gen
) (return 0.
)))
2319 (setq flag
(car p1
))
2324 (setq s1
(istrength p1
))
2325 (setq s2
(istrength p2
))
2328 ((isgreaterp s1 s2
) 1.
)
2329 ((isgreaterp s2 s1
) -
1.
)
2331 (setq *indicator nil
)
2333 ((and (poly? p1 var
) (poly? p2 var
))
2334 (setq p1
(m- p1 p2
))
2335 (cond ((zerop1 p1
) 0.
)
2336 (t (getsignl (hot-coef p1
)))))
2340 (list (m*t -
1 p2
))))
2341 (cond ((zerop2 s1
) 0.
)
2342 ((ratgreaterp s1
0.
) 1.
)
2346 (setq p1
(caddr p1
))
2347 (setq p2
(caddr p2
))
2348 (cond ((and (poly? p1 var
) (poly? p2 var
))
2349 (setq p1
(m- p1 p2
))
2350 (return (cond ((or (zerop1 p1
)
2351 (not (among var p1
)))
2353 (t (getsignl (hot-coef p1
))))))
2354 ((and (radicalp p1 var
) (radicalp p2 var
))
2357 (list (m*t -
1 p2
))))
2358 (return (cond ((eq s1
'$inf
) 1.
)
2359 ((eq s1
'$minf
) -
1.
)
2361 (cond ((ratgreaterp s1
0.
) 1.
)
2362 ((ratgreaterp 0. s1
) -
1.
)
2365 (t (return (cpa p1 p2 t
)))))
2367 (setq p1
(try-lhospital (asymredu p1
) (asymredu p2
) nil
))
2368 (return (cond ((zerop2 p1
) -
1.
)
2369 ((real-infinityp p1
) 1.
)
2372 ;;;EXPRESSIONS TO ISGREATERP ARE OF THE FOLLOWING FORMS
2373 ;;; ("VAR" POLY DEG)
2375 ;;; ("LOG" LOG(EXP))
2376 ;;; ("FACT" <A FACTORIAL EXPRESSION>)
2377 ;;; ("GEN" <ANY OTHER TYPE OF EXPRESSION>)
2379 (defun isgreaterp (a b
)
2382 (cond ((or (eq ta
'gen
)
2384 ((and (eq ta tb
) (eq ta
'var
))
2385 (ratgreaterp (caddr a
) (caddr b
)))
2386 ((and (eq ta tb
) (eq ta
'exp
))
2387 ;; Both are exponential order of infinity. Check the
2388 ;; exponents to determine which exponent is bigger.
2389 (eq (limit (m- `((%log
) ,(second a
)) `((%log
) ,(second b
)))
2392 ((member ta
(cdr (member tb
'(num log var exp fact gen
) :test
#'eq
)) :test
#'eq
)))))
2395 ;; Preprocess the list of products. Separate the terms that
2396 ;; exponentials and those that don't. Actually multiply the
2397 ;; exponential terms together to form a single term. Pass this and
2398 ;; the rest to ismax-core to find the max.
2399 (let (exp-terms non-exp-terms
)
2401 (if (eq 'exp
(car term
))
2402 (push term exp-terms
)
2403 (push term non-exp-terms
)))
2404 ;; Multiply the exp-terms together
2407 ;;(format t "exp-terms = ~A~%" exp-terms)
2408 (dolist (term exp-terms
)
2409 (setf product
(simplify (mul product
(second term
)))))
2410 ;;(format t "product = ~A~%" product)
2411 (setf product
`(exp ,($logcontract product
)))
2412 ;;(format t "product = ~A~%" product)
2413 (ismax-core (cons product non-exp-terms
)))
2416 (defun ismax-core (l)
2419 ((= (length l
) 1) (car l
)) ;If there is only 1 thing give it back.
2420 ((every #'(lambda (x)
2421 (not (eq (car x
) 'gen
))) l
)
2423 (do ((l1 (cdr l
) (cdr l1
))
2427 (cond ((isgreaterp temp-ans
(car l1
))
2428 (setq ans temp-ans
))
2429 ((isgreaterp (car l1
) temp-ans
)
2430 (setq temp-ans
(car l1
))
2431 (setq ans temp-ans
))
2432 (t (setq ans
())))))
2435 ;RETURNS LIST OF HIGH TERMS
2437 (cond ((atom all
) nil
)
2438 (t (do ((l (cdr all
) (cdr l
))
2440 (total 1) ; running total constant factor of hi-term
2441 (hi-terms (ncons (car all
)))
2443 ((null l
) (if (zerop2 total
) ; if high-order terms cancel each other
2444 all
; keep everything
2445 hi-terms
)) ; otherwise return list of high terms
2446 (setq compare
(limit (m// (car l
) hi-term
) var val
'think
))
2448 ((or (infinityp compare
)
2449 (and (eq compare
'$und
)
2450 (zerop2 (limit (m// hi-term
(car l
)) var val
'think
))))
2451 (setq total
1) ; have found new high term
2452 (setq hi-terms
(ncons (setq hi-term
(car l
)))))
2453 ((zerop2 compare
) nil
)
2454 ;; COMPARE IS IND, FINITE-VALUED, or und in both directions
2455 (t ; add to list of high terms
2456 (setq total
(m+ total compare
))
2457 (setq hi-terms
(append hi-terms
(ncons (car l
))))))))))
2461 (cond ((atom l
) (return nil
)))
2462 l1
(setq ans
(car l
))
2464 (cond ((null l
) (return ans
))
2465 ((ratgreaterp ans
(car l
)) (go l2
))
2470 (cond ((atom l
) (return nil
)))
2471 l1
(setq ans
(car l
))
2473 (cond ((null l
) (return ans
))
2474 ((ratgreaterp (car l
) ans
) (go l2
))
2482 ((or (mnump e
) (not (among var e
))) nil
)
2483 ((and (mexptp e
) (eq (cadr e
) var
))
2484 (push (caddr e
) nn
*))
2486 (t (mapc #'pofx
(cdr e
)))))
2489 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2490 (t (setq e
(subin (m+ var val
) e
))))
2492 (cond ((pofx e
) e
)))
2494 (defun gather (ind l
)
2496 (setq ind
(car ind
))
2497 loop
(cond ((null l
)
2498 (return (list ind
(m*l ans
))))
2499 ((equal (caar l
) ind
)
2500 (push (cadar l
) ans
)))
2504 ; returns rough class-of-growth of term
2505 (defun istrength (term)
2506 (cond ((mnump term
) (list 'num term
))
2507 ((atom term
) (cond ((eq term var
)
2509 (t (list 'num term
))))
2510 ((not (among var term
)) (list 'num term
))
2512 (let ((temp (ismax-core (mapcar #'istrength
(cdr term
)))))
2513 (cond ((not (null temp
)) temp
)
2516 (let ((temp (mapcar #'istrength
(cdr term
)))
2518 (setq temp1
(ismax temp
))
2519 (cond ((null temp1
) `(gen ,term
))
2520 ((eq (car temp1
) 'log
) `(log ,temp
))
2521 ((eq (car temp1
) 'var
) (add-up-deg temp
))
2524 (real-infinityp (limit term var val t
)))
2525 (let ((logterm (logred term
)))
2526 (cond ((and (among var
(caddr term
))
2527 (member (car (istrength logterm
))
2528 '(var exp fact
) :test
#'eq
)
2529 (real-infinityp (limit logterm var val t
)))
2530 (list 'exp
(m^
'$%e logterm
)))
2531 ((not (among var
(caddr term
)))
2532 (let ((temp (istrength (cadr term
))))
2533 (cond ((not (alike1 temp term
))
2534 (rplaca (cdr temp
) term
)
2535 (and (eq (car temp
) 'var
)
2537 (m* (caddr temp
) (caddr term
))))
2541 ((and (eq (caar term
) '%log
)
2542 (real-infinityp (limit term var val t
)))
2543 (let ((stren (istrength (cadr term
))))
2544 (cond ((member (car stren
) '(log var
) :test
#'eq
)
2546 ((and (eq (car stren
) 'exp
)
2547 (eq (caar (second stren
)) 'mexpt
))
2548 (istrength (logred (second stren
))))
2550 ((eq (caar term
) 'mfactorial
)
2552 ((let ((temp (hyperex term
)))
2553 (and (not (alike1 term temp
))
2555 (t (list 'gen term
))))
2557 ;; log reduce - returns log of s1
2559 (or (and (eq (cadr s1
) '$%e
) (caddr s1
))
2560 (m* (caddr s1
) `((%log
) ,(cadr s1
)))))
2562 (defun asymredu (rd)
2563 (cond ((atom rd
) rd
)
2565 ((not (among var rd
)) rd
)
2566 ((polyinx rd var t
))
2568 (cond ((eq (cadr rd
) var
) rd
)
2570 (factor ($expand
(m^
(polyinx (cadr rd
) var t
)
2574 (t (simplify (cons (list (caar rd
))
2575 (mapcar #'asymredu
(cdr rd
)))))))
2578 (let ((dn** ()) (nn** ()))
2579 (cond ((atom rd
) rd
)
2581 ((not (among var rd
)) rd
)
2585 (cond ((eq (cadr rd
) var
) rd
)
2586 (t (setq dn
** (caddr rd
))
2587 (setq nn
** (factor (cadr rd
)))
2588 (cond ((mtimesp nn
**)
2589 (m*l
(mapcar #'(lambda (j) (m^ j dn
**))
2592 (t (simplify (cons (ncons (caar rd
))
2593 (mapcar #'rdfact
(cdr rd
))))))))
2595 (defun cnv (expl val
)
2596 (mapcar #'(lambda (e)
2597 (maxima-substitute (cond ((eq val
'$zerob
)
2598 (m* -
1 (m^ var -
1)))
2603 (t (m^
(m+ var
(m* -
1 val
)) -
1.
)))
2608 (defun pwtaylor (exp var l terms
)
2609 (prog (coef ans c mc
)
2610 (cond ((=0 terms
) (return nil
)) ((=0 l
) (setq mc t
)))
2613 loop
(setq c
(1+ c
))
2614 (cond ((or (> c
10.
) (equal c terms
))
2616 (t (setq exp
(sdiff exp var
))))
2617 tag1
(setq coef
($radcan
(subin l exp
)))
2618 (cond ((=0 coef
) (setq terms
(1+ terms
)) (go loop
)))
2625 (m^
`((mfactorial) ,c
) -
1)
2626 (m^
(if mc var
(m+t
(m*t -
1 l
) var
)) c
)))))
2631 ((simplerd e
) (rdtay e
))
2632 (t (cons (list (caar e
))
2633 (mapcar #'rdsget
(cdr e
))))))
2636 (cond (limit-using-taylor ($ratdisrep
($taylor rd var val
1.
)))
2640 (prog (varlist p c e d $ratfac
)
2641 (setq varlist
(ncons var
))
2642 (setq p
(ratnumerator (cdr (ratrep* (cadr rd
)))))
2643 (cond ((< (length p
) 3.
) (return rd
)))
2646 (setq c
(m^ var
(m* d e
)))
2647 (setq d
($ratsimp
(varinvert (m* (pdis p
) (m^ var
(m- d
)))
2649 (setq d
(pwtaylor (m^ d e
) var
0.
3.
))
2650 (return (m* c
(varinvert d var
)))))
2652 (defun varinvert (e var
) (subin (m^t var -
1.
) e
))
2655 (prog ((varlist (list var
)))
2656 (return (let (($ratfac nil
))
2658 (pdegr (cadr (ratrep* p
)))))))
2660 (defun rat-no-ratfac (e)
2661 (let (($ratfac nil
))
2666 (defun rddeg (rd low
*)
2667 (cond ((or (mnump rd
)
2668 (not (among var rd
)))
2673 (m* (deg (cadr rd
)) (caddr rd
)))
2675 (addn (mapcar #'(lambda (j)
2679 (setq rd
(andmapcar #'(lambda (j) (rddeg j low
*))
2681 (cond (low* (ratmin rd
))
2685 (cond ((or (atom pf
) (not (eq (caadr (ratf var
)) (car pf
))))
2687 (low* (cadr (reverse pf
)))
2689 ;;There is some confusion here. We need to be aware of Branch cuts etc....
2690 ;;when doing this section of code. It is not very carefully done so there
2691 ;;are bugs still lurking. Another misfortune is that LIMIT or its inferiors
2692 ;;sometimes decides to change the limit VAL in midstream. This must be corrected
2693 ;;since LIMIT's interaction with the data base environment must be maintained.
2694 ;;I'm not sure that this code can ever be called with VAL other than $INF but
2695 ;;there is a hook in the first important cond clause to cathc them anyway.
2698 (let ((num-power (rddeg n nil
))
2699 (den-power (rddeg d nil
))
2700 (coef ()) (coef-sign ()) (power ()))
2701 (setq coef
(m// ($ratcoef
($expand n
) var num-power
)
2702 ($ratcoef
($expand d
) var den-power
)))
2703 (setq coef-sign
(getsignl coef
))
2704 (setq power
(m// num-power den-power
))
2705 (cond ((eq (ask-integer power
'$integer
) '$integer
)
2706 (cond ((eq (ask-integer power
'$even
) '$even
) '$even
)
2707 (t '$odd
)))) ;Can be extended from here.
2708 (cond ((or (eq val
'$minf
)
2711 (equal val
0)) ()) ;Can be extended to cover some these.
2712 ((ratgreaterp den-power num-power
)
2713 (cond ((equal coef-sign
1.
) '$zeroa
)
2714 ((equal coef-sign -
1) '$zerob
)
2715 ((equal coef-sign
0) 0)
2717 ((ratgreaterp num-power den-power
)
2718 (cond ((equal coef-sign
1.
) '$inf
)
2719 ((equal coef-sign -
1) '$minf
)
2720 ((equal coef-sign
0) nil
) ; should never be zero
2721 ((null coef-sign
) '$infinity
)))
2724 (defun radlim (e n d
)
2726 (cond ((eq val
'$infinity
) (throw 'limit nil
))
2728 (setq nl
(m* var -
1))
2729 (setq n
(subin nl n
))
2730 (setq d
(subin nl d
))
2731 (setq val
'$inf
))) ;This is the Culprit. Doesn't tell the DATABASE.
2732 (cond ((eq val
'$inf
)
2733 (setq nl
(asymredu n
))
2734 (setq dl
(asymredu d
))
2736 ((or (rptrouble n
) (rptrouble d
))
2737 (return (limit (m* (rdsget n
) (m^
(rdsget d
) -
1.
)) var val t
)))
2738 (t (return (asy nl dl
))))))
2739 (setq nl
(limit n var val t
))
2740 (setq dl
(limit d var val t
))
2741 (cond ((and (zerop2 nl
) (zerop2 dl
))
2742 (cond ((or (polyp n
) (polyp d
))
2743 (return (try-lhospital-quit n d t
)))
2744 (t (return (ser0 e n d val
)))))
2745 (t (return ($radcan
(ratrad (m// n d
) n d nl dl
)))))))
2747 (defun ratrad (e n d nl dl
)
2749 (cond ((equal nl
0) (return 0))
2752 (cond ((equal dl
0) (setq d1
'$infinity
)) ;No direction Info.
2755 ((equal (setq d
(behavior d var val
)) 1)
2757 ((equal d -
1) (setq d1
'$minf
))
2758 (t (throw 'limit nil
))))
2761 (cond ((equal (setq n
(behavior n var val
)) 1)
2763 ((equal n -
1) (setq n1
'$zerob
))
2765 (t (return ($radcan
(ridofab (subin val e
))))))
2766 (return (simplimtimes (list n1 d1
)))))
2768 ;;; Limit(log(XXX), var, 0, val), where val is either zerob (limit from below)
2769 ;;; or zeroa (limit from above).
2770 (defun simplimln (expr var val
)
2771 (let ((arglim (limit (cadr expr
) var val
'think
)) (dir))
2772 (cond ((eq arglim
'$inf
) '$inf
) ;log(inf) = inf
2773 ;;log(minf,infinity,zerob) = infinity & log(0) = infinity
2774 ((or (member arglim
'($minf $infinity $zerob
)))
2776 ((eq arglim
'$zeroa
) '$minf
) ;log(zeroa) = minf
2777 ;; log(ind)=und & log(und)=und
2778 ((member arglim
'($ind $und
)) '$und
)
2779 ;; log(1^(-)) = zerob, log(1^(+)) = zeroa & log(1)=0
2781 (if (or (eq val
'$zerob
) (eq var
'$zeroa
)) val
0))
2782 ;; Special case of arglim = 0
2784 (setq dir
(behavior (cadr expr
) var val
))
2785 (cond ((eql dir -
1) '$infinity
)
2786 ((eql dir
0) '$infinity
)
2787 ((eql dir
1) '$minf
)))
2788 ;; When arglim is off the negative real axis, use direct substitution
2789 ((off-negative-real-axisp arglim
)
2790 (ftake '%log arglim
))
2792 ;; We know that arglim is a negative real number, say xx.
2793 ;; When the imaginary part of (cadr expr) near var is negative,
2794 ;; return log(-x) - %i*pi; when the imaginary part of (cadr expr)
2795 ;; near var is positive return log(-x) + %i*pi; and when
2796 ;; we cannot determine the behavior of the imaginary part,
2797 ;; return a nounform. The value of val (either zeroa or zerob)
2798 ;; determines what is meant by "near" (smaller than var when
2799 ;; val is zerob and larger than var when val is zeroa).
2800 (setq dir
(behavior ($imagpart
(cadr expr
)) var val
))
2801 (cond ((or (eql dir
1) (eql dir -
1))
2802 (add (ftake '%log
(mul -
1 arglim
)) (mul dir
'$%i
'$%pi
)))
2803 (t (throw 'limit nil
))))))) ;do a nounform return
2804 (setf (get '%log
'simplim%function
) #'simplimln
)
2805 (setf (get '%plog
'simplim%function
) #'simplimln
)
2807 (defun simplim%limit
(e x pt
)
2808 (declare (ignore e x pt
))
2810 (setf (get '%limit
'simplim%function
) #'simplim%limit
)
2812 (defun simplim%unit_step
(e var val
)
2813 (let ((lim (limit (cadr e
) var val
'think
)))
2814 (cond ((eq lim
'$und
) '$und
)
2815 ((eq lim
'$ind
) '$ind
)
2816 ((eq lim
'$zerob
) 0)
2817 ((eq lim
'$zeroa
) 1)
2818 ((not (lenient-realp lim
)) (throw 'limit nil
)) ;catches infinity too
2819 ;; catches minf & inf cases
2820 ((eq t
(mgrp 0 lim
)) 0)
2821 ((eq t
(mgrp lim
0)) 1)
2823 (setf (get '$unit_step
'simplim%function
) #'simplim%unit_step
)
2825 (defun simplim%conjugate
(e var val
)
2826 (let ((lim (limit (cadr e
) var val
'think
)))
2827 (cond ((off-negative-real-axisp lim
)
2828 (ftake '$conjugate lim
))
2829 (t (throw 'limit nil
)))))
2830 (setf (get '$conjugate
'simplim%function
) #'simplim%conjugate
)
2832 (defun simplim%imagpart
(e var val
)
2833 (let ((lim (limit (cadr e
) var val
'think
)))
2834 (cond ((eq lim
'$und
) '$und
)
2836 ((eq lim
'$infinity
) (throw 'limit nil
))
2837 (t (mfuncall '$imagpart lim
)))))
2838 (setf (get '$imagpart
'simplim%function
) #'simplim%imagpart
)
2839 (setf (get '%imagpart
'simplim%function
) #'simplim%imagpart
)
2841 (defun simplim%realpart
(e var val
)
2842 (let ((lim (limit (cadr e
) var val
'think
)))
2843 (cond ((eq lim
'$und
) '$und
)
2844 ((eq lim
'$ind
) '$ind
)
2845 ((eq lim
'$infinity
) (throw 'limit nil
))
2846 (t (mfuncall '$realpart lim
)))))
2847 (setf (get '$realpart
'simplim%function
) #'simplim%realpart
)
2848 (setf (get '%realpart
'simplim%function
) #'simplim%realpart
)
2849 ;;; Limit of the Factorial function
2851 (defun simplimfact (expr var val
)
2852 (let* ((arglim (limit (cadr expr
) var val
'think
)) ; Limit of the argument.
2854 (cond ((eq arglim
'$inf
) '$inf
)
2855 ((member arglim
'($minf $infinity $und $ind
) :test
#'eq
) '$und
)
2857 ((and (or (maxima-integerp arglim
)
2858 (setq arg2
(integer-representation-p arglim
)))
2859 (eq ($sign arg2
) '$neg
))
2860 ;; A negative integer or float or bigfloat representation.
2861 (let ((dir (limit (add (cadr expr
) (mul -
1 arg2
)) var val
'think
))
2862 (even (mevenp arg2
)))
2863 (cond ((or (and even
2873 (t (throw 'limit nil
)))))
2875 ;; Call simplifier to get value at the limit of the argument.
2876 (simplify (list '(mfactorial) arglim
))))))
2877 (setf (get 'mfactorial
'simplim%function
) #'simplimfact
)
2879 (defun simplim%erf-%tanh
(fn arg
)
2880 (let ((arglim (limit arg var val
'think
))
2883 (cond ((eq arglim
'$inf
) 1)
2884 ((eq arglim
'$minf
) -
1)
2885 ((eq arglim
'$infinity
)
2886 (destructuring-bind (rpart . ipart
)
2888 (setq rlim
(limit rpart var origval
'think
))
2889 (cond ((eq fn
'%tanh
)
2890 (cond ((equal rlim
'$inf
) 1)
2891 ((equal rlim
'$minf
) -
1)))
2893 (setq ans
(limit (m* rpart
(m^t ipart -
1)) var origval
'think
))
2894 (setq ans
($asksign
(m+ `((mabs) ,ans
) -
1)))
2895 (cond ((or (eq ans
'$pos
) (eq ans
'$zero
))
2896 (cond ((eq rlim
'$inf
) 1)
2897 ((eq rlim
'$minf
) -
1)
2900 ((eq arglim
'$und
) '$und
)
2901 ((member arglim
'($zeroa $zerob $ind
) :test
#'eq
) arglim
)
2902 ;;;Ignore tanh(%pi/2*%I) and multiples of the argument.
2904 ;; erf (or tanh) of a known value is just erf(arglim).
2905 (simplify (list (ncons fn
) arglim
))))))
2907 (defun in-domain-of-atan (z)
2908 (setq z
(trisplit z
)) ; split z into real and imaginary parts
2909 (let ((x (car z
)) (y (cdr z
))) ;z = x+%i*y
2911 (eq t
(mnqp x
0)) ; x # 0
2913 (eq t
(mgrp -
1 y
)) ; y < -1
2914 (eq t
(mgrp y
1)))))) ; y > 1
2916 (defun simplim%atan
(e x pt
)
2917 (let ((lim (limit (cadr e
) x pt
'think
)))
2918 (cond ((or (eq lim
'$zeroa
) (eq lim
'$zerob
) (eq lim
0) (eq lim
'$ind
)) lim
)
2919 ((or (eq lim
'$und
) (eq lim
'$infinity
)) (throw 'limit nil
))
2920 ((eq lim
'$inf
) #$%pi
/2$
) ;atan(inf) -> %pi/2
2921 ((eq lim
'$minf
) #$-%pi
/2$
) ;atan(-inf) -> -%pi/2
2922 ((in-domain-of-atan (ridofab lim
)) ; direct substitution
2923 (ftake '%atan
(ridofab lim
)))
2924 (t (limit ($logarc e
) x pt
'think
)))))
2925 (setf (get '%atan
'simplim%function
) #'simplim%atan
)
2927 (defmvar extended-reals
2928 (append infinitesimals infinities
(list '$und
'$ind
)))
2930 ;; Most instances of atan2 are simplified to atan expressions, but this routine
2931 ;; handles tricky cases such as limit(atan2((x^2-2), x^3-2*x), x, sqrt(2), minus).
2932 ;; Taylor and Gruntz cannot handle the discontinuity at atan(0, -1)
2934 ;; When possible, we want to evaluate the limit of an atan2 expression using
2935 ;; direct substituion--that produces, I think, the least surprising values.
2937 ;; The general simplifier catches atan2(0,0) and it transforms atan2(minf or inf,X)
2938 ;; and atan2(X, minf or inf) into an atan expression, but it doesn't catch
2939 ;; atan2(X, zerob or zeroa) or atan2(zerob or zeroa, X). For the other extended
2940 ;; real (ind,und, or infinity) arguments, the general simplifer gives sign errors.
2942 (defun simplim%atan2
(e v pt
)
2943 (let ((y (second e
)) (x (third e
)) (xlim) (ylim) (xlim-z) (ylim-z) (q))
2944 (setq xlim
(limit x v pt
'think
))
2945 (setq ylim
(limit y v pt
'think
))
2946 (setq xlim-z
(ridofab xlim
)
2947 ylim-z
(ridofab ylim
))
2948 ;; For cases for which direct substitution fails, normalize
2949 ;; x & y and try again.
2950 (setq q
(cond ((eq xlim
'$inf
) x
)
2951 ((eq xlim
'$minf
) (mul -
1 x
))
2953 ((eq ylim
'$minf
) (mul -
1 y
))
2954 ((and (eq xlim
'$zerob
) (zerop2 ylim
)) (mul -
1 x
))
2955 ((and (eq xlim
'$zeroa
) (zerop2 ylim
)) x
)
2956 ((and (eq ylim
'$zerob
) (zerop2 xlim
)) (mul -
1 y
))
2957 ((and (eq ylim
'$zeroa
) (zerop2 xlim
)) y
)
2960 (when (not (eql q
1))
2963 (setq xlim
(limit x v pt
'think
))
2964 (setq ylim
(limit y v pt
'think
))
2965 (setq xlim-z
(ridofab xlim
)
2966 ylim-z
(ridofab ylim
)))
2969 ((and (eq '$zerob ylim
) (eq t
(mgrp 0 xlim
))) (mul -
1 '$%pi
))
2970 ((and (eq '$zerob ylim
) (eq t
(mgrp xlim
0))) 0)
2971 ((and (eq '$zeroa ylim
) (eq t
(mgrp 0 xlim
))) '$%pi
)
2972 ((and (eq '$zeroa ylim
) (eq t
(mgrp xlim
0))) 0)
2973 ((and (eql xlim
1) (eql ylim
'$inf
)) (div '$%pi
2))
2974 ((and (eql xlim -
1) (eql ylim
0)) '$ind
)
2976 ;; Use direct substitution when ylim-z # 0 or xlim-z > 0. We need
2977 ;; to check that xlim-z & ylim-z are real too.
2978 ((and (lenient-realp xlim-z
) (lenient-realp ylim-z
)
2979 (or (eq t
(mnqp ylim-z
0)) (eq t
(mgrp xlim-z
0))))
2980 (ftake '$atan2 ylim-z xlim-z
))
2982 (throw 'limit nil
)))))
2983 (setf (get '$atan2
'simplim%function
) #'simplim%atan2
)
2985 (defun simplimsch (sch arg
)
2986 (cond ((real-infinityp arg
)
2987 (cond ((eq sch
'%sinh
) arg
) (t '$inf
)))
2988 ((eq arg
'$infinity
) '$infinity
)
2989 ((eq arg
'$ind
) '$ind
)
2990 ((eq arg
'$und
) '$und
)
2991 ((and (eq sch
'%sinh
)
2992 (or (eq arg
'$zerob
) (eq arg
'$zeroa
)))
2994 (t (let (($exponentialize t
))
2995 (resimplify (list (ncons sch
) (ridofab arg
)))))))
2997 ;; simple limit of sin and cos
2998 (defun simplimsc (exp fn arg
)
2999 (cond ((member arg
'($inf $minf $ind
) :test
#'eq
) '$ind
)
3000 ((member arg
'($und $infinity
) :test
#'eq
)
3002 ((member arg
'($zeroa $zerob
) :test
#'eq
)
3003 (cond ((eq fn
'%sin
) arg
)
3004 (t (m+ 1 '$zerob
))))
3006 (simplify (list (ncons fn
) (ridofab arg
)))
3009 (defun simplim%tan
(arg)
3010 (let ((arg1 (ridofab (limit arg var val
'think
))))
3012 ((member arg1
'($inf $minf $infinity $ind $und
) :test
#'eq
) '$und
)
3014 (let ((c (trigred (pip arg1
))))
3015 (cond ((not (equal ($imagpart arg1
) 0)) '$infinity
)
3016 ((and (eq (caar c
) 'rat
)
3019 (setq arg1
(behavior arg var val
))
3020 (cond ((= arg1
1) '$inf
)
3021 ((= arg1 -
1) '$minf
)
3023 ((and (eq (caar c
) 'rat
)
3026 (setq arg1
(behavior arg var val
))
3027 (cond ((= arg1
1) '$minf
)
3030 (t (throw 'limit
())))))
3032 (setq arg1
(behavior arg var val
))
3033 (cond ((equal arg1
1) '$zeroa
)
3034 ((equal arg1 -
1) '$zerob
)
3036 (t (simp-%tan
(list '(%tan
) arg1
) 1 nil
)))))
3038 (defun simplim%asinh
(arg)
3039 (cond ((member arg
'($inf $minf $zeroa $zerob $ind $und
) :test
#'eq
)
3041 ((eq arg
'$infinity
) '$und
)
3042 (t (simplify (list '(%asinh
) (ridofab arg
))))))
3044 (defun simplim%acosh
(arg)
3045 (cond ((equal (ridofab arg
) 1) '$zeroa
)
3046 ((eq arg
'$inf
) arg
)
3047 ((eq arg
'$minf
) '$infinity
)
3048 ((member arg
'($und $ind $infinity
) :test
#'eq
) '$und
)
3049 (t (simplify (list '(%acosh
) (ridofab arg
))))))
3051 (defun simplim%atanh
(arg dir
)
3052 ;; Compute limit(atanh(x),x,arg). If ARG is +/-1, we need to take
3053 ;; into account which direction we're approaching ARG.
3054 (cond ((zerop2 arg
) arg
)
3055 ((member arg
'($ind $und $infinity $minf $inf
) :test
#'eq
)
3057 ((equal (setq arg
(ridofab arg
)) 1.
)
3058 ;; The limit at 1 should be complex infinity because atanh(x)
3059 ;; is complex for x > 1, but inf if we're approaching 1 from
3061 (if (eq dir
'$zerob
)
3065 ;; Same as above, except for the limit is at -1.
3066 (if (eq dir
'$zeroa
)
3069 (t (simplify (list '(%atanh
) arg
)))))
3071 (defun simplim%asin
(e x pt
)
3072 (let ((lim (limit (cadr e
) x pt
'think
)) (dir) (lim-sgn))
3073 (cond ((member lim
'($zeroa $zerob
)) lim
) ;asin(zeoroa/b) = zeroa/b
3074 ((in-domain-of-asin lim
) ;direct subsititution
3076 ((member lim
'($und $ind $inf $minf $infinity
)) ;boundary cases
3079 (setq e
(trisplit (cadr e
))) ;overwrite e!
3080 (setq dir
(behavior (cdr e
) x pt
))
3081 (setq lim-sgn
($csign lim
))
3084 (throw 'limit t
)) ;unable to find behavior of imaginary part
3085 ;; for the values of asin on the branch cuts, see DLMF 4.23.20 & 4.23.21
3086 ((or (eq '$pos lim-sgn
) (eq '$neg lim-sgn
))
3087 ;; continuous from above
3088 (if (eql dir
1) (ftake '%asin lim
) (sub '$%pi
(ftake '%asin lim
))))
3089 (t (throw 'limit t
))))))) ; unable to find sign of real part of lim.
3090 (setf (get '%asin
'simplim%function
) #'simplim%asin
)
3092 (defun simplim%acos
(e x pt
)
3093 (let ((lim (limit (cadr e
) x pt
'think
)) (dir) (lim-sgn))
3094 (cond ((in-domain-of-asin lim
) ;direct subsititution
3096 ((member lim
'($und $ind $inf $minf $infinity
)) ;boundary cases
3099 (setq e
(trisplit (cadr e
))) ;overwrite e!
3100 (setq dir
(behavior (cdr e
) x pt
))
3101 (setq lim-sgn
($csign lim
))
3104 (throw 'limit t
)) ;unable to find behavior of imaginary part
3105 ;; for the values of acos on the branch cuts, see DLMF 4.23.24 & 4.23.25
3106 ;; http://dlmf.nist.gov/4.23.E24
3107 ((or (eq '$pos lim-sgn
) (eq '$neg lim-sgn
))
3108 ;; continuous from above
3109 (if (eql dir
1) (ftake '%acos lim
) (sub (mul 2 '$%pi
) (ftake '%acos lim
))))
3110 (t (throw 'limit t
))))))) ; unable to find sign of real part of lim.
3111 (setf (get '%acos
'simplim%function
) #'simplim%acos
)
3113 ;; Limit of an %integrate expression. For a definite integral
3114 ;; integrate(ee,var,a,b), when ee is free of the limit variable
3115 (defun simplim%integrate
(e x pt
)
3116 (let* ((ee (second e
)) ;ee = integrand
3117 (var (third e
)) ;integration variable
3118 (a (fourth e
)) ;lower limit or nil if indefinite
3119 (b (fifth e
)) ;lower limit or nil if indefinite
3121 (cond ((and a b
($freeof x ee
) ($freeof x var
))
3122 (setq alim
(limit a x pt
'think
))
3123 (setq blim
(limit b x pt
'think
))
3124 (if (and (lenient-extended-realp alim
)
3125 (lenient-extended-realp blim
)
3126 (not (eq alim
'$infinity
))
3127 (not (eq blim
'$infinity
)))
3128 (ftake '%integrate ee var alim blim
)
3131 (throw 'limit t
)))))
3132 (setf (get '%integrate
'simplim%function
) #'simplim%integrate
)
3134 (defun subftake (op subarg arg
)
3135 (simplifya (subfunmake op subarg arg
) t
))
3137 (defun off-one-to-inf (z)
3138 (setq z
(trisplit z
)) ; split z into x+%i*y
3140 (eq t
(mnqp (cdr z
) 0)) ; y # 0
3141 (eq t
(mgrp 1 (car z
))))) ; x < 1
3143 (defun simplim%li
(expr x pt
)
3144 (let ((n (car (subfunsubs expr
))) (e (car (subfunargs expr
))))
3146 (setq e
(limit e x pt
'think
))
3147 (cond ((and (eq e
'$minf
) (eql n
2))
3149 ((and (eq e
'$inf
) (eql n
2))
3151 ((or (eql (ridofab e
) 1) (off-one-to-inf e
))
3152 ;; Limit of li[s](1) can be evaluated by just
3153 ;; substituting in 1.
3154 (subftake '$li
(list n
) (list e
)))
3155 (t (throw 'limit nil
))))
3156 ;; Claim ignorance when order depends on limit variable.
3157 (t (throw 'limit nil
)))))
3159 (setf (get '$li
'simplim%function
) #'simplim%li
)
3160 (setf (get '%li
'simplim%function
) #'simplim%li
)
3162 (defun simplim$psi
(order arg val
)
3163 (if (and (not (equal (length order
) 1))
3164 (not (equal (length arg
) 1)))
3166 (setq order
(car order
)
3168 (cond ((equal order
0)
3169 (destructuring-bind (rpart . ipart
)
3171 (cond ((not (equal ipart
0)) (throw 'limit
()))
3172 (t (setq rpart
(limit rpart var val
'think
))
3173 (cond ((eq rpart
'$zeroa
) '$minf
)
3174 ((eq rpart
'$zerob
) '$inf
)
3175 ((eq rpart
'$inf
) '$inf
)
3176 ((eq rpart
'$minf
) '$und
)
3177 ((equal (getsignl rpart
) -
1) (throw 'limit
()))
3178 (t (simplify (subfunmake '$psi
(list order
)
3179 (list rpart
)))))))))
3180 ((and (integerp order
) (> order
0)
3181 (equal (limit arg var val
'think
) '$inf
))
3182 (cond ((mevenp order
) '$zerob
)
3183 ((moddp order
) '$zeroa
)
3184 (t (throw 'limit
()))))
3185 (t (throw 'limit
()))))
3187 (defun simplim%inverse_jacobi_ns
(arg m
)
3188 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3190 `((%inverse_jacobi_ns
) ,arg
,m
)))
3192 (defun simplim%inverse_jacobi_nc
(arg m
)
3193 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3194 `((%elliptic_kc
) ,m
)
3195 `((%inverse_jacobi_nc
) ,arg
,m
)))
3197 (defun simplim%inverse_jacobi_sc
(arg m
)
3198 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3199 `((%elliptic_kc
) ,m
)
3200 `((%inverse_jacobi_sc
) ,arg
,m
)))
3202 (defun simplim%inverse_jacobi_dc
(arg m
)
3203 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3204 `((%elliptic_kc
) ,m
)
3205 `((%inverse_jacobi_dc
) ,arg
,m
)))
3207 (defun simplim%inverse_jacobi_cs
(arg m
)
3208 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3210 `((%inverse_jacobi_cs
) ,arg
,m
)))
3212 (defun simplim%inverse_jacobi_ds
(arg m
)
3213 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3215 `((%inverse_jacobi_ds
) ,arg
,m
)))
3217 (defun simplim%signum
(e x pt
)
3218 (let ((e (limit (cadr e
) x pt
'think
)) (sgn))
3219 (cond ((eq '$minf e
) -
1)
3221 ((eq '$infinity e
) '$und
)
3222 ((eq '$ind e
) '$ind
)
3227 (setq sgn
(mnqp e
0))
3228 (cond ((eq t sgn
) (ftake '%signum e
))
3229 (t (throw 'limit nil
))))))) ; don't know
3230 (setf (get '%signum
'simplim%function
) #'simplim%signum
)
3232 ;; more functions for limit to handle
3234 (defun lfibtophi (e)
3235 (cond ((not (involve e
'($fib
))) e
)
3236 ((eq (caar e
) '$fib
)
3237 (let ((lnorecurse t
))
3238 ($fibtophi
(list '($fib
) (lfibtophi (cadr e
))) lnorecurse
)))
3240 (mapcar #'lfibtophi
(cdr e
))))))
3242 ;;; FOLLOWING CODE MAKES $LDEFINT WORK
3244 (defmfun $ldefint
(exp var ll ul
&aux $logabs ans a1 a2
)
3245 (setq $logabs t ans
(sinint exp var
)
3246 a1
(toplevel-$limit ans var ul
'$minus
)
3247 a2
(toplevel-$limit ans var ll
'$plus
))
3248 (and (member a1
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3249 (setq a1
(nounlimit ans var ul
)))
3250 (and (member a2
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3251 (setq a2
(nounlimit ans var ll
)))
3252 ($expand
(m- a1 a2
)))
3254 (defun nounlimit (exp var val
)
3255 (setq exp
(restorelim exp
))
3256 (nconc (list '(%limit
) exp var
(ridofab val
))
3257 (cond ((eq val
'$zeroa
) '($plus
))
3258 ((eq val
'$zerob
) '($minus
)))))
3260 ;; substitute inside noun form of %derivative
3261 ;; for cases such as limit('diff(x+2,x), x, 1)
3262 ;; -> limit('diff(xx+3), xx, 0)
3264 ;; maxima-substitute with *atp* skips over %derivative
3266 ;; substitutes diff(f(realvar), realvar, n)
3267 ;; -> diff(f(var+val), var, n)
3268 (defun derivative-subst (exp val var realvar
)
3269 (cond ((atom exp
) exp
)
3270 ((eq '%derivative
(caar exp
))
3273 (cons ;; the function being differentiated
3274 (maxima-substitute (m+ val var
) realvar
(cadr exp
))
3275 (cons ;; the var of differentiation
3276 (maxima-substitute var realvar
(caddr exp
))
3277 (cdddr exp
))))) ;; the order of the derivative
3279 (mapcar (lambda (x) (derivative-subst x val var realvar
))
3282 ;;;Used by Defint also.
3284 (or (involve e
'(%sin %cos %tan
))
3285 (among '$%i
(%einvolve e
))))
3287 (defun %einvolve
(e)
3288 (when (among '$%e e
) (%einvolve01 e
)))
3290 (defun %einvolve01
(e)
3291 (cond ((atom e
) nil
)
3295 (among var
(caddr e
)))
3297 (t (some #'%einvolve
(cdr e
)))))
3299 (declare-top (unspecial *indicator nn
* dn
* exp var val origval taylored
3300 $tlimswitch logcombed lhp? lhcount $ratfac
))
3306 ;; "On Computing Limits in a Symbolic Manipulation System"
3307 ;; PhD Dissertation ETH Zurich 1996
3309 ;; The algorithm identifies the most rapidly varying (MRV) subexpression,
3310 ;; replaces it with a new variable w, rewrites the expression in terms
3311 ;; of the new variable, and then repeats.
3313 ;; The algorithm doesn't handle oscillating functions, so it can't do things like
3314 ;; limit(sin(x)/x, x, inf).
3316 ;; To handle limits involving functions like gamma(x) and erf(x), the
3317 ;; gruntz algorithm requires them to be written in terms of asymptotic
3318 ;; expansions, which maxima cannot currently do.
3320 ;; The algorithm assumes that everything is real, so it can't
3321 ;; currently handle limit((-2)^x, x, inf).
3323 ;; This is one of the methods used by maxima's $limit.
3324 ;; It is also directly available to the user as $gruntz.
3327 ;; most rapidly varying subexpression of expression exp with respect to limit variable var.
3328 ;; returns a list of subexpressions which are in the same MRV equivalence class.
3329 (defun mrv (exp var
)
3330 (cond ((freeof var exp
)
3335 (mrv-max (mrv (cadr exp
) var
)
3336 (mrv (m*l
(cddr exp
)) var
)
3339 (mrv-max (mrv (cadr exp
) var
)
3340 (mrv (m+l
(cddr exp
)) var
)
3343 (cond ((freeof var
(caddr exp
))
3344 (mrv (cadr exp
) var
))
3345 ((member (limitinf (logred exp
) var
) '($inf $minf
) :test
#'eq
)
3346 (mrv-max (list exp
) (mrv (caddr exp
) var
) var
))
3347 (t (mrv-max (mrv (cadr exp
) var
) (mrv (caddr exp
) var
) var
))))
3349 (mrv (cadr exp
) var
))
3350 ((equal (length (cdr exp
)) 1)
3351 (mrv (cadr exp
) var
))
3352 ((equal (length (cdr exp
)) 2)
3353 (mrv-max (mrv (cadr exp
) var
)
3354 (mrv (caddr exp
) var
)
3356 (t (tay-error "mrv not implemented" exp
))))
3358 ;; takes two lists of expresions, f and g, and limit variable var.
3359 ;; members in each list are assumed to be in same MRV equivalence
3360 ;; class. returns MRV set of the union of the inputs - either f or g
3361 ;; or the union of f and g.
3362 (defun mrv-max (f g var
)
3369 (return (union f g
))))
3370 (let ((c (mrv-compare (car f
) (car g
) var
)))
3376 (return (union f g
)))
3377 (t (merror "MRV-MAX: expected '>' '<' or '='; found: ~M" c
))))))
3379 (defun mrv-compare (a b var
)
3380 (let ((c (limitinf (m// `((%log
) ,a
) `((%log
) ,b
)) var
)))
3383 ((member c
'($inf $minf
) :test
#'eq
)
3387 ;; rewrite expression exp by replacing members of MRV set omega with
3388 ;; expressions in terms of new variable wsym. return cons pair of new
3389 ;; version of exp and the log of the new variable wsym.
3390 (defun mrv-rewrite (exp omega var wsym
)
3391 (setq omega
(stable-sort omega
(lambda (x y
) (> (length (mrv x var
))
3392 (length (mrv y var
))))));FIXME consider a total order function with #'sort
3393 (let* ((g (car (last omega
)))
3395 (sig (equal (mrv-sign logg var
) 1))
3396 (w (if sig
(m// 1 wsym
) wsym
))
3397 (logw (if sig
(m* -
1 logg
) logg
)))
3398 (mapcar (lambda (x y
)
3399 ;;(mtell "y:~M x:~M exp:~M~%" y x exp)
3400 (setq exp
(syntactic-substitute y x exp
)))
3402 (mapcar (lambda (f) ;; rewrite each element of omega
3403 (let* ((logf (logred f
))
3404 (c (mrv-leadterm (m// logf logg
) var nil
)))
3405 (cond ((not (equal (cadr c
) 0))
3406 (merror "MRV-REWRITE: expected leading term to be constant in ~M" c
)))
3407 ;;(mtell "logg: ~M logf: ~M~%" logg logf)
3410 (m* (car c
) logg
))))))
3414 ;;; if log w(x) = h(x), rewrite all subexpressions of the form
3415 ;;; log(f(x)) as log(w^-c f(x)) + c h(x) with c the unique constant
3416 ;;; such that w^-c f(x) is strictly less rapidly varying than w.
3417 (defun mrv-rewrite-logs (exp wsym logw
)
3418 (cond ((atom exp
) exp
)
3420 (not (freeof wsym exp
)))
3421 (let* ((f (cadr exp
))
3422 (c ($lopow
(calculate-series f wsym
)
3425 (m* (m^ wsym
(m- c
))
3426 (mrv-rewrite-logs f wsym logw
)))
3431 (mrv-rewrite-logs e wsym logw
))
3434 ;; returns list of two elements: coeff and exponent of leading term of exp,
3435 ;; after rewriting exp in term of its MRV set omega.
3436 (defun mrv-leadterm (exp var omega
)
3437 (prog ((new-omega ()))
3438 (cond ((freeof var exp
)
3439 (return (list exp
0))))
3440 (dolist (term omega
)
3441 (cond ((subexp exp term
)
3442 (push term new-omega
))))
3443 (setq omega new-omega
)
3445 (setq omega
(mrv exp var
))))
3446 (cond ((member var omega
:test
#'eq
)
3447 (let* ((omega-up (mrv-moveup omega var
))
3448 (e-up (car (mrv-moveup (list exp
) var
)))
3449 (mrv-leadterm-up (mrv-leadterm e-up var omega-up
)))
3450 (return (mrv-movedown mrv-leadterm-up var
)))))
3451 (destructuring-let* ((wsym (gensym "w"))
3454 ((f . logw
) (mrv-rewrite exp omega var wsym
))
3455 (series (calculate-series (mrv-rewrite-logs f wsym logw
)
3457 (setq series
(maxima-substitute logw
`((%log
) ,wsym
) series
))
3458 (setq lo
($lopow series wsym
))
3459 (when (or (not ($constantp lo
))
3460 (not (free series
'%derivative
)))
3461 ;; (mtell "series: ~M lo: ~M~%" series lo)
3462 (tay-error "error in series expansion" f
))
3463 (setq coef
($coeff series wsym lo
))
3464 (when (not (free coef wsym
))
3465 (tay-error "MRV-LEADTERM: failed to extract leading coefficient; obtained" coef
))
3466 ;;(mtell "exp: ~M f: ~M~%" exp f)
3467 ;;(mtell "series: ~M~%coeff: ~M~%pow: ~M~%" series coef lo)
3468 (return (list coef lo
)))))
3470 (defun mrv-moveup (l var
)
3471 (mapcar (lambda (exp)
3472 (simplify-log-of-exp
3473 (syntactic-substitute `((mexpt) $%e
,var
) var exp
)))
3476 (defun mrv-movedown (l var
)
3477 (mapcar (lambda (exp) (syntactic-substitute `((%log simp
) ,var
) var exp
))
3480 ;; test whether sub is a subexpression of exp
3481 (defun subexp (exp sub
)
3482 (let ((dummy (gensym)))
3483 (putprop dummy t
'internal
)
3484 (not (alike1 (maxima-substitute dummy
3489 ;; Generate $lhospitallim terms of taylor expansion.
3490 ;; Ideally we would use a lazy series representation that generates
3491 ;; more terms as higher order terms cancel.
3492 (defun calculate-series (exp var
)
3493 (let ((cntx ($supcontext
)))
3497 (mfuncall '$assume
(ftake 'mgreaterp var
0))
3498 (putprop var t
'internal
); keep var from appearing in questions to user
3499 ($taylor exp var
0 $lhospitallim
))
3500 (remprop var
'internal
)
3501 ($killcontext cntx
))))
3503 (defun mrv-sign (exp var
)
3504 (cond ((freeof var exp
)
3505 (let ((sign ($sign
($radcan exp
))))
3506 (cond ((eq sign
'$zero
)
3512 (t (tay-error " cannot determine mrv-sign" exp
)))))
3516 (* (mrv-sign (cadr exp
) var
)
3517 (mrv-sign (m*l
(cddr exp
)) var
)))
3519 (equal (mrv-sign (cadr exp
) var
) 1))
3522 (cond ((equal (mrv-sign (cadr exp
) var
) -
1)
3523 (tay-error " complex expression in gruntz limit" exp
)))
3524 (mrv-sign (m+ -
1 (cadr exp
)) var
))
3526 (mrv-sign (limitinf exp var
) var
))
3527 (t (tay-error " cannot determine mrv-sign" exp
))))
3529 ;; gruntz algorithm for limit of exp as var goes to positive infinity
3530 (defun limitinf (exp var
)
3531 (prog (($exptsubst nil
))
3532 (cond ((freeof var exp
)
3534 (destructuring-let* ((c0-e0 (mrv-leadterm exp var nil
))
3537 (sig (mrv-sign e0 var
)))
3538 (cond ((equal sig
1)
3541 (cond ((equal (mrv-sign c0 var
) 1)
3543 ((equal (mrv-sign c0 var
) -
1)
3547 ;; example: gruntz(n^n/(n^n+(n-1)^n), n, inf);
3548 (tay-error " infinite recursion in limitinf" exp
))
3549 (return (limitinf c0 var
)))))))
3551 ;; user-level function equivalent to $limit.
3552 ;; direction must be specified if limit point is not infinite
3553 ;; The arguments are checked and a failure of taylor is catched.
3555 (defmfun $gruntz
(expr var val
&rest rest
)
3557 (when (> (length rest
) 1)
3559 (intl:gettext
"gruntz: too many arguments; expected just 3 or 4")))
3560 (setq dir
(car rest
))
3561 (when (and (not (member val
'($inf $minf $zeroa $zerob
)))
3562 (not (member dir
'($plus $minus
))))
3564 (intl:gettext
"gruntz: direction must be 'plus' or 'minus'")))
3566 (catch 'taylor-catch
3567 (let ((silent-taylor-flag t
))
3568 (declare (special silent-taylor-flag
))
3569 (gruntz1 expr var val dir
))))
3570 (if (or (null ans
) (eq ans t
))
3572 `(($gruntz simp
) ,expr
,var
, val
,dir
)
3573 `(($gruntz simp
) ,expr
,var
,val
))
3576 ;; This function is for internal use in $limit.
3577 (defun gruntz1 (exp var val
&rest rest
)
3578 (cond ((> (length rest
) 1)
3579 (merror (intl:gettext
"gruntz: too many arguments; expected just 3 or 4"))))
3580 (let (($logexpand t
) ; gruntz needs $logexpand T
3581 (newvar (gensym "w"))
3583 (putprop newvar t
'internal
); keep var from appearing in questions to user
3584 (cond ((eq val
'$inf
)
3587 (setq exp
(maxima-substitute (m* -
1 newvar
) var exp
)))
3589 (setq exp
(maxima-substitute (m// 1 newvar
) var exp
)))
3591 (setq exp
(maxima-substitute (m// -
1 newvar
) var exp
)))
3593 (setq exp
(maxima-substitute (m+ val
(m// 1 newvar
)) var exp
)))
3595 (setq exp
(maxima-substitute (m+ val
(m// -
1 newvar
)) var exp
)))
3596 (t (merror (intl:gettext
"gruntz: direction must be 'plus' or 'minus'; found: ~M") dir
)))
3597 (limitinf exp newvar
)))
3599 ;; substitute y for x in exp
3600 ;; similar to maxima-substitute but does not simplify result
3601 (defun syntactic-substitute (y x exp
)
3602 (cond ((alike1 x exp
) y
)
3605 (mapcar (lambda (exp)
3606 (syntactic-substitute y x exp
))
3609 ;; log(exp(subexpr)) -> subexpr
3610 ;; without simplifying entire exp
3611 (defun simplify-log-of-exp (exp)
3612 (cond ((atom exp
) exp
)
3615 (eq '$%e
(cadadr exp
)))
3618 (mapcar #'simplify-log-of-exp
