1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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 origval
31 *indicator numer denom exp var val
33 $exponentialize lhp? lhcount
34 loginprod? a context limit-assumptions
35 limit-top integer-info old-integer-info
))
37 (defconstant +behavior-count
+ 4)
38 (defvar *behavior-count-now
*)
39 (defvar *getsignl-asksign-ok
* nil
)
41 (load-macsyma-macros rzmac
)
43 (defmvar simplimplus-problems
()
44 "A list of all problems in the stack of recursive calls to simplimplus.")
46 (defmvar limit-answers
()
47 "An association list for storing limit answers.")
49 (defmvar limit-using-taylor
()
50 "Is the current limit computation using taylor expansion?")
52 (defmvar preserve-direction
() "Makes `limit' return Direction info.")
54 (unless (boundp 'integer-info
) (setq integer-info
()))
56 ;; For limits toward infinity for the gruntz code, we assume that the limit
57 ;; variable exceeds *large-positive-number*. This value matches the value
58 ;; that the limit code uses for the same purpose.
59 (defmvar *large-positive-number
* (expt 10 8))
61 ;; Don't ask sign questions about $ind.
62 (putprop '$ind t
'internal
)
64 ;; This should be made to give more information about the error.
66 ;; (cond (errorsw (throw 'errorsw t))
67 ;; (t (merror "Discontinuity Encountered"))))
69 ;;(DEFUN PUTLIMVAL (E V)
70 ;; (let ((exp (cons '(%limit) (list e var val))))
71 ;; (cond ((not (assolike exp limit-answers))
72 ;; (setq limit-answers (cons (cons exp v) limit-answers))
76 (defun putlimval (e v
&aux exp
)
77 (setq exp
`((%limit
) ,e
,var
,val
))
78 (unless (assolike exp limit-answers
)
79 (push (cons exp v
) limit-answers
))
83 (let ((exp (cons '(%limit
) (list e var val
))))
84 (assolike exp limit-answers
)))
86 (defmacro limit-catch
(exp var val
)
88 (let ((ans (catch 'errorsw
89 (catch 'limit
(limit ,exp
,var
,val
'think
)))))
90 (if (or (null ans
) (eq ans t
))
94 (defmfun $limit
(&rest args
)
95 (let ((first-try (apply #'toplevel-$limit args
)))
96 (if (and (consp first-try
) (eq (caar first-try
) '%limit
))
97 (let ((*getsignl-asksign-ok
* t
))
98 (apply #'toplevel-$limit args
))
100 ;; The function toplevel-$limit sets $numer, $%enumer, %emode, and
101 ;; $%e_to_numlog to false. When any of these option variables are true,
102 ;; we resimplify the limit in the current context.
104 (when (or $numer $%enumer $%emode $%e_to_numlog
)
105 (setq first-try
(resimplify first-try
)))
108 (defun toplevel-$limit
(&rest args
)
109 (let ((limit-assumptions ())
110 (old-integer-info ())
117 (declare (special limit-assumptions old-integer-info
120 (setq old-integer-info integer-info
)
121 (setq integer-info
()))
124 (let ((exp1 ()) (lhcount $lhospitallim
) (*behavior-count-now
* 0)
125 (exp ()) (var ()) (val ()) (dr ())
126 (*indicator
()) (taylored ()) (origval ())
127 (logcombed ()) (lhp?
())
128 (varlist ()) (ans ()) (genvar ()) (loginprod?
())
129 (limit-answers ()) (limitp t
) (simplimplus-problems ())
130 (lenargs (length args
))
132 (declare (special lhcount
*behavior-count-now
* exp var val
*indicator
133 taylored origval logcombed lhp?
134 varlist genvar loginprod? limitp
))
136 (unless (or (= lenargs
3) (= lenargs
4) (= lenargs
1))
138 ;; Is it a LIST of Things?
139 (when (setq ans
(apply #'limit-list args
))
141 (setq exp1
(specrepcheck (first args
)))
142 (when (and (atom exp1
)
143 (member exp1
'(nil t
)))
144 ;; The expression is 'T or 'NIL. Return immediately.
147 (setq var
(setq genfoo
(gensym)) ; Use a gensym. Not foo.
150 (setq var
(second args
))
151 (when ($constantp var
)
152 (merror (intl:gettext
"limit: second argument must be a variable, not a constant; found: ~M") var
))
153 (unless (or ($subvarp var
) (atom var
))
154 (merror (intl:gettext
"limit: variable must be a symbol or subscripted symbol; found: ~M") var
))
155 (setq val
(infsimp (third args
)))
156 ;; infsimp converts -inf to minf. it also converts -infinity to
157 ;; infinity, although perhaps this should generate the error below.
158 (when (and (not (atom val
))
159 (some #'(lambda (x) (not (freeof x val
)))
161 (merror (intl:gettext
"limit: third argument must be a finite value or one of: inf, minf, infinity; found: ~M") val
))
162 (when (eq val
'$zeroa
) (setq dr
'$plus
))
163 (when (eq val
'$zerob
) (setq dr
'$minus
))))
165 (unless (member (fourth args
) '($plus $minus
) :test
#'eq
)
166 (merror (intl:gettext
"limit: direction must be either 'plus' or 'minus'; found: ~M") (fourth args
)))
167 (setq dr
(fourth args
))))
168 (if (and (atom var
) (not (among var val
)))
170 (let ((realvar var
)) ;; Var is funny so make it a gensym.
172 (setq exp
(maxima-substitute var realvar exp1
))
173 (putprop var realvar
'limitsub
)))
174 (unless (or $limsubst
(eq var genfoo
))
175 (when (limunknown exp
)
176 (return `((%limit
) ,@(cons exp1
(cdr args
))))))
177 (setq varlist
(ncons var
) genvar nil origval val
)
178 ;; Transform limits to minf to limits to inf by
179 ;; replacing var with -var everywhere.
180 (when (eq val
'$minf
)
183 exp
(subin (m* -
1 var
) exp
)))
185 ;; Transform the limit value.
186 (unless (infinityp val
)
188 (let ((*atp
* t
) (realvar var
))
189 ;; *atp* prevents substitution from applying to vars
190 ;; bound by %sum, %product, %integrate, %limit
192 (putprop var t
'internal
)
193 (setq exp
(derivative-subst exp val var realvar
))
194 (setq exp
(maxima-substitute (m+ val var
) realvar exp
))))
195 (setq val
(cond ((eq dr
'$plus
) '$zeroa
)
196 ((eq dr
'$minus
) '$zerob
)
200 ;; Make assumptions about limit var being very small or very large.
201 ;; Assumptions are forgotten upon exit.
202 (unless (= lenargs
1)
203 (limit-context var val dr
))
205 ;; Resimplify in light of new assumptions.
206 (setq exp
(resimplify
211 (limitsimp ($expand exp
1 0) var
)))))))
213 (if (not (or (real-epsilonp val
) ;; if direction of limit not specified
215 (setq ans
(both-side exp var val
)) ;; compute from both sides
216 (let ((d (catch 'mabs
(mabs-subst exp var val
))))
217 (cond ;; otherwise try to remove absolute value
218 ((eq d
'$und
) (return '$und
))
221 (setq ans
(limit-catch exp var val
));; and find limit from one side
225 (setq ans
(catch 'taylor-catch
226 (let ((silent-taylor-flag t
))
227 (declare (special silent-taylor-flag
))
228 (gruntz1 exp var val
)))))
230 ;; try taylor series expansion if simple limit didn't work
231 (if (and (null ans
) ;; if no limit found and
232 $tlimswitch
;; user says ok to use taylor and
233 (not limit-using-taylor
));; not already doing taylor
234 (let ((limit-using-taylor t
))
235 (declare (special limit-using-taylor
))
236 (setq ans
(limit-catch exp var val
))))))
239 (return (clean-limit-exp ans
))
240 (return (cons '(%limit
) args
))))) ;; failure: return nounform
241 (restore-assumptions))))
243 (defun clean-limit-exp (exp)
244 (setq exp
(restorelim exp
))
245 (if preserve-direction exp
(ridofab exp
)))
247 ;; Users who want limit to map over equality (mequal) will need to do that
249 (defun limit-list (exp1 &rest rest
)
250 (if (and (mbagp exp1
) (not (mequalp exp1
)))
251 `(,(car exp1
) ,@(mapcar #'(lambda (x) (apply #'toplevel-$limit
`(,x
,@rest
))) (cdr exp1
)))
254 (defun limit-context (var val direction
) ;Only works on entry!
256 (assume '((mgreaterp) lim-epsilon
0))
257 (assume '((mgreaterp) prin-inf
100000000))
258 (setq limit-assumptions
(make-limit-assumptions var val direction
))
263 (defun make-limit-assumptions (var val direction
)
264 (let ((new-assumptions))
265 (cond ((or (null var
) (null val
))
267 ((and (not (infinityp val
)) (null direction
))
270 `(,(assume `((mgreaterp) ,var
100000000)) ,@new-assumptions
))
272 `(,(assume `((mgreaterp) -
100000000 ,var
)) ,@new-assumptions
))
273 ((eq direction
'$plus
)
274 `(,(assume `((mgreaterp) ,var
0)) ,@new-assumptions
)) ;All limits around 0
275 ((eq direction
'$minus
)
276 `(,(assume `((mgreaterp) 0 ,var
)) ,@new-assumptions
))
280 (defun restore-assumptions ()
281 ;;;Hackery until assume and forget take reliable args. Nov. 9 1979.
283 (do ((assumption-list limit-assumptions
(cdr assumption-list
)))
284 ((null assumption-list
) t
)
285 (forget (car assumption-list
)))
286 (forget '((mgreaterp) lim-epsilon
0))
287 (forget '((mgreaterp) prin-inf
100000000))
288 (cond ((and (not (null integer-info
))
290 (do ((list integer-info
(cdr list
)))
292 (i-$remove
`(,(cadar list
) ,(caddar list
))))
293 (setq integer-info old-integer-info
))))
295 ;; The optional arg allows the caller to decide on the value of
296 ;; preserve-direction. Default is nil, since we immediately ridofab.
297 (defun both-side (exp var val
&optional
(preserve nil
))
298 (let* ((preserve-direction preserve
)
299 (la (toplevel-$limit exp var val
'$plus
)) lb
)
300 ; Immediately propagate an und without trying the
302 (when (eq la
'$und
) (return-from both-side
'$und
))
303 (setf lb
(toplevel-$limit exp var val
'$minus
))
304 ; Immediately propagate an und
305 (when (eq lb
'$und
) (return-from both-side
'$und
))
306 (let ((ra (ridofab la
))
308 (cond ((eq t
(meqp ra rb
))
312 ; Maxima does not consider equal(ind,ind) to be true, but
313 ; if both one-sided limits are ind then we want to call
314 ; the two-sided limit ind (e.g., limit(sin(1/x),x,0)).
316 ((or (not (free la
'%limit
))
317 (not (free lb
'%limit
)))
320 (let ((infa (infinityp la
))
321 (infb (infinityp lb
)))
322 (cond ((and infa infb
)
323 ; inf + minf => infinity
330 (defun limunknown (f)
331 (catch 'limunknown
(limunknown1 (specrepcheck f
))))
333 (defun limunknown1 (f)
334 (cond ((mapatom f
) nil
)
335 ;; special pass for fib(X). What else?
336 ((and (consp f
) (eq '$fib
(caar f
))) nil
)
337 ((or (not (safe-get (caar f
) 'operators
))
338 (member (caar f
) '(%sum %product mncexpt
) :test
#'eq
)
339 ;;Special function code here i.e. for li[2](x).
340 (and (eq (caar f
) 'mqapply
)
341 (not (get (subfunname f
) 'specsimp
))))
342 (if (not (free f var
)) (throw 'limunknown t
)))
343 (t (mapc #'limunknown1
(cdr f
)) nil
)))
346 (if (involve e
'(%gamma
)) (setq e
($makefact e
)))
347 (cond ((involve e
'(mfactorial))
348 (setq e
(simplify ($minfactorial e
))))
352 ;; or nil if sign unknown or complex
354 (let ((z (ridofab z
)))
355 (if (not (free z var
)) (setq z
(toplevel-$limit z var val
)))
356 (let ((*complexsign
* t
))
357 (let ((sign (if *getsignl-asksign-ok
* ($asksign z
) ($sign z
))))
358 (cond ((eq sign
'$pos
) 1)
360 ((eq sign
'$zero
) 0))))))
362 (defun restorelim (exp)
363 (cond ((null exp
) nil
)
364 ((atom exp
) (or (and (symbolp exp
) (get exp
'limitsub
)) exp
))
365 ((and (consp (car exp
)) (eq (caar exp
) 'mrat
))
367 (cons (restorelim (cadr exp
))
368 (restorelim (cddr exp
)))))
369 (t (cons (car exp
) (mapcar #'restorelim
(cdr exp
))))))
372 (defun mabs-subst (exp var val
) ; RETURNS EXP WITH MABS REMOVED, OR THROWS.
373 (let ((d (involve exp
'(mabs)))
377 ((not (and (equal ($imagpart
(let ((v (limit-catch d var val
)))
378 ;; The above call might
379 ;; throw 'limit, so we
380 ;; need to catch it. If
382 ;; limit without ABS, we
383 ;; assume the limit is
384 ;; undefined. Is this
385 ;; right? Anyway, this
386 ;; fixes Bug 1548643.
391 (equal ($imagpart var
) 0)))
392 (cond ((eq arglim
'$infinity
)
393 ;; Check for $infinity as limit of argument.
396 (throw 'mabs
'retn
))))
397 (t (do ((ans d
(involve exp
'(mabs))) (a () ()))
399 (setq a
(mabs-subst ans var val
))
400 (setq d
(limit a var val t
))
404 (setq d
(behavior a var val
))
405 (if (zerop1 d
) (throw 'mabs
'retn
))))
408 (cond ((or (eq d
'$zeroa
) (eq d
'$inf
)
410 ;; fails on limit(abs(sin(x))/sin(x), x, inf)
411 (eq ($sign d
) '$pos
))
412 (setq exp
(maxima-substitute a
`((mabs) ,ans
) exp
)))
413 ((or (eq d
'$zerob
) (eq d
'$minf
)
414 (eq ($sign d
) '$neg
))
415 (setq exp
(maxima-substitute (m* -
1 a
) `((mabs) ,ans
) exp
)))
416 (t (throw 'mabs
'retn
))))
417 (t (throw 'mabs
'retn
))))))))))
419 ;; Called on an expression that might contain $INF, $MINF, $ZEROA, $ZEROB. Tries
420 ;; to simplify it to sort out things like inf^inf or inf+1.
422 (simpinf-ic exp
(count-general-inf exp
)))
424 (defun count-general-inf (expr)
425 (count-atoms-matching
426 (lambda (x) (or (infinityp x
) (real-epsilonp x
))) expr
))
428 (defun count-atoms-matching (predicate expr
)
429 "Count the number of atoms in the Maxima expression EXPR matching PREDICATE,
430 ignoring dummy variables and array indices."
432 ((atom expr
) (if (funcall predicate expr
) 1 0))
433 ;; Don't count atoms that occur as a limit of %integrate, %sum, %product,
435 ((member (caar expr
) dummy-variable-operators
)
436 (count-atoms-matching predicate
(cadr expr
)))
437 ;; Ignore array indices
438 ((member 'array
(car expr
)) 0)
440 for arg in
(cdr expr
)
441 summing
(count-atoms-matching predicate arg
)))))
443 (defun simpinf-ic (exp &optional infinity-count
)
445 ;; A very slow identity transformation...
448 ;; If there's only one infinity, we replace it by a variable and take the
449 ;; limit as that variable goes to infinity. Use $gensym in case we can't
450 ;; compute the answer and the limit leaks out.
451 (1 (let* ((val (or (inf-typep exp
) (epsilon-typep exp
)))
453 (expr (subst var val exp
))
454 (limit (toplevel-$limit expr var val
)))
456 ;; Now we look to see whether the computed limit is any simpler than
457 ;; what we shoved in (which we'll define as "doesn't contain EXPR as a
458 ;; subtree"). If so, return it.
459 ((not (subtree-p expr limit
:test
#'equal
))
462 ;; Otherwise, return the original form: apparently, we can't compute
463 ;; the limit we needed, and it's uglier than what we started with.
466 ;; If more than one infinity, we have to be a bit more careful.
468 (let* ((arguments (mapcar 'simpinf
(cdr exp
)))
469 (new-expression (cons (list (caar exp
)) arguments
))
472 ;; If any of the arguments are undefined, we are too.
473 ((among '$und arguments
) '$und
)
474 ;; If we ended up with something indeterminate, we punt and just return
476 ((amongl '(%limit $ind
) arguments
) exp
)
478 ;; Exponentiation & multiplication
479 ((mexptp exp
) (simpinf-expt (first arguments
) (second arguments
)))
480 ((mtimesp exp
) (simpinf-times arguments
))
482 ;; Down to at most one infinity? We do this after exponentiation to
483 ;; avoid zeroa^zeroa => 0^0, which will raise an error rather than just
484 ;; returning und. We do it after multiplication to avoid zeroa * inf =>
486 ((<= (setf infinities-left
(count-general-inf new-expression
)) 1)
487 (simpinf-ic new-expression infinities-left
))
490 ((mplusp exp
) (simpinf-plus arguments
))
493 (t new-expression
))))))
495 (defun simpinf-times (arguments)
496 (declare (special exp var val
))
497 ;; When we have a product, we need to spot that zeroa * zerob = zerob, zeroa *
498 ;; inf = und etc. Note that (SIMPINF '$ZEROA) => 0, so a nonzero atom is not
499 ;; an infinitesimal. Moreover, we can assume that each of ARGUMENTS is either
500 ;; a number, computed successfully by the recursive SIMPINF call, or maybe a
501 ;; %LIMIT noun-form (in which case, we aren't going to be able to tell the
504 ((member 0 arguments
)
506 ((find-if #'infinityp arguments
) '$und
)
507 ((every #'atom arguments
) 0)
510 ((member '$infinity arguments
)
511 (if (every #'atom arguments
)
515 (t (simplimit (cons '(mtimes) arguments
) var val
))))
517 (defun simpinf-expt (base exponent
)
518 ;; In the comments below, zero* represents one of 0, zeroa, zerob.
520 ;; TODO: In some cases we give up too early. E.g. inf^(2 + 1/inf) => inf^2
521 ;; (which should simplify to inf)
531 ((0 $zeroa $zerob
) '$und
)
532 (t (list '(mexpt) base exponent
))))
533 ;; minf^inf = infinity <== Or should it be und?
536 ;; minf^foo = minf^foo
541 ((0 $zeroa $zerob
) '$und
)
542 (t (list '(mexpt) base exponent
))))
546 ;; zero*^foo = zero*^foo
551 ((0 $zeroa $zerob
) '$und
)
552 (t (list '(mexpt) base exponent
))))
553 ;; a^b where a is pretty much anything except for a naked
554 ;; inf,minf,zeroa,zerob or 0.
557 ;; When a isn't crazy, try a^b = e^(b log(a))
558 ((not (amongl (append *infinitesimals
* *infinities
*) base
))
559 (simpinf (m^
'$%e
(m* exponent
`((%log
) ,base
)))))
561 ;; No idea. Just return what we've found so far.
562 (t (list '(mexpt) base exponent
))))))
564 (defun simpinf-plus (arguments)
565 ;; We know that none of the arguments are infinitesimals, since SIMPINF never
566 ;; returns one of them. As such, we partition our arguments into infinities
567 ;; and everything else. The latter won't have any "hidden" infinities like
568 ;; lim(x,x,inf), since SIMPINF gave up on anything containing a %lim already.
569 (let ((bigs) (others))
570 (dolist (arg arguments
)
571 (cond ((infinityp arg
) (push arg bigs
))
572 (t (push arg others
))))
574 ;; inf + minf or the like
575 ((cdr (setf bigs
(delete-duplicates bigs
))) '$und
)
576 ;; inf + smaller + stuff
578 ;; I don't think this can happen, since SIMPINF goes back to the start if
579 ;; there are fewer than two infinities in the arguments, but let's be
581 (t (cons '(mplus) others
)))))
583 ;; Simplify expression with zeroa or zerob.
584 (defun simpab (small)
585 (cond ((null small
) ())
586 ((member small
'($zeroa $zerob $inf $minf $infinity
) :test
#'eq
) small
)
587 ((not (free small
'$ind
)) '$ind
) ;Not exactly right but not
588 ((not (free small
'$und
)) '$und
) ;causing trouble now.
589 ((mapatom small
) small
)
590 (t (let ((preserve-direction t
)
591 (new-small (subst (m^
'$inf -
1) '$zeroa
592 (subst (m^
'$minf -
1) '$zerob small
))))
593 (simpinf new-small
)))))
596 ;;;*I* INDICATES: T => USE LIMIT1,THINK, NIL => USE SIMPLIMIT.
597 (defun limit (exp var val
*i
*)
599 ((among '$und exp
) '$und
)
602 ((not (among var exp
))
603 (cond ((amongl '($inf $minf $infinity $ind
) exp
)
605 ((amongl '($zeroa $zerob
) exp
)
606 ;; Simplify expression with zeroa or zerob.
610 (t (putlimval exp
(cond ((and limit-using-taylor
613 (taylim exp var val
*i
*))
614 ((ratp exp var
) (ratlim exp
))
615 ((or (eq *i
* t
) (radicalp exp var
))
616 (limit1 exp var val
))
618 (cond ((or (mtimesp exp
) (mexptp exp
))
619 (limit1 exp var val
))
620 (t (simplimit exp var val
))))
621 (t (simplimit exp var val
)))))))
623 (defun limitsimp (exp var
)
624 (limitsimp-expt (sin-sq-cos-sq-sub exp
) var
))
625 ;;Hack for sin(x)^2+cos(x)^2.
627 ;; if var appears in base and power of expt,
628 ;; push var into power of of expt
629 (defun limitsimp-expt (exp var
)
630 (cond ((or (atom exp
)
632 (freeof var exp
)) exp
)
634 (not (freeof var
(cadr exp
)))
635 (not (freeof var
(caddr exp
))))
636 (m^
'$%e
(simplify `((%log
) ,exp
))))
637 (t (subst0 (cons (cons (caar exp
) ())
638 (mapcar #'(lambda (x)
639 (limitsimp-expt x var
))
643 (defun gather-args-of (e fn x
)
644 (cond (($mapatom e
) nil
)
645 ((and (eq fn
(caar e
)) (not (freeof x e
))) (cdr e
))
647 (reduce #'append
(mapcar #'(lambda (q) (gather-args-of q fn x
)) (cdr e
))))))
649 ;; When X depends on x, replace cos(X)^2 + sin(X)^2 by 1
650 (defun sin-sq-cos-sq-sub (e &optional
(x var
))
651 (let ((ccc nil
) (z) (ee))
652 (cond (($mapatom e
) e
)
654 (setq ccc
(gather-args-of e
'%cos x
))
657 (setq ee
(maxima-substitute z
(power (ftake '%cos g
) 2) e
))
658 (setq ee
(maxima-substitute (sub 1 z
) (power (ftake '%sin g
) 2) ee
))
659 (when (freeof z
(sratsimp ee
))
662 ;; maybe this isn't needed, but I think it's not wrong.
663 ((eq 'mqapply
(caar e
))
664 (subftake (caar (second e
))
665 (mapcar #'(lambda (q) (sin-sq-cos-sq-sub q x
)) (subfunsubs e
))
666 (mapcar #'(lambda (q) (sin-sq-cos-sq-sub q x
)) (subfunargs e
))))
667 ;; map sin-sq-cos-sq-sub over the args
669 (fapply (caar e
) (mapcar #'(lambda (q) (sin-sq-cos-sq-sub q x
)) (cdr e
)))))))
671 (defun expand-trigs (x var
)
674 ((and (or (eq (caar x
) '%sin
)
676 (not (free (cadr x
) var
)))
678 ((member 'array
(car x
))
679 ;; Some kind of array reference. Return it.
681 (t (simplify (cons (ncons (caar x
))
682 (mapcar #'(lambda (x)
683 (expand-trigs x var
))
688 (cond ((not (involve e
689 '(%cot %csc %binomial
690 %sec %coth %sech %csch
691 %acot %acsc %asec %acoth
693 %jacobi_ns %jacobi_nc %jacobi_cs
694 %jacobi_ds %jacobi_dc
)))
696 (t ($ratsimp
(tansc1 e
)))))
698 (defun tansc1 (e &aux tem
)
700 ((and (setq e
(cons (car e
) (mapcar 'tansc1
(cdr e
)))) ()))
701 ((setq tem
(assoc (caar e
) '((%cot . %tan
) (%coth . %tanh
)
702 (%sec . %cos
) (%sech . %cosh
)
703 (%csc . %sin
) (%csch . %sinh
)) :test
#'eq
))
704 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
)) -
1.
)))
705 ((setq tem
(assoc (caar e
) '((%jacobi_nc . %jacobi_cn
)
706 (%jacobi_ns . %jacobi_sn
)
707 (%jacobi_cs . %jacobi_sc
)
708 (%jacobi_ds . %jacobi_sd
)
709 (%jacobi_dc . %jacobi_cd
)) :test
#'eq
))
710 ;; Converts Jacobi elliptic function to its reciprocal
712 (tansc1 (m^
(list (ncons (cdr tem
)) (cadr e
) (third e
)) -
1.
)))
713 ((setq tem
(member (caar e
) '(%sinh %cosh %tanh
) :test
#'eq
))
714 (let (($exponentialize t
))
716 ((setq tem
(assoc (caar e
) '((%acsc . %asin
) (%asec . %acos
)
717 (%acot . %atan
) (%acsch . %asinh
)
718 (%asech . %acosh
) (%acoth . %atanh
)) :test
#'eq
))
719 (list (ncons (cdr tem
)) (m^t
(cadr e
) -
1.
)))
720 ((and (eq (caar e
) '%binomial
) (among var
(cdr e
)))
721 (m// `((mfactorial) ,(cadr e
))
722 (m* `((mfactorial) ,(m+t
(cadr e
) (m- (caddr e
))))
723 `((mfactorial) ,(caddr e
)))))
727 (cond ((not (involve ex
'(%sin %cos %tan %asin %acos %atan
728 %sinh %cosh %tanh %asinh %acosh %atanh
)))
734 ((eq (caar ex
) '%sinh
)
735 (m// (m+ (m^
'$%e
(cadr ex
)) (m- (m^
'$%e
(m- (cadr ex
)))))
737 ((eq (caar ex
) '%cosh
)
738 (m// (m+ (m^
'$%e
(cadr ex
)) (m^
'$%e
(m- (cadr ex
))))
740 ((and (member (caar ex
)
741 '(%sin %cos %tan %asin %acos %atan %sinh
742 %cosh %tanh %asinh %acosh %atanh
) :test
#'eq
)
745 (t (cons (car ex
) (mapcar #'hyperex0
(cdr ex
))))))
750 ;;Used by tlimit also.
751 (defun limit1 (exp var val
)
753 (let ((lhprogress? lhp?
)
756 (cond ((setq ans
(and (not (atom exp
)) (getlimval exp
)))
758 ((and (not (infinityp val
)) (setq ans
(simplimsubst val exp
)))
761 ;;;NUMDEN* => (values numerator denominator)
762 (multiple-value-bind (n dn
)
764 (cond ((not (among var dn
))
765 (return (simplimit (m// (simplimit n var val
) dn
) var val
)))
767 (return (simplimit (m* n
(simplimexpt dn -
1 (simplimit dn var val
) -
1)) var val
)))
769 (/#alike n
(car lhprogress?
))
770 (/#alike dn
(cdr lhprogress?
)))
771 (throw 'lhospital nil
)))
772 (return (limit2 n dn var val
))))))
777 (let ((deriv (sdiff (m// e f
) var
)))
779 ((=0 ($ratsimp deriv
)) t
)
782 ;; The standard /#alike function is possibly somewhat inefficient. Here is
783 ;; a possible replacement:
785 ;; If e/f is free of var (a special), return true. This code
786 ;; assumes that f is not zero. First, we test if e & f
787 ;; are alike--this check is relatively fast; second, we check
788 ;; if e/f is free of var.
790 ;;(defun /#alike (e f)
791 ;; (or (alike1 e f) (freeof var (sratsimp (div e f)))))
793 (defun limit2 (n dn var val
)
794 (prog (n1 d1 lim-sign gcp sheur-ans
)
795 (setq n
(hyperex n
) dn
(hyperex dn
))
796 ;;;Change to uniform limit call.
797 (cond ((infinityp val
)
798 (setq d1
(limit dn var val nil
))
799 (setq n1
(limit n var val nil
)))
800 (t (cond ((setq n1
(simplimsubst val n
)) nil
)
801 (t (setq n1
(limit n var val nil
))))
802 (cond ((setq d1
(simplimsubst val dn
)) nil
)
803 (t (setq d1
(limit dn var val nil
))))))
804 (cond ((or (null n1
) (null d1
)) (return nil
))
805 (t (setq n1
(sratsimp n1
) d1
(sratsimp d1
))))
806 (cond ((or (involve n
'(mfactorial)) (involve dn
'(mfactorial)))
807 (let ((ans (limfact2 n dn var val
)))
808 (cond (ans (return ans
))))))
809 (cond ((and (zerop2 n1
) (zerop2 d1
))
810 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
811 (return (colexpt n dn gcp
)))
812 ((and (real-epsilonp val
)
814 (not (free dn
'%log
)))
815 (return (liminv (m// n dn
))))
816 ((setq n1
(try-lhospital-quit n dn nil
))
818 ((and (zerop2 n1
) (not (member d1
'($ind $und
) :test
#'eq
))) (return 0))
820 (setq n1
(ridofab n1
))
821 (return (simplimtimes `(,n1
,(simplimexpt dn -
1 d1 -
1))))))
822 (setq n1
(ridofab n1
))
823 (setq d1
(ridofab d1
))
824 (cond ((or (eq d1
'$und
)
825 (and (eq n1
'$und
) (not (real-infinityp d1
))))
828 ;; At this point we have n1/$ind. Look if n1 is one of the
829 ;; infinities or zero.
830 (cond ((and (infinityp n1
) (eq ($sign dn
) '$pos
))
832 ((and (infinityp n1
) (eq ($sign dn
) '$neg
))
833 (return (simpinf (m* -
1 n1
))))
835 (or (eq ($sign dn
) '$pos
)
836 (eq ($sign dn
) '$neg
)))
839 ((eq n1
'$ind
) (return (cond ((infinityp d1
) 0)
842 ((and (real-infinityp d1
) (member n1
'($inf $und $minf
) :test
#'eq
))
843 (cond ((and (not (atom dn
)) (not (atom n
))
844 (cond ((not (equal (setq gcp
(gcpower n dn
)) 1))
845 (return (colexpt n dn gcp
)))
847 (or (involve dn
'(mfactorial %gamma
))
848 (involve n
'(mfactorial %gamma
))))
849 (return (limfact n dn
))))))
850 ((eq n1 d1
) (setq lim-sign
1) (go cp
))
851 (t (setq lim-sign -
1) (go cp
))))
852 ((and (infinityp d1
) (infinityp n1
))
853 (setq lim-sign
(if (or (eq d1
'$minf
) (eq n1
'$minf
)) -
1 1))
855 (t (return (simplimtimes `(,n1
,(m^ d1 -
1))))))
856 cp
(setq n
($expand n
) dn
($expand dn
))
858 (let ((new-n (m+l
(maxi (cdr n
)))))
859 (cond ((not (alike1 new-n n
))
860 (return (limit (m// new-n dn
) var val
'think
))))
864 (let ((new-dn (m+l
(maxi (cdr dn
)))))
865 (cond ((not (alike1 new-dn dn
))
866 (return (limit (m// n new-dn
) var val
'think
))))
869 (setq sheur-ans
(sheur0 n1 d1
))
870 (cond ((or (member sheur-ans
'($inf $zeroa
) :test
#'eq
)
871 (free sheur-ans var
))
872 (return (simplimtimes `(,lim-sign
,sheur-ans
))))
873 ((and (alike1 sheur-ans dn
)
875 ((member (setq n1
(cond ((expfactorp n1 d1
) (expfactor n1 d1 var
))
877 '($inf $zeroa
) :test
#'eq
)
879 ((not (null (setq n1
(cond ((expfactorp n dn
) (expfactor n dn var
))
882 ((and (alike1 sheur-ans dn
) (not (mplusp n
))))
883 ((not (alike1 sheur-ans
(m// n dn
)))
884 (return (simplimit (m// ($expand
(m// n sheur-ans
))
885 ($expand
(m// dn sheur-ans
)))
888 (cond ((and (not (and (eq val
'$inf
) (expp n
) (expp dn
)))
889 (setq n1
(try-lhospital-quit n dn nil
))
894 ;; Test whether both n and dn have form
895 ;; product of poly^poly
896 (defun expfactorp (n dn
)
897 (do ((llist (append (cond ((mtimesp n
) (cdr n
))
899 (cond ((mtimesp dn
) (cdr dn
))
902 (exp? t
) ;IS EVERY ELEMENT SO FAR
903 (factor nil
)) ;A POLY^POLY?
907 (setq factor
(car llist
))
908 (setq exp?
(or (polyinx factor var
())
910 (polyinx (cadr factor
) var
())
911 (polyinx (caddr factor
) var
()))))))
913 (defun expfactor (n dn var
) ;Attempts to evaluate limit by grouping
914 (prog (highest-deg) ; terms with similar exponents.
915 (let ((new-exp (exppoly n
))) ;exppoly unrats expon
916 (setq n
(car new-exp
) ;and rtns deg of expons
917 highest-deg
(cdr new-exp
)))
918 (cond ((null n
) (return nil
))) ;nil means expon is not
919 (let ((new-exp (exppoly dn
))) ;a rat func.
920 (setq dn
(car new-exp
)
921 highest-deg
(max highest-deg
(cdr new-exp
))))
923 (= highest-deg
0)) ; prevent infinite recursion
927 (degree highest-deg
(1- degree
))
934 (limit (m// numerator denominator
)
938 (let ((newnumer-factor (get-newexp&factors
942 (setq numerator
(car newnumer-factor
)
943 numfactors
(cdr newnumer-factor
)))
944 (let ((newdenom-factor (get-newexp&factors
948 (setq denominator
(car newdenom-factor
)
949 denfactors
(cdr newdenom-factor
)))
950 (setq answer
(simplimit (list '(mexpt)
952 (m// numfactors denfactors
))
953 (cond ((> degree
0) var
)
957 (cond ((member answer
'($ind $und
) :test
#'equal
)
958 ;; cannot handle limit(exp(x*%i)*x, x, inf);
960 ((member answer
'($inf $minf
) :test
#'equal
)
961 ;; 0, zeroa, zerob are passed through to next iteration
962 (return (simplimtimes (list (m// numerator denominator
) answer
)))))))))
964 (defun exppoly (exp) ;RETURNS EXPRESSION WITH UNRATTED EXPONENTS
968 (exp (cond ((mtimesp exp
)
972 ((null exp
) (cons new-exp highest-deg
))
973 (setq factor
(car exp
))
975 (m* (cond ((or (not (mexptp factor
))
976 (not (ratp (caddr factor
) var
)))
980 (ratdegree (caddr factor
))))
981 (m^
(cadr factor
) (unrat (caddr factor
)))))
984 (defun unrat (exp) ;RETURNS UNRATTED EXPRESION
985 (multiple-value-bind (n d
)
987 (let ((tem ($divide n d
)))
989 (m// (caddr tem
) d
)))))
991 (defun get-newexp&factors
(exp degree var
) ;RETURNS (CONS NEWEXP FACTORS)
992 (do ((terms (cond ((mtimesp exp
)(cdr exp
)) ; SUCH THAT
993 (t (ncons exp
))) ; NEWEXP*FACTORS^(VAR^DEGREE)
994 (cdr terms
)) ; IS EQUAL TO EXP.
1001 (setq factor
(car terms
))
1002 (cond ((not (mexptp factor
))
1004 (setq factors
(m* factor factors
)))
1005 (t (setq newexp
(m* factor newexp
)))))
1007 (= (ratdegree (caddr factor
))
1009 (setq factors
(m* (m^
(cadr factor
)
1010 (leading-coef (caddr factor
)))
1012 newexp
(m* (m^
(cadr factor
)
1014 (m* (leading-coef (caddr factor
))
1017 (t (setq newexp
(m* factor newexp
))))))
1019 (defun leading-coef (rat)
1020 (ratlim (m// rat
(m^ var
(ratdegree rat
)))))
1022 (defun ratdegree (rat)
1023 (multiple-value-bind (n d
)
1025 (- (deg n
) (deg d
))))
1027 (defun limfact2 (n d var val
)
1028 (let ((n1 (reflect0 n var val
))
1029 (d1 (reflect0 d var val
)))
1030 (cond ((and (alike1 n n1
)
1033 (t (limit (m// n1 d1
) var val
'think
)))))
1035 ;; takes expression and returns operator at front with all flags removed
1036 ;; except array flag.
1037 ;; array flag must match for alike1 to consider two things to be the same.
1038 ;; ((MTIMES SIMP) ... ) => (MTIMES)
1039 ;; ((PSI SIMP ARRAY) 0) => (PSI ARRAY)
1040 (defun operator-with-array-flag (exp)
1041 (cond ((member 'array
(car exp
) :test
#'eq
)
1042 (list (caar exp
) 'array
))
1043 (t (list (caar exp
)))))
1045 (defun reflect0 (exp var val
)
1046 (cond ((atom exp
) exp
)
1047 ((and (eq (caar exp
) 'mfactorial
)
1048 (let ((argval (limit (cadr exp
) var val
'think
)))
1049 (or (eq argval
'$minf
)
1050 (and (numberp argval
)
1052 (reflect (cadr exp
)))
1053 (t (cons (operator-with-array-flag exp
)
1056 (reflect0 term var val
)))
1059 (defun reflect (arg)
1062 (m^
(list (ncons 'mfactorial
)
1066 (m^
(list (ncons '%sin
)
1070 (defun limfact (n d
)
1072 (setq n
(stirling0 n
)
1074 (setq ans
(toplevel-$limit
(m// n d
) var
'$inf
))
1075 (cond ((and (atom ans
)
1076 (not (member ans
'(und ind
) :test
#'eq
))) ans
)
1077 ((eq (caar ans
) '%limit
) ())
1080 ;; substitute asymptotic approximations for gamma, factorial, and
1082 (defun stirling0 (e)
1084 ((and (setq e
(cons (car e
) (mapcar 'stirling0
(cdr e
))))
1086 ((and (eq (caar e
) '%gamma
)
1087 (eq (limit (cadr e
) var val
'think
) '$inf
))
1088 (stirling (cadr e
)))
1089 ((eq (caar e
) 'mfactorial
)
1090 (let ((n (limit (cadr e
) var val
'think
)))
1092 (m* (cadr e
) (stirling (cadr e
))))
1093 ((and (integerp n
) (< n
0))
1095 (div (power -
1 n
) (mul (ftake 'mfactorial n
) (add var n
))))
1097 ((and (eq (caar e
) 'mqapply
) ;; polylogarithm
1098 (eq (subfunname e
) '$li
)
1099 (let ((arglim (limit (car (subfunargs e
)) var val
'think
)))
1100 (or (eq arglim
'$inf
) (eq arglim
'$minf
)))
1101 (integerp (car (subfunsubs e
))))
1102 (li-asymptotic-expansion (m- (car (subfunsubs e
)) 1)
1103 (car (subfunsubs e
))
1104 (car (subfunargs e
))))
1108 (maxima-substitute x
'$z
1110 ((mexpt simp
) 2 ((rat simp
) 1 2))
1111 ((mexpt simp
) $%pi
((rat simp
) 1 2))
1112 ((mexpt simp
) $z
((mplus simp
) ((rat simp
) -
1 2) $z
))
1113 ((mexpt simp
) $%e
((mtimes simp
) -
1 $z
)))))
1115 (defun no-err-sub (v e
&aux ans
)
1116 (let ((errorsw t
) (*zexptsimp? t
)
1118 ;; Don't print any error messages
1120 ;; Should we just use IGNORE-ERRORS instead HANDLER-CASE here? I
1121 ;; (rtoy) am choosing the latter so that unexpected errors will
1122 ;; actually show up instead of being silently discarded.
1124 (setq ans
(catch 'errorsw
1126 (sratsimp (subin v e
)))))
1129 (cond ((null ans
) t
) ; Ratfun package returns NIL for failure.
1132 (defun extended-real-p (e)
1133 (member e
(list '$minf
'$zerob
'$zeroa
'$ind
'$und
'$inf
'$infinity
)))
1135 ;; Evaluate the limit of e as var (special) approaches v using direct substitution.
1136 ;; This function is the last of the chain of methods tried by limit. As such it
1137 ;; shouldn't call higher level functions such as simplimit or limit--doing so
1138 ;; would risk creating an infinite loop.
1140 ;; This function special cases sums, products, and powers. It declines to use
1141 ;; direct substitution on any expression whose main operator has a simplim%function
1142 ;; function. Generally, limits of functions that have a simplim%function should be
1143 ;; handled by those specialized functions, not by simplimsubst. Having said this,
1144 ;; it's OK for simplimsubst to special case an operator and do direct substitution
1145 ;; when its OK. The cases of log, cosine, and sine happen often enough that these
1146 ;; functions are special cased.
1148 ;; Possibly simplimsubst should decline to do direct substitution on functions
1149 ;; that have a limit function, for example inverse_jacobi_ns.
1151 ;; For all other kinds of expressions, this function assumes that if substitution
1152 ;; doesn't result in an error (division by zero, for example), the function
1153 ;; is continuous at the limit point. That's dodgy. This dodginess underscores the
1154 ;; need to define a simplim%function functions that are not continuous on their
1157 ;; This function sometimes receives an un-simplified expression. Maybe they should be
1158 ;; simplified, maybe not.
1160 ;; Locally setting $numer and $%enumer to nil keeps some limits, for example
1161 ;; limit(sin(x)/x,x,0) from returning 1.0 when numer is true. (Barton Willis)
1163 (defun simplimsubst (v e
)
1164 (let ((ans nil
) ($numer nil
) ($%enumer nil
) (ee))
1166 ;; When e is a number, return it.
1168 ;; When e is a mapatom, substitute v for var and return.
1169 (($mapatom e
) ($substitute v var e
))
1170 ;; Special case mexpt expressions. Decline direct substitution for
1172 ((and (mexptp e
) (not (extended-real-p v
)))
1173 (let ((x (simplimsubst v
(second e
)))
1174 (n (simplimsubst v
(third e
))))
1175 ;; Decline direct substitution (DS) for 0^negative. Also decline
1176 ;; DS when x is on the negative real axis and n isn't an integer.
1177 ;; Additionally, we require that DS is OK for both x & n.
1178 (if (and x n
(not (and (zerop2 x
) (eq t
(mgqp 0 n
)))) ; not 0^negative
1179 (or (off-negative-real-axisp x
) (integerp n
)))
1180 (ftake 'mexpt x n
) nil
)))
1181 ;; Special case product and sum expressions. Again, we decline direct
1182 ;; substitution for extended reals.
1183 ((and (or (mplusp e
) (mtimesp e
)) (not (extended-real-p v
)))
1184 (setq ee
(mapcar #'(lambda(q) (simplimsubst v q
)) (cdr e
)))
1185 (if (some #'(lambda (q) (eq q nil
)) ee
) nil
1186 (simplifya (cons (list (caar e
)) ee
) t
)))
1187 ;; Decline direct substitution for sums, products, and powers for
1188 ;; the extended real case.
1189 ((and (or (mexptp e
) (mplusp e
) (mtimesp e
)) (extended-real-p v
))
1191 ;; Special case log expressions--possibly the log case happens
1192 ;; often enough to make this worthwhile.
1193 ((and (consp e
) (consp (car e
)) (eq '%log
(caar e
)))
1194 (let ((w (simplimsubst v
(cadr e
))))
1195 (if (and w
(off-negative-real-axisp w
)) (ftake '%log w
) nil
)))
1196 ;; Special case %cos and %sin expressions--we could special case others.
1197 ;; The lenient-realp check declines direct substitution for complex arguments--
1198 ;; this check could be relaxed.
1199 ((and (consp e
) (consp (car e
)) (or (member (caar e
) (list '%cos
'%sin
))))
1200 (let ((op (caar e
)))
1201 (setq e
(simplimsubst v
(second e
)))
1202 (if (and e
(lenient-realp e
) (not (extended-real-p e
))) (ftake op e
) nil
)))
1203 ;; Don't use direct substitution on expressions whose main operator has
1204 ;; a simplim%function.
1205 ((and (consp e
) (consp (car e
)) (get (caar e
) 'simplim%function
)) nil
)
1206 ;; The function no-err-sub returns true when there is an error
1207 ((not (eq t
(setq ans
(no-err-sub (ridofab v
) e
))))
1208 ;; Previously the condition (zerop2 ans) was (=0 ($radcan ans)). In
1209 ;; December 2021, the testsuite + the share testsuite only gives ans = 0,
1210 ;; making the radcan unneeded.
1211 (cond ((and (member v
'($zeroa $zerob
) :test
#'eq
) (zerop2 ans
))
1212 (setq ans
(behavior e var v
))
1213 (cond ((eql ans
1) '$zeroa
)
1214 ((eql ans -
1) '$zerob
)
1215 (t nil
))) ;behavior can't find direction
1217 ;; direct substitution fails, so return nil.
1220 ;;;returns (cons numerator denominator)
1222 (let ((e (factor (simplify e
)))
1228 (mapc #'forq
(cdr e
)))
1234 (setq numer
(car numer
)))
1236 (setq numer
(m*l numer
))))
1240 (setq denom
(car denom
)))
1242 (setq denom
(m*l denom
))))
1243 (values (factor numer
) (factor denom
))))
1245 ;;;FACTOR OR QUOTIENT
1246 ;;;Setq's the special vars numer and denom from numden*
1248 (cond ((and (mexptp e
)
1249 (not (freeof var e
))
1250 (null (pos-neg-p (caddr e
))))
1251 (push (m^
(cadr e
) (m* -
1.
(caddr e
))) denom
))
1252 (t (push e numer
))))
1254 ;;;Predicate to tell whether an expression is pos,zero or neg as var -> val.
1255 ;;;returns T if pos,zero. () if negative or don't know.
1256 (defun pos-neg-p (exp)
1257 (let ((ans (limit exp var val
'think
)))
1258 (cond ((and (not (member ans
'($und $ind $infinity
) :test
#'eq
))
1259 (equal ($imagpart ans
) 0))
1260 (let ((sign (getsignl ans
)))
1261 (cond ((or (equal sign
1)
1264 ((equal sign -
1) nil
))))
1267 (declare-top (unspecial n dn
))
1270 (cond ((radicalp e var
) nil
)
1271 ((member (caar e
) '(%log %sin %cos %tan %sinh %cosh %tanh mfactorial
1272 %asin %acos %atan %asinh %acosh %atanh
) :test
#'eq
) nil
)
1274 ((do ((e (cdr e
) (cdr e
)))
1276 (and (expp (car e
)) (return t
))))))
1280 (radicalp (cadr e
) var
)
1281 (among var
(caddr e
))
1282 (radicalp (caddr e
) var
)))
1285 (defun gcpower (a b
)
1286 ($gcd
(getexp a
) (getexp b
)))
1289 (cond ((and (mexptp exp
)
1290 (free (caddr exp
) var
)
1291 (eq (ask-integer (caddr exp
) '$integer
) '$yes
))
1293 ((mtimesp exp
) (getexplist (cdr exp
)))
1296 (defun getexplist (list)
1297 (cond ((null (cdr list
))
1298 (getexp (car list
)))
1299 (t ($gcd
(getexp (car list
))
1300 (getexplist (cdr list
))))))
1302 (defun limroot (exp power
)
1303 (cond ((or (atom exp
) (not (member (caar exp
) '(mtimes mexpt
) :test
#'eq
)))
1304 (limroot (list '(mexpt) exp
1) power
)) ;This is strange-JIM.
1305 ((mexptp exp
) (m^
(cadr exp
)
1306 (sratsimp (m* (caddr exp
) (m^ power -
1.
)))))
1307 (t (m*l
(mapcar #'(lambda (x)
1311 ;;NUMERATOR AND DENOMINATOR HAVE EXPONENTS WITH GCD OF GCP.
1312 ;;; Used to call simplimit but some of the transformations used here
1313 ;;; were not stable w.r.t. the simplifier, so try keeping exponent separate
1316 (defun colexpt (n dn gcp
)
1317 (let ((bas (m* (limroot n gcp
) (limroot dn
(m* -
1 gcp
))))
1320 (setq baslim
(limit bas var val
'think
))
1321 (setq expolim
(limit expo var val
'think
))
1322 (simplimexpt bas expo baslim expolim
)))
1324 ;; this function is responsible for the following bug:
1325 ;; limit(x^2 + %i*x, x, inf) -> inf (should be infinity)
1327 (setq e
(sratsimp ($trigreduce e
)))
1328 (cond ((member val
'($inf $infinity
) :test
#'eq
)
1329 (setq e
(maxima-substitute (m^t
'x -
1) var e
)))
1331 (setq e
(maxima-substitute (m^t -
1 (m^t
'x -
1)) var e
)))
1333 (setq e
(maxima-substitute (m- 'x
) var e
)))
1335 (setq e
(maxima-substitute 'x var e
)))
1336 ((setq e
(maxima-substitute (m+t
'x val
) var e
))))
1337 (destructuring-let* ((e (let (($ratfac
()))
1338 ($rat
(sratsimp e
) 'x
)))
1346 (sratsimp (m// ($ratdisrep
`(,h
,(locoef n g
) .
1))
1347 ($ratdisrep
`(,h
,(locoef d g
) .
1))))))
1349 (cond ((not (member val
'($zerob $zeroa $inf $minf
) :test
#'eq
))
1351 ((not (equal ($imagpart e
) 0))
1353 ((null (setq e
(getsignl ($realpart e
))))
1355 ((equal e
1) '$zeroa
)
1356 ((equal e -
1) '$zerob
)
1359 ((not (member val
'($zerob $zeroa $infinity $inf $minf
) :test
#'eq
))
1361 ((eq val
'$infinity
) '$infinity
)
1362 ((not (equal ($imagpart e
) 0)) '$infinity
)
1363 ((null (setq e
(getsignl ($realpart e
))))
1366 ((equal e -
1) '$minf
)
1370 (if (or (atom n
) (not (eq (car n
) x
)))
1375 (if (or (atom n
) (not (eq (car n
) x
)))
1379 ;; This function tries to determine the increasing/decreasing
1380 ;; behavior of an expression exp with respect to some variable var.
1381 ;; val determines the mode:
1382 ;; - $zeroa: From "positive zero" towards the right
1383 ;; - $zerob: From "negative zero" towards the left
1384 ;; - $inf: From "positive infinity" towards the left
1385 ;; - $minf: From "negative infinity" towards the right
1386 ;; Return values are -1 (decreasing), +1 (increasing) or 0 (don't know).
1387 (defun behavior (exp var val
)
1388 ;; $inf/$minf is handled by substituting 1/var for var
1389 ;; and setting val to $zeroa/$zerob.
1390 (cond ((real-infinityp val
)
1392 ((eq val
'$inf
) '$zeroa
)
1393 ((eq val
'$minf
) '$zerob
))
1394 exp
(sratsimp (subin (m^ var -
1) exp
)))))
1396 ((not (member val
'($zeroa $zerob $inf $minf
)))
1398 ;; Shortcut for constants.
1401 ;; Shortcut for the variable itself.
1403 (if (member val
'($zeroa $minf
) :test
#'eq
) 1 -
1))
1404 ;; This limits the recursion depth of the function.
1405 ;; Before giving up, always try behavior-by-diff, which doesn't recurse.
1406 ((= *behavior-count-now
* +behavior-count
+)
1407 (behavior-by-diff exp var val
))
1409 (let ((*behavior-count-now
* (1+ *behavior-count-now
*)) pair sign ans
)
1411 ;; Pull out constant factors with known signs from a product:
1412 ;; behavior(c * f(x)) = signum(c) * behavior(f(x))
1414 (prog2 (setq pair
(partition exp var
1))
1415 (not (equal (car pair
) 1))))
1416 (setq sign
(getsignl (car pair
)))
1417 (if (not (fixnump sign
))
1419 (mul sign
(behavior (cdr pair
) var val
))))
1420 ;; Pull out constant terms from a sum:
1421 ;; behavior(c + f(x)) = behavior(f(x))
1423 (prog2 (setq pair
(partition exp var
0))
1424 (not (equal (car pair
) 0))))
1425 (behavior (cdr pair
) var val
))
1426 ;; Handle f(x)^c for the case f(0) = 0 and c > 0
1427 ;; (probably other cases could be handled, too)
1429 (free (caddr exp
) var
)
1430 (=0 (no-err-sub 0 exp
))
1431 (equal (getsignl (caddr exp
)) 1)
1432 (not (equal 0 (setq ans
(behavior-expt (cadr exp
) (caddr exp
))))))
1435 ;; behavior(1 / f(x)) = -behavior(f(x))
1436 ((equal ($num exp
) 1)
1437 (- (behavior ($denom exp
) var val
)))
1438 ;; Handle c^f(x) for c > 1:
1439 ;; behavior(c^f(x)) = behavior(f(x))
1441 (free (cadr exp
) var
)
1442 (equal 1 (getsignl (m- (cadr exp
) 1)))
1443 (not (equal 0 (setq ans
(behavior (caddr exp
) var val
)))))
1445 ;; Handle c^f(x) for 0 < c < 1:
1446 ;; behavior(c^f(x)) = -behavior(f(x))
1448 (free (cadr exp
) var
)
1449 (equal 1 (getsignl (cadr exp
)))
1450 (equal -
1 (getsignl (m- (cadr exp
) 1)))
1451 (not (equal 0 (setq ans
(behavior (caddr exp
) var val
)))))
1453 ;; atan, erf, sinh and tanh are monotonically increasing,
1454 ;; so their behavior is equal to the behavior of their arguments.
1455 ;; behavior(monotonically_increasing(f(x))) = behavior(f(x))
1456 ((member (caar exp
) '(%atan %erf %sinh %tanh
) :test
#'eq
)
1457 (behavior (cadr exp
) var val
))
1458 ;; Similarly, acot is monotonically decreasing left and right of x = 0.
1459 ;; The singularity doesn't matter for our purposes.
1460 ;; behavior(monotonically_decreasing(f(x))) = -behavior(f(x))
1461 ((eq (caar exp
) '%acot
)
1462 (- (behavior (cadr exp
) var val
)))
1463 ;; Note: More functions could be added here.
1464 ;; Ideally, use properties defined on the functions.
1465 ;; Handle log(f(x)) for f(0) = 0:
1466 ;; If behavior(f(x)) = 1, then behavior(log(f(x))) = 1
1468 (=0 (no-err-sub 0 (cadr exp
)))
1469 (equal 1 (behavior (cadr exp
) var val
)))
1471 ;; Try to determine the behavior by taking the derivative.
1472 ((not (equal 0 (setq ans
(behavior-by-diff exp var val
))))
1474 ;; If exp is a sum and all terms have the same behavior, return that.
1475 ;; The sum of increasing functions is increasing.
1476 ;; The sum of decreasing functions is decreasing.
1478 (not (equal 0 (setq ans
(behavior-all-same exp var val
)))))
1482 (defun behavior-all-same (exp var val
)
1483 (let* ((exps (cdr exp
)) (first-behavior (behavior (first exps
) var val
)))
1484 (if (or (equal first-behavior
0)
1485 (not (every #'(lambda (exp) (equal (behavior exp var val
) first-behavior
))
1490 (defun behavior-expt (bas expo
)
1491 (let ((behavior (behavior bas var val
)))
1492 (cond ((= behavior
1) 1)
1494 ((eq (ask-integer expo
'$integer
) '$yes
)
1495 (cond ((eq (ask-integer expo
'$even
) '$yes
) 1)
1498 (cond ((evenp (cadr expo
)) 1)
1499 ((oddp (caddr expo
)) behavior
)
1503 (defun behavior-by-diff (exp var val
)
1505 (exp (sratsimp (sdiff exp var
)) (sratsimp (sdiff exp var
)))
1507 (ans ())) ; This do wins by a return.
1508 ((> ct
1) 0) ; This loop used to run up to 5 times,
1509 ;; but the size of some expressions would blow up.
1510 (setq ans
(no-err-sub 0 exp
)) ;Why not do an EVENFN and ODDFN
1514 ((=0 ans
) ()) ;Do it again.
1515 (t (setq ans
(getsignl ans
))
1516 (cond (n (return ans
))
1518 (return (if (eq val
'$zeroa
) 1 -
1)))
1520 (return (if (eq val
'$zeroa
) -
1 1)))
1523 (defun try-lhospital (n d ind
)
1524 ;;Make one catch for the whole bunch of lhospital trials.
1525 (let ((ans (lhospital-catch n d ind
)))
1526 (cond ((null ans
) ())
1527 ((not (free-infp ans
)) (simpinf ans
))
1528 ((not (free-epsilonp ans
)) (simpab ans
))
1531 (defun try-lhospital-quit (n d ind
)
1532 (let ((ans (or (lhospital-catch n d ind
)
1533 (lhospital-catch (m^ d -
1) (m^ n -
1) ind
))))
1534 (cond ((null ans
) (throw 'limit t
))
1535 ((not (free-infp ans
)) (simpinf ans
))
1536 ((not (free-epsilonp ans
)) (simpab ans
))
1539 (defun lhospital-catch (n d ind
)
1540 (cond ((> 0 lhcount
)
1541 (setq lhcount $lhospitallim
)
1542 (throw 'lhospital nil
))
1543 ((equal lhcount $lhospitallim
)
1544 (let ((lhcount (m+ lhcount -
1)))
1545 (catch 'lhospital
(lhospital n d ind
))))
1546 (t (setq lhcount
(m+ lhcount -
1))
1547 (prog1 (lhospital n d ind
)
1548 (setq lhcount
(m+ lhcount
1))))))
1549 ;;If this succeeds then raise LHCOUNT.
1552 (defun lhospital (n d ind
)
1553 (declare (special val lhp?
))
1555 (setq n
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr n
)))))
1557 (setq d
(m*l
(mapcar #'(lambda (term) (lhsimp term var val
)) (cdr d
)))))
1558 (multiple-value-bind (n d
)
1560 (let (const nconst dconst
)
1561 (setq lhp?
(and (null ind
) (cons n d
)))
1562 (multiple-value-setq (nconst n
) (var-or-const n
))
1563 (multiple-value-setq (dconst d
) (var-or-const d
))
1565 (setq n
(stirling0 n
)) ;; replace factorial and %gamma
1566 (setq d
(stirling0 d
)) ;; with approximations
1568 (setq n
(sdiff n var
) ;; take derivatives for l'hospital
1571 (if (or (not (free n
'%derivative
)) (not (free d
'%derivative
)))
1572 (throw 'lhospital
()))
1573 (setq n
(expand-trigs (tansc n
) var
))
1574 (setq d
(expand-trigs (tansc d
) var
))
1576 (multiple-value-setq (const n d
) (remove-singularities n d
))
1577 (setq const
(m* const
(m// nconst dconst
)))
1578 (simpinf (let ((ans (if ind
1579 (limit2 n d var val
)
1580 (limit-numden n d val
))))
1581 ;; When the limit function returns, it's possible that it will return NIL
1582 ;; (gave up without finding a limit). It's also possible that it will
1583 ;; return something containing UND. We treat that as a failure too.
1584 (when (and ans
(freeof '$und ans
))
1585 (sratsimp (m* const ans
))))))))
1587 ;; Try to compute the limit of a quotient NUM/DEN, trying to massage the input
1588 ;; into a convenient form for LIMIT on the way.
1589 (defun limit-numden (n d val
)
1591 ;; For general arguments, the best approach seems to be to use
1592 ;; sratsimp to simplify the quotient as much as we can, then
1593 ;; $multthru, which splits it up into a sum (presumably useful
1594 ;; because limit(a+b) = limit(a) + limit(b) if the limits exist, and
1595 ;; the right hand side might be easier to calculate)
1597 ($multthru
(sratsimp (m// n d
))))
1599 ;; If we've already got a sum in the numerator, it seems to be
1600 ;; better not to recombine it. Call LIMIT on the whole lot, though,
1601 ;; because terms with infinite limits might cancel to give a finite
1604 (m+l
(mapcar #'(lambda (x)
1605 (sratsimp (m// x d
)))
1608 (limit expr var val
'think
)))
1610 ;; Heuristics for picking the right way to express a LHOSPITAL problem.
1611 (defun lhop-numden (num denom
)
1612 (declare (special var
))
1613 (cond ((let ((log-num (involve num
'(%log
))))
1614 (cond ((null log-num
) ())
1615 ((lessthan (num-of-logs (factor (sratsimp (sdiff (m^ num -
1) var
))))
1616 (num-of-logs (factor (sratsimp (sdiff num var
)))))
1617 (psetq num
(m^ denom -
1) denom
(m^ num -
1))
1620 ((let ((log-denom (involve denom
'(%log
))))
1621 (cond ((null log-denom
) ())
1622 ((lessthan (num-of-logs (sratsimp (sdiff (m^ denom -
1) var
)))
1623 (num-of-logs (sratsimp (sdiff denom var
))))
1624 (psetq denom
(m^ num -
1) num
(m^ denom -
1))
1627 ((let ((exp-num (%einvolve num
)))
1629 (cond ((%e-right-placep exp-num
)
1631 (t (psetq num
(m^ denom -
1)
1632 denom
(m^ num -
1)) t
)))
1634 ((let ((exp-den (%einvolve denom
)))
1636 (cond ((%e-right-placep exp-den
)
1638 (t (psetq num
(m^ denom -
1)
1639 denom
(m^ num -
1)) t
)))
1641 ((let ((scnum (involve num
'(%sin
))))
1642 (cond (scnum (cond ((trig-right-placep '%sin scnum
) t
)
1643 (t (psetq num
(m^ denom -
1)
1644 denom
(m^ num -
1)) t
)))
1646 ((let ((scden (involve denom
'(%sin
))))
1647 (cond (scden (cond ((trig-right-placep '%sin scden
) t
)
1648 (t (psetq num
(m^ denom -
1)
1649 denom
(m^ num -
1)) t
)))
1651 ((let ((scnum (involve num
'(%asin %acos %atan
))))
1652 ;; If the numerator contains an inverse trig and the
1653 ;; denominator or reciprocal of denominator is polynomial,
1654 ;; leave everything as is. If the inverse trig is moved to
1655 ;; the denominator, things get messy, even if the numerator
1656 ;; becomes a polynomial. This is not perfect.
1657 (cond ((and scnum
(or (polyinx denom var
())
1658 (polyinx (m^ denom -
1) var
())))
1661 ((or (oscip num
) (oscip denom
)))
1663 (psetq num
(m^ denom -
1) denom
(m^ num -
1))))
1666 ;;i don't know what to do here for some cases, may have to be refined.
1667 (defun num-of-logs (exp)
1668 (cond ((mapatom exp
) 0)
1669 ((equal (caar exp
) '%log
)
1670 (m+ 1 (num-of-log-l (cdr exp
))))
1671 ((and (mexptp exp
) (mnump (caddr exp
)))
1672 (m* (simplify `((mabs) ,(caddr exp
)))
1673 (num-of-logs (cadr exp
))))
1674 (t (num-of-log-l (cdr exp
)))))
1676 (defun num-of-log-l (llist)
1677 (do ((temp llist
(cdr temp
)) (ans 0))
1679 (setq ans
(m+ ans
(num-of-logs (car temp
))))))
1681 (defun %e-right-placep
(%e-arg
)
1682 (let ((%e-arg-diff
(sdiff %e-arg var
)))
1684 ((free %e-arg-diff var
)) ;simple cases
1685 ((or (and (mexptp denom
)
1686 (equal (cadr denom
) -
1))
1687 (polyinx (m^ denom -
1) var
())) ())
1688 ((let ((%e-arg-diff-lim
(ridofab (limit %e-arg-diff var val
'think
)))
1689 (%e-arg-exp-lim
(ridofab (limit (m^
'$%e %e-arg
) var val
'think
))))
1692 (format t
"%e-arg-dif-lim = ~A~%" %e-arg-diff-lim
)
1693 (format t
"%e-arg-exp-lim = ~A~%" %e-arg-exp-lim
))
1694 (cond ((equal %e-arg-diff-lim %e-arg-exp-lim
)
1696 ((and (mnump %e-arg-diff-lim
) (mnump %e-arg-exp-lim
))
1698 ((and (mnump %e-arg-diff-lim
) (infinityp %e-arg-exp-lim
))
1699 ;; This is meant to make maxima handle bug 1469411
1700 ;; correctly. Undoubtedly, this needs work.
1704 (defun trig-right-placep (trig-type arg
)
1705 (let ((arglim (ridofab (limit arg var val
'think
)))
1706 (triglim (ridofab (limit `((,trig-type
) ,arg
) var val
'think
))))
1707 (cond ((and (equal arglim
0) (equal triglim
0)) t
)
1708 ((and (infinityp arglim
) (infinityp triglim
)) t
)
1711 ;;Takes a numerator and a denominator. If they tries all combinations of
1712 ;;products to try and make a simpler set of subproblems for LHOSPITAL.
1713 (defun remove-singularities (numer denom
)
1714 (cond ((or (null numer
) (null denom
)
1715 (atom numer
) (atom denom
)
1716 (not (mtimesp numer
)) ;Leave this here for a while.
1717 (not (mtimesp denom
)))
1718 (values 1 numer denom
))
1721 (multiple-value-bind (num-consts num-vars
)
1722 (var-or-const numer
)
1723 (multiple-value-bind (denom-consts denom-vars
)
1724 (var-or-const denom
)
1725 (if (not (mtimesp num-vars
))
1726 (setq num-vars
(list num-vars
))
1727 (setq num-vars
(cdr num-vars
)))
1728 (if (not (mtimesp denom-vars
))
1729 (setq denom-vars
(list denom-vars
))
1730 (setq denom-vars
(cdr denom-vars
)))
1731 (do ((nl num-vars
(cdr nl
))
1732 (num-list (copy-list num-vars
))
1733 (den-list denom-vars den-list-temp
)
1734 (den-list-temp (copy-list denom-vars
)))
1735 ((null nl
) (values (m* const
(m// num-consts denom-consts
))
1737 (m*l den-list-temp
)))
1738 (do ((dl den-list
(cdr dl
)))
1740 (if (or (%einvolve
(car nl
)) (%einvolve
(car nl
)))
1742 (let ((lim (catch 'limit
(simpinf (simpab (limit (m// (car nl
) (car dl
))
1743 var val
'think
))))))
1744 (cond ((or (eq lim t
)
1746 (equal (ridofab lim
) 0)
1748 (not (free lim
'$inf
))
1749 (not (free lim
'$minf
))
1750 (not (free lim
'$infinity
))
1751 (not (free lim
'$ind
))
1752 (not (free lim
'$und
)))
1755 (setq const
(m* lim const
))
1756 (setq num-list
(delete (car nl
) num-list
:count
1 :test
#'equal
))
1757 (setq den-list-temp
(delete (car dl
) den-list-temp
:count
1 :test
#'equal
))
1758 (return t
)))))))))))))
1760 ;; separate terms that contain var from constant terms
1761 ;; returns (const-terms . var-terms)
1762 (defun var-or-const (expr)
1763 (setq expr
($factor expr
))
1771 (do ((l (cdr expr
) (cdr l
))
1774 ((null l
) (values const varl
))
1775 (if (free (car l
) var
)
1776 (setq const
(m* (car l
) const
))
1777 (setq varl
(m* (car l
) varl
)))))
1781 ;; if term goes to non-zero constant, replace with constant
1782 (defun lhsimp (term var val
)
1783 (cond (($mapatom term
) term
)
1785 (let ((term-value (ridofab (limit term var val
'think
))))
1786 (if (and (not (member term-value
'($inf $minf $und $ind $infinity
)))
1787 (eq t
(mnqp term-value
0))) term-value term
)))))
1789 (defun bylog (expo bas
)
1792 (try-lhospital-quit (simplify `((%log
) ,(tansc bas
)))
1797 (defun simplimexpt (bas expo bl el
)
1798 (cond ((or (eq bl
'$und
) (eq el
'$und
)) '$und
)
1800 (cond ((eq el
'$inf
) (if (eq bl
'$zeroa
) bl
0))
1801 ((eq el
'$minf
) (if (eq bl
'$zeroa
) '$inf
'$infinity
))
1802 ((eq el
'$ind
) '$ind
)
1803 ((eq el
'$infinity
) '$und
)
1804 ((zerop2 el
) (bylog expo bas
))
1805 (t (cond ((equal (getsignl el
) -
1)
1806 (cond ((eq bl
'$zeroa
) '$inf
)
1808 (cond ((even1 el
) '$inf
)
1809 ((eq (ask-integer el
'$integer
) '$yes
)
1810 (if (eq (ask-integer el
'$even
) '$yes
)
1812 '$minf
)))) ;Gotta be ODD.
1813 (t (setq bas
(behavior bas var val
))
1814 (cond ((equal bas
1) '$inf
)
1815 ((equal bas -
1) '$minf
)
1816 (t (throw 'limit t
))))))
1818 (member bl
'($zeroa $zerob
) :test
#'eq
))
1819 (cond ((even1 el
) '$zeroa
)
1820 ((and (eq bl
'$zerob
)
1822 (evenp (caddr el
))) 0)
1824 ((equal (getsignl el
) 1)
1825 (if (eq bl
'$zeroa
) bl
0))
1826 ((equal (getsignl el
) 0)
1828 (t (throw 'limit t
))))))
1830 (cond ((zerop2 el
) (bylog expo bas
))
1832 ((eq el
'$inf
) '$infinity
)
1833 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1834 ((equal (setq el
(getsignl el
)) 1) '$infinity
)
1837 (t (throw 'limit t
))))
1839 (cond ((eq el
'$inf
) '$inf
)
1840 ((equal el
'$minf
) 0)
1841 ((zerop2 el
) (bylog expo bas
))
1842 ((member el
'($infinity $ind
) :test
#'eq
) '$und
)
1843 (t (cond ((eql 0 (getsignl el
)) 1)
1844 ((ratgreaterp 0 el
) '$zeroa
)
1845 ((ratgreaterp el
0) '$inf
)
1846 (t (throw 'limit t
))))))
1848 (cond ((zerop2 el
) (bylog expo bas
))
1849 ((eq el
'$inf
) '$und
)
1850 ((equal el
'$minf
) 0)
1851 ;;;Why not generalize this. We can ask about the number. -Jim 2/23/81
1852 ((mnump el
) (cond ((mnegp el
)
1855 (if (eq (ask-integer el
'$integer
) '$yes
)
1856 (if (eq (ask-integer el
'$even
) '$yes
)
1860 (t (cond ((even1 el
) '$inf
)
1861 ((eq (ask-integer el
'$integer
) '$yes
)
1862 (if (eq (ask-integer el
'$even
) '$yes
)
1866 (loginprod?
(throw 'lip?
'lip
!))
1868 ((equal (simplify (ratdisrep (ridofab bl
))) 1)
1869 (if (infinityp el
) (bylog expo bas
) 1))
1870 ((and (equal (ridofab bl
) -
1)
1871 (infinityp el
)) '$ind
) ;LB
1872 ((eq bl
'$ind
) (cond ((or (zerop2 el
) (infinityp el
)) '$und
)
1873 ((not (equal (getsignl el
) -
1)) '$ind
)
1875 ((eq el
'$inf
) (cond ((abeq1 bl
)
1876 (if (equal (getsignl bl
) 1) 1 '$ind
))
1878 (if (equal (getsignl bl
) 1) '$zeroa
0))
1879 ((equal (getsignl (m1- bl
)) 1) '$inf
)
1880 ((equal (getsignl (m1- `((mabs) ,bl
))) 1) '$infinity
)
1881 (t (throw 'limit t
))))
1882 ((eq el
'$minf
) (cond ((abeq1 bl
)
1883 (if (equal (getsignl bl
) 1) 1 '$ind
))
1885 (if (equal (getsignl bl
) 1) '$zeroa
0))
1886 ((ratgreaterp 0 bl
) '$infinity
)
1889 (if (equal val
'$infinity
)
1890 '$und
;Not enough info to do anything.
1891 (destructuring-bind (real-el . imag-el
)
1893 (setq real-el
(limit real-el var origval nil
))
1894 (cond ((eq real-el
'$minf
)
1896 ((and (eq real-el
'$inf
)
1897 (not (equal (ridofab (limit imag-el var origval nil
)) 0)))
1899 ((eq real-el
'$infinity
)
1900 (throw 'limit t
)) ;; don't really know real component
1904 ((eq el
'$ind
) '$ind
)
1909 (cond ((numberp x
) (and (integerp x
) (evenp x
)))
1910 ((and (mnump x
) (evenp (cadr x
))))))
1912 ;; is absolute value less than one?
1916 (and (ratgreaterp 1. bl
) (ratgreaterp bl -
1.
)))
1917 (t (equal (getsignl (m1- `((mabs) ,bl
))) -
1.
))))
1919 ;; is absolute value equal to one?
1923 (or (equal 1. bl
) (equal bl -
1.
)))
1924 (t (equal (getsignl (m1- `((mabs) ,bl
))) 0))))
1926 (defun simplimit (exp var val
&aux op
)
1929 ((or (atom exp
) (mnump exp
)) exp
)
1930 ;; Lookup and dispatch a simplim%function from the property list
1931 ((setq op
(safe-get (mop exp
) 'simplim%function
))
1932 (funcall op exp var val
))
1934 ;; And do the same for subscripted:
1935 ((and (consp exp
) (consp (car exp
)) (eq (caar exp
) 'mqapply
)
1936 (setq op
(safe-get (subfunname exp
) 'simplim%function
)))
1937 (funcall op exp var val
))
1939 ;; Without the call to rootscontract,
1940 ;; limit(((-4)*x^2-10*x+24)/((4*x+8)^(1/3)+2),x,-4)
1941 ;; returns zero (should be 66). This bug is due to the fact that
1942 ;; 4^(1/3)*2^(1/3)-2 does not simplify to zero. Although the call
1943 ;; to rootscontract takes care of this case, almost surely there are
1944 ;; many other limit problems that need more than rootscontract.
1945 ((mplusp exp
) (let (($rootsconmode nil
)) ($rootscontract
(simplimplus exp
))))
1946 ((mtimesp exp
) (simplimtimes (cdr exp
)))
1947 ((mexptp exp
) (simplimexpt (cadr exp
) (caddr exp
)
1948 (limit (cadr exp
) var val
'think
)
1949 (limit (caddr exp
) var val
'think
)))
1950 ((member (caar exp
) '(%sin %cos
) :test
#'eq
)
1951 (simplimsc exp
(caar exp
) (limit (cadr exp
) var val
'think
)))
1952 ((eq (caar exp
) '%tan
) (simplim%tan
(cadr exp
)))
1953 ((member (caar exp
) '(%sinh %cosh
) :test
#'eq
)
1954 (simplimsch (caar exp
) (limit (cadr exp
) var val
'think
)))
1955 ((member (caar exp
) '(%erf %tanh
) :test
#'eq
)
1956 (simplim%erf-%tanh
(caar exp
) (cadr exp
)))
1957 ((eq (caar exp
) '%atanh
)
1958 (simplim%atanh
(limit (cadr exp
) var val
'think
) val
))
1959 ((eq (caar exp
) '%acosh
)
1960 (simplim%acosh
(limit (cadr exp
) var val
'think
)))
1961 ((eq (caar exp
) '%asinh
)
1962 (simplim%asinh
(limit (cadr exp
) var val
'think
)))
1963 ((eq (caar exp
) '%inverse_jacobi_ns
)
1964 (simplim%inverse_jacobi_ns
(limit (cadr exp
) var val
'think
) (third exp
)))
1965 ((eq (caar exp
) '%inverse_jacobi_nc
)
1966 (simplim%inverse_jacobi_nc
(limit (cadr exp
) var val
'think
) (third exp
)))
1967 ((eq (caar exp
) '%inverse_jacobi_sc
)
1968 (simplim%inverse_jacobi_sc
(limit (cadr exp
) var val
'think
) (third exp
)))
1969 ((eq (caar exp
) '%inverse_jacobi_cs
)
1970 (simplim%inverse_jacobi_cs
(limit (cadr exp
) var val
'think
) (third exp
)))
1971 ((eq (caar exp
) '%inverse_jacobi_dc
)
1972 (simplim%inverse_jacobi_dc
(limit (cadr exp
) var val
'think
) (third exp
)))
1973 ((eq (caar exp
) '%inverse_jacobi_ds
)
1974 (simplim%inverse_jacobi_ds
(limit (cadr exp
) var val
'think
) (third exp
)))
1975 ((and (eq (caar exp
) 'mqapply
)
1976 (eq (subfunname exp
) '$psi
))
1977 (simplim$psi
(subfunsubs exp
) (subfunargs exp
) val
))
1978 ((and (eq (caar exp
) var
)
1979 (member 'array
(car exp
) :test
#'eq
)
1980 (every #'(lambda (sub-exp)
1983 exp
) ;LIMIT(B[I],B,INF); -> B[I]
1984 ;; When limsubst is true, limit(f(n+1)/f(n),n,inf) = 1. The user documentation
1985 ;; warns against setting limsubst to true.
1987 (simplify (cons (operator-with-array-flag exp
)
1988 (mapcar #'(lambda (a)
1989 (limit a var val
'think
))
1991 (throw 'limit t
)))))
1994 (setq e
(resimplify (subst (m// 1 var
) var e
)))
1995 (let ((new-val (cond ((eq val
'$zeroa
) '$inf
)
1996 ((eq val
'$zerob
) '$minf
))))
1997 (if new-val
(let ((preserve-direction t
))
1998 (toplevel-$limit e var new-val
)) (throw 'limit t
))))
2000 (defun simplimtimes (exp)
2001 ;; The following test
2002 ;; handles (-1)^x * 2^x => (-2)^x => $infinity
2003 ;; wants to avoid (-1)^x * 2^x => $ind * $inf => $und
2005 (and (expfactorp (cons '(mtimes) exp
) 1)
2006 (expfactor (cons '(mtimes) exp
) 1 var
))))
2007 (when try
(return-from simplimtimes try
)))
2009 (let ((prod 1) (num 1) (denom 1)
2010 (zf nil
) (ind-flag nil
) (inf-type nil
)
2011 (constant-zero nil
) (constant-infty nil
))
2013 (let* ((loginprod?
(involve term
'(%log
)))
2014 (y (catch 'lip?
(limit term var val
'think
))))
2016 ;; limit failed due to log in product
2018 (return-from simplimtimes
(liminv (cons '(mtimes simp
) exp
))))
2020 ;; If the limit is infinitesimal or zero
2022 (setf num
(m* num term
)
2023 constant-zero
(or constant-zero
(not (among var term
))))
2026 (unless zf
(setf zf
1)))
2028 (setf zf
(* -
1 (or zf
1))))))
2030 ;; If the limit is not some form of infinity or
2031 ;; undefined/indeterminate.
2032 ((not (member y
'($inf $minf $infinity $ind $und
) :test
#'eq
))
2033 (setq prod
(m* prod y
)))
2035 ((eq y
'$und
) (return-from simplimtimes
'$und
))
2036 ((eq y
'$ind
) (setq ind-flag t
))
2038 ;; Some form of infinity
2040 (setf denom
(m* denom term
)
2041 constant-infty
(or constant-infty
(not (among var term
))))
2042 (unless (eq inf-type
'$infinity
)
2044 ((eq y
'$infinity
) (setq inf-type
'$infinity
))
2045 ((null inf-type
) (setf inf-type y
))
2046 ;; minf * minf or inf * inf
2047 ((eq y inf-type
) (setf inf-type
'$inf
))
2049 (t (setf inf-type
'$minf
))))))))
2052 ;; If there are zeros and infinities among the terms that are free of
2053 ;; VAR, then we have an expression like "inf * zeroa * f(x)" or
2054 ;; similar. That gives an undefined result. Note that we don't
2055 ;; necessarily have something undefined if only the zeros have a term
2056 ;; free of VAR. For example "zeroa * exp(-1/x) * 1/x" as x -> 0. And
2057 ;; similarly for the infinities.
2058 ((and constant-zero constant-infty
) '$und
)
2060 ;; If num=denom=1, we didn't find any explicit infinities or zeros, so we
2061 ;; either return the simplified product or ind
2062 ((and (eql num
1) (eql denom
1))
2063 (if ind-flag
'$ind prod
))
2064 ;; If denom=1 (and so num != 1), we have some form of zero
2068 (let ((sign (getsignl prod
)))
2069 (if (or (not sign
) (eq sign
'complex
))
2075 ;; If num=1 (and so denom != 1), we have some form of infinity
2077 (let ((sign ($csign prod
)))
2080 ((eq sign
'$pos
) inf-type
)
2081 ((eq sign
'$neg
) (case inf-type
2085 ((member sign
'($complex $imaginary
)) '$infinity
)
2086 ; sign is '$zero, $pnz, $pz, etc
2087 (t (throw 'limit t
)))))
2088 ;; Both zeros and infinities
2090 ;; All bets off if there are some infinities or some zeros, but it
2091 ;; needn't be undefined (see above)
2092 (when (or constant-zero constant-infty
) (throw 'limit t
))
2094 (let ((ans (limit2 num
(m^ denom -
1) var val
)))
2096 (simplimtimes (list prod ans
))
2097 (throw 'limit t
)))))))
2099 ;;;PUT CODE HERE TO ELIMINATE FAKE SINGULARITIES??
2101 (defun simplimplus (exp)
2102 (cond ((memalike exp simplimplus-problems
)
2105 (progn (push exp simplimplus-problems
)
2106 (let ((ans (catch 'limit
(simplimplus1 exp
))))
2107 (cond ((or (eq ans
())
2109 (among '%limit ans
))
2110 (let ((new-exp (sratsimp exp
)))
2111 (cond ((not (alike1 exp new-exp
))
2113 (limit new-exp var val
'think
))))
2114 (cond ((or (eq ans
())
2119 (pop simplimplus-problems
)))))
2121 (defun simplimplus1 (exp)
2122 (prog (sum y infl infinityl minfl indl undl
)
2124 (do ((exp (cdr exp
) (cdr exp
)) (f))
2125 ((or y
(null exp
)) nil
)
2126 (setq f
(limit (car exp
) var val
'think
))
2129 ((eq f
'$und
) (push (car exp
) undl
))
2130 ((not (member f
'($inf $minf $infinity $ind
) :test
#'eq
))
2131 (setq sum
(m+ sum f
)))
2132 ((eq f
'$ind
) (push (car exp
) indl
))
2133 (infinityl (throw 'limit t
))
2134 ;;;Don't know what to do with an '$infinity and an $inf or $minf
2135 ((eq f
'$inf
) (push (car exp
) infl
))
2136 ((eq f
'$minf
) (push (car exp
) minfl
))
2138 (cond ((or infl minfl
)
2140 (t (push (car exp
) infinityl
))))))
2142 (if (or infl minfl indl infinityl
)
2143 (setq infinityl
(append undl infinityl
)); x^2 + x*sin(x)
2144 (return '$und
))) ; 1 + x*sin(x)
2145 ((not (or infl minfl indl infinityl
))
2146 (return (cond ((atom sum
) sum
)
2147 ((or (not (free sum
'$zeroa
))
2148 (not (free sum
'$zerob
)))
2151 (t (cond ((null infinityl
)
2152 (cond (infl (cond ((null minfl
) (return '$inf
))
2154 (minfl (return '$minf
))
2155 ((> (length indl
) 1)
2156 ;; At this point we have a sum of '$ind. We factor
2157 ;; the sum and try again. This way we get the limit
2158 ;; of expressions like (a-b)*ind, where (a-b)--> 0.
2159 (cond ((not (alike1 (setq y
($factorsum exp
)) exp
))
2160 (return (limit y var val
'think
)))
2163 (t (return '$ind
))))
2164 (t (setq infl
(append infl infinityl
))))))
2166 oon
(setq y
(m+l
(append minfl infl
)))
2167 (cond ((alike1 exp
(setq y
(sratsimp (hyperex y
))))
2168 (cond ((not (infinityp val
))
2169 (setq infl
(cnv infl val
)) ;THIS IS HORRIBLE!!!!
2170 (setq minfl
(cnv minfl val
))))
2172 (cond ((every #'(lambda (j) (radicalp j var
))
2173 (append infl minfl
))
2174 (setq y
(rheur infl minfl
)))
2175 (t (setq y
(sheur infl minfl
))))))
2176 (t (setq y
(limit y var val
'think
))))
2177 (cond ((or (eq y
())
2178 (eq y t
)) (return ()))
2179 ((infinityp y
) (return y
))
2180 (t (return (m+ sum y
))))))
2182 ;; Limit n/d, using heuristics on the order of growth.
2185 (cond ((and (free n var
)
2188 (t (setq n
(cpa n d nil
))
2190 (cond ((oscip orig-n
) '$und
)
2192 ((equal n -
1) '$zeroa
)
2193 ((equal n
0) (m// orig-n d
)))))))
2196 ;;;L1 is a list of INF's and L2 is a list of MINF's. Added together
2197 ;;;it is indeterminate.
2198 (defun sheur (l1 l2
)
2199 (let ((term (sheur1 l1 l2
)))
2200 (cond ((equal term
'$inf
) '$inf
)
2201 ((equal term
'$minf
) '$minf
)
2202 (t (let ((new-num (m+l
(mapcar #'(lambda (num-term)
2203 (m// num-term
(car l1
)))
2205 (cond ((limit2 new-num
(m// 1 (car l1
)) var val
))))))))
2208 (cond ((atom exp
) nil
)
2209 (t (setq exp
(nformat exp
))
2210 (cond ((and (eq (caar exp
) 'mquotient
)
2211 (among var
(caddr exp
)))
2214 (defun zerop2 (z) (=0 (ridofab z
)))
2216 (defun raise (a) (m+ a
'$zeroa
))
2218 (defun lower (a) (m+ a
'$zerob
))
2220 (defun sincoshk (exp1 l sc
)
2221 (cond ((equal l
1) (lower l
))
2222 ((equal l -
1) (raise l
))
2224 ((member val
'($zeroa $zerob
) :test
#'eq
) (spangside exp1 l
))
2227 (defun spangside (e l
)
2228 (setq e
(behavior e var val
))
2229 (cond ((equal e
1) (raise l
))
2230 ((equal e -
1) (lower l
))
2233 ;; get rid of zeroa and zerob
2235 (if (among '$zeroa e
) (setq e
(maxima-substitute 0 '$zeroa e
)))
2236 (if (among '$zerob e
) (setq e
(maxima-substitute 0 '$zerob e
)))
2240 ;; returns true if exp is a polynomial raised to a numeric power
2241 (defun simplerd (exp)
2243 (mnump (caddr exp
)) ;; exponent must be a number - no variables
2244 (polyp (cadr exp
))))
2246 (defun branch1 (exp val
)
2247 (cond ((polyp exp
) nil
)
2248 ((simplerd exp
) (zerop2 (subin val
(cadr exp
))))
2250 (loop for v on
(cdr exp
)
2251 when
(branch1 (car v
) val
)
2254 (defun branch (exp val
)
2255 (cond ((polyp exp
) nil
)
2256 ((or (simplerd exp
) (mtimesp exp
))
2259 (every #'(lambda (j) (branch j val
)) (the list
(cdr exp
))))))
2261 (defun ser0 (e n d val
)
2262 (cond ((and (branch n val
) (branch d val
))
2266 ;;;NN* gets set by POFX, called by SER1, to get a list of exponents.
2267 (setq nn
* (ratmin nn
*))
2268 (setq n
(sratsimp (m^ n
(m^ var nn
*))))
2269 (setq d
(sratsimp (m^ d
(m^ var nn
*))))
2270 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2273 (t (try-lhospital-quit n d nil
))))
2275 (defun rheur (l1 l2
)
2277 (setq m1
(mapcar (function asymredu
) l1
))
2278 (setq m2
(mapcar (function asymredu
) l2
))
2279 (setq ans
(m+l
(append m1 m2
)))
2280 (cond ((rptrouble (m+l
(append l1 l2
)))
2281 (return (limit (simplify (rdsget (m+l
(append l1 l2
))))
2285 ((mplusp ans
) (return (sheur m1 m2
)))
2286 (t (return (limit ans var val t
))))))
2288 (defun rptrouble (rp)
2289 (not (equal (rddeg rp nil
) (rddeg (asymredu rp
) nil
))))
2291 (defun radicalp (exp var
)
2292 (cond ((polyinx exp var
()))
2293 ((mexptp exp
) (cond ((equal (caddr exp
) -
1.
)
2294 (radicalp (cadr exp
) var
))
2296 ((member (caar exp
) '(mplus mtimes
) :test
#'eq
)
2297 (every #'(lambda (j) (radicalp j var
))
2300 (defun involve (e nn
*)
2301 (declare (special var
))
2302 (cond ((atom e
) nil
)
2304 ((and (member (caar e
) nn
* :test
#'eq
) (among var
(cdr e
))) (cadr e
))
2305 (t (some #'(lambda (j) (involve j nn
*)) (cdr e
)))))
2307 (defun notinvolve (exp nn
*)
2308 (cond ((atom exp
) t
)
2310 ((member (caar exp
) nn
* :test
#'eq
) (not (among var
(cdr exp
))))
2311 ((every #'(lambda (j) (notinvolve j nn
*))
2314 (defun sheur1 (l1 l2
)
2316 (setq l1
(mapcar #'stirling0 l1
))
2317 (setq l2
(mapcar #'stirling0 l2
))
2318 (setq l1
(m+l
(maxi l1
)))
2319 (setq l2
(m+l
(maxi l2
)))
2320 (setq ans
(cpa l1 l2 t
))
2321 (return (cond ((=0 ans
) (m+ l1 l2
))
2322 ((equal ans
1.
) '$inf
)
2325 (defun zero-lim (cpa-list)
2326 (do ((l cpa-list
(cdr l
)))
2328 (and (eq (caar l
) 'gen
)
2329 (zerop2 (limit (cadar l
) var val
'think
))
2332 ;; Compare order of growth for R1 and R2. The result is 0, -1, +1
2333 ;; depending on the relative order of growth. 0 is returned if R1 and
2334 ;; R2 have the same growth; -1 if R1 grows much more slowly than R2;
2335 ;; +1 if R1 grows much more quickly than R2.
2336 (defun cpa (r1 r2 flag
)
2339 (cond ((alike1 t1 t2
) 0.
)
2341 (cond ((free t2 var
) 0.
)
2342 (t (let ((lim-ans (limit1 t2 var val
)))
2343 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2346 (let ((lim-ans (limit1 t1 var val
)))
2347 (cond ((not (member lim-ans
'($inf $minf $und $ind
) :test
#'eq
)) 0.
)
2350 ;; Make T1 and T2 be a list of terms that are multiplied
2352 (cond ((mtimesp t1
) (setq t1
(cdr t1
)))
2353 (t (setq t1
(list t1
))))
2354 (cond ((mtimesp t2
) (setq t2
(cdr t2
)))
2355 (t (setq t2
(list t2
))))
2356 ;; Find the strengths of each term of T1 and T2
2357 (setq t1
(mapcar (function istrength
) t1
))
2358 (setq t2
(mapcar (function istrength
) t2
))
2359 ;; Compute the max of the strengths of the terms.
2360 (let ((ans (ismax t1
))
2362 (cond ((or (null ans
) (null d
))
2363 ;;(eq (car ans) 'gen) (eq (car d) 'gen))
2364 ;; ismax couldn't find highest term; give up
2367 (if (eq (car ans
) 'var
) (setq ans
(add-up-deg t1
)))
2368 (if (eq (car d
) 'var
) (setq d
(add-up-deg t2
)))
2369 ;; Can't just just compare dominating terms if there are
2370 ;; indeterm-inates present; e.g. X-X^2*LOG(1+1/X). So
2372 (cond ((or (zero-lim t1
)
2374 (cpa-indeterm ans d t1 t2 flag
))
2375 ((isgreaterp ans d
) 1.
)
2376 ((isgreaterp d ans
) -
1.
)
2379 (defun cpa-indeterm (ans d t1 t2 flag
)
2380 (cond ((not (eq (car ans
) 'var
))
2381 (setq ans
(gather ans t1
) d
(gather d t2
))))
2382 (let ((*indicator
(and (eq (car ans
) 'exp
)
2385 (setq test
(cpa1 ans d
))
2386 (cond ((and (zerop1 test
)
2387 (or (equal ($radcan
(m// (cadr ans
) (cadr d
))) 1.
)
2388 (and (polyp (cadr ans
))
2390 (equal (limit (m// (cadr ans
) (cadr d
)) var val
'think
)
2392 (let ((new-term1 (m// t1
(cadr ans
)))
2393 (new-term2 (m// t2
(cadr d
))))
2394 (cpa new-term1 new-term2 flag
)))
2397 (defun add-up-deg (strengthl)
2398 (do ((stl strengthl
(cdr stl
))
2401 ((null stl
) (list 'var
(m*l poxl
) (m+l degl
)))
2402 (cond ((eq (caar stl
) 'var
)
2403 (push (cadar stl
) poxl
)
2404 (push (caddar stl
) degl
)))))
2408 (cond ((eq (car p1
) 'gen
) (return 0.
)))
2409 (setq flag
(car p1
))
2414 (setq s1
(istrength p1
))
2415 (setq s2
(istrength p2
))
2418 ((isgreaterp s1 s2
) 1.
)
2419 ((isgreaterp s2 s1
) -
1.
)
2421 (setq *indicator nil
)
2423 ((and (poly? p1 var
) (poly? p2 var
))
2424 (setq p1
(m- p1 p2
))
2425 (cond ((zerop1 p1
) 0.
)
2426 (t (getsignl (hot-coef p1
)))))
2430 (list (m*t -
1 p2
))))
2431 (cond ((zerop2 s1
) 0.
)
2432 ((ratgreaterp s1
0.
) 1.
)
2436 (setq p1
(caddr p1
))
2437 (setq p2
(caddr p2
))
2438 (cond ((and (poly? p1 var
) (poly? p2 var
))
2439 (setq p1
(m- p1 p2
))
2440 (return (cond ((or (zerop1 p1
)
2441 (not (among var p1
)))
2443 (t (getsignl (hot-coef p1
))))))
2444 ((and (radicalp p1 var
) (radicalp p2 var
))
2447 (list (m*t -
1 p2
))))
2448 (return (cond ((eq s1
'$inf
) 1.
)
2449 ((eq s1
'$minf
) -
1.
)
2451 (cond ((ratgreaterp s1
0.
) 1.
)
2452 ((ratgreaterp 0. s1
) -
1.
)
2455 (t (return (cpa p1 p2 t
)))))
2457 (setq p1
(try-lhospital (asymredu p1
) (asymredu p2
) nil
))
2458 (return (cond ((zerop2 p1
) -
1.
)
2459 ((real-infinityp p1
) 1.
)
2462 ;;;EXPRESSIONS TO ISGREATERP ARE OF THE FOLLOWING FORMS
2463 ;;; ("VAR" POLY DEG)
2465 ;;; ("LOG" LOG(EXP))
2466 ;;; ("FACT" <A FACTORIAL EXPRESSION>)
2467 ;;; ("GEN" <ANY OTHER TYPE OF EXPRESSION>)
2469 (defun isgreaterp (a b
)
2472 (cond ((or (eq ta
'gen
)
2474 ((and (eq ta tb
) (eq ta
'var
))
2475 (ratgreaterp (caddr a
) (caddr b
)))
2476 ((and (eq ta tb
) (eq ta
'exp
))
2477 ;; Both are exponential order of infinity. Check the
2478 ;; exponents to determine which exponent is bigger.
2479 (eq (limit (m- `((%log
) ,(second a
)) `((%log
) ,(second b
)))
2482 ((member ta
(cdr (member tb
'(num log var exp fact gen
) :test
#'eq
)) :test
#'eq
)))))
2485 ;; Preprocess the list of products. Separate the terms that
2486 ;; exponentials and those that don't. Actually multiply the
2487 ;; exponential terms together to form a single term. Pass this and
2488 ;; the rest to ismax-core to find the max.
2489 (let (exp-terms non-exp-terms
)
2491 (if (eq 'exp
(car term
))
2492 (push term exp-terms
)
2493 (push term non-exp-terms
)))
2494 ;; Multiply the exp-terms together
2497 ;;(format t "exp-terms = ~A~%" exp-terms)
2498 (dolist (term exp-terms
)
2499 (setf product
(simplify (mul product
(second term
)))))
2500 ;;(format t "product = ~A~%" product)
2501 (setf product
`(exp ,($logcontract product
)))
2502 ;;(format t "product = ~A~%" product)
2503 (ismax-core (cons product non-exp-terms
)))
2506 (defun ismax-core (l)
2509 ((= (length l
) 1) (car l
)) ;If there is only 1 thing give it back.
2510 ((every #'(lambda (x)
2511 (not (eq (car x
) 'gen
))) l
)
2513 (do ((l1 (cdr l
) (cdr l1
))
2517 (cond ((isgreaterp temp-ans
(car l1
))
2518 (setq ans temp-ans
))
2519 ((isgreaterp (car l1
) temp-ans
)
2520 (setq temp-ans
(car l1
))
2521 (setq ans temp-ans
))
2522 (t (setq ans
())))))
2525 ;RETURNS LIST OF HIGH TERMS
2527 (cond ((atom all
) nil
)
2528 (t (do ((l (cdr all
) (cdr l
))
2530 (total 1) ; running total constant factor of hi-term
2531 (hi-terms (ncons (car all
)))
2533 ((null l
) (if (zerop2 total
) ; if high-order terms cancel each other
2534 all
; keep everything
2535 hi-terms
)) ; otherwise return list of high terms
2536 (setq compare
(limit (m// (car l
) hi-term
) var val
'think
))
2538 ((or (infinityp compare
)
2539 (and (eq compare
'$und
)
2540 (zerop2 (limit (m// hi-term
(car l
)) var val
'think
))))
2541 (setq total
1) ; have found new high term
2542 (setq hi-terms
(ncons (setq hi-term
(car l
)))))
2543 ((zerop2 compare
) nil
)
2544 ;; COMPARE IS IND, FINITE-VALUED, or und in both directions
2545 (t ; add to list of high terms
2546 (setq total
(m+ total compare
))
2547 (setq hi-terms
(append hi-terms
(ncons (car l
))))))))))
2551 (cond ((atom l
) (return nil
)))
2552 l1
(setq ans
(car l
))
2554 (cond ((null l
) (return ans
))
2555 ((ratgreaterp ans
(car l
)) (go l2
))
2560 (cond ((atom l
) (return nil
)))
2561 l1
(setq ans
(car l
))
2563 (cond ((null l
) (return ans
))
2564 ((ratgreaterp (car l
) ans
) (go l2
))
2572 ((or (mnump e
) (not (among var e
))) nil
)
2573 ((and (mexptp e
) (eq (cadr e
) var
))
2574 (push (caddr e
) nn
*))
2576 (t (mapc #'pofx
(cdr e
)))))
2579 (cond ((member val
'($zeroa $zerob
) :test
#'eq
) nil
)
2580 (t (setq e
(subin (m+ var val
) e
))))
2582 (cond ((pofx e
) e
)))
2584 (defun gather (ind l
)
2586 (setq ind
(car ind
))
2587 loop
(cond ((null l
)
2588 (return (list ind
(m*l ans
))))
2589 ((equal (caar l
) ind
)
2590 (push (cadar l
) ans
)))
2594 ; returns rough class-of-growth of term
2595 (defun istrength (term)
2596 (cond ((mnump term
) (list 'num term
))
2597 ((atom term
) (cond ((eq term var
)
2599 (t (list 'num term
))))
2600 ((not (among var term
)) (list 'num term
))
2602 (let ((temp (ismax-core (mapcar #'istrength
(cdr term
)))))
2603 (cond ((not (null temp
)) temp
)
2606 (let ((temp (mapcar #'istrength
(cdr term
)))
2608 (setq temp1
(ismax temp
))
2609 (cond ((null temp1
) `(gen ,term
))
2610 ((eq (car temp1
) 'log
) `(log ,temp
))
2611 ((eq (car temp1
) 'var
) (add-up-deg temp
))
2614 (real-infinityp (limit term var val t
)))
2615 (let ((logterm (logred term
)))
2616 (cond ((and (among var
(caddr term
))
2617 (member (car (istrength logterm
))
2618 '(var exp fact
) :test
#'eq
)
2619 (real-infinityp (limit logterm var val t
)))
2620 (list 'exp
(m^
'$%e logterm
)))
2621 ((not (among var
(caddr term
)))
2622 (let ((temp (istrength (cadr term
))))
2623 (cond ((not (alike1 temp term
))
2624 (rplaca (cdr temp
) term
)
2625 (and (eq (car temp
) 'var
)
2627 (m* (caddr temp
) (caddr term
))))
2631 ((and (eq (caar term
) '%log
)
2632 (real-infinityp (limit term var val t
)))
2633 (let ((stren (istrength (cadr term
))))
2634 (cond ((member (car stren
) '(log var
) :test
#'eq
)
2636 ((and (eq (car stren
) 'exp
)
2637 (eq (caar (second stren
)) 'mexpt
))
2638 (istrength (logred (second stren
))))
2640 ((eq (caar term
) 'mfactorial
)
2642 ((let ((temp (hyperex term
)))
2643 (and (not (alike1 term temp
))
2645 (t (list 'gen term
))))
2647 ;; log reduce - returns log of s1
2649 (or (and (eq (cadr s1
) '$%e
) (caddr s1
))
2650 (m* (caddr s1
) `((%log
) ,(cadr s1
)))))
2652 (defun asymredu (rd)
2653 (cond ((atom rd
) rd
)
2655 ((not (among var rd
)) rd
)
2656 ((polyinx rd var t
))
2658 (cond ((eq (cadr rd
) var
) rd
)
2660 (factor ($expand
(m^
(polyinx (cadr rd
) var t
)
2664 (t (simplify (cons (list (caar rd
))
2665 (mapcar #'asymredu
(cdr rd
)))))))
2668 (let ((dn** ()) (nn** ()))
2669 (cond ((atom rd
) rd
)
2671 ((not (among var rd
)) rd
)
2675 (cond ((eq (cadr rd
) var
) rd
)
2676 (t (setq dn
** (caddr rd
))
2677 (setq nn
** (factor (cadr rd
)))
2678 (cond ((mtimesp nn
**)
2679 (m*l
(mapcar #'(lambda (j) (m^ j dn
**))
2682 (t (simplify (cons (ncons (caar rd
))
2683 (mapcar #'rdfact
(cdr rd
))))))))
2685 (defun cnv (expl val
)
2686 (mapcar #'(lambda (e)
2687 (maxima-substitute (cond ((eq val
'$zerob
)
2688 (m* -
1 (m^ var -
1)))
2693 (t (m^
(m+ var
(m* -
1 val
)) -
1.
)))
2698 (defun pwtaylor (exp var l terms
)
2699 (prog (coef ans c mc
)
2700 (cond ((=0 terms
) (return nil
)) ((=0 l
) (setq mc t
)))
2703 loop
(setq c
(1+ c
))
2704 (cond ((or (> c
10.
) (equal c terms
))
2706 (t (setq exp
(sdiff exp var
))))
2707 tag1
(setq coef
($radcan
(subin l exp
)))
2708 (cond ((=0 coef
) (setq terms
(1+ terms
)) (go loop
)))
2715 (m^
`((mfactorial) ,c
) -
1)
2716 (m^
(if mc var
(m+t
(m*t -
1 l
) var
)) c
)))))
2721 ((simplerd e
) (rdtay e
))
2722 (t (cons (list (caar e
))
2723 (mapcar #'rdsget
(cdr e
))))))
2726 (cond (limit-using-taylor ($ratdisrep
($taylor rd var val
1.
)))
2730 (prog (varlist p c e d $ratfac
)
2731 (setq varlist
(ncons var
))
2732 (setq p
(ratnumerator (cdr (ratrep* (cadr rd
)))))
2733 (cond ((< (length p
) 3.
) (return rd
)))
2736 (setq c
(m^ var
(m* d e
)))
2737 (setq d
($ratsimp
(varinvert (m* (pdis p
) (m^ var
(m- d
)))
2739 (setq d
(pwtaylor (m^ d e
) var
0.
3.
))
2740 (return (m* c
(varinvert d var
)))))
2742 (defun varinvert (e var
) (subin (m^t var -
1.
) e
))
2745 (prog ((varlist (list var
)))
2746 (return (let (($ratfac nil
))
2748 (pdegr (cadr (ratrep* p
)))))))
2750 (defun rat-no-ratfac (e)
2751 (let (($ratfac nil
))
2756 (defun rddeg (rd low
*)
2757 (cond ((or (mnump rd
)
2758 (not (among var rd
)))
2763 (m* (deg (cadr rd
)) (caddr rd
)))
2765 (addn (mapcar #'(lambda (j)
2769 (setq rd
(andmapcar #'(lambda (j) (rddeg j low
*))
2771 (cond (low* (ratmin rd
))
2775 (cond ((or (atom pf
) (not (eq (caadr (ratf var
)) (car pf
))))
2777 (low* (cadr (reverse pf
)))
2779 ;;There is some confusion here. We need to be aware of Branch cuts etc....
2780 ;;when doing this section of code. It is not very carefully done so there
2781 ;;are bugs still lurking. Another misfortune is that LIMIT or its inferiors
2782 ;;sometimes decides to change the limit VAL in midstream. This must be corrected
2783 ;;since LIMIT's interaction with the data base environment must be maintained.
2784 ;;I'm not sure that this code can ever be called with VAL other than $INF but
2785 ;;there is a hook in the first important cond clause to cathc them anyway.
2788 (let ((num-power (rddeg n nil
))
2789 (den-power (rddeg d nil
))
2790 (coef ()) (coef-sign ()) (power ()))
2791 (setq coef
(m// ($ratcoef
($expand n
) var num-power
)
2792 ($ratcoef
($expand d
) var den-power
)))
2793 (setq coef-sign
(getsignl coef
))
2794 (setq power
(m// num-power den-power
))
2795 (cond ((eq (ask-integer power
'$integer
) '$integer
)
2796 (cond ((eq (ask-integer power
'$even
) '$even
) '$even
)
2797 (t '$odd
)))) ;Can be extended from here.
2798 (cond ((or (eq val
'$minf
)
2801 (equal val
0)) ()) ;Can be extended to cover some these.
2802 ((ratgreaterp den-power num-power
)
2803 (cond ((equal coef-sign
1.
) '$zeroa
)
2804 ((equal coef-sign -
1) '$zerob
)
2805 ((equal coef-sign
0) 0)
2807 ((ratgreaterp num-power den-power
)
2808 (cond ((equal coef-sign
1.
) '$inf
)
2809 ((equal coef-sign -
1) '$minf
)
2810 ((equal coef-sign
0) nil
) ; should never be zero
2811 ((null coef-sign
) '$infinity
)))
2814 (defun radlim (e n d
)
2816 (cond ((eq val
'$infinity
) (throw 'limit nil
))
2818 (setq nl
(m* var -
1))
2819 (setq n
(subin nl n
))
2820 (setq d
(subin nl d
))
2821 (setq val
'$inf
))) ;This is the Culprit. Doesn't tell the DATABASE.
2822 (cond ((eq val
'$inf
)
2823 (setq nl
(asymredu n
))
2824 (setq dl
(asymredu d
))
2826 ((or (rptrouble n
) (rptrouble d
))
2827 (return (limit (m* (rdsget n
) (m^
(rdsget d
) -
1.
)) var val t
)))
2828 (t (return (asy nl dl
))))))
2829 (setq nl
(limit n var val t
))
2830 (setq dl
(limit d var val t
))
2831 (cond ((and (zerop2 nl
) (zerop2 dl
))
2832 (cond ((or (polyp n
) (polyp d
))
2833 (return (try-lhospital-quit n d t
)))
2834 (t (return (ser0 e n d val
)))))
2835 (t (return ($radcan
(ratrad (m// n d
) n d nl dl
)))))))
2837 (defun ratrad (e n d nl dl
)
2839 (cond ((equal nl
0) (return 0))
2842 (cond ((equal dl
0) (setq d1
'$infinity
)) ;No direction Info.
2845 ((equal (setq d
(behavior d var val
)) 1)
2847 ((equal d -
1) (setq d1
'$minf
))
2848 (t (throw 'limit nil
))))
2851 (cond ((equal (setq n
(behavior n var val
)) 1)
2853 ((equal n -
1) (setq n1
'$zerob
))
2855 (t (return ($radcan
(ridofab (subin val e
))))))
2856 (return (simplimtimes (list n1 d1
)))))
2858 ;;; Limit(log(XXX), var, 0, val), where val is either zerob (limit from below)
2859 ;;; or zeroa (limit from above).
2860 (defun simplimln (expr var val
)
2861 (let ((arglim (limit (cadr expr
) var val
'think
)) (dir))
2862 (cond ((eq arglim
'$inf
) '$inf
) ;log(inf) = inf
2863 ;;log(minf,infinity,zerob) = infinity & log(0) = infinity
2864 ((or (member arglim
'($minf $infinity $zerob
)))
2866 ((eq arglim
'$zeroa
) '$minf
) ;log(zeroa) = minf
2867 ;; log(ind)=und & log(und)=und
2868 ((member arglim
'($ind $und
)) '$und
)
2869 ;; log(1^(-)) = zerob, log(1^(+)) = zeroa & log(1)=0
2871 (if (or (eq val
'$zerob
) (eq var
'$zeroa
)) val
0))
2872 ;; Special case of arglim = 0
2874 (setq dir
(behavior (cadr expr
) var val
))
2875 (cond ((eql dir -
1) '$infinity
)
2876 ((eql dir
0) '$infinity
)
2877 ((eql dir
1) '$minf
)))
2878 ;; When arglim is off the negative real axis, use direct substitution
2879 ((off-negative-real-axisp arglim
)
2880 (ftake '%log arglim
))
2882 ;; We know that arglim is a negative real number, say xx.
2883 ;; When the imaginary part of (cadr expr) near var is negative,
2884 ;; return log(-x) - %i*pi; when the imaginary part of (cadr expr)
2885 ;; near var is positive return log(-x) + %i*pi; and when
2886 ;; we cannot determine the behavior of the imaginary part,
2887 ;; return a nounform. The value of val (either zeroa or zerob)
2888 ;; determines what is meant by "near" (smaller than var when
2889 ;; val is zerob and larger than var when val is zeroa).
2890 (setq dir
(behavior ($imagpart
(cadr expr
)) var val
))
2891 (cond ((or (eql dir
1) (eql dir -
1))
2892 (add (ftake '%log
(mul -
1 arglim
)) (mul dir
'$%i
'$%pi
)))
2893 (t (throw 'limit nil
))))))) ;do a nounform return
2894 (setf (get '%log
'simplim%function
) 'simplimln
)
2895 (setf (get '%plog
'simplim%function
) 'simplimln
)
2897 (defun simplim%limit
(e x pt
)
2898 (declare (ignore e x pt
))
2900 (setf (get '%limit
'simplim%function
) 'simplim%limit
)
2902 (defun simplim%unit_step
(e var val
)
2903 (let ((lim (limit (cadr e
) var val
'think
)))
2904 (cond ((eq lim
'$und
) '$und
)
2905 ((eq lim
'$ind
) '$ind
)
2906 ((eq lim
'$zerob
) 0)
2907 ((eq lim
'$zeroa
) 1)
2908 ((not (lenient-realp lim
)) (throw 'limit nil
)) ;catches infinity too
2909 ;; catches minf & inf cases
2910 ((eq t
(mgrp 0 lim
)) 0)
2911 ((eq t
(mgrp lim
0)) 1)
2913 (setf (get '$unit_step
'simplim%function
) 'simplim%unit_step
)
2915 (defun simplim%conjugate
(e var val
)
2916 (let ((lim (limit (cadr e
) var val
'think
)))
2917 (cond ((off-negative-real-axisp lim
)
2918 (ftake '$conjugate lim
))
2919 (t (throw 'limit nil
)))))
2920 (setf (get '$conjugate
'simplim%function
) 'simplim%conjugate
)
2922 (defun simplim%imagpart
(e var val
)
2923 (let ((lim (limit (cadr e
) var val
'think
)))
2924 (cond ((eq lim
'$und
) '$und
)
2926 ((eq lim
'$infinity
) (throw 'limit nil
))
2927 (t (mfuncall '$imagpart lim
)))))
2928 (setf (get '$imagpart
'simplim%function
) 'simplim%imagpart
)
2929 (setf (get '%imagpart
'simplim%function
) 'simplim%imagpart
)
2931 (defun simplim%realpart
(e var val
)
2932 (let ((lim (limit (cadr e
) var val
'think
)))
2933 (cond ((eq lim
'$und
) '$und
)
2934 ((eq lim
'$ind
) '$ind
)
2935 ((eq lim
'$infinity
) (throw 'limit nil
))
2936 (t (mfuncall '$realpart lim
)))))
2937 (setf (get '$realpart
'simplim%function
) 'simplim%realpart
)
2938 (setf (get '%realpart
'simplim%function
) 'simplim%realpart
)
2939 ;;; Limit of the Factorial function
2941 (defun simplimfact (expr var val
)
2942 (let* ((arglim (limit (cadr expr
) var val
'think
)) ; Limit of the argument.
2944 (cond ((eq arglim
'$inf
) '$inf
)
2945 ((member arglim
'($minf $infinity $und $ind
) :test
#'eq
) '$und
)
2947 ((and (or (maxima-integerp arglim
)
2948 (setq arg2
(integer-representation-p arglim
)))
2949 (eq ($sign arg2
) '$neg
))
2950 ;; A negative integer or float or bigfloat representation.
2951 (let ((dir (limit (add (cadr expr
) (mul -
1 arg2
)) var val
'think
))
2952 (even (mevenp arg2
)))
2953 (cond ((or (and even
2963 (t (throw 'limit nil
)))))
2965 ;; Call simplifier to get value at the limit of the argument.
2966 (simplify (list '(mfactorial) arglim
))))))
2967 (setf (get 'mfactorial
'simplim%function
) 'simplimfact
)
2969 (defun simplim%erf-%tanh
(fn arg
)
2970 (let ((arglim (limit arg var val
'think
))
2973 (cond ((eq arglim
'$inf
) 1)
2974 ((eq arglim
'$minf
) -
1)
2975 ((eq arglim
'$infinity
)
2976 (destructuring-bind (rpart . ipart
)
2978 (setq rlim
(limit rpart var origval
'think
))
2979 (cond ((eq fn
'%tanh
)
2980 (cond ((equal rlim
'$inf
) 1)
2981 ((equal rlim
'$minf
) -
1)))
2983 (setq ans
(limit (m* rpart
(m^t ipart -
1)) var origval
'think
))
2984 (setq ans
($asksign
(m+ `((mabs) ,ans
) -
1)))
2985 (cond ((or (eq ans
'$pos
) (eq ans
'$zero
))
2986 (cond ((eq rlim
'$inf
) 1)
2987 ((eq rlim
'$minf
) -
1)
2990 ((eq arglim
'$und
) '$und
)
2991 ((member arglim
'($zeroa $zerob $ind
) :test
#'eq
) arglim
)
2992 ;;;Ignore tanh(%pi/2*%I) and multiples of the argument.
2994 ;; erf (or tanh) of a known value is just erf(arglim).
2995 (simplify (list (ncons fn
) arglim
))))))
2997 (defun in-domain-of-atan (z)
2998 (setq z
(trisplit z
)) ; split z into real and imaginary parts
2999 (let ((x (car z
)) (y (cdr z
))) ;z = x+%i*y
3002 (eq t
(meqp x
0)) ;Re(z) = 0
3003 (or (eq t
(mgqp -
1 y
)) ;-1 >= Im(z)
3004 (eq t
(mgqp y
1))))))) ; Im(z) >= 1
3006 (defun simplim%atan
(e x pt
)
3007 (let ((lim (limit (cadr e
) x pt
'think
)))
3008 (cond ((or (eq lim
'$zeroa
) (eq lim
'$zerob
) (eq lim
0) (eq lim
'$ind
)) lim
)
3009 ((or (eq lim
'$und
) (eq lim
'$infinity
)) (throw 'limit nil
))
3010 ((eq lim
'$inf
) #$%pi
/2$
) ;atan(inf) -> %pi/2
3011 ((eq lim
'$minf
) #$-%pi
/2$
) ;atan(-inf) -> -%pi/2
3012 ((in-domain-of-atan (ridofab lim
)) ; direct substitution
3013 (ftake '%atan
(ridofab lim
)))
3014 (t (limit ($logarc e
) x pt
'think
)))))
3015 (setf (get '%atan
'simplim%function
) 'simplim%atan
)
3017 (defmvar extended-reals
3018 (append *infinitesimals
* *infinities
* (list '$und
'$ind
)))
3020 ;; Most instances of atan2 are simplified to atan expressions, but this routine
3021 ;; handles tricky cases such as limit(atan2((x^2-2), x^3-2*x), x, sqrt(2), minus).
3022 ;; Taylor and Gruntz cannot handle the discontinuity at atan(0, -1)
3024 ;; When possible, we want to evaluate the limit of an atan2 expression using
3025 ;; direct substitution--that produces, I think, the least surprising values.
3027 ;; The general simplifier catches atan2(0,0) and it transforms atan2(minf or inf,X)
3028 ;; and atan2(X, minf or inf) into an atan expression, but it doesn't catch
3029 ;; atan2(X, zerob or zeroa) or atan2(zerob or zeroa, X). For the other extended
3030 ;; real (ind,und, or infinity) arguments, the general simplifier gives sign errors.
3032 (defun simplim%atan2
(e v pt
)
3033 (let ((y (second e
)) (x (third e
)) (xlim) (ylim) (xlim-z) (ylim-z) (q))
3034 (setq xlim
(limit x v pt
'think
))
3035 (setq ylim
(limit y v pt
'think
))
3036 (setq xlim-z
(ridofab xlim
)
3037 ylim-z
(ridofab ylim
))
3038 ;; For cases for which direct substitution fails, normalize
3039 ;; x & y and try again.
3040 (setq q
(cond ((eq xlim
'$inf
) x
)
3041 ((eq xlim
'$minf
) (mul -
1 x
))
3043 ((eq ylim
'$minf
) (mul -
1 y
))
3044 ((and (eq xlim
'$zerob
) (zerop2 ylim
)) (mul -
1 x
))
3045 ((and (eq xlim
'$zeroa
) (zerop2 ylim
)) x
)
3046 ((and (eq ylim
'$zerob
) (zerop2 xlim
)) (mul -
1 y
))
3047 ((and (eq ylim
'$zeroa
) (zerop2 xlim
)) y
)
3050 (when (not (eql q
1))
3053 (setq xlim
(limit x v pt
'think
))
3054 (setq ylim
(limit y v pt
'think
))
3055 (setq xlim-z
(ridofab xlim
)
3056 ylim-z
(ridofab ylim
)))
3059 ((and (eq '$zerob ylim
) (eq t
(mgrp 0 xlim
))) (mul -
1 '$%pi
))
3060 ((and (eq '$zerob ylim
) (eq t
(mgrp xlim
0))) 0)
3061 ((and (eq '$zeroa ylim
) (eq t
(mgrp 0 xlim
))) '$%pi
)
3062 ((and (eq '$zeroa ylim
) (eq t
(mgrp xlim
0))) 0)
3063 ((and (eql xlim
1) (eql ylim
'$inf
)) (div '$%pi
2))
3064 ((and (eql xlim -
1) (eql ylim
0)) '$ind
)
3066 ;; Use direct substitution when ylim-z # 0 or xlim-z > 0. We need
3067 ;; to check that xlim-z & ylim-z are real too.
3068 ((and (lenient-realp xlim-z
) (lenient-realp ylim-z
)
3069 (or (eq t
(mnqp ylim-z
0)) (eq t
(mgrp xlim-z
0))))
3070 (ftake '$atan2 ylim-z xlim-z
))
3072 (throw 'limit nil
)))))
3073 (setf (get '$atan2
'simplim%function
) 'simplim%atan2
)
3075 (defun simplimsch (sch arg
)
3076 (cond ((real-infinityp arg
)
3077 (cond ((eq sch
'%sinh
) arg
) (t '$inf
)))
3078 ((eq arg
'$infinity
) '$infinity
)
3079 ((eq arg
'$ind
) '$ind
)
3080 ((eq arg
'$und
) '$und
)
3081 ((and (eq sch
'%sinh
)
3082 (or (eq arg
'$zerob
) (eq arg
'$zeroa
)))
3084 (t (let (($exponentialize t
))
3085 (resimplify (list (ncons sch
) (ridofab arg
)))))))
3087 ;; simple limit of sin and cos
3088 (defun simplimsc (exp fn arg
)
3089 (cond ((member arg
'($inf $minf $ind
) :test
#'eq
) '$ind
)
3090 ((member arg
'($und $infinity
) :test
#'eq
)
3092 ((member arg
'($zeroa $zerob
) :test
#'eq
)
3093 (cond ((eq fn
'%sin
) arg
)
3094 (t (m+ 1 '$zerob
))))
3096 (simplify (list (ncons fn
) (ridofab arg
)))
3099 (defun simplim%tan
(arg)
3100 (let ((arglim (limit arg var val
'think
)))
3102 ((member arglim
'($inf $minf $ind
) :test
#'eq
)
3104 ((member arglim
'($und $infinity
) :test
#'eq
)
3107 ;; Write the limit of the argument as c*%pi + rest.
3109 ((c (or (pip arglim
) 0))
3110 (rest (sratsimp (m- arglim
(m* '$%pi c
))))
3112 ;; Check if tan(x) has a zero or pole at x=arglim.
3113 ;; zero: tan(n*%pi + 0*)
3114 ;; pole: tan((2*n+1)*%pi/2 + 0*)
3115 ;; 0* can be $zeroa, $zerob or 0.
3116 (if (and (member rest
'(0 $zeroa $zerob
) :test
#'equal
)
3117 (or (setq hit-zero
(integerp c
))
3118 (and (ratnump c
) (equal (caddr c
) 2))))
3119 ;; This is a zero or a pole.
3120 ;; Determine on which side of the zero/pole we are.
3121 ;; If rest is $zeroa or $zerob, use that.
3122 ;; Otherwise (rest = 0), try to determine the side
3123 ;; using the behavior of the argument.
3125 ((side (cond ((eq rest
'$zeroa
) 1)
3126 ((eq rest
'$zerob
) -
1)
3127 (t (behavior arg var val
)))))
3129 ;; For a zero, if we don't know the side, just return 0.
3131 ((equal side
1) '$zeroa
)
3132 ((equal side -
1) '$zerob
)
3134 ;; For a pole, we need to know the side.
3135 ;; Otherwise, we can't determine the limit.
3137 ((equal side
1) '$minf
)
3138 ((equal side -
1) '$inf
)
3139 (t (throw 'limit t
)))))
3140 ;; No zero or pole - substitute in the limit of the argument.
3141 (take '(%tan
) (ridofab arglim
))))))))
3143 (defun simplim%asinh
(arg)
3144 (cond ((member arg
'($inf $minf $zeroa $zerob $ind $und
) :test
#'eq
)
3146 ((eq arg
'$infinity
) '$und
)
3147 (t (simplify (list '(%asinh
) (ridofab arg
))))))
3149 (defun simplim%acosh
(arg)
3150 (cond ((equal (ridofab arg
) 1) '$zeroa
)
3151 ((eq arg
'$inf
) arg
)
3152 ((eq arg
'$minf
) '$infinity
)
3153 ((member arg
'($und $ind $infinity
) :test
#'eq
) '$und
)
3154 (t (simplify (list '(%acosh
) (ridofab arg
))))))
3156 (defun simplim%atanh
(arg dir
)
3157 ;; Compute limit(atanh(x),x,arg). If ARG is +/-1, we need to take
3158 ;; into account which direction we're approaching ARG.
3159 (cond ((zerop2 arg
) arg
)
3160 ((member arg
'($ind $und $infinity $minf $inf
) :test
#'eq
)
3162 ((equal (setq arg
(ridofab arg
)) 1.
)
3163 ;; The limit at 1 should be complex infinity because atanh(x)
3164 ;; is complex for x > 1, but inf if we're approaching 1 from
3166 (if (eq dir
'$zerob
)
3170 ;; Same as above, except for the limit is at -1.
3171 (if (eq dir
'$zeroa
)
3174 (t (simplify (list '(%atanh
) arg
)))))
3176 (defun simplim%asin
(e x pt
)
3177 (let ((lim (limit (cadr e
) x pt
'think
)) (dir) (lim-sgn))
3178 (cond ((member lim
'($zeroa $zerob
)) lim
) ;asin(zeoroa/b) = zeroa/b
3179 ((member lim
'($minf
'$inf
'$infinity
)) '$infinity
)
3180 ((eq lim
'$ind
) '$ind
) ;asin(ind)=ind
3181 ((eq lim
'$und
) '$und
) ;asin(und)=und
3182 ((in-domain-of-asin lim
) ;direct substitution
3185 (setq e
(trisplit (cadr e
))) ;overwrite e!
3186 (setq dir
(behavior (cdr e
) x pt
))
3187 (setq lim-sgn
($csign
(car e
))) ;lim-sgn = sign limit(Re(e))
3190 (throw 'limit t
)) ;unable to find behavior of imaginary part
3192 ;; For the values of asin on the branch cuts, see DLMF 4.23.20 & 4.23.21
3193 ;; Diagram of the values of asin just above and below the branch cuts
3195 ;; asin(x) pi - asin(x)
3196 ;;................ -1 ....0.... 1 ...............
3197 ;; -pi - asin(x) asin(x)
3199 ;; Let's start in northwest and rotate counterclockwise:
3200 ((and (eq '$neg lim-sgn
) (eq dir
1))
3202 ((and (eq '$pos lim-sgn
) (eq dir
1))
3203 (sub '$%pi
(ftake '%asin lim
)))
3204 ((and (eq '$pos lim-sgn
) (eq dir -
1))
3206 ((and (eq '$neg lim-sgn
) (eq dir -
1))
3207 (sub (mul -
1 '$%pi
) (ftake '%asin lim
)))
3208 (t (throw 'limit t
))))))) ; unable to find sign of real part of lim.
3209 (setf (get '%asin
'simplim%function
) 'simplim%asin
)
3211 (defun simplim%acos
(e x pt
)
3212 (let ((lim (limit (cadr e
) x pt
'think
)) (dir) (lim-sgn))
3213 (cond ((in-domain-of-asin lim
) ;direct substitution
3215 ((member lim
'($und $ind $inf $minf $infinity
)) ;boundary cases
3218 (setq e
(trisplit (cadr e
))) ;overwrite e!
3219 (setq dir
(behavior (cdr e
) x pt
))
3220 (setq lim-sgn
($csign lim
))
3223 (throw 'limit t
)) ;unable to find behavior of imaginary part
3224 ;; for the values of acos on the branch cuts, see DLMF 4.23.24 & 4.23.25
3225 ;; http://dlmf.nist.gov/4.23.E24
3226 ((or (eq '$pos lim-sgn
) (eq '$neg lim-sgn
))
3227 ;; continuous from above
3228 (if (eql dir
1) (ftake '%acos lim
) (sub (mul 2 '$%pi
) (ftake '%acos lim
))))
3229 (t (throw 'limit t
))))))) ; unable to find sign of real part of lim.
3230 (setf (get '%acos
'simplim%function
) 'simplim%acos
)
3232 ;; Limit of an %integrate expression. For a definite integral
3233 ;; integrate(ee,var,a,b), when ee is free of the limit variable
3234 (defun simplim%integrate
(e x pt
)
3235 (let* ((ee (second e
)) ;ee = integrand
3236 (var (third e
)) ;integration variable
3237 (a (fourth e
)) ;lower limit or nil if indefinite
3238 (b (fifth e
)) ;lower limit or nil if indefinite
3240 (cond ((and a b
($freeof x ee
) ($freeof x var
))
3241 (setq alim
(limit a x pt
'think
))
3242 (setq blim
(limit b x pt
'think
))
3243 (if (and (lenient-extended-realp alim
)
3244 (lenient-extended-realp blim
)
3245 (not (eq alim
'$infinity
))
3246 (not (eq blim
'$infinity
)))
3247 (ftake '%integrate ee var alim blim
)
3250 (throw 'limit t
)))))
3251 (setf (get '%integrate
'simplim%function
) 'simplim%integrate
)
3253 (defun subftake (op subarg arg
)
3254 (simplifya (subfunmake op subarg arg
) t
))
3256 (defun off-one-to-inf (z)
3257 (setq z
(trisplit z
)) ; split z into x+%i*y
3259 (eq t
(mnqp (cdr z
) 0)) ; y # 0
3260 (eq t
(mgrp 1 (car z
))))) ; x < 1
3262 (defun simplim%li
(expr x pt
)
3263 (let ((n (car (subfunsubs expr
))) (e (car (subfunargs expr
))))
3265 (setq e
(limit e x pt
'think
))
3266 (cond ((and (eq e
'$minf
) (integerp n
) (>= n
2))
3268 ((and (eq e
'$inf
) (integerp n
) (>= n
2))
3270 ((or (eql (ridofab e
) 1) (and (not (extended-real-p e
)) (off-one-to-inf e
)))
3271 ;; Limit of li[s](1) can be evaluated by just
3272 ;; substituting in 1.
3273 ;; Same for li[s](x) when x is < 1.
3274 (subftake '$li
(list n
) (list e
)))
3275 (t (throw 'limit nil
))))
3276 ;; Claim ignorance when order depends on limit variable.
3277 (t (throw 'limit nil
)))))
3279 (setf (get '$li
'simplim%function
) 'simplim%li
)
3280 (setf (get '%li
'simplim%function
) 'simplim%li
)
3282 (defun simplim$psi
(order arg val
)
3283 (if (and (not (equal (length order
) 1))
3284 (not (equal (length arg
) 1)))
3286 (setq order
(car order
)
3288 (cond ((equal order
0)
3289 (destructuring-bind (rpart . ipart
)
3291 (cond ((not (equal ipart
0)) (throw 'limit
()))
3292 (t (setq rpart
(limit rpart var val
'think
))
3293 (cond ((eq rpart
'$zeroa
) '$minf
)
3294 ((eq rpart
'$zerob
) '$inf
)
3295 ((eq rpart
'$inf
) '$inf
)
3296 ((eq rpart
'$minf
) '$und
)
3297 ((equal (getsignl rpart
) -
1) (throw 'limit
()))
3298 (t (simplify (subfunmake '$psi
(list order
)
3299 (list rpart
)))))))))
3300 ((and (integerp order
) (> order
0)
3301 (equal (limit arg var val
'think
) '$inf
))
3302 (cond ((mevenp order
) '$zerob
)
3303 ((moddp order
) '$zeroa
)
3304 (t (throw 'limit
()))))
3305 (t (throw 'limit
()))))
3307 (defun simplim%inverse_jacobi_ns
(arg m
)
3308 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3310 `((%inverse_jacobi_ns
) ,arg
,m
)))
3312 (defun simplim%inverse_jacobi_nc
(arg m
)
3313 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3314 `((%elliptic_kc
) ,m
)
3315 `((%inverse_jacobi_nc
) ,arg
,m
)))
3317 (defun simplim%inverse_jacobi_sc
(arg m
)
3318 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3319 `((%elliptic_kc
) ,m
)
3320 `((%inverse_jacobi_sc
) ,arg
,m
)))
3322 (defun simplim%inverse_jacobi_dc
(arg m
)
3323 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3324 `((%elliptic_kc
) ,m
)
3325 `((%inverse_jacobi_dc
) ,arg
,m
)))
3327 (defun simplim%inverse_jacobi_cs
(arg m
)
3328 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3330 `((%inverse_jacobi_cs
) ,arg
,m
)))
3332 (defun simplim%inverse_jacobi_ds
(arg m
)
3333 (if (or (eq arg
'$inf
) (eq arg
'$minf
))
3335 `((%inverse_jacobi_ds
) ,arg
,m
)))
3337 (defun simplim%signum
(e x pt
)
3338 (let ((e (limit (cadr e
) x pt
'think
)) (sgn))
3339 (cond ((eq '$minf e
) -
1)
3341 ((eq '$infinity e
) '$und
)
3342 ((eq '$ind e
) '$ind
)
3347 (setq sgn
(mnqp e
0))
3348 (cond ((eq t sgn
) (ftake '%signum e
))
3349 (t (throw 'limit nil
))))))) ; don't know
3350 (setf (get '%signum
'simplim%function
) 'simplim%signum
)
3352 ;; more functions for limit to handle
3354 (defun lfibtophi (e)
3355 (cond ((not (involve e
'($fib
))) e
)
3356 ((eq (caar e
) '$fib
)
3357 (let ((lnorecurse t
))
3358 ($fibtophi
(list '($fib
) (lfibtophi (cadr e
))) lnorecurse
)))
3360 (mapcar #'lfibtophi
(cdr e
))))))
3362 ;;; FOLLOWING CODE MAKES $LDEFINT WORK
3364 (defmfun $ldefint
(exp var ll ul
&aux $logabs ans a1 a2
)
3365 (setq $logabs t ans
(sinint exp var
)
3366 a1
(toplevel-$limit ans var ul
'$minus
)
3367 a2
(toplevel-$limit ans var ll
'$plus
))
3368 (and (member a1
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3369 (setq a1
(nounlimit ans var ul
)))
3370 (and (member a2
'($inf $minf $infinity $und $ind
) :test
#'eq
)
3371 (setq a2
(nounlimit ans var ll
)))
3372 ($expand
(m- a1 a2
)))
3374 (defun nounlimit (exp var val
)
3375 (setq exp
(restorelim exp
))
3376 (nconc (list '(%limit
) exp var
(ridofab val
))
3377 (cond ((eq val
'$zeroa
) '($plus
))
3378 ((eq val
'$zerob
) '($minus
)))))
3380 ;; substitute inside noun form of %derivative
3381 ;; for cases such as limit('diff(x+2,x), x, 1)
3382 ;; -> limit('diff(xx+3), xx, 0)
3384 ;; maxima-substitute with *atp* skips over %derivative
3386 ;; substitutes diff(f(realvar), realvar, n)
3387 ;; -> diff(f(var+val), var, n)
3388 (defun derivative-subst (exp val var realvar
)
3389 (cond ((atom exp
) exp
)
3390 ((eq '%derivative
(caar exp
))
3393 (cons ;; the function being differentiated
3394 (maxima-substitute (m+ val var
) realvar
(cadr exp
))
3395 (cons ;; the var of differentiation
3396 (maxima-substitute var realvar
(caddr exp
))
3397 (cdddr exp
))))) ;; the order of the derivative
3399 (mapcar (lambda (x) (derivative-subst x val var realvar
))
3402 ;;;Used by Defint also.
3404 (or (involve e
'(%sin %cos %tan
))
3405 (among '$%i
(%einvolve e
))))
3407 (defun %einvolve
(e)
3408 (when (among '$%e e
) (%einvolve01 e
)))
3410 (defun %einvolve01
(e)
3411 (cond ((atom e
) nil
)
3415 (among var
(caddr e
)))
3417 (t (some #'%einvolve
(cdr e
)))))
3419 (declare-top (unspecial *indicator exp var val origval taylored
3420 $tlimswitch logcombed lhp? lhcount
))
3426 ;; "On Computing Limits in a Symbolic Manipulation System"
3427 ;; PhD Dissertation ETH Zurich 1996
3429 ;; The algorithm identifies the most rapidly varying (MRV) subexpression,
3430 ;; replaces it with a new variable w, rewrites the expression in terms
3431 ;; of the new variable, and then repeats.
3433 ;; The algorithm doesn't handle oscillating functions, so it can't do things like
3434 ;; limit(sin(x)/x, x, inf).
3436 ;; To handle limits involving functions like gamma(x) and erf(x), the
3437 ;; gruntz algorithm requires them to be written in terms of asymptotic
3438 ;; expansions, which maxima cannot currently do.
3440 ;; The algorithm assumes that everything is real, so it can't
3441 ;; currently handle limit((-2)^x, x, inf).
3443 ;; This is one of the methods used by maxima's $limit.
3444 ;; It is also directly available to the user as $gruntz.
3447 ;; most rapidly varying subexpression of expression exp with respect to limit variable var.
3448 ;; returns a list of subexpressions which are in the same MRV equivalence class.
3449 (defun mrv (exp var
)
3450 (cond ((freeof var exp
)
3455 (mrv-max (mrv (cadr exp
) var
)
3456 (mrv (m*l
(cddr exp
)) var
)
3459 (mrv-max (mrv (cadr exp
) var
)
3460 (mrv (m+l
(cddr exp
)) var
)
3463 (cond ((freeof var
(caddr exp
))
3464 (mrv (cadr exp
) var
))
3465 ((member (limitinf (logred exp
) var
) '($inf $minf
) :test
#'eq
)
3466 (mrv-max (list exp
) (mrv (caddr exp
) var
) var
))
3467 (t (mrv-max (mrv (cadr exp
) var
) (mrv (caddr exp
) var
) var
))))
3469 (mrv (cadr exp
) var
))
3470 ((equal (length (cdr exp
)) 1)
3471 (mrv (cadr exp
) var
))
3472 ((equal (length (cdr exp
)) 2)
3473 (mrv-max (mrv (cadr exp
) var
)
3474 (mrv (caddr exp
) var
)
3476 (t (tay-error "mrv not implemented" exp
))))
3478 ;; takes two lists of expressions, f and g, and limit variable var.
3479 ;; members in each list are assumed to be in same MRV equivalence
3480 ;; class. returns MRV set of the union of the inputs - either f or g
3481 ;; or the union of f and g.
3482 (defun mrv-max (f g var
)
3489 (return (union f g
))))
3490 (let ((c (mrv-compare (car f
) (car g
) var
)))
3496 (return (union f g
)))
3497 (t (merror "MRV-MAX: expected '>' '<' or '='; found: ~M" c
))))))
3499 (defun mrv-compare (a b var
)
3500 (let ((c (limitinf (m// `((%log
) ,a
) `((%log
) ,b
)) var
)))
3503 ((member c
'($inf $minf
) :test
#'eq
)
3507 ;; rewrite expression exp by replacing members of MRV set omega with
3508 ;; expressions in terms of new variable wsym. return cons pair of new
3509 ;; version of exp and the log of the new variable wsym.
3510 (defun mrv-rewrite (exp omega var wsym
)
3511 (setq omega
(stable-sort omega
(lambda (x y
) (> (length (mrv x var
))
3512 (length (mrv y var
))))));FIXME consider a total order function with #'sort
3513 (let* ((g (car (last omega
)))
3515 (sig (equal (mrv-sign logg var
) 1))
3516 (w (if sig
(m// 1 wsym
) wsym
))
3517 (logw (if sig
(m* -
1 logg
) logg
)))
3518 (mapcar (lambda (x y
)
3519 ;;(mtell "y:~M x:~M exp:~M~%" y x exp)
3520 (setq exp
(syntactic-substitute y x exp
)))
3522 (mapcar (lambda (f) ;; rewrite each element of omega
3523 (let* ((logf (logred f
))
3524 (c (mrv-leadterm (m// logf logg
) var nil
)))
3525 (cond ((not (equal (cadr c
) 0))
3526 (merror "MRV-REWRITE: expected leading term to be constant in ~M" c
)))
3527 ;;(mtell "logg: ~M logf: ~M~%" logg logf)
3530 (m* (car c
) logg
))))))
3534 ;;; if log w(x) = h(x), rewrite all subexpressions of the form
3535 ;;; log(f(x)) as log(w^-c f(x)) + c h(x) with c the unique constant
3536 ;;; such that w^-c f(x) is strictly less rapidly varying than w.
3537 (defun mrv-rewrite-logs (exp wsym logw
)
3538 (cond ((atom exp
) exp
)
3540 (not (freeof wsym exp
)))
3541 (let* ((f (cadr exp
))
3542 (c ($lopow
(calculate-series f wsym
)
3545 (m* (m^ wsym
(m- c
))
3546 (mrv-rewrite-logs f wsym logw
)))
3551 (mrv-rewrite-logs e wsym logw
))
3554 ;; returns list of two elements: coeff and exponent of leading term of exp,
3555 ;; after rewriting exp in term of its MRV set omega.
3556 (defun mrv-leadterm (exp var omega
)
3557 (prog ((new-omega ()))
3558 (cond ((freeof var exp
)
3559 (return (list exp
0))))
3560 (dolist (term omega
)
3561 (cond ((subexp exp term
)
3562 (push term new-omega
))))
3563 (setq omega new-omega
)
3565 (setq omega
(mrv exp var
))))
3566 (cond ((member var omega
:test
#'eq
)
3567 (let* ((omega-up (mrv-moveup omega var
))
3568 (e-up (car (mrv-moveup (list exp
) var
)))
3569 (mrv-leadterm-up (mrv-leadterm e-up var omega-up
)))
3570 (return (mrv-movedown mrv-leadterm-up var
)))))
3571 (destructuring-let* ((wsym (gensym "w"))
3574 ((f . logw
) (mrv-rewrite exp omega var wsym
))
3575 (series (calculate-series (mrv-rewrite-logs f wsym logw
)
3577 (setq series
(maxima-substitute logw
`((%log
) ,wsym
) series
))
3578 (setq lo
($lopow series wsym
))
3579 (when (or (not ($constantp lo
))
3580 (not (free series
'%derivative
)))
3581 ;; (mtell "series: ~M lo: ~M~%" series lo)
3582 (tay-error "error in series expansion" f
))
3583 (setq coef
($coeff series wsym lo
))
3584 (when (not (free coef wsym
))
3585 (tay-error "MRV-LEADTERM: failed to extract leading coefficient; obtained" coef
))
3586 ;;(mtell "exp: ~M f: ~M~%" exp f)
3587 ;;(mtell "series: ~M~%coeff: ~M~%pow: ~M~%" series coef lo)
3588 (return (list coef lo
)))))
3590 (defun mrv-moveup (l var
)
3591 (mapcar (lambda (exp)
3592 (simplify-log-of-exp
3593 (syntactic-substitute `((mexpt) $%e
,var
) var exp
)))
3596 (defun mrv-movedown (l var
)
3597 (mapcar (lambda (exp) (syntactic-substitute `((%log simp
) ,var
) var exp
))
3600 ;; test whether sub is a subexpression of exp
3601 (defun subexp (exp sub
)
3602 (let ((dummy (gensym)))
3603 (putprop dummy t
'internal
)
3604 (not (alike1 (maxima-substitute dummy
3609 ;; Generate $lhospitallim terms of taylor expansion.
3610 ;; Ideally we would use a lazy series representation that generates
3611 ;; more terms as higher order terms cancel.
3612 (defun calculate-series (exp var
)
3613 (let ((cntx ($supcontext
)))
3617 (mfuncall '$assume
(ftake 'mgreaterp var
0))
3618 (putprop var t
'internal
); keep var from appearing in questions to user
3619 ($taylor exp var
0 $lhospitallim
))
3620 (remprop var
'internal
)
3621 ($killcontext cntx
))))
3623 (defun mrv-sign (exp var
)
3624 (cond ((freeof var exp
)
3625 (let ((sign ($sign
($radcan exp
))))
3626 (cond ((eq sign
'$zero
)
3632 (t (tay-error " cannot determine mrv-sign" exp
)))))
3636 (* (mrv-sign (cadr exp
) var
)
3637 (mrv-sign (m*l
(cddr exp
)) var
)))
3639 (equal (mrv-sign (cadr exp
) var
) 1))
3642 (cond ((equal (mrv-sign (cadr exp
) var
) -
1)
3643 (tay-error " complex expression in gruntz limit" exp
)))
3644 (mrv-sign (m+ -
1 (cadr exp
)) var
))
3646 (mrv-sign (limitinf exp var
) var
))
3647 (t (tay-error " cannot determine mrv-sign" exp
))))
3649 ;; gruntz algorithm for limit of exp as var goes to positive infinity
3650 (defun limitinf (exp var
)
3651 (prog (($exptsubst nil
))
3652 (cond ((freeof var exp
)
3654 (destructuring-let* ((c0-e0 (mrv-leadterm exp var nil
))
3657 (sig (mrv-sign e0 var
)))
3658 (cond ((equal sig
1)
3661 (cond ((equal (mrv-sign c0 var
) 1)
3663 ((equal (mrv-sign c0 var
) -
1)
3667 ;; example: gruntz(n^n/(n^n+(n-1)^n), n, inf);
3668 (tay-error " infinite recursion in limitinf" exp
))
3669 (return (limitinf c0 var
)))))))
3671 ;; user-level function equivalent to $limit.
3672 ;; direction must be specified if limit point is not infinite
3673 ;; The arguments are checked and a failure of taylor is catched.
3675 (defmfun $gruntz
(expr var val
&rest rest
)
3677 (when (> (length rest
) 1)
3679 (intl:gettext
"gruntz: too many arguments; expected just 3 or 4")))
3680 (setq dir
(car rest
))
3681 (when (and (not (member val
'($inf $minf $zeroa $zerob
)))
3682 (not (member dir
'($plus $minus
))))
3684 (intl:gettext
"gruntz: direction must be 'plus' or 'minus'")))
3686 (catch 'taylor-catch
3687 (let ((silent-taylor-flag t
))
3688 (gruntz1 expr var val dir
))))
3689 (if (or (null ans
) (eq ans t
))
3691 `(($gruntz simp
) ,expr
,var
, val
,dir
)
3692 `(($gruntz simp
) ,expr
,var
,val
))
3695 ;; Additional simplifications for the limit function. Specifically:
3696 ;; (a) replace every mapatom that is declared to be zero by zero
3697 ;; (b) dispatch radcan on expressions of the form (positive integer)^XXX
3698 ;; The mechanism (a) isn't perfect--if a+b is declared to zero, it doesn't
3699 ;; simplify a+b+c to c, for example.
3701 ;; For efficiency, the functionality of factosimp, tansc, lfibtophi,
3702 ;; and limitsimp should be incorporated into extra-simp.
3703 (defun extra-simp (e)
3704 (cond (($subvarp e
) e
) ;return e
3705 ((extended-real-p e
) e
) ;we don't want to call sign on ind, so catch this
3706 (($mapatom e
) ;if e is declared zero, return 0; otherwise e
3707 (if (eq '$zero
($csign e
)) 0 e
))
3708 ;; dispatch radcan on (positive integer)^Y
3709 ((and (mexptp e
) (integerp (cadr e
)) (> (cadr e
) 0))
3710 ($radcan
(ftake 'mexpt
(cadr e
) (extra-simp (caddr e
)))))
3711 (($subvarp
(mop e
)) ;subscripted function
3714 (mapcar #'extra-simp
(subfunsubs e
))
3715 (mapcar #'extra-simp
(subfunargs e
))))
3717 (simplifya (cons (list (caar e
)) (mapcar #'extra-simp
(cdr e
))) t
))))
3719 ;; This function is for internal use in $limit.
3721 ;; The function gruntz1 standardizes the limit point to inf and the limit variable
3722 ;; to a gensym. Since the limit point is possibly altered by this function, we
3723 ;; need to make the appropriate assumptions on the limit variable. This is done
3725 (defun gruntz1 (exp var val
&rest rest
)
3726 (cond ((> (length rest
) 1)
3727 (merror (intl:gettext
"gruntz: too many arguments; expected just 3 or 4"))))
3728 (let (($logexpand t
) ; gruntz needs $logexpand T
3729 (newvar (gensym "w"))
3731 (putprop newvar t
'internal
); keep var from appearing in questions to user
3732 (cond ((eq val
'$inf
)
3733 (setq exp
(maxima-substitute newvar var exp
)))
3735 (setq exp
(maxima-substitute (m* -
1 newvar
) var exp
)))
3737 (setq exp
(maxima-substitute (m// 1 newvar
) var exp
)))
3739 (setq exp
(maxima-substitute (m// -
1 newvar
) var exp
)))
3741 (setq exp
(maxima-substitute (m+ val
(m// 1 newvar
)) var exp
)))
3743 (setq exp
(maxima-substitute (m+ val
(m// -
1 newvar
)) var exp
)))
3744 (t (merror (intl:gettext
"gruntz: direction must be 'plus' or 'minus'; found: ~M") dir
)))
3745 (let ((cx ($supcontext
)))
3748 (mfuncall '$assume
(ftake 'mlessp
*large-positive-number
* newvar
)) ; *large-positive-number* < newvar
3749 (mfuncall '$assume
(ftake 'mlessp
0 'lim-epsilon
)) ; 0 < lim-epsilon
3750 (mfuncall '$assume
(ftake 'mlessp
*large-positive-number
* 'prin-inf
)) ; *large-positive-number* < prin-inf
3751 (mfuncall '$activate cx
) ;not sure this is needed, but OK
3752 (setq exp
(resimplify exp
)) ;simplify in new context
3753 (setq exp
(extra-simp (sratsimp exp
))) ;additional simplifications
3754 (limitinf exp newvar
)) ;compute & return limit
3755 ($killcontext cx
))))) ;kill context & forget all new facts.
3757 ;; substitute y for x in exp
3758 ;; similar to maxima-substitute but does not simplify result
3759 (defun syntactic-substitute (y x exp
)
3760 (cond ((alike1 x exp
) y
)
3763 (mapcar (lambda (exp)
3764 (syntactic-substitute y x exp
))
3767 ;; log(exp(subexpr)) -> subexpr
3768 ;; without simplifying entire exp
3769 (defun simplify-log-of-exp (exp)
3770 (cond ((atom exp
) exp
)
3773 (eq '$%e
(cadadr exp
)))
3776 (mapcar #'simplify-log-of-exp
