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 (when (eq la
'$und
) (return-from both-side
'$und
))
303 (setf lb
(toplevel-$limit exp var val
'$minus
))
304 (let ((ra (ridofab la
))
306 (cond ((eq t
(meqp ra rb
))
310 ; Maxima does not consider equal(ind,ind) to be true, but
311 ; if both one-sided limits are ind then we want to call
312 ; the two-sided limit ind (e.g., limit(sin(1/x),x,0)).
314 ((or (not (free la
'%limit
))
315 (not (free lb
'%limit
)))
317 ((and (infinityp la
) (infinityp lb
))
318 ; inf + minf => infinity
323 (defun limunknown (f)
324 (catch 'limunknown
(limunknown1 (specrepcheck f
))))
326 (defun limunknown1 (f)
327 (cond ((mapatom f
) nil
)
328 ((or (not (safe-get (caar f
) 'operators
))
329 (member (caar f
) '(%sum %product mncexpt
) :test
#'eq
)
330 ;;Special function code here i.e. for li[2](x).
331 (and (eq (caar f
) 'mqapply
)
332 (not (get (subfunname f
) 'specsimp
))))
333 (if (not (free f var
)) (throw 'limunknown t
)))
334 (t (mapc #'limunknown1
(cdr f
)) nil
)))
337 (if (involve e
'(%gamma
)) (setq e
($makefact e
)))
338 (cond ((involve e
'(mfactorial))
339 (setq e
(simplify ($minfactorial e
))))
343 ;; or nil if sign unknown or complex
345 (let ((z (ridofab z
)))
346 (if (not (free z var
)) (setq z
(toplevel-$limit z var val
)))
347 (let ((*complexsign
* t
))
348 (let ((sign (if *getsignl-asksign-ok
* ($asksign z
) ($sign z
))))
349 (cond ((eq sign
'$pos
) 1)
351 ((eq sign
'$zero
) 0))))))
353 (defun restorelim (exp)
354 (cond ((null exp
) nil
)
355 ((atom exp
) (or (and (symbolp exp
) (get exp
'limitsub
)) exp
))
356 ((and (consp (car exp
)) (eq (caar exp
) 'mrat
))
358 (cons (restorelim (cadr exp
))
359 (restorelim (cddr exp
)))))
360 (t (cons (car exp
) (mapcar #'restorelim
(cdr exp
))))))
363 (defun mabs-subst (exp var val
) ; RETURNS EXP WITH MABS REMOVED, OR THROWS.
364 (let ((d (involve exp
'(mabs)))
368 ((not (and (equal ($imagpart
(let ((v (limit-catch d var val
)))
369 ;; The above call might
370 ;; throw 'limit, so we
371 ;; need to catch it. If
373 ;; limit without ABS, we
374 ;; assume the limit is
375 ;; undefined. Is this
376 ;; right? Anyway, this
377 ;; fixes Bug 1548643.
382 (equal ($imagpart var
) 0)))
383 (cond ((eq arglim
'$infinity
)
384 ;; Check for $infinity as limit of argument.
387 (throw 'mabs
'retn
))))
388 (t (do ((ans d
(involve exp
'(mabs))) (a () ()))
390 (setq a
(mabs-subst ans var val
))
391 (setq d
(limit a var val t
))
395 (setq d
(behavior a var val
))
396 (if (zerop1 d
) (throw 'mabs
'retn
))))
399 (cond ((or (eq d
'$zeroa
) (eq d
'$inf
)
401 ;; fails on limit(abs(sin(x))/sin(x), x, inf)
402 (eq ($sign d
) '$pos
))
403 (setq exp
(maxima-substitute a
`((mabs) ,ans
) exp
)))
404 ((or (eq d
'$zerob
) (eq d
'$minf
)
405 (eq ($sign d
) '$neg
))
406 (setq exp
(maxima-substitute (m* -
1 a
) `((mabs) ,ans
) exp
)))
407 (t (throw 'mabs
'retn
))))
408 (t (throw 'mabs
'retn
))))))))))
410 ;; Called on an expression that might contain $INF, $MINF, $ZEROA, $ZEROB. Tries
411 ;; to simplify it to sort out things like inf^inf or inf+1.
413 (simpinf-ic exp
(count-general-inf exp
)))
415 (defun count-general-inf (expr)
416 (count-atoms-matching
417 (lambda (x) (or (infinityp x
) (real-epsilonp x
))) expr
))
419 (defun count-atoms-matching (predicate expr
)
420 "Count the number of atoms in the Maxima expression EXPR matching PREDICATE,
421 ignoring dummy variables and array indices."
423 ((atom expr
) (if (funcall predicate expr
) 1 0))
424 ;; Don't count atoms that occur as a limit of %integrate, %sum, %product,
426 ((member (caar expr
) dummy-variable-operators
)
427 (count-atoms-matching predicate
(cadr expr
)))
428 ;; Ignore array indices
429 ((member 'array
(car expr
)) 0)
431 for arg in
(cdr expr
)
432 summing
(count-atoms-matching predicate arg
)))))
434 (defun simpinf-ic (exp &optional infinity-count
)
436 ;; A very slow identity transformation...
439 ;; If there's only one infinity, we replace it by a variable and take the
440 ;; limit as that variable goes to infinity. Use $gensym in case we can't
441 ;; compute the answer and the limit leaks out.
442 (1 (let* ((val (or (inf-typep exp
) (epsilon-typep exp
)))
444 (expr (subst var val exp
))
445 (limit (toplevel-$limit expr var val
)))
447 ;; Now we look to see whether the computed limit is any simpler than
448 ;; what we shoved in (which we'll define as "doesn't contain EXPR as a
449 ;; subtree"). If so, return it.
450 ((not (subtree-p expr limit
:test
#'equal
))
453 ;; Otherwise, return the original form: apparently, we can't compute
454 ;; the limit we needed, and it's uglier than what we started with.
457 ;; If more than one infinity, we have to be a bit more careful.
459 (let* ((arguments (mapcar 'simpinf
(cdr exp
)))
460 (new-expression (cons (list (caar exp
)) arguments
))
463 ;; If any of the arguments are undefined, we are too.
464 ((among '$und arguments
) '$und
)
465 ;; If we ended up with something indeterminate, we punt and just return
467 ((amongl '(%limit $ind
) arguments
) exp
)
469 ;; Exponentiation & multiplication
470 ((mexptp exp
) (simpinf-expt (first arguments
) (second arguments
)))
471 ((mtimesp exp
) (simpinf-times arguments
))
473 ;; Down to at most one infinity? We do this after exponentiation to
474 ;; avoid zeroa^zeroa => 0^0, which will raise an error rather than just
475 ;; returning und. We do it after multiplication to avoid zeroa * inf =>
477 ((<= (setf infinities-left
(count-general-inf new-expression
)) 1)
478 (simpinf-ic new-expression infinities-left
))
481 ((mplusp exp
) (simpinf-plus arguments
))
484 (t new-expression
))))))
486 (defun simpinf-times (arguments)
487 (declare (special exp var val
))
488 ;; When we have a product, we need to spot that zeroa * zerob = zerob, zeroa *
489 ;; inf = und etc. Note that (SIMPINF '$ZEROA) => 0, so a nonzero atom is not
490 ;; an infinitesimal. Moreover, we can assume that each of ARGUMENTS is either
491 ;; a number, computed successfully by the recursive SIMPINF call, or maybe a
492 ;; %LIMIT noun-form (in which case, we aren't going to be able to tell the
495 ((member 0 arguments
)
497 ((find-if #'infinityp arguments
) '$und
)
498 ((every #'atom arguments
) 0)
501 ((member '$infinity arguments
)
502 (if (every #'atom arguments
)
506 (t (simplimit (cons '(mtimes) arguments
) var val
))))
508 (defun simpinf-expt (base exponent
)
509 ;; In the comments below, zero* represents one of 0, zeroa, zerob.
511 ;; TODO: In some cases we give up too early. E.g. inf^(2 + 1/inf) => inf^2
512 ;; (which should simplify to inf)
522 ((0 $zeroa $zerob
) '$und
)
523 (t (list '(mexpt) base exponent
))))
524 ;; minf^inf = infinity <== Or should it be und?
527 ;; minf^foo = minf^foo
532 ((0 $zeroa $zerob
) '$und
)
533 (t (list '(mexpt) base exponent
))))
537 ;; zero*^foo = zero*^foo
542 ((0 $zeroa $zerob
) '$und
)
543 (t (list '(mexpt) base exponent
))))
544 ;; a^b where a is pretty much anything except for a naked
545 ;; inf,minf,zeroa,zerob or 0.
548 ;; When a isn't crazy, try a^b = e^(b log(a))
549 ((not (amongl (append infinitesimals infinities
) base
))
550 (simpinf (m^
'$%e
(m* exponent
`((%log
) ,base
)))))
552 ;; No idea. Just return what we've found so far.
553 (t (list '(mexpt) base exponent
))))))
555 (defun simpinf-plus (arguments)
556 ;; We know that none of the arguments are infinitesimals, since SIMPINF never
557 ;; returns one of them. As such, we partition our arguments into infinities
558 ;; and everything else. The latter won't have any "hidden" infinities like
559 ;; lim(x,x,inf), since SIMPINF gave up on anything containing a %lim already.
560 (let ((bigs) (others))
561 (dolist (arg arguments
)
562 (cond ((infinityp arg
) (push arg bigs
))
563 (t (push arg others
))))
565 ;; inf + minf or the like
566 ((cdr (setf bigs
(delete-duplicates bigs
))) '$und
)
567 ;; inf + smaller + stuff
569 ;; I don't think this can happen, since SIMPINF goes back to the start if
570 ;; there are fewer than two infinities in the arguments, but let's be
572 (t (cons '(mplus) others
)))))
574 ;; Simplify expression with zeroa or zerob.
575 (defun simpab (small)
576 (cond ((null small
) ())
577 ((member small
'($zeroa $zerob $inf $minf $infinity
) :test
#'eq
) small
)
578 ((not (free small
'$ind
)) '$ind
) ;Not exactly right but not
579 ((not (free small
'$und
)) '$und
) ;causing trouble now.
580 ((mapatom small
) small
)
581 (t (let ((preserve-direction t
)
582 (new-small (subst (m^
'$inf -
1) '$zeroa
583 (subst (m^
'$minf -
1) '$zerob small
))))
584 (simpinf new-small
)))))
587 ;;;*I* INDICATES: T => USE LIMIT1,THINK, NIL => USE SIMPLIMIT.
588 (defun limit (exp var val
*i
*)
590 ((among '$und exp
) '$und
)
593 ((not (among var exp
))
594 (cond ((amongl '($inf $minf $infinity $ind
) exp
)
596 ((amongl '($zeroa $zerob
) exp
)
597 ;; Simplify expression with zeroa or zerob.
601 (t (putlimval exp
(cond ((and limit-using-taylor
604 (taylim exp var val
*i
*))
605 ((ratp exp var
) (ratlim exp
))
606 ((or (eq *i
* t
) (radicalp exp var
))
607 (limit1 exp var val
))
609 (cond ((or (mtimesp exp
) (mexptp exp
))
610 (limit1 exp var val
))
611 (t (simplimit exp var val
))))
612 (t (simplimit exp var val
)))))))
614 (defun limitsimp (exp var
)
615 (limitsimp-expt (sin-sq-cos-sq-sub exp
) var
))
616 ;;Hack for sin(x)^2+cos(x)^2.
618 ;; if var appears in base and power of expt,
619 ;; push var into power of of expt
620 (defun limitsimp-expt (exp var
)
621 (cond ((or (atom exp
)
623 (freeof var exp
)) exp
)
625 (not (freeof var
(cadr exp
)))
626 (not (freeof var
(caddr exp
))))
627 (m^
'$%e
(simplify `((%log
) ,exp
))))
628 (t (subst0 (cons (cons (caar exp
) ())
629 (mapcar #'(lambda (x)
630 (limitsimp-expt x var
))
634 (defun sin-sq-cos-sq-sub (exp) ;Hack ... Hack
635 (let ((arg (involve exp
'(%sin %cos
))))
638 (t (let ((new-exp ($substitute
(m+t
1 (m- (m^t
`((%sin simp
) ,arg
) 2)))
639 (m^t
`((%cos simp
) ,arg
) 2)
641 (m+t
1 (m- (m^t
`((%cos simp
) ,arg
) 2)))
642 (m^t
`((%sin simp
) ,arg
) 2)
644 (cond ((not (involve new-exp
'(%sin %cos
))) new-exp
)
647 (defun expand-trigs (x var
)
650 ((and (or (eq (caar x
) '%sin
)
652 (not (free (cadr x
) var
)))
654 ((member 'array
(car x
))
655 ;; Some kind of array reference. Return it.
657 (t (simplify (cons (ncons (caar x
))
658 (mapcar #'(lambda (x)
659 (expand-trigs x var
))
664 (cond ((not (involve e
665 '(%cot %csc %binomial
666 %sec %coth %sech %csch
667 %acot %acsc %asec %acoth
669 %jacobi_ns %jacobi_nc %jacobi_cs
670 %jacobi_ds %jacobi_dc
)))
672 (t ($ratsimp
(tansc1 e
)))))
674 (defun tansc1 (e &aux tem
)
676 ((and (setq e
(cons (car e
) (mapcar 'tansc1
(cdr e
)))) ()))
677 ((setq tem
(assoc (caar e
) '((%cot . %tan
) (%coth . %tanh
)
678 (%sec . %cos
) (%sech . %cosh
)
679 (%csc . %sin
) (%csch . %sinh
)) :test
#'eq
))
680 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
)) -
1.
)))
681 ((setq tem
(assoc (caar e
) '((%jacobi_nc . %jacobi_cn
)
682 (%jacobi_ns . %jacobi_sn
)
683 (%jacobi_cs . %jacobi_sc
)
684 (%jacobi_ds . %jacobi_sd
)
685 (%jacobi_dc . %jacobi_cd
)) :test
#'eq
))
686 ;; Converts Jacobi elliptic function to its reciprocal
688 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
) (third e
)) -
1.
)))
689 ((setq tem
(member (caar e
) '(%sinh %cosh %tanh
) :test
#'eq
))
690 (let (($exponentialize t
))
692 ((setq tem
(assoc (caar e
) '((%acsc . %asin
) (%asec . %acos
)
693 (%acot . %atan
) (%acsch . %asinh
)
694 (%asech . %acosh
) (%acoth . %atanh
)) :test
#'eq
))
695 (list (ncons (cdr tem
)) (m^t
(cadr e
) -
1.
)))
696 ((and (eq (caar e
) '%binomial
) (among var
(cdr e
)))
697 (m// `((mfactorial) ,(cadr e
))
698 (m* `((mfactorial) ,(m+t
(cadr e
) (m- (caddr e
))))
699 `((mfactorial) ,(caddr e
)))))
703 (cond ((not (involve ex
'(%sin %cos %tan %asin %acos %atan
704 %sinh %cosh %tanh %asinh %acosh %atanh
)))
710 ((eq (caar ex
) '%sinh
)
711 (m// (m+ (m^
'$%e
(cadr ex
)) (m- (m^
'$%e
(m- (cadr ex
)))))
713 ((eq (caar ex
) '%cosh
)
714 (m// (m+ (m^
'$%e
(cadr ex
)) (m^
'$%e
(m- (cadr ex
))))
716 ((and (member (caar ex
)
717 '(%sin %cos %tan %asin %acos %atan %sinh
718 %cosh %tanh %asinh %acosh %atanh
) :test
#'eq
)
721 (t (cons (car ex
) (mapcar #'hyperex0
(cdr ex
))))))
726 ;;Used by tlimit also.
727 (defun limit1 (exp var val
)
729 (let ((lhprogress? lhp?
)
732 (cond ((setq ans
(and (not (atom exp
)) (getlimval exp
)))
734 ((and (not (infinityp val
)) (setq ans
(simplimsubst val exp
)))
737 ;;;NUMDEN* => (values numerator denominator)
738 (multiple-value-bind (n dn
)
740 (cond ((not (among var dn
))
741 (return (simplimit (m// (simplimit n var val
) dn
) var val
)))
743 (return (simplimit (m* n
(simplimexpt dn -
1 (simplimit dn var val
) -
1)) var val
)))
745 (/#alike n
(car lhprogress?
))
746 (/#alike dn
(cdr lhprogress?
)))
747 (throw 'lhospital nil
)))
748 (return (limit2 n dn var val
))))))
753 (let ((deriv (sdiff (m// e f
) var
)))
755 ((=0 ($ratsimp deriv
)) t
)
758 ;; The standard /#alike function is possibly somewhat inefficient. Here is
759 ;; a possible replacement:
761 ;; If e/f is free of var (a special), return true. This code
762 ;; assumes that f is not zero. First, we test if e & f
763 ;; are alike--this check is relatively fast; second, we check
764 ;; if e/f is free of var.
766 ;;(defun /#alike (e f)
767 ;; (or (alike1 e f) (freeof var (sratsimp (div e f)))))
769 (defun limit2 (n dn var val
)
770 (prog (n1 d1 lim-sign gcp sheur-ans
)
771 (setq n
(hyperex n
) dn
(hyperex dn
))
772 ;;;Change to uniform limit call.
773 (cond ((infinityp val
)
774 (setq d1
(limit dn var val nil
))
775 (setq n1
(limit n var val nil
)))
776 (t (cond ((setq n1
(simplimsubst val n
)) nil
)
777 (t (setq n1
(limit n var val nil
))))
778 (cond ((setq d1
(simplimsubst val dn
)) nil
)
779 (t (setq d1
(limit dn var val nil
))))))
780 (cond ((or (null n1
) (null d1
)) (return nil
))
781 (t (setq n1
(sratsimp n1
) d1
(sratsimp d1
))))
782 (cond ((or (involve n
'(mfactorial)) (involve dn
'(mfactorial)))
783 (let ((ans (limfact2 n dn var val
)))
784 (cond (ans (return ans
))))))
785 (cond ((and (zerop2 n1
) (zerop2 d1
))
786 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
787 (return (colexpt n dn gcp
)))
788 ((and (real-epsilonp val
)
790 (not (free dn
'%log
)))
791 (return (liminv (m// n dn
))))
792 ((setq n1
(try-lhospital-quit n dn nil
))
794 ((and (zerop2 n1
) (not (member d1
'($ind $und
) :test
#'eq
))) (return 0))
796 (setq n1
(ridofab n1
))
797 (return (simplimtimes `(,n1
,(simplimexpt dn -
1 d1 -
1))))))
798 (setq n1
(ridofab n1
))
799 (setq d1
(ridofab d1
))
800 (cond ((or (eq d1
'$und
)
801 (and (eq n1
'$und
) (not (real-infinityp d1
))))
804 ;; At this point we have n1/$ind. Look if n1 is one of the
805 ;; infinities or zero.
806 (cond ((and (infinityp n1
) (eq ($sign dn
) '$pos
))
808 ((and (infinityp n1
) (eq ($sign dn
) '$neg
))
809 (return (simpinf (m* -
1 n1
))))
811 (or (eq ($sign dn
) '$pos
)
812 (eq ($sign dn
) '$neg
)))
815 ((eq n1
'$ind
) (return (cond ((infinityp d1
) 0)
818 ((and (real-infinityp d1
) (member n1
'($inf $und $minf
) :test
#'eq
))
819 (cond ((and (not (atom dn
)) (not (atom n
))
820 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
821 (return (colexpt n dn gcp
)))
823 (or (involve dn
'(mfactorial %gamma
))
824 (involve n
'(mfactorial %gamma
))))
825 (return (limfact n dn
))))))
826 ((eq n1 d1
) (setq lim-sign
1) (go cp
))
827 (t (setq lim-sign -
1) (go cp
))))
828 ((and (infinityp d1
) (infinityp n1
))
829 (setq lim-sign
(if (or (eq d1
'$minf
) (eq n1
'$minf
)) -
1 1))
831 (t (return (simplimtimes `(,n1
,(m^ d1 -
1))))))
832 cp
(setq n
($expand n
) dn
($expand dn
))
834 (let ((new-n (m+l
(maxi (cdr n
)))))
835 (cond ((not (alike1 new-n n
))
836 (return (limit (m// new-n dn
) var val
'think
))))
840 (let ((new-dn (m+l
(maxi (cdr dn
)))))
841 (cond ((not (alike1 new-dn dn
))
842 (return (limit (m// n new-dn
) var val
'think
))))
845 (setq sheur-ans
(sheur0 n1 d1
))
846 (cond ((or (member sheur-ans
'($inf $zeroa
) :test
#'eq
)
847 (free sheur-ans var
))
848 (return (simplimtimes `(,lim-sign
,sheur-ans
))))
849 ((and (alike1 sheur-ans dn
)
851 ((member (setq n1
(cond ((expfactorp n1 d1
) (expfactor n1 d1 var
))
853 '($inf $zeroa
) :test
#'eq
)
855 ((not (null (setq n1
(cond ((expfactorp n dn
) (expfactor n dn var
))
858 ((and (alike1 sheur-ans dn
) (not (mplusp n
))))
859 ((not (alike1 sheur-ans
(m// n dn
)))
860 (return (simplimit (m// ($expand
(m// n sheur-ans
))
861 ($expand
(m// dn sheur-ans
)))
864 (cond ((and (not (and (eq val
'$inf
) (expp n
) (expp dn
)))
865 (setq n1
(try-lhospital-quit n dn nil
))
870 ;; Test whether both n and dn have form
871 ;; product of poly^poly
872 (defun expfactorp (n dn
)
873 (do ((llist (append (cond ((mtimesp n
) (cdr n
))
875 (cond ((mtimesp dn
) (cdr dn
))
878 (exp? t
) ;IS EVERY ELEMENT SO FAR
879 (factor nil
)) ;A POLY^POLY?
883 (setq factor
(car llist
))
884 (setq exp?
(or (polyinx factor var
())
886 (polyinx (cadr factor
) var
())
887 (polyinx (caddr factor
) var
()))))))
889 (defun expfactor (n dn var
) ;Attempts to evaluate limit by grouping
890 (prog (highest-deg) ; terms with similar exponents.
891 (let ((new-exp (exppoly n
))) ;exppoly unrats expon
892 (setq n
(car new-exp
) ;and rtns deg of expons
893 highest-deg
(cdr new-exp
)))
894 (cond ((null n
) (return nil
))) ;nil means expon is not
895 (let ((new-exp (exppoly dn
))) ;a rat func.
896 (setq dn
(car new-exp
)
897 highest-deg
(max highest-deg
(cdr new-exp
))))
899 (= highest-deg
0)) ; prevent infinite recursion
903 (degree highest-deg
(1- degree
))
910 (limit (m// numerator denominator
)
914 (let ((newnumer-factor (get-newexp&factors
918 (setq numerator
(car newnumer-factor
)
919 numfactors
(cdr newnumer-factor
)))
920 (let ((newdenom-factor (get-newexp&factors
924 (setq denominator
(car newdenom-factor
)
925 denfactors
(cdr newdenom-factor
)))
926 (setq answer
(simplimit (list '(mexpt)
928 (m// numfactors denfactors
))
929 (cond ((> degree
0) var
)
933 (cond ((member answer
'($ind $und
) :test
#'equal
)
934 ;; cannot handle limit(exp(x*%i)*x, x, inf);
936 ((member answer
'($inf $minf
) :test
#'equal
)
937 ;; 0, zeroa, zerob are passed through to next iteration
938 (return (simplimtimes (list (m// numerator denominator
) answer
)))))))))
940 (defun exppoly (exp) ;RETURNS EXPRESSION WITH UNRATTED EXPONENTS
944 (exp (cond ((mtimesp exp
)
948 ((null exp
) (cons new-exp highest-deg
))
949 (setq factor
(car exp
))
951 (m* (cond ((or (not (mexptp factor
))
952 (not (ratp (caddr factor
) var
)))
956 (ratdegree (caddr factor
))))
957 (m^
(cadr factor
) (unrat (caddr factor
)))))
960 (defun unrat (exp) ;RETURNS UNRATTED EXPRESION
961 (multiple-value-bind (n d
)
963 (let ((tem ($divide n d
)))
965 (m// (caddr tem
) d
)))))
967 (defun get-newexp&factors
(exp degree var
) ;RETURNS (CONS NEWEXP FACTORS)
968 (do ((terms (cond ((mtimesp exp
)(cdr exp
)) ; SUCH THAT
969 (t (ncons exp
))) ; NEWEXP*FACTORS^(VAR^DEGREE)
970 (cdr terms
)) ; IS EQUAL TO EXP.
977 (setq factor
(car terms
))
978 (cond ((not (mexptp factor
))
980 (setq factors
(m* factor factors
)))
981 (t (setq newexp
(m* factor newexp
)))))
983 (= (ratdegree (caddr factor
))
985 (setq factors
(m* (m^
(cadr factor
)
986 (leading-coef (caddr factor
)))
988 newexp
(m* (m^
(cadr factor
)
990 (m* (leading-coef (caddr factor
))
993 (t (setq newexp
(m* factor newexp
))))))
995 (defun leading-coef (rat)
996 (ratlim (m// rat
(m^ var
(ratdegree rat
)))))
998 (defun ratdegree (rat)
999 (multiple-value-bind (n d
)
1001 (- (deg n
) (deg d
))))
1003 (defun limfact2 (n d var val
)
1004 (let ((n1 (reflect0 n var val
))
1005 (d1 (reflect0 d var val
)))
1006 (cond ((and (alike1 n n1
)
1009 (t (limit (m// n1 d1
) var val
'think
)))))
1011 ;; takes expression and returns operator at front with all flags removed
1012 ;; except array flag.
1013 ;; array flag must match for alike1 to consider two things to be the same.
1014 ;; ((MTIMES SIMP) ... ) => (MTIMES)
1015 ;; ((PSI SIMP ARRAY) 0) => (PSI ARRAY)
1016 (defun operator-with-array-flag (exp)
1017 (cond ((member 'array
(car exp
) :test
#'eq
)
1018 (list (caar exp
) 'array
))
1019 (t (list (caar exp
)))))
1021 (defun reflect0 (exp var val
)
1022 (cond ((atom exp
) exp
)
1023 ((and (eq (caar exp
) 'mfactorial
)
1024 (let ((argval (limit (cadr exp
) var val
'think
)))
1025 (or (eq argval
'$minf
)
1026 (and (numberp argval
)
1028 (reflect (cadr exp
)))
1029 (t (cons (operator-with-array-flag exp
)
1032 (reflect0 term var val
)))
1035 (defun reflect (arg)
1038 (m^
(list (ncons 'mfactorial
)
1042 (m^
(list (ncons '%sin
)
1046 (defun limfact (n d
)
1048 (setq n
(stirling0 n
)
1050 (setq ans
(toplevel-$limit
(m// n d
) var
'$inf
))
1051 (cond ((and (atom ans
)
1052 (not (member ans
'(und ind
) :test
#'eq
))) ans
)
1053 ((eq (caar ans
) '%limit
) ())
1056 ;; substitute asymptotic approximations for gamma, factorial, and
1058 (defun stirling0 (e)
1060 ((and (setq e
(cons (car e
) (mapcar 'stirling0
(cdr e
))))
1062 ((and (eq (caar e
) '%gamma
)
1063 (eq (limit (cadr e
) var val
'think
) '$inf
))
1064 (stirling (cadr e
)))
1065 ((eq (caar e
) 'mfactorial
)
1066 (let ((n (limit (cadr e
) var val
'think
)))
1068 (m* (cadr e
) (stirling (cadr e
))))
1069 ((and (integerp n
) (< n
0))
1071 (div (power -
1 n
) (mul (ftake 'mfactorial n
) (add var n
))))
1073 ((and (eq (caar e
) 'mqapply
) ;; polylogarithm
1074 (eq (subfunname e
) '$li
)
1075 (integerp (car (subfunsubs e
))))
1076 (li-asymptotic-expansion (m- (car (subfunsubs e
)) 1)
1077 (car (subfunsubs e
))
1078 (car (subfunargs e
))))
1082 (maxima-substitute x
'$z
1084 ((mexpt simp
) 2 ((rat simp
) 1 2))
1085 ((mexpt simp
) $%pi
((rat simp
) 1 2))
1086 ((mexpt simp
) $z
((mplus simp
) ((rat simp
) -
1 2) $z
))
1087 ((mexpt simp
) $%e
((mtimes simp
) -
1 $z
)))))
1089 (defun no-err-sub (v e
&aux ans
)
1090 (let ((errorsw t
) (*zexptsimp? t
)
1092 ;; Don't print any error messages
1094 (declare (special errcatch
))
1095 ;; Should we just use IGNORE-ERRORS instead HANDLER-CASE here? I
1096 ;; (rtoy) am choosing the latter so that unexpected errors will
1097 ;; actually show up instead of being silently discarded.
1099 (setq ans
(catch 'errorsw
1101 (sratsimp (subin v e
)))))
1104 (cond ((null ans
) t
) ; Ratfun package returns NIL for failure.
1107 (defun extended-real-p (e)
1108 (member e
(list '$minf
'$zerob
'$zeroa
'$ind
'$und
'$inf
'$infinity
)))
1110 ;; Evaluate the limit of e as var (special) approaches v using direct substitution.
1111 ;; This function is the last of the chain of methods tried by limit. As such it
1112 ;; shouldn't call higher level functions such as simplimit or limit--doing so
1113 ;; would risk creating an infinite loop.
1115 ;; This function special cases sums, products, and powers. It declines to use
1116 ;; direct subsitution on any expression whose main operator has a simplim%function
1117 ;; function. Generally, limits of functions that have a simplim%function should be
1118 ;; handled by those specialized functions, not by simplimsubst. Having said this,
1119 ;; it's OK for simplimsubst to special case an operator and do direct substitution
1120 ;; when its OK. The cases of log, cosine, and sine happen often enough that these
1121 ;; functions are special cased.
1123 ;; Possibly simplimsubst should decline to do direct substituion on functions
1124 ;; that have a limit function, for example inverse_jacobi_ns.
1126 ;; For all other kinds of expressions, this function assumes that if substiution
1127 ;; doesn't result in an error (division by zero, for example), the function
1128 ;; is continuous at the limit point. That's dodgy. This dodginess underscores the
1129 ;; need to define a simplim%function functions that are not continous on their
1132 ;; This function sometimes recieves unsimplifed expressions. Maybe they should be
1133 ;; simplified, maybe not.
1135 ;; Locally setting $numer and $%enumer to nil keeps some limits, for example
1136 ;; limit(sin(x)/x,x,0) from returning 1.0 when numer is true. (Barton Willis)
1138 (defun simplimsubst (v e
)
1139 (let ((ans nil
) ($numer nil
) ($%enumer nil
) (ee))
1141 ;; When e is a number, return it.
1143 ;; When e is a mapatom, substitute v for var and return.
1144 (($mapatom e
) ($substitute v var e
))
1145 ;; Special case mexpt expressions. Decline direct subsitution for
1147 ((and (mexptp e
) (not (extended-real-p v
)))
1148 (let ((x (simplimsubst v
(second e
)))
1149 (n (simplimsubst v
(third e
))))
1150 ;; Decline direct substitution (DS) for 0^negative. Also decline
1151 ;; DS when x is on the negative real axis and n isn't an integer.
1152 ;; Additionally, we require that DS is OK for both x & n.
1153 (if (and x n
(not (and (zerop2 x
) (eq t
(mgqp 0 n
)))) ; not 0^negative
1154 (or (off-negative-real-axisp x
) (integerp n
)))
1155 (ftake 'mexpt x n
) nil
)))
1156 ;; Special case product and sum expressions. Again, we decline direct
1157 ;; subsitution for extended reals.
1158 ((and (or (mplusp e
) (mtimesp e
)) (not (extended-real-p v
)))
1159 (setq ee
(mapcar #'(lambda(q) (simplimsubst v q
)) (cdr e
)))
1160 (if (some #'(lambda (q) (eq q nil
)) ee
) nil
1161 (simplifya (cons (list (caar e
)) ee
) t
)))
1162 ;; Decline direct substitution for sums, products, and powers for
1163 ;; the extended real case.
1164 ((and (or (mexptp e
) (mplusp e
) (mtimesp e
)) (extended-real-p v
))
1166 ;; Special case log expressions--possibly the log case happens
1167 ;; often enough to make this worthwile.
1168 ((and (consp e
) (consp (car e
)) (eq '%log
(caar e
)))
1169 (let ((w (simplimsubst v
(cadr e
))))
1170 (if (and w
(off-negative-real-axisp w
)) (ftake '%log w
) nil
)))
1171 ;; Special case %cos and %sin expressions--we could special case others.
1172 ;; The lenient-realp check declines direct substitution for complex arguments--
1173 ;; this check could be relaxed.
1174 ((and (consp e
) (consp (car e
)) (or (member (caar e
) (list '%cos
'%sin
))))
1175 (let ((op (caar e
)))
1176 (setq e
(simplimsubst v
(second e
)))
1177 (if (and e
(lenient-realp e
) (not (extended-real-p e
))) (ftake op e
) nil
)))
1178 ;; Don't use direct substitution on expressions whose main operator has
1179 ;; a simplim%function.
1180 ((and (consp e
) (consp (car e
)) (get (caar e
) 'simplim%function
)) nil
)
1181 ;; The function no-err-sub returns true when there is an error
1182 ((not (eq t
(setq ans
(no-err-sub (ridofab v
) e
))))
1183 ;; Previously the condition (zerop2 ans) was (=0 ($radcan ans)). In
1184 ;; December 2021, the testsuite + the share testsuite only gives ans = 0,
1185 ;; making the radcan unneeded.
1186 (cond ((and (member v
'($zeroa $zerob
) :test
#'eq
) (zerop2 ans
))
1187 (setq ans
(behavior e var v
))
1188 (cond ((eql ans
1) '$zeroa
)
1189 ((eql ans -
1) '$zerob
)
1190 (t nil
))) ;behavior can't find direction
1192 ;; direct substitution fails, so return nil.
1195 ;;;returns (cons numerator denominator)
1197 (let ((e (factor (simplify e
)))
1203 (mapc #'forq
(cdr e
)))
1209 (setq numer
(car numer
)))
1211 (setq numer
(m*l numer
))))
1215 (setq denom
(car denom
)))
1217 (setq denom
(m*l denom
))))
1218 (values (factor numer
) (factor denom
))))
1220 ;;;FACTOR OR QUOTIENT
1221 ;;;Setq's the special vars numer and denom from numden*
1223 (cond ((and (mexptp e
)
1224 (not (freeof var e
))
1225 (null (pos-neg-p (caddr e
))))
1226 (push (m^
(cadr e
) (m* -
1.
(caddr e
))) denom
))
1227 (t (push e numer
))))
1229 ;;;Predicate to tell whether an expression is pos,zero or neg as var -> val.
1230 ;;;returns T if pos,zero. () if negative or don't know.
1231 (defun pos-neg-p (exp)
1232 (let ((ans (limit exp var val
'think
)))
1233 (cond ((and (not (member ans
'($und $ind $infinity
) :test
#'eq
))
1234 (equal ($imagpart ans
) 0))
1235 (let ((sign (getsignl ans
)))
1236 (cond ((or (equal sign
1)
1239 ((equal sign -
1) nil
))))
1242 (declare-top (unspecial n dn
))
1245 (cond ((radicalp e var
) nil
)
1246 ((member (caar e
) '(%log %sin %cos %tan %sinh %cosh %tanh mfactorial
1247 %asin %acos %atan %asinh %acosh %atanh
) :test
#'eq
) nil
)
1249 ((do ((e (cdr e
) (cdr e
)))
1251 (and (expp (car e
)) (return t
))))))
1255 (radicalp (cadr e
) var
)
1256 (among var
(caddr e
))
1257 (radicalp (caddr e
) var
)))
1260 (defun gcpower (a b
)
1261 ($gcd
(getexp a
) (getexp b
)))
1264 (cond ((and (mexptp exp
)
1265 (free (caddr exp
) var
)
1266 (eq (ask-integer (caddr exp
) '$integer
) '$yes
))
1268 ((mtimesp exp
) (getexplist (cdr exp
)))
1271 (defun getexplist (list)
1272 (cond ((null (cdr list
))
1273 (getexp (car list
)))
1274 (t ($gcd
(getexp (car list
))
1275 (getexplist (cdr list
))))))
1277 (defun limroot (exp power
)
1278 (cond ((or (atom exp
) (not (member (caar exp
) '(mtimes mexpt
) :test
#'eq
)))
1279 (limroot (list '(mexpt) exp
1) power
)) ;This is strange-JIM.
1280 ((mexptp exp
) (m^
(cadr exp
)
1281 (sratsimp (m* (caddr exp
) (m^ power -
1.
)))))
1282 (t (m*l
(mapcar #'(lambda (x)
1286 ;;NUMERATOR AND DENOMINATOR HAVE EXPONENTS WITH GCD OF GCP.
1287 ;;; Used to call simplimit but some of the transformations used here
1288 ;;; were not stable w.r.t. the simplifier, so try keeping exponent separate
1291 (defun colexpt (n dn gcp
)
1292 (let ((bas (m* (limroot n gcp
) (limroot dn
(m* -
1 gcp
))))
1295 (setq baslim
(limit bas var val
'think
))
1296 (setq expolim
(limit expo var val
'think
))
1297 (simplimexpt bas expo baslim expolim
)))
1299 ;; this function is responsible for the following bug:
1300 ;; limit(x^2 + %i*x, x, inf) -> inf (should be infinity)
1302 (cond ((member val
'($inf $infinity
) :test
#'eq
)
1303 (setq e
(maxima-substitute (m^t
'x -
1) var e
)))
1305 (setq e
(maxima-substitute (m^t -
1 (m^t
'x -
1)) var e
)))
1307 (setq e
(maxima-substitute (m- 'x
) var e
)))
1309 (setq e
(maxima-substitute 'x var e
)))
1310 ((setq e
(maxima-substitute (m+t
'x val
) var e
))))
1311 (destructuring-let* ((e (let (($ratfac
()))
1312 ($rat
(sratsimp e
) 'x
)))
1320 (sratsimp (m// ($ratdisrep
`(,h
,(locoef n g
) .
1))
1321 ($ratdisrep
`(,h
,(locoef d g
) .
1))))))
1323 (cond ((not (member val
'($zerob $zeroa $inf $minf
) :test
#'eq
))
1325 ((not (equal ($imagpart e
) 0))
1327 ((null (setq e
(getsignl ($realpart e
))))
1329 ((equal e
1) '$zeroa
)
1330 ((equal e -
1) '$zerob
)
1333 ((not (member val
'($zerob $zeroa $infinity $inf $minf
) :test
#'eq
))
1335 ((eq val
'$infinity
) '$infinity
)
1336 ((not (equal ($imagpart e
) 0)) '$infinity
)
1337 ((null (setq e
(getsignl ($realpart e
))))
1340 ((equal e -
1) '$minf
)
1344 (if (or (atom n
) (not (eq (car n
) x
)))
1349 (if (or (atom n
) (not (eq (car n
) x
)))
1353 (defun behavior (exp var val
) ; returns either -1, 0, 1.
1354 (if (= *behavior-count-now
* +behavior-count
+)
1356 (let ((*behavior-count-now
* (1+ *behavior-count-now
*)) pair sign
)
1357 (cond ((real-infinityp val
)
1358 (setq val
(cond ((eq val
'$inf
) '$zeroa
)
1359 ((eq val
'$minf
) '$zerob
)))
1360 (setq exp
(sratsimp (subin (m^ var -
1) exp
)))))
1361 (cond ((eq val
'$infinity
) 0) ; Needs more hacking for complex.
1363 (prog2 (setq pair
(partition exp var
1))
1364 (not (mtimesp (cdr pair
)))))
1365 (setq sign
(getsignl (car pair
)))
1366 (if (not (fixnump sign
))
1368 (mul sign
(behavior (cdr pair
) var val
))))
1369 ((and (=0 (no-err-sub (ridofab val
) exp
))
1371 (free (caddr exp
) var
)
1372 (equal (getsignl (caddr exp
)) 1))
1373 (let ((bas (cadr exp
)) (expo (caddr exp
)))
1374 (behavior-expt bas expo
)))
1375 (t (behavior-by-diff exp var val
))))))
1377 (defun behavior-expt (bas expo
)
1378 (let ((behavior (behavior bas var val
)))
1379 (cond ((= behavior
1) 1)
1381 ((eq (ask-integer expo
'$integer
) '$yes
)
1382 (cond ((eq (ask-integer expo
'$even
) '$yes
) 1)
1385 (cond ((evenp (cadr expo
)) 1)
1386 ((oddp (caddr expo
)) behavior
)
1390 (defun behavior-by-diff (exp var val
)
1391 (cond ((not (or (eq val
'$zeroa
) (eq val
'$zerob
))) 0)
1392 (t (let ((old-val val
) (old-exp exp
))
1393 (setq val
(ridofab val
))
1395 (exp (sratsimp (sdiff exp var
)) (sratsimp (sdiff exp var
)))
1397 (ans ())) ; This do wins by a return.
1398 ((> ct
1) 0) ; This loop used to run up to 5 times,
1399 ;; but the size of some expressions would blow up.
1400 (setq ans
(no-err-sub val exp
)) ;Why not do an EVENFN and ODDFN
1403 (return (behavior-numden old-exp var old-val
)))
1404 ((=0 ans
) ()) ;Do it again.
1405 (t (setq ans
(getsignl ans
))
1406 (cond (n (return ans
))
1408 (return (if (eq old-val
'$zeroa
) 1 -
1)))
1410 (return (if (eq old-val
'$zeroa
) -
1 1)))
1411 (t (return 0))))))))))
1413 (defun behavior-numden (exp var val
)
1414 (let ((num ($num exp
)) (denom ($denom exp
)))
1415 (cond ((equal denom
1) 0) ;Could be hacked more from here.
1416 (t (let ((num-behav (behavior num var val
))
1417 (denom-behav (behavior denom var val
)))
1418 (cond ((or (= num-behav
0) (= denom-behav
0)) 0)
1419 ((= num-behav denom-behav
) 1)
1422 (defun try-lhospital (n d ind
)
1423 ;;Make one catch for the whole bunch of lhospital trials.
1424 (let ((ans (lhospital-catch n d ind
)))
1425 (cond ((null ans
) ())
1426 ((not (free-infp ans
)) (simpinf ans
))
1427 ((not (free-epsilonp ans
)) (simpab ans
))
1430 (defun try-lhospital-quit (n d ind
)
1431 (let ((ans (or (lhospital-catch n d ind
)
1432 (lhospital-catch (m^ d -
1) (m^ n -
1) ind
))))
1433 (cond ((null ans
) (throw 'limit t
))
1434 ((not (free-infp ans
)) (simpinf ans
))
1435 ((not (free-epsilonp ans
)) (simpab ans
))
1438 (defun lhospital-catch (n d ind
)
1439 (cond ((> 0 lhcount
)
1440 (setq lhcount $lhospitallim
)
1441 (throw 'lhospital nil
))
1442 ((equal lhcount $lhospitallim
)
1443 (let ((lhcount (m+ lhcount -
1)))
1444 (catch 'lhospital
(lhospital n d ind
))))
1445 (t (setq lhcount
(m+ lhcount -
1))
1446 (prog1 (lhospital n d ind
)
1447 (setq lhcount
(m+ lhcount
1))))))
1448 ;;If this succeeds then raise LHCOUNT.
1451 (defun lhospital (n d ind
)
1452 (declare (special val lhp?
))
1454 (setq n
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr n
)))))
1456 (setq d
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr d
)))))
1457 (multiple-value-bind (n d
)
1459 (let (const nconst dconst
)
1460 (setq lhp?
(and (null ind
) (cons n d
)))
1461 (multiple-value-setq (nconst n
) (var-or-const n
))
1462 (multiple-value-setq (dconst d
) (var-or-const d
))
1464 (setq n
(stirling0 n
)) ;; replace factorial and %gamma
1465 (setq d
(stirling0 d
)) ;; with approximations
1467 (setq n
(sdiff n var
) ;; take derivatives for l'hospital
1470 (if (or (not (free n
'%derivative
)) (not (free d
'%derivative
)))
1471 (throw 'lhospital
()))
1472 (setq n
(expand-trigs (tansc n
) var
))
1473 (setq d
(expand-trigs (tansc d
) var
))
1475 (multiple-value-setq (const n d
) (remove-singularities n d
))
1476 (setq const
(m* const
(m// nconst dconst
)))
1477 (simpinf (let ((ans (if ind
1478 (limit2 n d var val
)
1479 (limit-numden n d val
))))
1480 ;; When the limit function returns, it's possible that it will return NIL
1481 ;; (gave up without finding a limit). It's also possible that it will
1482 ;; return something containing UND. We treat that as a failure too.
1483 (when (and ans
(freeof '$und ans
))
1484 (m* const ans
)))))))
1486 ;; Try to compute the limit of a quotient NUM/DEN, trying to massage the input
1487 ;; into a convenient form for LIMIT on the way.
1488 (defun limit-numden (n d val
)
1490 ;; For general arguments, the best approach seems to be to use
1491 ;; sratsimp to simplify the quotient as much as we can, then
1492 ;; $multthru, which splits it up into a sum (presumably useful
1493 ;; because limit(a+b) = limit(a) + limit(b) if the limits exist, and
1494 ;; the right hand side might be easier to calculate)
1496 ($multthru
(sratsimp (m// n d
))))
1498 ;; If we've already got a sum in the numerator, it seems to be
1499 ;; better not to recombine it. Call LIMIT on the whole lot, though,
1500 ;; because terms with infinite limits might cancel to give a finite
1503 (m+l
(mapcar #'(lambda (x)
1504 (sratsimp (m// x d
)))
1507 (limit expr var val
'think
)))
1509 ;; Heuristics for picking the right way to express a LHOSPITAL problem.
1510 (defun lhop-numden (num denom
)
1511 (declare (special var
))
1512 (cond ((let ((log-num (involve num
'(%log
))))
1513 (cond ((null log-num
) ())
1514 ((lessthan (num-of-logs (factor (sratsimp (sdiff (m^ num -
1) var
))))
1515 (num-of-logs (factor (sratsimp (sdiff num var
)))))
1516 (psetq num
(m^ denom -
1) denom
(m^ num -
1))
1519 ((let ((log-denom (involve denom
'(%log
))))
1520 (cond ((null log-denom
) ())
1521 ((lessthan (num-of-logs (sratsimp (sdiff (m^ denom -
1) var
)))
1522 (num-of-logs (sratsimp (sdiff denom var
))))
1523 (psetq denom
(m^ num -
1) num
(m^ denom -
1))
1526 ((let ((exp-num (%einvolve num
)))
1528 (cond ((%e-right-placep exp-num
)
1530 (t (psetq num
(m^ denom -
1)
1531 denom
(m^ num -
1)) t
)))
1533 ((let ((exp-den (%einvolve denom
)))
1535 (cond ((%e-right-placep exp-den
)
1537 (t (psetq num
(m^ denom -
1)
1538 denom
(m^ num -
1)) t
)))
1540 ((let ((scnum (involve num
'(%sin
))))
1541 (cond (scnum (cond ((trig-right-placep '%sin scnum
) t
)
1542 (t (psetq num
(m^ denom -
1)
1543 denom
(m^ num -
1)) t
)))
1545 ((let ((scden (involve denom
'(%sin
))))
1546 (cond (scden (cond ((trig-right-placep '%sin scden
) t
)
1547 (t (psetq num
(m^ denom -
1)
1548 denom
(m^ num -
1)) t
)))
1550 ((let ((scnum (involve num
'(%asin %acos %atan
))))
1551 ;; If the numerator contains an inverse trig and the
1552 ;; denominator or reciprocal of denominator is polynomial,
1553 ;; leave everything as is. If the inverse trig is moved to
1554 ;; the denominator, things get messy, even if the numerator
1555 ;; becomes a polynomial. This is not perfect.
1556 (cond ((and scnum
(or (polyinx denom var
())
1557 (polyinx (m^ denom -
1) var
())))
1560 ((or (oscip num
) (oscip denom
)))
1562 (psetq num
(m^ denom -
1) denom
(m^ num -
1))))
1565 ;;i don't know what to do here for some cases, may have to be refined.
1566 (defun num-of-logs (exp)
1567 (cond ((mapatom exp
) 0)
1568 ((equal (caar exp
) '%log
)
1569 (m+ 1 (num-of-log-l (cdr exp
))))
1570 ((and (mexptp exp
) (mnump (caddr exp
)))
1571 (m* (simplify `((mabs) ,(caddr exp
)))
1572 (num-of-logs (cadr exp
))))
1573 (t (num-of-log-l (cdr exp
)))))
1575 (defun num-of-log-l (llist)
1576 (do ((temp llist
(cdr temp
)) (ans 0))
1578 (setq ans
(m+ ans
(num-of-logs (car temp
))))))
1580 (defun %e-right-placep
(%e-arg
)
1581 (let ((%e-arg-diff
(sdiff %e-arg var
)))
1583 ((free %e-arg-diff var
)) ;simple cases
1584 ((or (and (mexptp denom
)
1585 (equal (cadr denom
) -
1))
1586 (polyinx (m^ denom -
1) var
())) ())
1587 ((let ((%e-arg-diff-lim
(ridofab (limit %e-arg-diff var val
'think
)))
1588 (%e-arg-exp-lim
(ridofab (limit (m^
'$%e %e-arg
) var val
'think
))))
1591 (format t
"%e-arg-dif-lim = ~A~%" %e-arg-diff-lim
)
1592 (format t
"%e-arg-exp-lim = ~A~%" %e-arg-exp-lim
))
1593 (cond ((equal %e-arg-diff-lim %e-arg-exp-lim
)
1595 ((and (mnump %e-arg-diff-lim
) (mnump %e-arg-exp-lim
))
1597 ((and (mnump %e-arg-diff-lim
) (infinityp %e-arg-exp-lim
))
1598 ;; This is meant to make maxima handle bug 1469411
1599 ;; correctly. Undoubtedly, this needs work.
1603 (defun trig-right-placep (trig-type arg
)
1604 (let ((arglim (ridofab (limit arg var val
'think
)))
1605 (triglim (ridofab (limit `((,trig-type
) ,arg
) var val
'think
))))
1606 (cond ((and (equal arglim
0) (equal triglim
0)) t
)
1607 ((and (infinityp arglim
) (infinityp triglim
)) t
)
1610 ;;Takes a numerator and a denominator. If they tries all combinations of
1611 ;;products to try and make a simpler set of subproblems for LHOSPITAL.
1612 (defun remove-singularities (numer denom
)
1613 (cond ((or (null numer
) (null denom
)
1614 (atom numer
) (atom denom
)
1615 (not (mtimesp numer
)) ;Leave this here for a while.
1616 (not (mtimesp denom
)))
1617 (values 1 numer denom
))
1620 (multiple-value-bind (num-consts num-vars
)
1621 (var-or-const numer
)
1622 (multiple-value-bind (denom-consts denom-vars
)
1623 (var-or-const denom
)
1624 (if (not (mtimesp num-vars
))
1625 (setq num-vars
(list num-vars
))
1626 (setq num-vars
(cdr num-vars
)))
1627 (if (not (mtimesp denom-vars
))
1628 (setq denom-vars
(list denom-vars
))
1629 (setq denom-vars
(cdr denom-vars
)))
1630 (do ((nl num-vars
(cdr nl
))
1631 (num-list (copy-list num-vars
))
1632 (den-list denom-vars den-list-temp
)
1633 (den-list-temp (copy-list denom-vars
)))
1634 ((null nl
) (values (m* const
(m// num-consts denom-consts
))
1636 (m*l den-list-temp
)))
1637 (do ((dl den-list
(cdr dl
)))
1639 (if (or (%einvolve
(car nl
)) (%einvolve
(car nl
)))
1641 (let ((lim (catch 'limit
(simpinf (simpab (limit (m// (car nl
) (car dl
))
1642 var val
'think
))))))
1643 (cond ((or (eq lim t
)
1645 (equal (ridofab lim
) 0)
1647 (not (free lim
'$inf
))
1648 (not (free lim
'$minf
))
1649 (not (free lim
'$infinity
))
1650 (not (free lim
'$ind
))
1651 (not (free lim
'$und
)))
1654 (setq const
(m* lim const
))
1655 (setq num-list
(delete (car nl
) num-list
:count
1 :test
#'equal
))
1656 (setq den-list-temp
(delete (car dl
) den-list-temp
:count
1 :test
#'equal
))
1657 (return t
)))))))))))))
1659 ;; separate terms that contain var from constant terms
1660 ;; returns (const-terms . var-terms)
1661 (defun var-or-const (expr)
1662 (setq expr
($factor expr
))
1670 (do ((l (cdr expr
) (cdr l
))
1673 ((null l
) (values const varl
))
1674 (if (free (car l
) var
)
1675 (setq const
(m* (car l
) const
))
1676 (setq varl
(m* (car l
) varl
)))))
1680 ;; if term goes to non-zero constant, replace with constant
1681 (defun lhsimp (term var val
)
1682 (cond ((atom term
) term
)
1684 (let ((term-value (ridofab (limit term var val
'think
))))
1685 (cond ((not (member term-value
1686 '($inf $minf $und $ind $infinity
0)))
1690 (defun bylog (expo bas
)
1693 (try-lhospital-quit (simplify `((%log
) ,(tansc bas
)))
1698 (defun simplimexpt (bas expo bl el
)
1699 (cond ((or (eq bl
'$und
) (eq el
'$und
)) '$und
)
1701 (cond ((eq el
'$inf
) (if (eq bl
'$zeroa
) bl
0))
1702 ((eq el
'$minf
) (if (eq bl
'$zeroa
) '$inf
'$infinity
))
1703 ((eq el
'$ind
) '$ind
)
1704 ((eq el
'$infinity
) '$und
)
1705 ((zerop2 el
) (bylog expo bas
))
1706 (t (cond ((equal (getsignl el
) -
1)
1707 (cond ((eq bl
'$zeroa
) '$inf
)
1709 (cond ((even1 el
) '$inf
)
1710 ((eq (ask-integer el
'$integer
) '$yes
)
1711 (if (eq (ask-integer el
'$even
) '$yes
)
1713 '$minf
)))) ;Gotta be ODD.
1714 (t (setq bas
(behavior bas var val
))
1715 (cond ((equal bas
1) '$inf
)
1716 ((equal bas -
1) '$minf
)
1717 (t (throw 'limit t
))))))
1719 (member bl
'($zeroa $zerob
) :test
#'eq
))
1720 (cond ((even1 el
) '$zeroa
)
1721 ((and (eq bl
'$zerob
)
1723 (evenp (caddr el
))) 0)
1725 ((and (equal (getsignl el
) 1)
1726 (eq bl
'$zeroa
)) bl
)
1727 ((equal (getsignl el
) 0)
1731 (cond ((zerop2 el
) (bylog expo bas
))
1733 ((eq el
'$inf
) '$infinity
)
1734 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1735 ((equal (setq el
(getsignl el
)) 1) '$infinity
)
1738 (t (throw 'limit t
))))
1740 (cond ((eq el
'$inf
) '$inf
)
1741 ((equal el
'$minf
) 0)
1742 ((zerop2 el
) (bylog expo bas
))
1743 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1744 (t (cond ((eql 0 (getsignl el
)) 1)
1745 ((ratgreaterp 0 el
) '$zeroa
)
1746 ((ratgreaterp el
0) '$inf
)
1747 (t (throw 'limit t
))))))
1749 (cond ((zerop2 el
) (bylog expo bas
))
1750 ((eq el
'$inf
) '$und
)
1751 ((equal el
'$minf
) 0)
1752 ;;;Why not generalize this. We can ask about the number. -Jim 2/23/81
1753 ((mnump el
) (cond ((mnegp el
)
1756 (if (eq (ask-integer el
'$integer
) '$yes
)
1757 (if (eq (ask-integer el
'$even
) '$yes
)
1761 (t (cond ((even1 el
) '$inf
)
1762 ((eq (ask-integer el
'$integer
) '$yes
)
1763 (if (eq (ask-integer el
'$even
) '$yes
)
1767 (loginprod?
(throw 'lip?
'lip
!))
1769 ((equal (simplify (ratdisrep (ridofab bl
))) 1)
1770 (if (infinityp el
) (bylog expo bas
) 1))
1771 ((and (equal (ridofab bl
) -
1)
1772 (infinityp el
)) '$ind
) ;LB
1773 ((eq bl
'$ind
) (cond ((or (zerop2 el
) (infinityp el
)) '$und
)
1774 ((not (equal (getsignl el
) -
1)) '$ind
)
1776 ((eq el
'$inf
) (cond ((abeq1 bl
)
1777 (if (equal (getsignl bl
) 1) 1 '$ind
))
1779 (if (equal (getsignl bl
) 1) '$zeroa
0))
1780 ((equal (getsignl (m1- bl
)) 1) '$inf
)
1781 ((equal (getsignl (m1- `((mabs) ,bl
))) 1) '$infinity
)
1782 (t (throw 'limit t
))))
1783 ((eq el
'$minf
) (cond ((abeq1 bl
)
1784 (if (equal (getsignl bl
) 1) 1 '$ind
))
1786 (if (equal (getsignl bl
) 1) '$zeroa
0))
1787 ((ratgreaterp 0 bl
) '$infinity
)
1790 (if (equal val
'$infinity
)
1791 '$und
;Not enough info to do anything.
1792 (destructuring-bind (real-el . imag-el
)
1794 (setq real-el
(limit real-el var origval nil
))
1795 (cond ((eq real-el
'$minf
)
1797 ((and (eq real-el
'$inf
)
1798 (not (equal (ridofab (limit imag-el var origval nil
)) 0)))
1800 ((eq real-el
'$infinity
)
1801 (throw 'limit t
)) ;; don't really know real component
1805 ((eq el
'$ind
) '$ind
)
1810 (cond ((numberp x
) (and (integerp x
) (evenp x
)))
1811 ((and (mnump x
) (evenp (cadr x
))))))
1813 ;; is absolute value less than one?
1817 (and (ratgreaterp 1. bl
) (ratgreaterp bl -
1.
)))
1818 (t (equal (getsignl (m1- `((mabs) ,bl
))) -
1.
))))
1820 ;; is absolute value equal to one?
1824 (or (equal 1. bl
) (equal bl -
1.
)))
1825 (t (equal (getsignl (m1- `((mabs) ,bl
))) 0))))
1827 (defun simplimit (exp var val
&aux op
)
1830 ((or (atom exp
) (mnump exp
)) exp
)
1831 ;; Lookup and dispatch a simplim%function from the property list
1832 ((setq op
(safe-get (mop exp
) 'simplim%function
))
1833 (funcall op exp var val
))
1835 ;; And do the same for subscripted:
1836 ((and (consp exp
) (consp (car exp
)) (eq (caar exp
) 'mqapply
)
1837 (setq op
(safe-get (subfunname exp
) 'simplim%function
)))
1838 (funcall op exp var val
))
1840 ;; Without the call to rootscontract,
1841 ;; limit(((-4)*x^2-10*x+24)/((4*x+8)^(1/3)+2),x,-4)
1842 ;; returns zero (should be 66). This bug is due to the fact that
1843 ;; 4^(1/3)*2^(1/3)-2 does not simplify to zero. Although the call
1844 ;; to rootscontract takes care of this case, almost surely there are
1845 ;; many other limit problems that need more than rootscontract.
1846 ((mplusp exp
) (let (($rootsconmode nil
)) ($rootscontract
(simplimplus exp
))))
1847 ((mtimesp exp
) (simplimtimes (cdr exp
)))
1848 ((mexptp exp
) (simplimexpt (cadr exp
) (caddr exp
)
1849 (limit (cadr exp
) var val
'think
)
1850 (limit (caddr exp
) var val
'think
)))
1851 ((member (caar exp
) '(%sin %cos
) :test
#'eq
)
1852 (simplimsc exp
(caar exp
) (limit (cadr exp
) var val
'think
)))
1853 ((eq (caar exp
) '%tan
) (simplim%tan
(cadr exp
)))
1854 ((member (caar exp
) '(%sinh %cosh
) :test
#'eq
)
1855 (simplimsch (caar exp
) (limit (cadr exp
) var val
'think
)))
1856 ((member (caar exp
) '(%erf %tanh
) :test
#'eq
)
1857 (simplim%erf-%tanh
(caar exp
) (cadr exp
)))
1858 ((eq (caar exp
) '%atanh
)
1859 (simplim%atanh
(limit (cadr exp
) var val
'think
) val
))
1860 ((eq (caar exp
) '%acosh
)
1861 (simplim%acosh
(limit (cadr exp
) var val
'think
)))
1862 ((eq (caar exp
) '%asinh
)
1863 (simplim%asinh
(limit (cadr exp
) var val
'think
)))
1864 ((eq (caar exp
) '%inverse_jacobi_ns
)
1865 (simplim%inverse_jacobi_ns
(limit (cadr exp
) var val
'think
) (third exp
)))
1866 ((eq (caar exp
) '%inverse_jacobi_nc
)
1867 (simplim%inverse_jacobi_nc
(limit (cadr exp
) var val
'think
) (third exp
)))
1868 ((eq (caar exp
) '%inverse_jacobi_sc
)
1869 (simplim%inverse_jacobi_sc
(limit (cadr exp
) var val
'think
) (third exp
)))
1870 ((eq (caar exp
) '%inverse_jacobi_cs
)
1871 (simplim%inverse_jacobi_cs
(limit (cadr exp
) var val
'think
) (third exp
)))
1872 ((eq (caar exp
) '%inverse_jacobi_dc
)
1873 (simplim%inverse_jacobi_dc
(limit (cadr exp
) var val
'think
) (third exp
)))
1874 ((eq (caar exp
) '%inverse_jacobi_ds
)
1875 (simplim%inverse_jacobi_ds
(limit (cadr exp
) var val
'think
) (third exp
)))
1876 ((and (eq (caar exp
) 'mqapply
)
1877 (eq (subfunname exp
) '$psi
))
1878 (simplim$psi
(subfunsubs exp
) (subfunargs exp
) val
))
1879 ((and (eq (caar exp
) var
)
1880 (member 'array
(car exp
) :test
#'eq
)
1881 (every #'(lambda (sub-exp)
1884 exp
) ;LIMIT(B[I],B,INF); -> B[I]
1885 ;; When limsubst is true, limit(f(n+1)/f(n),n,inf) = 1. The user documentation
1886 ;; warns against setting limsubst to true.
1888 (simplify (cons (operator-with-array-flag exp
)
1889 (mapcar #'(lambda (a)
1890 (limit a var val
'think
))
1892 (throw 'limit t
)))))
1895 (setq e
(resimplify (subst (m// 1 var
) var e
)))
1896 (let ((new-val (cond ((eq val
'$zeroa
) '$inf
)
1897 ((eq val
'$zerob
) '$minf
))))
1898 (if new-val
(let ((preserve-direction t
))
1899 (toplevel-$limit e var new-val
)) (throw 'limit t
))))
1901 (defun simplimtimes (exp)
1902 ;; The following test
1903 ;; handles (-1)^x * 2^x => (-2)^x => $infinity
1904 ;; wants to avoid (-1)^x * 2^x => $ind * $inf => $und
1906 (and (expfactorp (cons '(mtimes) exp
) 1)
1907 (expfactor (cons '(mtimes) exp
) 1 var
))))
1908 (when try
(return-from simplimtimes try
)))
1910 (let ((prod 1) (num 1) (denom 1)
1911 (zf nil
) (ind-flag nil
) (inf-type nil
)
1912 (constant-zero nil
) (constant-infty nil
))
1914 (let* ((loginprod?
(involve term
'(%log
)))
1915 (y (catch 'lip?
(limit term var val
'think
))))
1917 ;; limit failed due to log in product
1919 (return-from simplimtimes
(liminv (cons '(mtimes simp
) exp
))))
1921 ;; If the limit is infinitesimal or zero
1923 (setf num
(m* num term
)
1924 constant-zero
(or constant-zero
(not (among var term
))))
1927 (unless zf
(setf zf
1)))
1929 (setf zf
(* -
1 (or zf
1))))))
1931 ;; If the limit is not some form of infinity or
1932 ;; undefined/indeterminate.
1933 ((not (member y
'($inf $minf $infinity $ind $und
) :test
#'eq
))
1934 (setq prod
(m* prod y
)))
1936 ((eq y
'$und
) (return-from simplimtimes
'$und
))
1937 ((eq y
'$ind
) (setq ind-flag t
))
1939 ;; Some form of infinity
1941 (setf denom
(m* denom term
)
1942 constant-infty
(or constant-infty
(not (among var term
))))
1943 (unless (eq inf-type
'$infinity
)
1945 ((eq y
'$infinity
) (setq inf-type
'$infinity
))
1946 ((null inf-type
) (setf inf-type y
))
1947 ;; minf * minf or inf * inf
1948 ((eq y inf-type
) (setf inf-type
'$inf
))
1950 (t (setf inf-type
'$minf
))))))))
1953 ;; If there are zeros and infinities among the terms that are free of
1954 ;; VAR, then we have an expression like "inf * zeroa * f(x)" or
1955 ;; similar. That gives an undefined result. Note that we don't
1956 ;; necessarily have something undefined if only the zeros have a term
1957 ;; free of VAR. For example "zeroa * exp(-1/x) * 1/x" as x -> 0. And
1958 ;; similarly for the infinities.
1959 ((and constant-zero constant-infty
) '$und
)
1961 ;; If num=denom=1, we didn't find any explicit infinities or zeros, so we
1962 ;; either return the simplified product or ind
1963 ((and (eql num
1) (eql denom
1))
1964 (if ind-flag
'$ind prod
))
1965 ;; If denom=1 (and so num != 1), we have some form of zero
1969 (let ((sign (getsignl prod
)))
1970 (if (or (not sign
) (eq sign
'complex
))
1976 ;; If num=1 (and so denom != 1), we have some form of infinity
1978 (let ((sign ($csign prod
)))
1981 ((eq sign
'$pos
) inf-type
)
1982 ((eq sign
'$neg
) (case inf-type
1986 ((member sign
'($complex $imaginary
)) '$infinity
)
1987 ; sign is '$zero, $pnz, $pz, etc
1988 (t (throw 'limit t
)))))
1989 ;; Both zeros and infinities
1991 ;; All bets off if there are some infinities or some zeros, but it
1992 ;; needn't be undefined (see above)
1993 (when (or constant-zero constant-infty
) (throw 'limit t
))
1995 (let ((ans (limit2 num
(m^ denom -
1) var val
)))
1997 (simplimtimes (list prod ans
))
1998 (throw 'limit t
)))))))
2000 ;;;PUT CODE HERE TO ELIMINATE FAKE SINGULARITIES??
2002 (defun simplimplus (exp)
2003 (cond ((memalike exp simplimplus-problems
)
2006 (progn (push exp simplimplus-problems
)
2007 (let ((ans (catch 'limit
(simplimplus1 exp
))))
2008 (cond ((or (eq ans
())
2010 (among '%limit ans
))
2011 (let ((new-exp (sratsimp exp
)))
2012 (cond ((not (alike1 exp new-exp
))
2014 (limit new-exp var val
'think
))))
2015 (cond ((or (eq ans
())
2020 (pop simplimplus-problems
)))))
2022 (defun simplimplus1 (exp)
2023 (prog (sum y infl infinityl minfl indl
)
2025 (do ((exp (cdr exp
) (cdr exp
)) (f))
2026 ((or y
(null exp
)) nil
)
2027 (setq f
(limit (car exp
) var val
'think
))
2030 ((eq f
'$und
) (setq y t
))
2031 ((not (member f
'($inf $minf $infinity $ind
) :test
#'eq
))
2032 (setq sum
(m+ sum f
)))
2033 ((eq f
'$ind
) (push (car exp
) indl
))
2034 (infinityl (throw 'limit t
))
2035 ;;;Don't know what to do with an '$infinity and an $inf or $minf
2036 ((eq f
'$inf
) (push (car exp
) infl
))
2037 ((eq f
'$minf
) (push (car exp
) minfl
))
2039 (cond ((or infl minfl
)
2041 (t (push (car exp
) infinityl
))))))
2042 (cond ((not (or infl minfl indl infinityl
))
2043 (return (cond ((atom sum
) sum
)
2044 ((or (not (free sum
'$zeroa
))
2045 (not (free sum
'$zerob
)))
2048 (t (cond ((null infinityl
)
2049 (cond (infl (cond ((null minfl
) (return '$inf
))
2051 (minfl (return '$minf
))
2052 ((> (length indl
) 1)
2053 ;; At this point we have a sum of '$ind. We factor
2054 ;; the sum and try again. This way we get the limit
2055 ;; of expressions like (a-b)*ind, where (a-b)--> 0.
2056 (cond ((not (alike1 (setq y
($factorsum exp
)) exp
))
2057 (return (limit y var val
'think
)))
2060 (t (return '$ind
))))
2061 (t (setq infl
(append infl infinityl
))))))
2063 oon
(setq y
(m+l
(append minfl infl
)))
2064 (cond ((alike1 exp
(setq y
(sratsimp (hyperex y
))))
2065 (cond ((not (infinityp val
))
2066 (setq infl
(cnv infl val
)) ;THIS IS HORRIBLE!!!!
2067 (setq minfl
(cnv minfl val
))))
2069 (cond ((every #'(lambda (j) (radicalp j var
))
2070 (append infl minfl
))
2071 (setq y
(rheur infl minfl
)))
2072 (t (setq y
(sheur infl minfl
))))))
2073 (t (setq y
(limit y var val
'think
))))
2074 (cond ((or (eq y
())
2075 (eq y t
)) (return ()))
2076 ((infinityp y
) (return y
))
2077 (t (return (m+ sum y
))))))
2079 ;; Limit n/d, using heuristics on the order of growth.
2082 (cond ((and (free n var
)
2085 (t (setq n
(cpa n d nil
))
2087 (cond ((oscip orig-n
) '$und
)
2089 ((equal n -
1) '$zeroa
)
2090 ((equal n
0) (m// orig-n d
)))))))
2093 ;;;L1 is a list of INF's and L2 is a list of MINF's. Added together
2094 ;;;it is indeterminate.
2095 (defun sheur (l1 l2
)
2096 (let ((term (sheur1 l1 l2
)))
2097 (cond ((equal term
'$inf
) '$inf
)
2098 ((equal term
'$minf
) '$minf
)
2099 (t (let ((new-num (m+l
(mapcar #'(lambda (num-term)
2100 (m// num-term
(car l1
)))
2102 (cond ((limit2 new-num
(m// 1 (car l1
)) var val
))))))))
2105 (cond ((atom exp
) nil
)
2106 (t (setq exp
(nformat exp
))
2107 (cond ((and (eq (caar exp
) 'mquotient
)
2108 (among var
(caddr exp
)))
2111 (defun zerop2 (z) (=0 (ridofab z
)))
2113 (defun raise (a) (m+ a
'$zeroa
))
2115 (defun lower (a) (m+ a
'$zerob
))
2117 (defun sincoshk (exp1 l sc
)
2118 (cond ((equal l
1) (lower l
))
2119 ((equal l -
1) (raise l
))
2121 ((member val
'($zeroa $zerob
) :test
#'eq
) (spangside exp1 l
))
2124 (defun spangside (e l
)
2125 (setq e
(behavior e var val
))
2126 (cond ((equal e
1) (raise l
))
2127 ((equal e -
1) (lower l
))
2130 ;; get rid of zeroa and zerob
2132 (if (among '$zeroa e
) (setq e
(maxima-substitute 0 '$zeroa e
)))
2133 (if (among '$zerob e
) (setq e
(maxima-substitute 0 '$zerob e
)))
2137 ;; returns true if exp is a polynomial raised to a numeric power
2138 (defun simplerd (exp)
2140 (mnump (caddr exp
)) ;; exponent must be a number - no variables
2141 (polyp (cadr exp
))))
2143 (defun branch1 (exp val
)
2144 (cond ((polyp exp
) nil
)
2145 ((simplerd exp
) (zerop2 (subin val
(cadr exp
))))
2147 (loop for v on
(cdr exp
)
2148 when
(branch1 (car v
) val
)
2151 (defun branch (exp val
)
2152 (cond ((polyp exp
) nil
)
2153 ((or (simplerd exp
) (mtimesp exp
))
2156 (every #'(lambda (j) (branch j val
)) (the list
(cdr exp
))))))
2158 (defun ser0 (e n d val
)
2159 (cond ((and (branch n val
) (branch d val
))
2163 ;;;NN* gets set by POFX, called by SER1, to get a list of exponents.
2164 (setq nn
* (ratmin nn
*))
2165 (setq n
(sratsimp (m^ n
(m^ var nn
*))))
2166 (setq d
(sratsimp (m^ d
(m^ var nn
*))))
2167 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2170 (t (try-lhospital-quit n d nil
))))
2172 (defun rheur (l1 l2
)
2174 (setq m1
(mapcar (function asymredu
) l1
))
2175 (setq m2
(mapcar (function asymredu
) l2
))
2176 (setq ans
(m+l
(append m1 m2
)))
2177 (cond ((rptrouble (m+l
(append l1 l2
)))
2178 (return (limit (simplify (rdsget (m+l
(append l1 l2
))))
2182 ((mplusp ans
) (return (sheur m1 m2
)))
2183 (t (return (limit ans var val t
))))))
2185 (defun rptrouble (rp)
2186 (not (equal (rddeg rp nil
) (rddeg (asymredu rp
) nil
))))
2188 (defun radicalp (exp var
)
2189 (cond ((polyinx exp var
()))
2190 ((mexptp exp
) (cond ((equal (caddr exp
) -
1.
)
2191 (radicalp (cadr exp
) var
))
2193 ((member (caar exp
) '(mplus mtimes
) :test
#'eq
)
2194 (every #'(lambda (j) (radicalp j var
))
2197 (defun involve (e nn
*)
2198 (declare (special var
))
2199 (cond ((atom e
) nil
)
2201 ((and (member (caar e
) nn
* :test
#'eq
) (among var
(cdr e
))) (cadr e
))
2202 (t (some #'(lambda (j) (involve j nn
*)) (cdr e
)))))
2204 (defun notinvolve (exp nn
*)
2205 (cond ((atom exp
) t
)
2207 ((member (caar exp
) nn
* :test
#'eq
) (not (among var
(cdr exp
))))
2208 ((every #'(lambda (j) (notinvolve j nn
*))
2211 (defun sheur1 (l1 l2
)
2213 (setq l1
(mapcar #'stirling0 l1
))
2214 (setq l2
(mapcar #'stirling0 l2
))
2215 (setq l1
(m+l
(maxi l1
)))
2216 (setq l2
(m+l
(maxi l2
)))
2217 (setq ans
(cpa l1 l2 t
))
2218 (return (cond ((=0 ans
) (m+ l1 l2
))
2219 ((equal ans
1.
) '$inf
)
2222 (defun zero-lim (cpa-list)
2223 (do ((l cpa-list
(cdr l
)))
2225 (and (eq (caar l
) 'gen
)
2226 (zerop2 (limit (cadar l
) var val
'think
))
2229 ;; Compare order of growth for R1 and R2. The result is 0, -1, +1
2230 ;; depending on the relative order of growth. 0 is returned if R1 and
2231 ;; R2 have the same growth; -1 if R1 grows much more slowly than R2;
2232 ;; +1 if R1 grows much more quickly than R2.
2233 (defun cpa (r1 r2 flag
)
2236 (cond ((alike1 t1 t2
) 0.
)
2238 (cond ((free t2 var
) 0.
)
2239 (t (let ((lim-ans (limit1 t2 var val
)))
2240 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2243 (let ((lim-ans (limit1 t1 var val
)))
2244 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2247 ;; Make T1 and T2 be a list of terms that are multiplied
2249 (cond ((mtimesp t1
) (setq t1
(cdr t1
)))
2250 (t (setq t1
(list t1
))))
2251 (cond ((mtimesp t2
) (setq t2
(cdr t2
)))
2252 (t (setq t2
(list t2
))))
2253 ;; Find the strengths of each term of T1 and T2
2254 (setq t1
(mapcar (function istrength
) t1
))
2255 (setq t2
(mapcar (function istrength
) t2
))
2256 ;; Compute the max of the strengths of the terms.
2257 (let ((ans (ismax t1
))
2259 (cond ((or (null ans
) (null d
))
2260 ;;(eq (car ans) 'gen) (eq (car d) 'gen))
2261 ;; ismax couldn't find highest term; give up
2264 (if (eq (car ans
) 'var
) (setq ans
(add-up-deg t1
)))
2265 (if (eq (car d
) 'var
) (setq d
(add-up-deg t2
)))
2266 ;; Can't just just compare dominating terms if there are
2267 ;; indeterm-inates present; e.g. X-X^2*LOG(1+1/X). So
2269 (cond ((or (zero-lim t1
)
2271 (cpa-indeterm ans d t1 t2 flag
))
2272 ((isgreaterp ans d
) 1.
)
2273 ((isgreaterp d ans
) -
1.
)
2276 (defun cpa-indeterm (ans d t1 t2 flag
)
2277 (cond ((not (eq (car ans
) 'var
))
2278 (setq ans
(gather ans t1
) d
(gather d t2
))))
2279 (let ((*indicator
(and (eq (car ans
) 'exp
)
2282 (setq test
(cpa1 ans d
))
2283 (cond ((and (zerop1 test
)
2284 (or (equal ($radcan
(m// (cadr ans
) (cadr d
))) 1.
)
2285 (and (polyp (cadr ans
))
2287 (equal (limit (m// (cadr ans
) (cadr d
)) var val
'think
)
2289 (let ((new-term1 (m// t1
(cadr ans
)))
2290 (new-term2 (m// t2
(cadr d
))))
2291 (cpa new-term1 new-term2 flag
)))
2294 (defun add-up-deg (strengthl)
2295 (do ((stl strengthl
(cdr stl
))
2298 ((null stl
) (list 'var
(m*l poxl
) (m+l degl
)))
2299 (cond ((eq (caar stl
) 'var
)
2300 (push (cadar stl
) poxl
)
2301 (push (caddar stl
) degl
)))))
2305 (cond ((eq (car p1
) 'gen
) (return 0.
)))
2306 (setq flag
(car p1
))
2311 (setq s1
(istrength p1
))
2312 (setq s2
(istrength p2
))
2315 ((isgreaterp s1 s2
) 1.
)
2316 ((isgreaterp s2 s1
) -
1.
)
2318 (setq *indicator nil
)
2320 ((and (poly? p1 var
) (poly? p2 var
))
2321 (setq p1
(m- p1 p2
))
2322 (cond ((zerop1 p1
) 0.
)
2323 (t (getsignl (hot-coef p1
)))))
2327 (list (m*t -
1 p2
))))
2328 (cond ((zerop2 s1
) 0.
)
2329 ((ratgreaterp s1
0.
) 1.
)
2333 (setq p1
(caddr p1
))
2334 (setq p2
(caddr p2
))
2335 (cond ((and (poly? p1 var
) (poly? p2 var
))
2336 (setq p1
(m- p1 p2
))
2337 (return (cond ((or (zerop1 p1
)
2338 (not (among var p1
)))
2340 (t (getsignl (hot-coef p1
))))))
2341 ((and (radicalp p1 var
) (radicalp p2 var
))
2344 (list (m*t -
1 p2
))))
2345 (return (cond ((eq s1
'$inf
) 1.
)
2346 ((eq s1
'$minf
) -
1.
)
2348 (cond ((ratgreaterp s1
0.
) 1.
)
2349 ((ratgreaterp 0. s1
) -
1.
)
2352 (t (return (cpa p1 p2 t
)))))
2354 (setq p1
(try-lhospital (asymredu p1
) (asymredu p2
) nil
))
2355 (return (cond ((zerop2 p1
) -
1.
)
2356 ((real-infinityp p1
) 1.
)
2359 ;;;EXPRESSIONS TO ISGREATERP ARE OF THE FOLLOWING FORMS
2360 ;;; ("VAR" POLY DEG)
2362 ;;; ("LOG" LOG(EXP))
2363 ;;; ("FACT" <A FACTORIAL EXPRESSION>)
2364 ;;; ("GEN" <ANY OTHER TYPE OF EXPRESSION>)
2366 (defun isgreaterp (a b
)
2369 (cond ((or (eq ta
'gen
)
2371 ((and (eq ta tb
) (eq ta
'var
))
2372 (ratgreaterp (caddr a
) (caddr b
)))
2373 ((and (eq ta tb
) (eq ta
'exp
))
2374 ;; Both are exponential order of infinity. Check the
2375 ;; exponents to determine which exponent is bigger.
2376 (eq (limit (m- `((%log
) ,(second a
)) `((%log
) ,(second b
)))
2379 ((member ta
(cdr (member tb
'(num log var exp fact gen
) :test
#'eq
)) :test
#'eq
)))))
2382 ;; Preprocess the list of products. Separate the terms that
2383 ;; exponentials and those that don't. Actually multiply the
2384 ;; exponential terms together to form a single term. Pass this and
2385 ;; the rest to ismax-core to find the max.
2386 (let (exp-terms non-exp-terms
)
2388 (if (eq 'exp
(car term
))
2389 (push term exp-terms
)
2390 (push term non-exp-terms
)))
2391 ;; Multiply the exp-terms together
2394 ;;(format t "exp-terms = ~A~%" exp-terms)
2395 (dolist (term exp-terms
)
2396 (setf product
(simplify (mul product
(second term
)))))
2397 ;;(format t "product = ~A~%" product)
2398 (setf product
`(exp ,($logcontract product
)))
2399 ;;(format t "product = ~A~%" product)
2400 (ismax-core (cons product non-exp-terms
)))
2403 (defun ismax-core (l)
2406 ((= (length l
) 1) (car l
)) ;If there is only 1 thing give it back.
2407 ((every #'(lambda (x)
2408 (not (eq (car x
) 'gen
))) l
)
2410 (do ((l1 (cdr l
) (cdr l1
))
2414 (cond ((isgreaterp temp-ans
(car l1
))
2415 (setq ans temp-ans
))
2416 ((isgreaterp (car l1
) temp-ans
)
2417 (setq temp-ans
(car l1
))
2418 (setq ans temp-ans
))
2419 (t (setq ans
())))))
2422 ;RETURNS LIST OF HIGH TERMS
2424 (cond ((atom all
) nil
)
2425 (t (do ((l (cdr all
) (cdr l
))
2427 (total 1) ; running total constant factor of hi-term
2428 (hi-terms (ncons (car all
)))
2430 ((null l
) (if (zerop2 total
) ; if high-order terms cancel each other
2431 all
; keep everything
2432 hi-terms
)) ; otherwise return list of high terms
2433 (setq compare
(limit (m// (car l
) hi-term
) var val
'think
))
2435 ((or (infinityp compare
)
2436 (and (eq compare
'$und
)
2437 (zerop2 (limit (m// hi-term
(car l
)) var val
'think
))))
2438 (setq total
1) ; have found new high term
2439 (setq hi-terms
(ncons (setq hi-term
(car l
)))))
2440 ((zerop2 compare
) nil
)
2441 ;; COMPARE IS IND, FINITE-VALUED, or und in both directions
2442 (t ; add to list of high terms
2443 (setq total
(m+ total compare
))
2444 (setq hi-terms
(append hi-terms
(ncons (car l
))))))))))
2448 (cond ((atom l
) (return nil
)))
2449 l1
(setq ans
(car l
))
2451 (cond ((null l
) (return ans
))
2452 ((ratgreaterp ans
(car l
)) (go l2
))
2457 (cond ((atom l
) (return nil
)))
2458 l1
(setq ans
(car l
))
2460 (cond ((null l
) (return ans
))
2461 ((ratgreaterp (car l
) ans
) (go l2
))
2469 ((or (mnump e
) (not (among var e
))) nil
)
2470 ((and (mexptp e
) (eq (cadr e
) var
))
2471 (push (caddr e
) nn
*))
2473 (t (mapc #'pofx
(cdr e
)))))
2476 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2477 (t (setq e
(subin (m+ var val
) e
))))
2479 (cond ((pofx e
) e
)))
2481 (defun gather (ind l
)
2483 (setq ind
(car ind
))
2484 loop
(cond ((null l
)
2485 (return (list ind
(m*l ans
))))
2486 ((equal (caar l
) ind
)
2487 (push (cadar l
) ans
)))
2491 ; returns rough class-of-growth of term
2492 (defun istrength (term)
2493 (cond ((mnump term
) (list 'num term
))
2494 ((atom term
) (cond ((eq term var
)
2496 (t (list 'num term
))))
2497 ((not (among var term
)) (list 'num term
))
2499 (let ((temp (ismax-core (mapcar #'istrength
(cdr term
)))))
2500 (cond ((not (null temp
)) temp
)
2503 (let ((temp (mapcar #'istrength
(cdr term
)))
2505 (setq temp1
(ismax temp
))
2506 (cond ((null temp1
) `(gen ,term
))
2507 ((eq (car temp1
) 'log
) `(log ,temp
))
2508 ((eq (car temp1
) 'var
) (add-up-deg temp
))
2511 (real-infinityp (limit term var val t
)))
2512 (let ((logterm (logred term
)))
2513 (cond ((and (among var
(caddr term
))
2514 (member (car (istrength logterm
))
2515 '(var exp fact
) :test
#'eq
)
2516 (real-infinityp (limit logterm var val t
)))
2517 (list 'exp
(m^
'$%e logterm
)))
2518 ((not (among var
(caddr term
)))
2519 (let ((temp (istrength (cadr term
))))
2520 (cond ((not (alike1 temp term
))
2521 (rplaca (cdr temp
) term
)
2522 (and (eq (car temp
) 'var
)
2524 (m* (caddr temp
) (caddr term
))))
2528 ((and (eq (caar term
) '%log
)
2529 (real-infinityp (limit term var val t
)))
2530 (let ((stren (istrength (cadr term
))))
2531 (cond ((member (car stren
) '(log var
) :test
#'eq
)
2533 ((and (eq (car stren
) 'exp
)
2534 (eq (caar (second stren
)) 'mexpt
))
2535 (istrength (logred (second stren
))))
2537 ((eq (caar term
) 'mfactorial
)
2539 ((let ((temp (hyperex term
)))
2540 (and (not (alike1 term temp
))
2542 (t (list 'gen term
))))
2544 ;; log reduce - returns log of s1
2546 (or (and (eq (cadr s1
) '$%e
) (caddr s1
))
2547 (m* (caddr s1
) `((%log
) ,(cadr s1
)))))
2549 (defun asymredu (rd)
2550 (cond ((atom rd
) rd
)
2552 ((not (among var rd
)) rd
)
2553 ((polyinx rd var t
))
2555 (cond ((eq (cadr rd
) var
) rd
)
2557 (factor ($expand
(m^
(polyinx (cadr rd
) var t
)
2561 (t (simplify (cons (list (caar rd
))
2562 (mapcar #'asymredu
(cdr rd
)))))))
2565 (let ((dn** ()) (nn** ()))
2566 (cond ((atom rd
) rd
)
2568 ((not (among var rd
)) rd
)
2572 (cond ((eq (cadr rd
) var
) rd
)
2573 (t (setq dn
** (caddr rd
))
2574 (setq nn
** (factor (cadr rd
)))
2575 (cond ((mtimesp nn
**)
2576 (m*l
(mapcar #'(lambda (j) (m^ j dn
**))
2579 (t (simplify (cons (ncons (caar rd
))
2580 (mapcar #'rdfact
(cdr rd
))))))))
2582 (defun cnv (expl val
)
2583 (mapcar #'(lambda (e)
2584 (maxima-substitute (cond ((eq val
'$zerob
)
2585 (m* -
1 (m^ var -
1)))
2590 (t (m^
(m+ var
(m* -
1 val
)) -
1.
)))
2595 (defun pwtaylor (exp var l terms
)
2596 (prog (coef ans c mc
)
2597 (cond ((=0 terms
) (return nil
)) ((=0 l
) (setq mc t
)))
2600 loop
(setq c
(1+ c
))
2601 (cond ((or (> c
10.
) (equal c terms
))
2603 (t (setq exp
(sdiff exp var
))))
2604 tag1
(setq coef
($radcan
(subin l exp
)))
2605 (cond ((=0 coef
) (setq terms
(1+ terms
)) (go loop
)))
2612 (m^
`((mfactorial) ,c
) -
1)
2613 (m^
(if mc var
(m+t
(m*t -
1 l
) var
)) c
)))))
2618 ((simplerd e
) (rdtay e
))
2619 (t (cons (list (caar e
))
2620 (mapcar #'rdsget
(cdr e
))))))
2623 (cond (limit-using-taylor ($ratdisrep
($taylor rd var val
1.
)))
2627 (prog (varlist p c e d $ratfac
)
2628 (setq varlist
(ncons var
))
2629 (setq p
(ratnumerator (cdr (ratrep* (cadr rd
)))))
2630 (cond ((< (length p
) 3.
) (return rd
)))
2633 (setq c
(m^ var
(m* d e
)))
2634 (setq d
($ratsimp
(varinvert (m* (pdis p
) (m^ var
(m- d
)))
2636 (setq d
(pwtaylor (m^ d e
) var
0.
3.
))
2637 (return (m* c
(varinvert d var
)))))
2639 (defun varinvert (e var
) (subin (m^t var -
1.
) e
))
2642 (prog ((varlist (list var
)))
2643 (return (let (($ratfac nil
))
2645 (pdegr (cadr (ratrep* p
)))))))
2647 (defun rat-no-ratfac (e)
2648 (let (($ratfac nil
))
2653 (defun rddeg (rd low
*)
2654 (cond ((or (mnump rd
)
2655 (not (among var rd
)))
2660 (m* (deg (cadr rd
)) (caddr rd
)))
2662 (addn (mapcar #'(lambda (j)
2666 (setq rd
(andmapcar #'(lambda (j) (rddeg j low
*))
2668 (cond (low* (ratmin rd
))
2672 (cond ((or (atom pf
) (not (eq (caadr (ratf var
)) (car pf
))))
2674 (low* (cadr (reverse pf
)))
2676 ;;There is some confusion here. We need to be aware of Branch cuts etc....
2677 ;;when doing this section of code. It is not very carefully done so there
2678 ;;are bugs still lurking. Another misfortune is that LIMIT or its inferiors
2679 ;;sometimes decides to change the limit VAL in midstream. This must be corrected
2680 ;;since LIMIT's interaction with the data base environment must be maintained.
2681 ;;I'm not sure that this code can ever be called with VAL other than $INF but
2682 ;;there is a hook in the first important cond clause to cathc them anyway.
2685 (let ((num-power (rddeg n nil
))
2686 (den-power (rddeg d nil
))
2687 (coef ()) (coef-sign ()) (power ()))
2688 (setq coef
(m// ($ratcoef
($expand n
) var num-power
)
2689 ($ratcoef
($expand d
) var den-power
)))
2690 (setq coef-sign
(getsignl coef
))
2691 (setq power
(m// num-power den-power
))
2692 (cond ((eq (ask-integer power
'$integer
) '$integer
)
2693 (cond ((eq (ask-integer power
'$even
) '$even
) '$even
)
2694 (t '$odd
)))) ;Can be extended from here.
2695 (cond ((or (eq val
'$minf
)
2698 (equal val
0)) ()) ;Can be extended to cover some these.
2699 ((ratgreaterp den-power num-power
)
2700 (cond ((equal coef-sign
1.
) '$zeroa
)
2701 ((equal coef-sign -
1) '$zerob
)
2702 ((equal coef-sign
0) 0)
2704 ((ratgreaterp num-power den-power
)
2705 (cond ((equal coef-sign
1.
) '$inf
)
2706 ((equal coef-sign -
1) '$minf
)
2707 ((equal coef-sign
0) nil
) ; should never be zero
2708 ((null coef-sign
) '$infinity
)))
2711 (defun radlim (e n d
)
2713 (cond ((eq val
'$infinity
) (throw 'limit nil
))
2715 (setq nl
(m* var -
1))
2716 (setq n
(subin nl n
))
2717 (setq d
(subin nl d
))
2718 (setq val
'$inf
))) ;This is the Culprit. Doesn't tell the DATABASE.
2719 (cond ((eq val
'$inf
)
2720 (setq nl
(asymredu n
))
2721 (setq dl
(asymredu d
))
2723 ((or (rptrouble n
) (rptrouble d
))
2724 (return (limit (m* (rdsget n
) (m^
(rdsget d
) -
1.
)) var val t
)))
2725 (t (return (asy nl dl
))))))
2726 (setq nl
(limit n var val t
))
2727 (setq dl
(limit d var val t
))
2728 (cond ((and (zerop2 nl
) (zerop2 dl
))
2729 (cond ((or (polyp n
) (polyp d
))
2730 (return (try-lhospital-quit n d t
)))
2731 (t (return (ser0 e n d val
)))))
2732 (t (return ($radcan
(ratrad (m// n d
) n d nl dl
)))))))
2734 (defun ratrad (e n d nl dl
)
2736 (cond ((equal nl
0) (return 0))
2739 (cond ((equal dl
0) (setq d1
'$infinity
)) ;No direction Info.
2742 ((equal (setq d
(behavior d var val
)) 1)
2744 ((equal d -
1) (setq d1
'$minf
))
2745 (t (throw 'limit nil
))))
2748 (cond ((equal (setq n
(behavior n var val
)) 1)
2750 ((equal n -
1) (setq n1
'$zerob
))
2752 (t (return ($radcan
(ridofab (subin val e
))))))
2753 (return (simplimtimes (list n1 d1
)))))
2755 ;;; Limit(log(XXX), var, 0, val), where val is either zerob (limit from below)
2756 ;;; or zeroa (limit from above).
2757 (defun simplimln (expr var val
)
2758 (let ((arglim (limit (cadr expr
) var val
'think
)) (dir))
2759 (cond ((eq arglim
'$inf
) '$inf
) ;log(inf) = inf
2760 ;;log(minf,infinity,zerob) = infinity & log(0) = infinity
2761 ((or (member arglim
'($minf $infinity $zerob
)))
2763 ((eq arglim
'$zeroa
) '$minf
) ;log(zeroa) = minf
2764 ;; log(ind)=und & log(und)=und
2765 ((member arglim
'($ind $und
)) '$und
)
2766 ;; log(1^(-)) = zerob, log(1^(+)) = zeroa & log(1)=0
2768 (if (or (eq val
'$zerob
) (eq var
'$zeroa
)) val
0))
2769 ;; Special case of arglim = 0
2771 (setq dir
(behavior (cadr expr
) var val
))
2772 (cond ((eql dir -
1) '$infinity
)
2773 ((eql dir
0) '$infinity
)
2774 ((eql dir
1) '$minf
)))
2775 ;; When arglim is off the negative real axis, use direct substitution
2776 ((off-negative-real-axisp arglim
)
2777 (ftake '%log arglim
))
2779 ;; We know that arglim is a negative real number, say xx.
2780 ;; When the imaginary part of (cadr expr) near var is negative,
2781 ;; return log(-x) - %i*pi; when the imaginary part of (cadr expr)
2782 ;; near var is positive return log(-x) + %i*pi; and when
2783 ;; we cannot determine the behavior of the imaginary part,
2784 ;; return a nounform. The value of val (either zeroa or zerob)
2785 ;; determines what is meant by "near" (smaller than var when
2786 ;; val is zerob and larger than var when val is zeroa).
2787 (setq dir
(behavior ($imagpart
(cadr expr
)) var val
))
2788 (cond ((or (eql dir
1) (eql dir -
1))
2789 (add (ftake '%log
(mul -
1 arglim
)) (mul dir
'$%i
'$%pi
)))
2790 (t (throw 'limit nil
))))))) ;do a nounform return
2791 (setf (get '%log
'simplim%function
) #'simplimln
)
2792 (setf (get '%plog
'simplim%function
) #'simplimln
)
2794 (defun simplim%limit
(e x pt
)
2795 (declare (ignore e x pt
))
2797 (setf (get '%limit
'simplim%function
) #'simplim%limit
)
2799 (defun simplim%unit_step
(e var val
)
2800 (let ((lim (limit (cadr e
) var val
'think
)))
2801 (cond ((eq lim
'$und
) '$und
)
2802 ((eq lim
'$ind
) '$ind
)
2803 ((eq lim
'$zerob
) 0)
2804 ((eq lim
'$zeroa
) 1)
2805 ((not (lenient-realp lim
)) (throw 'limit nil
)) ;catches infinity too
2806 ;; catches minf & inf cases
2807 ((eq t
(mgrp 0 lim
)) 0)
2808 ((eq t
(mgrp lim
0)) 1)
2810 (setf (get '$unit_step
'simplim%function
) #'simplim%unit_step
)
2812 (defun simplim%conjugate
(e var val
)
2813 (let ((lim (limit (cadr e
) var val
'think
)))
2814 (cond ((off-negative-real-axisp lim
)
2815 (ftake '$conjugate lim
))
2816 (t (throw 'limit nil
)))))
2817 (setf (get '$conjugate
'simplim%function
) #'simplim%conjugate
)
2819 (defun simplim%imagpart
(e var val
)
2820 (let ((lim (limit (cadr e
) var val
'think
)))
2821 (cond ((eq lim
'$und
) '$und
)
2823 ((eq lim
'$infinity
) (throw 'limit nil
))
2824 (t (mfuncall '$imagpart lim
)))))
2825 (setf (get '$imagpart
'simplim%function
) #'simplim%imagpart
)
2826 (setf (get '%imagpart
'simplim%function
) #'simplim%imagpart
)
2828 (defun simplim%realpart
(e var val
)
2829 (let ((lim (limit (cadr e
) var val
'think
)))
2830 (cond ((eq lim
'$und
) '$und
)
2831 ((eq lim
'$ind
) '$ind
)
2832 ((eq lim
'$infinity
) (throw 'limit nil
))
2833 (t (mfuncall '$realpart lim
)))))
2834 (setf (get '$realpart
'simplim%function
) #'simplim%realpart
)
2835 (setf (get '%realpart
'simplim%function
) #'simplim%realpart
)
2836 ;;; Limit of the Factorial function
2838 (defun simplimfact (expr var val
)
2839 (let* ((arglim (limit (cadr expr
) var val
'think
)) ; Limit of the argument.
2841 (cond ((eq arglim
'$inf
) '$inf
)
2842 ((member arglim
'($minf $infinity $und $ind
) :test
#'eq
) '$und
)
2844 ((and (or (maxima-integerp arglim
)
2845 (setq arg2
(integer-representation-p arglim
)))
2846 (eq ($sign arg2
) '$neg
))
2847 ;; A negative integer or float or bigfloat representation.
2848 (let ((dir (limit (add (cadr expr
) (mul -
1 arg2
)) var val
'think
))
2849 (even (mevenp arg2
)))
2850 (cond ((or (and even
2860 (t (throw 'limit nil
)))))
2862 ;; Call simplifier to get value at the limit of the argument.
2863 (simplify (list '(mfactorial) arglim
))))))
2864 (setf (get 'mfactorial
'simplim%function
) #'simplimfact
)
2866 (defun simplim%erf-%tanh
(fn arg
)
2867 (let ((arglim (limit arg var val
'think
))
2870 (cond ((eq arglim
'$inf
) 1)
2871 ((eq arglim
'$minf
) -
1)
2872 ((eq arglim
'$infinity
)
2873 (destructuring-bind (rpart . ipart
)
2875 (setq rlim
(limit rpart var origval
'think
))
2876 (cond ((eq fn
'%tanh
)
2877 (cond ((equal rlim
'$inf
) 1)
2878 ((equal rlim
'$minf
) -
1)))
2880 (setq ans
(limit (m* rpart
(m^t ipart -
1)) var origval
'think
))
2881 (setq ans
($asksign
(m+ `((mabs) ,ans
) -
1)))
2882 (cond ((or (eq ans
'$pos
) (eq ans
'$zero
))
2883 (cond ((eq rlim
'$inf
) 1)
2884 ((eq rlim
'$minf
) -
1)
2887 ((eq arglim
'$und
) '$und
)
2888 ((member arglim
'($zeroa $zerob $ind
) :test
#'eq
) arglim
)
2889 ;;;Ignore tanh(%pi/2*%I) and multiples of the argument.
2891 ;; erf (or tanh) of a known value is just erf(arglim).
2892 (simplify (list (ncons fn
) arglim
))))))
2894 (defun in-domain-of-atan (z)
2895 (setq z
(trisplit z
)) ; split z into real and imaginary parts
2896 (let ((x (car z
)) (y (cdr z
))) ;z = x+%i*y
2898 (eq t
(mnqp x
0)) ; x # 0
2900 (eq t
(mgrp -
1 y
)) ; y < -1
2901 (eq t
(mgrp y
1)))))) ; y > 1
2903 (defun simplim%atan
(e x pt
)
2904 (let ((lim (limit (cadr e
) x pt
'think
)))
2905 (cond ((or (eq lim
'$zeroa
) (eq lim
'$zerob
) (eq lim
0) (eq lim
'$ind
)) lim
)
2906 ((or (eq lim
'$und
) (eq lim
'$infinity
)) (throw 'limit nil
))
2907 ((eq lim
'$inf
) #$%pi
/2$
) ;atan(inf) -> %pi/2
2908 ((eq lim
'$minf
) #$-%pi
/2$
) ;atan(-inf) -> -%pi/2
2909 ((in-domain-of-atan (ridofab lim
)) ; direct substitution
2910 (ftake '%atan
(ridofab lim
)))
2911 (t (limit ($logarc e
) x pt
'think
)))))
2912 (setf (get '%atan
'simplim%function
) #'simplim%atan
)
2914 (defmvar extended-reals
2915 (append infinitesimals infinities
(list '$und
'$ind
)))
2917 ;; Most instances of atan2 are simplified to atan expressions, but this routine
2918 ;; handles tricky cases such as limit(atan2((x^2-2), x^3-2*x), x, sqrt(2), minus).
2919 ;; Taylor and Gruntz cannot handle the discontinuity at atan(0, -1)
2921 ;; When possible, we want to evaluate the limit of an atan2 expression using
2922 ;; direct substituion--that produces, I think, the least surprising values.
2924 ;; The general simplifier catches atan2(0,0) and it transforms atan2(minf or inf,X)
2925 ;; and atan2(X, minf or inf) into an atan expression, but it doesn't catch
2926 ;; atan2(X, zerob or zeroa) or atan2(zerob or zeroa, X). For the other extended
2927 ;; real (ind,und, or infinity) arguments, the general simplifer gives sign errors.
2929 (defun simplim%atan2
(e v pt
)
2930 (let ((y (second e
)) (x (third e
)) (xlim) (ylim) (xlim-z) (ylim-z) (q))
2931 (setq xlim
(limit x v pt
'think
))
2932 (setq ylim
(limit y v pt
'think
))
2933 (setq xlim-z
(ridofab xlim
)
2934 ylim-z
(ridofab ylim
))
2935 ;; For cases for which direct substitution fails, normalize
2936 ;; x & y and try again.
2937 (setq q
(cond ((eq xlim
'$inf
) x
)
2938 ((eq xlim
'$minf
) (mul -
1 x
))
2940 ((eq ylim
'$minf
) (mul -
1 y
))
2941 ((and (eq xlim
'$zerob
) (zerop2 ylim
)) (mul -
1 x
))
2942 ((and (eq xlim
'$zeroa
) (zerop2 ylim
)) x
)
2943 ((and (eq ylim
'$zerob
) (zerop2 xlim
)) (mul -
1 y
))
2944 ((and (eq ylim
'$zeroa
) (zerop2 xlim
)) y
)
2947 (when (not (eql q
1))
2950 (setq xlim
(limit x v pt
'think
))
2951 (setq ylim
(limit y v pt
'think
))
2952 (setq xlim-z
(ridofab xlim
)
2953 ylim-z
(ridofab ylim
)))
2956 ((and (eq '$zerob ylim
) (eq t
(mgrp 0 xlim
))) (mul -
1 '$%pi
))
2957 ((and (eq '$zerob ylim
) (eq t
(mgrp xlim
0))) 0)
2958 ((and (eq '$zeroa ylim
) (eq t
(mgrp 0 xlim
))) '$%pi
)
2959 ((and (eq '$zeroa ylim
) (eq t
(mgrp xlim
0))) 0)
2960 ((and (eql xlim
1) (eql ylim
'$inf
)) (div '$%pi
2))
2961 ((and (eql xlim -
1) (eql ylim
0)) '$ind
)
2963 ;; Use direct substitution when ylim-z # 0 or xlim-z > 0. We need
2964 ;; to check that xlim-z & ylim-z are real too.
2965 ((and (lenient-realp xlim-z
) (lenient-realp ylim-z
)
2966 (or (eq t
(mnqp ylim-z
0)) (eq t
(mgrp xlim-z
0))))
2967 (ftake '$atan2 ylim-z xlim-z
))
2969 (throw 'limit nil
)))))
2970 (setf (get '$atan2
'simplim%function
) #'simplim%atan2
)
2972 (defun simplimsch (sch arg
)
2973 (cond ((real-infinityp arg
)
2974 (cond ((eq sch
'%sinh
) arg
) (t '$inf
)))
2975 ((eq arg
'$infinity
) '$infinity
)
2976 ((eq arg
'$ind
) '$ind
)
2977 ((eq arg
'$und
) '$und
)
2978 ((and (eq sch
'%sinh
)
2979 (or (eq arg
'$zerob
) (eq arg
'$zeroa
)))
2981 (t (let (($exponentialize t
))
2982 (resimplify (list (ncons sch
) (ridofab arg
)))))))
2984 ;; simple limit of sin and cos
2985 (defun simplimsc (exp fn arg
)
2986 (cond ((member arg
'($inf $minf $ind
) :test
#'eq
) '$ind
)
2987 ((member arg
'($und $infinity
) :test
#'eq
)
2989 ((member arg
'($zeroa $zerob
) :test
#'eq
)
2990 (cond ((eq fn
'%sin
) arg
)
2991 (t (m+ 1 '$zerob
))))
2993 (simplify (list (ncons fn
) (ridofab arg
)))
2996 (defun simplim%tan
(arg)
2997 (let ((arg1 (ridofab (limit arg var val
'think
))))
2999 ((member arg1
'($inf $minf $infinity $ind $und
) :test
#'eq
) '$und
)
3001 (let ((c (trigred (pip arg1
))))
3002 (cond ((not (equal ($imagpart arg1
) 0)) '$infinity
)
3003 ((and (eq (caar c
) 'rat
)
3006 (setq arg1
(behavior arg var val
))
3007 (cond ((= arg1
1) '$inf
)
3008 ((= arg1 -
1) '$minf
)
3010 ((and (eq (caar c
) 'rat
)
3013 (setq arg1
(behavior arg var val
))
3014 (cond ((= arg1
1) '$minf
)
3017 (t (throw 'limit
())))))
3019 (setq arg1
(behavior arg var val
))
3020 (cond ((equal arg1
1) '$zeroa
)
3021 ((equal arg1 -
1) '$zerob
)
3023 (t (simp-%tan
(list '(%tan
) arg1
) 1 nil
)))))
3025 (defun simplim%asinh
(arg)
3026 (cond ((member arg
'($inf $minf $zeroa $zerob $ind $und
) :test
#'eq
)
3028 ((eq arg
'$infinity
) '$und
)
3029 (t (simplify (list '(%asinh
) (ridofab arg
))))))
3031 (defun simplim%acosh
(arg)
3032 (cond ((equal (ridofab arg
) 1) '$zeroa
)
3033 ((eq arg
'$inf
) arg
)
3034 ((eq arg
'$minf
) '$infinity
)
3035 ((member arg
'($und $ind $infinity
) :test
#'eq
) '$und
)
3036 (t (simplify (list '(%acosh
) (ridofab arg
))))))
3038 (defun simplim%atanh
(arg dir
)
3039 ;; Compute limit(atanh(x),x,arg). If ARG is +/-1, we need to take
3040 ;; into account which direction we're approaching ARG.
3041 (cond ((zerop2 arg
) arg
)
3042 ((member arg
'($ind $und $infinity $minf $inf
) :test
#'eq
)
3044 ((equal (setq arg
(ridofab arg
)) 1.
)
3045 ;; The limit at 1 should be complex infinity because atanh(x)
3046 ;; is complex for x > 1, but inf if we're approaching 1 from
3048 (if (eq dir
'$zerob
)
3052 ;; Same as above, except for the limit is at -1.
3053 (if (eq dir
'$zeroa
)
3056 (t (simplify (list '(%atanh
) arg
)))))
3058 (defun simplim%asin
(e x pt
)
3059 (let ((lim (limit (cadr e
) x pt
'think
)) (dir) (lim-sgn))
3060 (cond ((member lim
'($zeroa $zerob
)) lim
) ;asin(zeoroa/b) = zeroa/b
3061 ((in-domain-of-asin lim
) ;direct subsititution
3063 ((member lim
'($und $ind $inf $minf $infinity
)) ;boundary cases
3066 (setq e
(trisplit (cadr e
))) ;overwrite e!
3067 (setq dir
(behavior (cdr e
) x pt
))
3068 (setq lim-sgn
($csign lim
))
3071 (throw 'limit t
)) ;unable to find behavior of imaginary part
3072 ;; for the values of asin on the branch cuts, see DLMF 4.23.20 & 4.23.21
3073 ((or (eq '$pos lim-sgn
) (eq '$neg lim-sgn
))
3074 ;; continuous from above
3075 (if (eql dir
1) (ftake '%asin lim
) (sub '$%pi
(ftake '%asin lim
))))
3076 (t (throw 'limit t
))))))) ; unable to find sign of real part of lim.
3077 (setf (get '%asin
'simplim%function
) #'simplim%asin
)
3079 (defun simplim%acos
(e x pt
)
3080 (let ((lim (limit (cadr e
) x pt
'think
)) (dir) (lim-sgn))
3081 (cond ((in-domain-of-asin lim
) ;direct subsititution
3083 ((member lim
'($und $ind $inf $minf $infinity
)) ;boundary cases
3086 (setq e
(trisplit (cadr e
))) ;overwrite e!
3087 (setq dir
(behavior (cdr e
) x pt
))
3088 (setq lim-sgn
($csign lim
))
3091 (throw 'limit t
)) ;unable to find behavior of imaginary part
3092 ;; for the values of acos on the branch cuts, see DLMF 4.23.24 & 4.23.25
3093 ;; http://dlmf.nist.gov/4.23.E24
3094 ((or (eq '$pos lim-sgn
) (eq '$neg lim-sgn
))
3095 ;; continuous from above
3096 (if (eql dir
1) (ftake '%acos lim
) (sub (mul 2 '$%pi
) (ftake '%acos lim
))))
3097 (t (throw 'limit t
))))))) ; unable to find sign of real part of lim.
3098 (setf (get '%acos
'simplim%function
) #'simplim%acos
)
3100 ;; Limit of an %integrate expression. For a definite integral
3101 ;; integrate(ee,var,a,b), when ee is free of the limit variable
3102 (defun simplim%integrate
(e x pt
)
3103 (let* ((ee (second e
)) ;ee = integrand
3104 (var (third e
)) ;integration variable
3105 (a (fourth e
)) ;lower limit or nil if indefinite
3106 (b (fifth e
)) ;lower limit or nil if indefinite
3108 (cond ((and a b
($freeof x ee
) ($freeof x var
))
3109 (setq alim
(limit a x pt
'think
))
3110 (setq blim
(limit b x pt
'think
))
3111 (if (and (lenient-extended-realp alim
)
3112 (lenient-extended-realp blim
)
3113 (not (eq alim
'$infinity
))
3114 (not (eq blim
'$infinity
)))
3115 (ftake '%integrate ee var alim blim
)
3118 (throw 'limit t
)))))
3119 (setf (get '%integrate
'simplim%function
) #'simplim%integrate
)
3121 (defun subftake (op subarg arg
)
3122 (simplifya (subfunmake op subarg arg
) t
))
3124 (defun off-one-to-inf (z)
3125 (setq z
(trisplit z
)) ; split z into x+%i*y
3127 (eq t
(mnqp (cdr z
) 0)) ; y # 0
3128 (eq t
(mgrp 1 (car z
))))) ; x < 1
3130 (defun simplim%li
(expr x pt
)
3131 (let ((n (car (subfunsubs expr
))) (e (car (subfunargs expr
))))
3133 (setq e
(limit e x pt
'think
))
3134 (cond ((and (eq e
'$minf
) (eql n
2)) '$minf
)
3135 ((and (eq e
'$inf
) (eql n
2)) '$infinity
)
3136 ((or (eql e
1) (off-one-to-inf e
))
3137 (subftake '$li
(list n
) (list e
)))
3138 (t (throw 'limit nil
))))
3139 ;; Claim ignorance when order depends on limit variable.
3140 (t (throw 'limit nil
)))))
3141 (setf (get '$li
'simplim%function
) #'simplim%li
)
3142 (setf (get '%li
'simplim%function
) #'simplim%li
)
3144 (defun simplim$psi
(order arg val
)
3145 (if (and (not (equal (length order
) 1))
3146 (not (equal (length arg
) 1)))
3148 (setq order
(car order
)
3150 (cond ((equal order
0)
3151 (destructuring-bind (rpart . ipart
)
3153 (cond ((not (equal ipart
0)) (throw 'limit
()))
3154 (t (setq rpart
(limit rpart var val
'think
))
3155 (cond ((eq rpart
'$zeroa
) '$minf
)
3156 ((eq rpart
'$zerob
) '$inf
)
3157 ((eq rpart
'$inf
) '$inf
)
3158 ((eq rpart
'$minf
) '$und
)
3159 ((equal (getsignl rpart
) -
1) (throw 'limit
()))
3160 (t (simplify (subfunmake '$psi
(list order
)
3161 (list rpart
)))))))))
3162 ((and (integerp order
) (> order
0)
3163 (equal (limit arg var val
'think
) '$inf
))
3164 (cond ((mevenp order
) '$zerob
)
3165 ((moddp order
) '$zeroa
)
3166 (t (throw 'limit
()))))
3167 (t (throw 'limit
()))))
3169 (defun simplim%inverse_jacobi_ns
(arg m
)
3170 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3172 `((%inverse_jacobi_ns
) ,arg
,m
)))
3174 (defun simplim%inverse_jacobi_nc
(arg m
)
3175 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3176 `((%elliptic_kc
) ,m
)
3177 `((%inverse_jacobi_nc
) ,arg
,m
)))
3179 (defun simplim%inverse_jacobi_sc
(arg m
)
3180 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3181 `((%elliptic_kc
) ,m
)
3182 `((%inverse_jacobi_sc
) ,arg
,m
)))
3184 (defun simplim%inverse_jacobi_dc
(arg m
)
3185 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3186 `((%elliptic_kc
) ,m
)
3187 `((%inverse_jacobi_dc
) ,arg
,m
)))
3189 (defun simplim%inverse_jacobi_cs
(arg m
)
3190 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3192 `((%inverse_jacobi_cs
) ,arg
,m
)))
3194 (defun simplim%inverse_jacobi_ds
(arg m
)
3195 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3197 `((%inverse_jacobi_ds
) ,arg
,m
)))
3199 (defun simplim%signum
(e x pt
)
3200 (let ((e (limit (cadr e
) x pt
'think
)) (sgn))
3201 (cond ((eq '$minf e
) -
1)
3203 ((eq '$infinity e
) '$und
)
3204 ((eq '$ind e
) '$ind
)
3209 (setq sgn
(mnqp e
0))
3210 (cond ((eq t sgn
) (ftake '%signum e
))
3211 (t (throw 'limit nil
))))))) ; don't know
3212 (setf (get '%signum
'simplim%function
) #'simplim%signum
)
3214 ;; more functions for limit to handle
3216 (defun lfibtophi (e)
3217 (cond ((not (involve e
'($fib
))) e
)
3218 ((eq (caar e
) '$fib
)
3219 (let ((lnorecurse t
))
3220 ($fibtophi
(list '($fib
) (lfibtophi (cadr e
))) lnorecurse
)))
3222 (mapcar #'lfibtophi
(cdr e
))))))
3224 ;;; FOLLOWING CODE MAKES $LDEFINT WORK
3226 (defmfun $ldefint
(exp var ll ul
&aux $logabs ans a1 a2
)
3227 (setq $logabs t ans
(sinint exp var
)
3228 a1
(toplevel-$limit ans var ul
'$minus
)
3229 a2
(toplevel-$limit ans var ll
'$plus
))
3230 (and (member a1
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3231 (setq a1
(nounlimit ans var ul
)))
3232 (and (member a2
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3233 (setq a2
(nounlimit ans var ll
)))
3234 ($expand
(m- a1 a2
)))
3236 (defun nounlimit (exp var val
)
3237 (setq exp
(restorelim exp
))
3238 (nconc (list '(%limit
) exp var
(ridofab val
))
3239 (cond ((eq val
'$zeroa
) '($plus
))
3240 ((eq val
'$zerob
) '($minus
)))))
3242 ;; substitute inside noun form of %derivative
3243 ;; for cases such as limit('diff(x+2,x), x, 1)
3244 ;; -> limit('diff(xx+3), xx, 0)
3246 ;; maxima-substitute with *atp* skips over %derivative
3248 ;; substitutes diff(f(realvar), realvar, n)
3249 ;; -> diff(f(var+val), var, n)
3250 (defun derivative-subst (exp val var realvar
)
3251 (cond ((atom exp
) exp
)
3252 ((eq '%derivative
(caar exp
))
3255 (cons ;; the function being differentiated
3256 (maxima-substitute (m+ val var
) realvar
(cadr exp
))
3257 (cons ;; the var of differentiation
3258 (maxima-substitute var realvar
(caddr exp
))
3259 (cdddr exp
))))) ;; the order of the derivative
3261 (mapcar (lambda (x) (derivative-subst x val var realvar
))
3264 ;;;Used by Defint also.
3266 (or (involve e
'(%sin %cos %tan
))
3267 (among '$%i
(%einvolve e
))))
3269 (defun %einvolve
(e)
3270 (when (among '$%e e
) (%einvolve01 e
)))
3272 (defun %einvolve01
(e)
3273 (cond ((atom e
) nil
)
3277 (among var
(caddr e
)))
3279 (t (some #'%einvolve
(cdr e
)))))
3281 (declare-top (unspecial *indicator nn
* dn
* exp var val origval taylored
3282 $tlimswitch logcombed lhp? lhcount $ratfac
))
3288 ;; "On Computing Limits in a Symbolic Manipulation System"
3289 ;; PhD Dissertation ETH Zurich 1996
3291 ;; The algorithm identifies the most rapidly varying (MRV) subexpression,
3292 ;; replaces it with a new variable w, rewrites the expression in terms
3293 ;; of the new variable, and then repeats.
3295 ;; The algorithm doesn't handle oscillating functions, so it can't do things like
3296 ;; limit(sin(x)/x, x, inf).
3298 ;; To handle limits involving functions like gamma(x) and erf(x), the
3299 ;; gruntz algorithm requires them to be written in terms of asymptotic
3300 ;; expansions, which maxima cannot currently do.
3302 ;; The algorithm assumes that everything is real, so it can't
3303 ;; currently handle limit((-2)^x, x, inf).
3305 ;; This is one of the methods used by maxima's $limit.
3306 ;; It is also directly available to the user as $gruntz.
3309 ;; most rapidly varying subexpression of expression exp with respect to limit variable var.
3310 ;; returns a list of subexpressions which are in the same MRV equivalence class.
3311 (defun mrv (exp var
)
3312 (cond ((freeof var exp
)
3317 (mrv-max (mrv (cadr exp
) var
)
3318 (mrv (m*l
(cddr exp
)) var
)
3321 (mrv-max (mrv (cadr exp
) var
)
3322 (mrv (m+l
(cddr exp
)) var
)
3325 (cond ((freeof var
(caddr exp
))
3326 (mrv (cadr exp
) var
))
3327 ((member (limitinf (logred exp
) var
) '($inf $minf
) :test
#'eq
)
3328 (mrv-max (list exp
) (mrv (caddr exp
) var
) var
))
3329 (t (mrv-max (mrv (cadr exp
) var
) (mrv (caddr exp
) var
) var
))))
3331 (mrv (cadr exp
) var
))
3332 ((equal (length (cdr exp
)) 1)
3333 (mrv (cadr exp
) var
))
3334 ((equal (length (cdr exp
)) 2)
3335 (mrv-max (mrv (cadr exp
) var
)
3336 (mrv (caddr exp
) var
)
3338 (t (tay-error "mrv not implemented" exp
))))
3340 ;; takes two lists of expresions, f and g, and limit variable var.
3341 ;; members in each list are assumed to be in same MRV equivalence
3342 ;; class. returns MRV set of the union of the inputs - either f or g
3343 ;; or the union of f and g.
3344 (defun mrv-max (f g var
)
3351 (return (union f g
))))
3352 (let ((c (mrv-compare (car f
) (car g
) var
)))
3358 (return (union f g
)))
3359 (t (merror "MRV-MAX: expected '>' '<' or '='; found: ~M" c
))))))
3361 (defun mrv-compare (a b var
)
3362 (let ((c (limitinf (m// `((%log
) ,a
) `((%log
) ,b
)) var
)))
3365 ((member c
'($inf $minf
) :test
#'eq
)
3369 ;; rewrite expression exp by replacing members of MRV set omega with
3370 ;; expressions in terms of new variable wsym. return cons pair of new
3371 ;; version of exp and the log of the new variable wsym.
3372 (defun mrv-rewrite (exp omega var wsym
)
3373 (setq omega
(stable-sort omega
(lambda (x y
) (> (length (mrv x var
))
3374 (length (mrv y var
))))));FIXME consider a total order function with #'sort
3375 (let* ((g (car (last omega
)))
3377 (sig (equal (mrv-sign logg var
) 1))
3378 (w (if sig
(m// 1 wsym
) wsym
))
3379 (logw (if sig
(m* -
1 logg
) logg
)))
3380 (mapcar (lambda (x y
)
3381 ;;(mtell "y:~M x:~M exp:~M~%" y x exp)
3382 (setq exp
(syntactic-substitute y x exp
)))
3384 (mapcar (lambda (f) ;; rewrite each element of omega
3385 (let* ((logf (logred f
))
3386 (c (mrv-leadterm (m// logf logg
) var nil
)))
3387 (cond ((not (equal (cadr c
) 0))
3388 (merror "MRV-REWRITE: expected leading term to be constant in ~M" c
)))
3389 ;;(mtell "logg: ~M logf: ~M~%" logg logf)
3392 (m* (car c
) logg
))))))
3396 ;;; if log w(x) = h(x), rewrite all subexpressions of the form
3397 ;;; log(f(x)) as log(w^-c f(x)) + c h(x) with c the unique constant
3398 ;;; such that w^-c f(x) is strictly less rapidly varying than w.
3399 (defun mrv-rewrite-logs (exp wsym logw
)
3400 (cond ((atom exp
) exp
)
3402 (not (freeof wsym exp
)))
3403 (let* ((f (cadr exp
))
3404 (c ($lopow
(calculate-series f wsym
)
3407 (m* (m^ wsym
(m- c
))
3408 (mrv-rewrite-logs f wsym logw
)))
3413 (mrv-rewrite-logs e wsym logw
))
3416 ;; returns list of two elements: coeff and exponent of leading term of exp,
3417 ;; after rewriting exp in term of its MRV set omega.
3418 (defun mrv-leadterm (exp var omega
)
3419 (prog ((new-omega ()))
3420 (cond ((freeof var exp
)
3421 (return (list exp
0))))
3422 (dolist (term omega
)
3423 (cond ((subexp exp term
)
3424 (push term new-omega
))))
3425 (setq omega new-omega
)
3427 (setq omega
(mrv exp var
))))
3428 (cond ((member var omega
:test
#'eq
)
3429 (let* ((omega-up (mrv-moveup omega var
))
3430 (e-up (car (mrv-moveup (list exp
) var
)))
3431 (mrv-leadterm-up (mrv-leadterm e-up var omega-up
)))
3432 (return (mrv-movedown mrv-leadterm-up var
)))))
3433 (destructuring-let* ((wsym (gensym "w"))
3436 ((f . logw
) (mrv-rewrite exp omega var wsym
))
3437 (series (calculate-series (mrv-rewrite-logs f wsym logw
)
3439 (setq series
(maxima-substitute logw
`((%log
) ,wsym
) series
))
3440 (setq lo
($lopow series wsym
))
3441 (when (or (not ($constantp lo
))
3442 (not (free series
'%derivative
)))
3443 ;; (mtell "series: ~M lo: ~M~%" series lo)
3444 (tay-error "error in series expansion" f
))
3445 (setq coef
($coeff series wsym lo
))
3446 (when (not (free coef wsym
))
3447 (tay-error "MRV-LEADTERM: failed to extract leading coefficient; obtained" coef
))
3448 ;;(mtell "exp: ~M f: ~M~%" exp f)
3449 ;;(mtell "series: ~M~%coeff: ~M~%pow: ~M~%" series coef lo)
3450 (return (list coef lo
)))))
3452 (defun mrv-moveup (l var
)
3453 (mapcar (lambda (exp)
3454 (simplify-log-of-exp
3455 (syntactic-substitute `((mexpt) $%e
,var
) var exp
)))
3458 (defun mrv-movedown (l var
)
3459 (mapcar (lambda (exp) (syntactic-substitute `((%log simp
) ,var
) var exp
))
3462 ;; test whether sub is a subexpression of exp
3463 (defun subexp (exp sub
)
3464 (let ((dummy (gensym)))
3465 (putprop dummy t
'internal
)
3466 (not (alike1 (maxima-substitute dummy
3471 ;; Generate $lhospitallim terms of taylor expansion.
3472 ;; Ideally we would use a lazy series representation that generates
3473 ;; more terms as higher order terms cancel.
3474 (defun calculate-series (exp var
)
3475 (let ((cntx ($supcontext
)))
3479 (mfuncall '$assume
(ftake 'mgreaterp var
0))
3480 (putprop var t
'internal
); keep var from appearing in questions to user
3481 ($taylor exp var
0 $lhospitallim
))
3482 (remprop var
'internal
)
3483 ($killcontext cntx
))))
3485 (defun mrv-sign (exp var
)
3486 (cond ((freeof var exp
)
3487 (let ((sign ($sign
($radcan exp
))))
3488 (cond ((eq sign
'$zero
)
3494 (t (tay-error " cannot determine mrv-sign" exp
)))))
3498 (* (mrv-sign (cadr exp
) var
)
3499 (mrv-sign (m*l
(cddr exp
)) var
)))
3501 (equal (mrv-sign (cadr exp
) var
) 1))
3504 (cond ((equal (mrv-sign (cadr exp
) var
) -
1)
3505 (tay-error " complex expression in gruntz limit" exp
)))
3506 (mrv-sign (m+ -
1 (cadr exp
)) var
))
3508 (mrv-sign (limitinf exp var
) var
))
3509 (t (tay-error " cannot determine mrv-sign" exp
))))
3511 ;; gruntz algorithm for limit of exp as var goes to positive infinity
3512 (defun limitinf (exp var
)
3513 (prog (($exptsubst nil
))
3514 (cond ((freeof var exp
)
3516 (destructuring-let* ((c0-e0 (mrv-leadterm exp var nil
))
3519 (sig (mrv-sign e0 var
)))
3520 (cond ((equal sig
1)
3523 (cond ((equal (mrv-sign c0 var
) 1)
3525 ((equal (mrv-sign c0 var
) -
1)
3529 ;; example: gruntz(n^n/(n^n+(n-1)^n), n, inf);
3530 (tay-error " infinite recursion in limitinf" exp
))
3531 (return (limitinf c0 var
)))))))
3533 ;; user-level function equivalent to $limit.
3534 ;; direction must be specified if limit point is not infinite
3535 ;; The arguments are checked and a failure of taylor is catched.
3537 (defmfun $gruntz
(expr var val
&rest rest
)
3539 (when (> (length rest
) 1)
3541 (intl:gettext
"gruntz: too many arguments; expected just 3 or 4")))
3542 (setq dir
(car rest
))
3543 (when (and (not (member val
'($inf $minf $zeroa $zerob
)))
3544 (not (member dir
'($plus $minus
))))
3546 (intl:gettext
"gruntz: direction must be 'plus' or 'minus'")))
3548 (catch 'taylor-catch
3549 (let ((silent-taylor-flag t
))
3550 (declare (special silent-taylor-flag
))
3551 (gruntz1 expr var val dir
))))
3552 (if (or (null ans
) (eq ans t
))
3554 `(($gruntz simp
) ,expr
,var
, val
,dir
)
3555 `(($gruntz simp
) ,expr
,var
,val
))
3558 ;; This function is for internal use in $limit.
3559 (defun gruntz1 (exp var val
&rest rest
)
3560 (cond ((> (length rest
) 1)
3561 (merror (intl:gettext
"gruntz: too many arguments; expected just 3 or 4"))))
3562 (let (($logexpand t
) ; gruntz needs $logexpand T
3563 (newvar (gensym "w"))
3565 (putprop newvar t
'internal
); keep var from appearing in questions to user
3566 (cond ((eq val
'$inf
)
3569 (setq exp
(maxima-substitute (m* -
1 newvar
) var exp
)))
3571 (setq exp
(maxima-substitute (m// 1 newvar
) var exp
)))
3573 (setq exp
(maxima-substitute (m// -
1 newvar
) var exp
)))
3575 (setq exp
(maxima-substitute (m+ val
(m// 1 newvar
)) var exp
)))
3577 (setq exp
(maxima-substitute (m+ val
(m// -
1 newvar
)) var exp
)))
3578 (t (merror (intl:gettext
"gruntz: direction must be 'plus' or 'minus'; found: ~M") dir
)))
3579 (limitinf exp newvar
)))
3581 ;; substitute y for x in exp
3582 ;; similar to maxima-substitute but does not simplify result
3583 (defun syntactic-substitute (y x exp
)
3584 (cond ((alike1 x exp
) y
)
3587 (mapcar (lambda (exp)
3588 (syntactic-substitute y x exp
))
3591 ;; log(exp(subexpr)) -> subexpr
3592 ;; without simplifying entire exp
3593 (defun simplify-log-of-exp (exp)
3594 (cond ((atom exp
) exp
)
3597 (eq '$%e
(cadadr exp
)))
3600 (mapcar #'simplify-log-of-exp
