1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 4;;; ;;;;; 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;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module limit) 14 15 16;;; ************************************************************** 17;;; ** ** 18;;; ** LIMIT PACKAGE ** 19;;; ** ** 20;;; ************************************************************** 21 22;;; I believe a large portion of this file is described in Paul 23;;; Wang's thesis, "Evaluation of Definite Integrals by Symbolic 24;;; Integration," MIT/LCS/TR-92, Oct. 1971. This can be found at 25;;; http://www.lcs.mit.edu/publications/specpub.php?id=660, but some 26;;; important pages are black. 27 28;;; TOP LEVEL FUNCTION(S): $LIMIT $LDEFINT 29 30(declare-top (special errorsw origval $lhospitallim low* 31 *indicator half%pi nn* dn* numer denom exp var val varlist 32 *zexptsimp? $tlimswitch $logarc taylored logcombed 33 $exponentialize lhp? lhcount $ratfac genvar 34 loginprod? $limsubst $logabs a context limit-assumptions 35 limit-top limitp integer-info old-integer-info $keepfloat $logexpand)) 36 37(defconstant +behavior-count+ 4) 38(defvar *behavior-count-now*) 39(defvar *getsignl-asksign-ok* nil) 40 41(load-macsyma-macros rzmac) 42 43(defmvar infinities '($inf $minf $infinity) 44 "The types of infinities recognized by Maxima. 45 INFINITY is complex infinity") 46 47(defmvar real-infinities '($inf $minf) 48 "The real infinities, `inf' is positive infinity, `minf' negative infinity") 49 50(defmvar infinitesimals '($zeroa $zerob) 51 "The infinitesimals recognized by Maxima. ZEROA zero from above, 52 ZEROB zero from below") 53 54(defmvar simplimplus-problems () 55 "A list of all problems in the stack of recursive calls to simplimplus.") 56 57(defmvar limit-answers () 58 "An association list for storing limit answers.") 59 60(defmvar limit-using-taylor () 61 "Is the current limit computation using taylor expansion?") 62 63(defmvar preserve-direction () "Makes `limit' return Direction info.") 64 65(unless (boundp 'integer-info) (setq integer-info ())) 66 67;; This should be made to give more information about the error. 68;;(DEFun DISCONT () 69;; (cond (errorsw (throw 'errorsw t)) 70;; (t (merror "Discontinuity Encountered")))) 71 72;;(DEFUN PUTLIMVAL (E V) 73;; (let ((exp (cons '(%limit) (list e var val)))) 74;; (cond ((not (assolike exp limit-answers)) 75;; (setq limit-answers (cons (cons exp v) limit-answers)) 76;; v) 77;; (t ())))) 78 79(defun putlimval (e v &aux exp) 80 (setq exp `((%limit) ,e ,var ,val)) 81 (unless (assolike exp limit-answers) 82 (push (cons exp v) limit-answers)) 83 v) 84 85(defun getlimval (e) 86 (let ((exp (cons '(%limit) (list e var val)))) 87 (assolike exp limit-answers))) 88 89(defmacro limit-catch (exp var val) 90 `(let ((errorsw t)) 91 (let ((ans (catch 'errorsw 92 (catch 'limit (limit ,exp ,var ,val 'think))))) 93 (if (or (null ans) (eq ans t)) 94 () 95 ans)))) 96 97(defmfun $limit (&rest args) 98 (let ((first-try (apply #'toplevel-$limit args))) 99 (if (and (consp first-try) (eq (caar first-try) '%limit)) 100 (let ((*getsignl-asksign-ok* t)) 101 (apply #'toplevel-$limit args)) 102 first-try))) 103 104(defun toplevel-$limit (&rest args) 105 (let ((limit-assumptions ()) 106 (old-integer-info ()) 107 ($keepfloat t) 108 (limit-top t)) 109 (declare (special limit-assumptions old-integer-info 110 $keepfloat limit-top)) 111 (unless limitp 112 (setq old-integer-info integer-info) 113 (setq integer-info ())) 114 115 (unwind-protect 116 (let ((exp1 ()) (lhcount $lhospitallim) (*behavior-count-now* 0) 117 (exp ()) (var ()) (val ()) (dr ()) 118 (*indicator ()) (taylored ()) (origval ()) 119 (logcombed ()) (lhp? ()) 120 (varlist ()) (ans ()) (genvar ()) (loginprod? ()) 121 (limit-answers ()) (limitp t) (simplimplus-problems ()) 122 (lenargs (length args)) 123 (genfoo ())) 124 (declare (special lhcount *behavior-count-now* exp var val *indicator 125 taylored origval logcombed lhp? 126 varlist genvar loginprod? limitp)) 127 (prog () 128 (unless (or (= lenargs 3) (= lenargs 4) (= lenargs 1)) 129 (wna-err '$limit)) 130 ;; Is it a LIST of Things? 131 (when (setq ans (apply #'limit-list args)) 132 (return ans)) 133 (setq exp1 (specrepcheck (first args))) 134 (when (and (atom exp1) 135 (member exp1 '(nil t))) 136 ;; The expression is 'T or 'NIL. Return immediately. 137 (return exp1)) 138 (cond ((= lenargs 1) 139 (setq var (setq genfoo (gensym)) ; Use a gensym. Not foo. 140 val 0)) 141 (t 142 (setq var (second args)) 143 (when ($constantp var) 144 (merror (intl:gettext "limit: second argument must be a variable, not a constant; found: ~M") var)) 145 (unless (or ($subvarp var) (atom var)) 146 (merror (intl:gettext "limit: variable must be a symbol or subscripted symbol; found: ~M") var)) 147 (setq val (infsimp (third args))) 148 ;; infsimp converts -inf to minf. it also converts -infinity to 149 ;; infinity, although perhaps this should generate the error below. 150 (when (and (not (atom val)) 151 (some #'(lambda (x) (not (freeof x val))) 152 infinities)) 153 (merror (intl:gettext "limit: third argument must be a finite value or one of: inf, minf, infinity; found: ~M") val)) 154 (when (eq val '$zeroa) (setq dr '$plus)) 155 (when (eq val '$zerob) (setq dr '$minus)))) 156 (cond ((= lenargs 4) 157 (unless (member (fourth args) '($plus $minus) :test #'eq) 158 (merror (intl:gettext "limit: direction must be either 'plus' or 'minus'; found: ~M") (fourth args))) 159 (setq dr (fourth args)))) 160 (if (and (atom var) (not (among var val))) 161 (setq exp exp1) 162 (let ((realvar var)) ;; Var is funny so make it a gensym. 163 (setq var (gensym)) 164 (setq exp (maxima-substitute var realvar exp1)) 165 (putprop var realvar 'limitsub))) 166 (unless (or $limsubst (eq var genfoo)) 167 (when (limunknown exp) 168 (return `((%limit) ,@(cons exp1 (cdr args)))))) 169 (setq varlist (ncons var) genvar nil origval val) 170 ;; Transform limits to minf to limits to inf by 171 ;; replacing var with -var everywhere. 172 (when (eq val '$minf) 173 (setq val '$inf 174 origval '$inf 175 exp (subin (m* -1 var) exp))) 176 177 ;; Hide noun form of %derivative, %integrate. 178 (setq exp (hide exp)) 179 180 ;; Transform the limit value. 181 (unless (infinityp val) 182 (unless (zerop2 val) 183 (let ((*atp* t) (realvar var)) 184 ;; *atp* prevents substitution from applying to vars 185 ;; bound by %sum, %product, %integrate, %limit 186 (setq var (gensym)) 187 (putprop var t 'internal) 188 (setq exp (maxima-substitute (m+ val var) realvar exp)))) 189 (setq val (cond ((eq dr '$plus) '$zeroa) 190 ((eq dr '$minus) '$zerob) 191 (t 0))) 192 (setq origval 0)) 193 194 ;; Make assumptions about limit var being very small or very large. 195 ;; Assumptions are forgotten upon exit. 196 (unless (= lenargs 1) 197 (limit-context var val dr)) 198 199 ;; Resimplify in light of new assumptions. 200 (setq exp (resimplify 201 (factosimp 202 (tansc 203 (lfibtophi 204 (limitsimp ($expand exp 1 0) var)))))) 205 206 (if (not (or (real-epsilonp val) ;; if direction of limit not specified 207 (infinityp val))) 208 (setq ans (both-side exp var val)) ;; compute from both sides 209 (let ((d (catch 'mabs (mabs-subst exp var val)))) 210 (cond ;; otherwise try to remove absolute value 211 ((eq d '$und) (return '$und)) 212 ((eq d 'retn) ) 213 (t (setq exp d))) 214 (setq ans (limit-catch exp var val));; and find limit from one side 215 216 ;; try gruntz 217 (if (not ans) 218 (setq ans (catch 'taylor-catch 219 (let ((silent-taylor-flag t)) 220 (declare (special silent-taylor-flag)) 221 (gruntz1 exp var val))))) 222 223 ;; try taylor series expansion if simple limit didn't work 224 (if (and (null ans) ;; if no limit found and 225 $tlimswitch ;; user says ok to use taylor and 226 (not limit-using-taylor));; not already doing taylor 227 (let ((limit-using-taylor t)) 228 (declare (special limit-using-taylor)) 229 (setq ans (limit-catch exp var val)))))) 230 231 (if ans 232 (return (clean-limit-exp ans)) 233 (return (cons '(%limit) args))))) ;; failure: return nounform 234 (restore-assumptions)))) 235 236(defun clean-limit-exp (exp) 237 (setq exp (restorelim exp)) 238 (if preserve-direction exp (ridofab exp))) 239 240(defun limit-list (exp1 &rest rest) 241 (if (mbagp exp1) 242 `(,(car exp1) ,@(mapcar #'(lambda (x) (apply #'toplevel-$limit `(,x ,@rest))) (cdr exp1))) 243 ())) 244 245(defun limit-context (var val direction) ;Only works on entry! 246 (cond (limit-top 247 (assume '((mgreaterp) lim-epsilon 0)) 248 (assume '((mgreaterp) prin-inf 100000000)) 249 (setq limit-assumptions (make-limit-assumptions var val direction)) 250 (setq limit-top ())) 251 (t ())) 252 limit-assumptions) 253 254(defun make-limit-assumptions (var val direction) 255 (let ((new-assumptions)) 256 (cond ((or (null var) (null val)) 257 ()) 258 ((and (not (infinityp val)) (null direction)) 259 ()) 260 ((eq val '$inf) 261 `(,(assume `((mgreaterp) ,var 100000000)) ,@new-assumptions)) 262 ((eq val '$minf) 263 `(,(assume `((mgreaterp) -100000000 ,var)) ,@new-assumptions)) 264 ((eq direction '$plus) 265 `(,(assume `((mgreaterp) ,var 0)) ,@new-assumptions)) ;All limits around 0 266 ((eq direction '$minus) 267 `(,(assume `((mgreaterp) 0 ,var)) ,@new-assumptions)) 268 (t 269 ())))) 270 271(defun restore-assumptions () 272;;;Hackery until assume and forget take reliable args. Nov. 9 1979. 273;;;JIM. 274 (do ((assumption-list limit-assumptions (cdr assumption-list))) 275 ((null assumption-list) t) 276 (forget (car assumption-list))) 277 (forget '((mgreaterp) lim-epsilon 0)) 278 (forget '((mgreaterp) prin-inf 100000000)) 279 (cond ((and (not (null integer-info)) 280 (not limitp)) 281 (do ((list integer-info (cdr list))) 282 ((null list) t) 283 (i-$remove `(,(cadar list) ,(caddar list)))) 284 (setq integer-info old-integer-info)))) 285 286;; The optional arg allows the caller to decide on the value of 287;; preserve-direction. Default is T, like it used to be. 288(defun both-side (exp var val &optional (preserve t)) 289 (let* ((preserve-direction preserve) 290 (la (toplevel-$limit exp var val '$plus)) lb) 291 (when (eq la '$und) (return-from both-side '$und)) 292 (setf lb (toplevel-$limit exp var val '$minus)) 293 (let ((ra (ridofab la)) 294 (rb (ridofab lb))) 295 (cond ((eq t (meqp ra rb)) 296 ra) 297 ((and (eq ra '$ind) 298 (eq rb '$ind)) 299 ; Maxima does not consider equal(ind,ind) to be true, but 300 ; if both one-sided limits are ind then we want to call 301 ; the two-sided limit ind (e.g., limit(sin(1/x),x,0)). 302 '$ind) 303 ((or (not (free la '%limit)) 304 (not (free lb '%limit))) 305 ()) 306 ((and (infinityp la) (infinityp lb)) 307 ; inf + minf => infinity 308 '$infinity) 309 (t 310 '$und))))) 311 312(defun limunknown (f) 313 (catch 'limunknown (limunknown1 (specrepcheck f)))) 314 315(defun limunknown1 (f) 316 (cond ((mapatom f) nil) 317 ((or (not (safe-get (caar f) 'operators)) 318 (member (caar f) '(%sum %product mncexpt) :test #'eq) 319 ;;Special function code here i.e. for li[2](x). 320 (and (eq (caar f) 'mqapply) 321 (not (get (subfunname f) 'specsimp)))) 322 (if (not (free f var)) (throw 'limunknown t))) 323 (t (mapc #'limunknown1 (cdr f)) nil))) 324 325(defun factosimp(e) 326 (if (involve e '(%gamma)) (setq e ($makefact e))) 327 (cond ((involve e '(mfactorial)) 328 (setq e (simplify ($minfactorial e)))) 329 (t e))) 330 331;; returns 1, 0, -1 332;; or nil if sign unknown or complex 333(defun getsignl (z) 334 (let ((z (ridofab z))) 335 (if (not (free z var)) (setq z (toplevel-$limit z var val))) 336 (let ((*complexsign* t)) 337 (let ((sign (if *getsignl-asksign-ok* ($asksign z) ($sign z)))) 338 (cond ((eq sign '$pos) 1) 339 ((eq sign '$neg) -1) 340 ((eq sign '$zero) 0)))))) 341 342(defun restorelim (exp) 343 (cond ((null exp) nil) 344 ((atom exp) (or (and (symbolp exp) (get exp 'limitsub)) exp)) 345 ((and (consp (car exp)) (eq (caar exp) 'mrat)) 346 (cons (car exp) 347 (cons (restorelim (cadr exp)) 348 (restorelim (cddr exp))))) 349 (t (cons (car exp) (mapcar #'restorelim (cdr exp)))))) 350 351 352(defun mabs-subst (exp var val) ; RETURNS EXP WITH MABS REMOVED, OR THROWS. 353 (let ((d (involve exp '(mabs))) 354 arglim) 355 (cond ((null d) exp) 356 (t (cond 357 ((not (and (equal ($imagpart (let ((v (limit-catch d var val))) 358 ;; The above call might 359 ;; throw 'limit, so we 360 ;; need to catch it. If 361 ;; we can't find the 362 ;; limit without ABS, we 363 ;; assume the limit is 364 ;; undefined. Is this 365 ;; right? Anyway, this 366 ;; fixes Bug 1548643. 367 (unless v 368 (throw 'mabs '$und)) 369 (setq arglim v))) 370 0) 371 (equal ($imagpart var) 0))) 372 (cond ((eq arglim '$infinity) 373 ;; Check for $infinity as limit of argument. 374 '$inf) 375 (t 376 (throw 'mabs 'retn)))) 377 (t (do ((ans d (involve exp '(mabs))) (a () ())) 378 ((null ans) exp) 379 (setq a (mabs-subst ans var val)) 380 (setq d (limit a var val t)) 381 (cond 382 ((and a d) 383 (cond ((zerop1 d) 384 (setq d (behavior a var val)) 385 (if (zerop1 d) (throw 'mabs 'retn)))) 386 (if (eq d '$und) 387 (throw 'mabs d)) 388 (cond ((or (eq d '$zeroa) (eq d '$inf) 389 (eq d '$ind) 390 ;; fails on limit(abs(sin(x))/sin(x), x, inf) 391 (eq ($sign d) '$pos)) 392 (setq exp (maxima-substitute a `((mabs) ,ans) exp))) 393 ((or (eq d '$zerob) (eq d '$minf) 394 (eq ($sign d) '$neg)) 395 (setq exp (maxima-substitute (m* -1 a) `((mabs) ,ans) exp))) 396 (t (throw 'mabs 'retn)))) 397 (t (throw 'mabs 'retn)))))))))) 398 399;; Called on an expression that might contain $INF, $MINF, $ZEROA, $ZEROB. Tries 400;; to simplify it to sort out things like inf^inf or inf+1. 401(defun simpinf (exp) 402 (simpinf-ic exp (count-general-inf exp))) 403 404(defun count-general-inf (expr) 405 (count-atoms-matching 406 (lambda (x) (or (infinityp x) (real-epsilonp x))) expr)) 407 408(defun count-atoms-matching (predicate expr) 409 "Count the number of atoms in the Maxima expression EXPR matching PREDICATE, 410ignoring dummy variables and array indices." 411 (cond 412 ((atom expr) (if (funcall predicate expr) 1 0)) 413 ;; Don't count atoms that occur as a limit of %integrate, %sum, %product, 414 ;; %limit etc. 415 ((member (caar expr) dummy-variable-operators) 416 (count-atoms-matching predicate (cadr expr))) 417 ;; Ignore array indices 418 ((member 'array (car expr)) 0) 419 (t (loop 420 for arg in (cdr expr) 421 summing (count-atoms-matching predicate arg))))) 422 423(defun simpinf-ic (exp &optional infinity-count) 424 (case infinity-count 425 ;; A very slow identity transformation... 426 (0 exp) 427 428 ;; If there's only one infinity, we replace it by a variable and take the 429 ;; limit as that variable goes to infinity. Use $gensym in case we can't 430 ;; compute the answer and the limit leaks out. 431 (1 (let* ((val (or (inf-typep exp) (epsilon-typep exp))) 432 (var ($gensym)) 433 (expr (subst var val exp)) 434 (limit (toplevel-$limit expr var val))) 435 (cond 436 ;; Now we look to see whether the computed limit is any simpler than 437 ;; what we shoved in (which we'll define as "doesn't contain EXPR as a 438 ;; subtree"). If so, return it. 439 ((not (subtree-p expr limit :test #'equal)) 440 limit) 441 442 ;; Otherwise, return the original form: apparently, we can't compute 443 ;; the limit we needed, and it's uglier than what we started with. 444 (t exp)))) 445 446 ;; If more than one infinity, we have to be a bit more careful. 447 (otherwise 448 (let* ((arguments (mapcar 'simpinf (cdr exp))) 449 (new-expression (cons (list (caar exp)) arguments)) 450 infinities-left) 451 (cond 452 ;; If any of the arguments are undefined, we are too. 453 ((among '$und arguments) '$und) 454 ;; If we ended up with something indeterminate, we punt and just return 455 ;; the input. 456 ((amongl '(%limit $ind) arguments) exp) 457 458 ;; Exponentiation & multiplication 459 ((mexptp exp) (simpinf-expt (first arguments) (second arguments))) 460 ((mtimesp exp) (simpinf-times arguments)) 461 462 ;; Down to at most one infinity? We do this after exponentiation to 463 ;; avoid zeroa^zeroa => 0^0, which will raise an error rather than just 464 ;; returning und. We do it after multiplication to avoid zeroa * inf => 465 ;; 0 * inf => 0. 466 ((<= (setf infinities-left (count-general-inf new-expression)) 1) 467 (simpinf-ic new-expression infinities-left)) 468 469 ;; Addition 470 ((mplusp exp) (simpinf-plus arguments)) 471 472 ;; Give up! 473 (t new-expression)))))) 474 475(defun simpinf-times (arguments) 476 (declare (special exp var val)) 477 ;; When we have a product, we need to spot that zeroa * zerob = zerob, zeroa * 478 ;; inf = und etc. Note that (SIMPINF '$ZEROA) => 0, so a nonzero atom is not 479 ;; an infinitesimal. Moreover, we can assume that each of ARGUMENTS is either 480 ;; a number, computed successfully by the recursive SIMPINF call, or maybe a 481 ;; %LIMIT noun-form (in which case, we aren't going to be able to tell the 482 ;; answer). 483 (cond 484 ((member 0 arguments) 485 (cond 486 ((find-if #'infinityp arguments) '$und) 487 ((every #'atom arguments) 0) 488 (t exp))) 489 490 ((member '$infinity arguments) 491 (if (every #'atom arguments) 492 '$infinity 493 exp)) 494 495 (t (simplimit (cons '(mtimes) arguments) var val)))) 496 497(defun simpinf-expt (base exponent) 498 ;; In the comments below, zero* represents one of 0, zeroa, zerob. 499 ;; 500 ;; TODO: In some cases we give up too early. E.g. inf^(2 + 1/inf) => inf^2 501 ;; (which should simplify to inf) 502 (case base 503 ;; inf^inf = inf 504 ;; inf^minf = 0 505 ;; inf^zero* = und 506 ;; inf^foo = inf^foo 507 ($inf 508 (case exponent 509 ($inf '$inf) 510 ($minf 0) 511 ((0 $zeroa $zerob) '$und) 512 (t (list '(mexpt) base exponent)))) 513 ;; minf^inf = infinity <== Or should it be und? 514 ;; minf^minf = 0 515 ;; minf^zero* = und 516 ;; minf^foo = minf^foo 517 ($minf 518 (case exponent 519 ($inf '$infinity) 520 ($minf 0) 521 ((0 $zeroa $zerob) '$und) 522 (t (list '(mexpt) base exponent)))) 523 ;; zero*^inf = 0 524 ;; zero*^minf = und 525 ;; zero*^zero* = und 526 ;; zero*^foo = zero*^foo 527 ((0 $zeroa $zerob) 528 (case exponent 529 ($inf 0) 530 ($minf '$und) 531 ((0 $zeroa $zerob) '$und) 532 (t (list '(mexpt) base exponent)))) 533 ;; a^b where a is pretty much anything except for a naked 534 ;; inf,minf,zeroa,zerob or 0. 535 (t 536 (cond 537 ;; When a isn't crazy, try a^b = e^(b log(a)) 538 ((not (amongl (append infinitesimals infinities) base)) 539 (simpinf (m^ '$%e (m* exponent `((%log) ,base))))) 540 541 ;; No idea. Just return what we've found so far. 542 (t (list '(mexpt) base exponent)))))) 543 544(defun simpinf-plus (arguments) 545 ;; We know that none of the arguments are infinitesimals, since SIMPINF never 546 ;; returns one of them. As such, we partition our arguments into infinities 547 ;; and everything else. The latter won't have any "hidden" infinities like 548 ;; lim(x,x,inf), since SIMPINF gave up on anything containing a %lim already. 549 (let ((bigs) (others)) 550 (dolist (arg arguments) 551 (cond ((infinityp arg) (push arg bigs)) 552 (t (push arg others)))) 553 (cond 554 ;; inf + minf or the like 555 ((cdr (setf bigs (delete-duplicates bigs))) '$und) 556 ;; inf + smaller + stuff 557 (bigs (car bigs)) 558 ;; I don't think this can happen, since SIMPINF goes back to the start if 559 ;; there are fewer than two infinities in the arguments, but let's be 560 ;; careful. 561 (t (cons '(mplus) others))))) 562 563;; Simplify expression with zeroa or zerob. 564(defun simpab (small) 565 (cond ((null small) ()) 566 ((member small '($zeroa $zerob $inf $minf $infinity) :test #'eq) small) 567 ((not (free small '$ind)) '$ind) ;Not exactly right but not 568 ((not (free small '$und)) '$und) ;causing trouble now. 569 ((mapatom small) small) 570 (t (let ((preserve-direction t) 571 (new-small (subst (m^ '$inf -1) '$zeroa 572 (subst (m^ '$minf -1) '$zerob small)))) 573 (simpinf new-small))))) 574 575 576;;;*I* INDICATES: T => USE LIMIT1,THINK, NIL => USE SIMPLIMIT. 577(defun limit (exp var val *i*) 578 (cond 579 ((among '$und exp) '$und) 580 ((eq var exp) val) 581 ((atom exp) exp) 582 ((not (among var exp)) 583 (cond ((amongl '($inf $minf $infinity $ind) exp) 584 (simpinf exp)) 585 ((amongl '($zeroa $zerob) exp) 586 ;; Simplify expression with zeroa or zerob. 587 (simpab exp)) 588 (t exp))) 589 ((getlimval exp)) 590 (t (putlimval exp (cond ((and limit-using-taylor 591 (null taylored) 592 (tlimp exp)) 593 (taylim exp var val *i*)) 594 ((ratp exp var) (ratlim exp)) 595 ((or (eq *i* t) (radicalp exp var)) 596 (limit1 exp var val)) 597 ((eq *i* 'think) 598 (cond ((or (mtimesp exp) (mexptp exp)) 599 (limit1 exp var val)) 600 (t (simplimit exp var val)))) 601 (t (simplimit exp var val))))))) 602 603(defun limitsimp (exp var) 604 (limitsimp-expt (sin-sq-cos-sq-sub exp) var)) 605;;Hack for sin(x)^2+cos(x)^2. 606 607;; if var appears in base and power of expt, 608;; push var into power of of expt 609(defun limitsimp-expt (exp var) 610 (cond ((or (atom exp) 611 (mnump exp) 612 (freeof var exp)) exp) 613 ((and (mexptp exp) 614 (not (freeof var (cadr exp))) 615 (not (freeof var (caddr exp)))) 616 (m^ '$%e (simplify `((%log) ,exp)))) 617 (t (subst0 (cons (cons (caar exp) ()) 618 (mapcar #'(lambda (x) 619 (limitsimp-expt x var)) 620 (cdr exp))) 621 exp)))) 622 623(defun sin-sq-cos-sq-sub (exp) ;Hack ... Hack 624 (let ((arg (involve exp '(%sin %cos)))) 625 (cond 626 ((null arg) exp) 627 (t (let ((new-exp ($substitute (m+t 1 (m- (m^t `((%sin simp) ,arg) 2))) 628 (m^t `((%cos simp) ,arg) 2) 629 ($substitute 630 (m+t 1 (m- (m^t `((%cos simp) ,arg) 2))) 631 (m^t `((%sin simp) ,arg) 2) 632 exp)))) 633 (cond ((not (involve new-exp '(%sin %cos))) new-exp) 634 (t exp))))))) 635 636(defun expand-trigs (x var) 637 (cond ((atom x) x) 638 ((mnump x) x) 639 ((and (or (eq (caar x) '%sin) 640 (eq (caar x) '%cos)) 641 (not (free (cadr x) var))) 642 ($trigexpand x)) 643 ((member 'array (car x)) 644 ;; Some kind of array reference. Return it. 645 x) 646 (t (simplify (cons (ncons (caar x)) 647 (mapcar #'(lambda (x) 648 (expand-trigs x var)) 649 (cdr x))))))) 650 651 652(defun tansc (e) 653 (cond ((not (involve e 654 '(%cot %csc %binomial 655 %sec %coth %sech %csch 656 %acot %acsc %asec %acoth 657 %asech %acsch 658 %jacobi_ns %jacobi_nc %jacobi_cs 659 %jacobi_ds %jacobi_dc))) 660 e) 661 (t ($ratsimp (tansc1 e))))) 662 663(defun tansc1 (e &aux tem) 664 (cond ((atom e) e) 665 ((and (setq e (cons (car e) (mapcar 'tansc1 (cdr e)))) ())) 666 ((setq tem (assoc (caar e) '((%cot . %tan) (%coth . %tanh) 667 (%sec . %cos) (%sech . %cosh) 668 (%csc . %sin) (%csch . %sinh)) :test #'eq)) 669 (tansc1 (m^ (list (ncons (cdr tem)) (cadr e)) -1.))) 670 ((setq tem (assoc (caar e) '((%jacobi_nc . %jacobi_cn) 671 (%jacobi_ns . %jacobi_sn) 672 (%jacobi_cs . %jacobi_sc) 673 (%jacobi_ds . %jacobi_sd) 674 (%jacobi_dc . %jacobi_cd)) :test #'eq)) 675 ;; Converts Jacobi elliptic function to its reciprocal 676 ;; function. 677 (tansc1 (m^ (list (ncons (cdr tem)) (cadr e) (third e)) -1.))) 678 ((setq tem (member (caar e) '(%sinh %cosh %tanh) :test #'eq)) 679 (let (($exponentialize t)) 680 (resimplify e))) 681 ((setq tem (assoc (caar e) '((%acsc . %asin) (%asec . %acos) 682 (%acot . %atan) (%acsch . %asinh) 683 (%asech . %acosh) (%acoth . %atanh)) :test #'eq)) 684 (list (ncons (cdr tem)) (m^t (cadr e) -1.))) 685 ((and (eq (caar e) '%binomial) (among var (cdr e))) 686 (m// `((mfactorial) ,(cadr e)) 687 (m* `((mfactorial) ,(m+t (cadr e) (m- (caddr e)))) 688 `((mfactorial) ,(caddr e))))) 689 (t e))) 690 691(defun hyperex (ex) 692 (cond ((not (involve ex '(%sin %cos %tan %asin %acos %atan 693 %sinh %cosh %tanh %asinh %acosh %atanh))) 694 ex) 695 (t (hyperex0 ex)))) 696 697(defun hyperex0 (ex) 698 (cond ((atom ex) ex) 699 ((eq (caar ex) '%sinh) 700 (m// (m+ (m^ '$%e (cadr ex)) (m- (m^ '$%e (m- (cadr ex))))) 701 2)) 702 ((eq (caar ex) '%cosh) 703 (m// (m+ (m^ '$%e (cadr ex)) (m^ '$%e (m- (cadr ex)))) 704 2)) 705 ((and (member (caar ex) 706 '(%sin %cos %tan %asin %acos %atan %sinh 707 %cosh %tanh %asinh %acosh %atanh) :test #'eq) 708 (among var ex)) 709 (hyperex1 ex)) 710 (t (cons (car ex) (mapcar #'hyperex0 (cdr ex)))))) 711 712(defun hyperex1 (ex) 713 (resimplify ex)) 714 715;;Used by tlimit also. 716(defun limit1 (exp var val) 717 (prog () 718 (let ((lhprogress? lhp?) 719 (lhp? ()) 720 (ans ())) 721 (cond ((setq ans (and (not (atom exp)) (getlimval exp))) 722 (return ans)) 723 ((and (not (infinityp val)) (setq ans (simplimsubst val exp))) 724 (return ans)) 725 (t nil)) 726;;;NUMDEN* => (values numerator denominator) 727 (multiple-value-bind (n dn) 728 (numden* exp) 729 (cond ((not (among var dn)) 730 (return (simplimit (m// (simplimit n var val) dn) var val))) 731 ((not (among var n)) 732 (return (simplimit (m* n (simplimexpt dn -1 (simplimit dn var val) -1)) var val))) 733 ((and lhprogress? 734 (/#alike n (car lhprogress?)) 735 (/#alike dn (cdr lhprogress?))) 736 (throw 'lhospital nil))) 737 (return (limit2 n dn var val)))))) 738 739(defun /#alike (e f) 740 (if (alike1 e f) 741 t 742 (let ((deriv (sdiff (m// e f) var))) 743 (cond ((=0 deriv) t) 744 ((=0 ($ratsimp deriv)) t) 745 (t nil))))) 746 747(defun limit2 (n dn var val) 748 (prog (n1 d1 lim-sign gcp sheur-ans) 749 (setq n (hyperex n) dn (hyperex dn)) 750;;;Change to uniform limit call. 751 (cond ((infinityp val) 752 (setq d1 (limit dn var val nil)) 753 (setq n1 (limit n var val nil))) 754 (t (cond ((setq n1 (simplimsubst val n)) nil) 755 (t (setq n1 (limit n var val nil)))) 756 (cond ((setq d1 (simplimsubst val dn)) nil) 757 (t (setq d1 (limit dn var val nil)))))) 758 (cond ((or (null n1) (null d1)) (return nil)) 759 (t (setq n1 (sratsimp n1) d1 (sratsimp d1)))) 760 (cond ((or (involve n '(mfactorial)) (involve dn '(mfactorial))) 761 (let ((ans (limfact2 n dn var val))) 762 (cond (ans (return ans)))))) 763 (cond ((and (zerop2 n1) (zerop2 d1)) 764 (cond ((not (equal (setq gcp (gcpower n dn)) 1)) 765 (return (colexpt n dn gcp))) 766 ((and (real-epsilonp val) 767 (not (free n '%log)) 768 (not (free dn '%log))) 769 (return (liminv (m// n dn)))) 770 ((setq n1 (try-lhospital-quit n dn nil)) 771 (return n1)))) 772 ((and (zerop2 n1) (not (member d1 '($ind $und) :test #'eq))) (return 0)) 773 ((zerop2 d1) 774 (setq n1 (ridofab n1)) 775 (return (simplimtimes `(,n1 ,(simplimexpt dn -1 d1 -1)))))) 776 (setq n1 (ridofab n1)) 777 (setq d1 (ridofab d1)) 778 (cond ((or (eq d1 '$und) 779 (and (eq n1 '$und) (not (real-infinityp d1)))) 780 (return '$und)) 781 ((eq d1 '$ind) 782 ;; At this point we have n1/$ind. Look if n1 is one of the 783 ;; infinities or zero. 784 (cond ((and (infinityp n1) (eq ($sign dn) '$pos)) 785 (return n1)) 786 ((and (infinityp n1) (eq ($sign dn) '$neg)) 787 (return (simpinf (m* -1 n1)))) 788 ((and (not (eq n1 '$ind)) 789 (eq ($csign n1) '$zero)) 790 (return 0)) 791 (t (return '$und)))) 792 ((eq n1 '$ind) (return (cond ((infinityp d1) 0) 793 ((equal d1 0) '$und) 794 (t '$ind)))) ;SET LB 795 ((and (real-infinityp d1) (member n1 '($inf $und $minf) :test #'eq)) 796 (cond ((and (not (atom dn)) (not (atom n)) 797 (cond ((not (equal (setq gcp (gcpower n dn)) 1)) 798 (return (colexpt n dn gcp))) 799 ((and (eq '$inf val) 800 (or (involve dn '(mfactorial %gamma)) 801 (involve n '(mfactorial %gamma)))) 802 (return (limfact n dn)))))) 803 ((eq n1 d1) (setq lim-sign 1) (go cp)) 804 (t (setq lim-sign -1) (go cp)))) 805 ((and (infinityp d1) (infinityp n1)) 806 (setq lim-sign (if (or (eq d1 '$minf) (eq n1 '$minf)) -1 1)) 807 (go cp)) 808 (t (return (simplimtimes `(,n1 ,(m^ d1 -1)))))) 809 cp (setq n ($expand n) dn ($expand dn)) 810 (cond ((mplusp n) 811 (let ((new-n (m+l (maxi (cdr n))))) 812 (cond ((not (alike1 new-n n)) 813 (return (limit (m// new-n dn) var val 'think)))) 814 (setq n1 new-n))) 815 (t (setq n1 n))) 816 (cond ((mplusp dn) 817 (let ((new-dn (m+l (maxi (cdr dn))))) 818 (cond ((not (alike1 new-dn dn)) 819 (return (limit (m// n new-dn) var val 'think)))) 820 (setq d1 new-dn))) 821 (t (setq d1 dn))) 822 (setq sheur-ans (sheur0 n1 d1)) 823 (cond ((or (member sheur-ans '($inf $zeroa) :test #'eq) 824 (free sheur-ans var)) 825 (return (simplimtimes `(,lim-sign ,sheur-ans)))) 826 ((and (alike1 sheur-ans dn) 827 (not (mplusp n)))) 828 ((member (setq n1 (cond ((expfactorp n1 d1) (expfactor n1 d1 var)) 829 (t ()))) 830 '($inf $zeroa) :test #'eq) 831 (return n1)) 832 ((not (null (setq n1 (cond ((expfactorp n dn) (expfactor n dn var)) 833 (t ()))))) 834 (return n1)) 835 ((and (alike1 sheur-ans dn) (not (mplusp n)))) 836 ((not (alike1 sheur-ans (m// n dn))) 837 (return (simplimit (m// ($expand (m// n sheur-ans)) 838 ($expand (m// dn sheur-ans))) 839 var 840 val)))) 841 (cond ((and (not (and (eq val '$inf) (expp n) (expp dn))) 842 (setq n1 (try-lhospital-quit n dn nil)) 843 (not (eq n1 '$und))) 844 (return n1))) 845 (throw 'limit t))) 846 847;; Test whether both n and dn have form 848;; product of poly^poly 849(defun expfactorp (n dn) 850 (do ((llist (append (cond ((mtimesp n) (cdr n)) 851 (t (ncons n))) 852 (cond ((mtimesp dn) (cdr dn)) 853 (t (ncons dn)))) 854 (cdr llist)) 855 (exp? t) ;IS EVERY ELEMENT SO FAR 856 (factor nil)) ;A POLY^POLY? 857 ((or (null llist) 858 (not exp?)) 859 exp?) 860 (setq factor (car llist)) 861 (setq exp? (or (polyinx factor var ()) 862 (and (mexptp factor) 863 (polyinx (cadr factor) var ()) 864 (polyinx (caddr factor) var ())))))) 865 866(defun expfactor (n dn var) ;Attempts to evaluate limit by grouping 867 (prog (highest-deg) ; terms with similar exponents. 868 (let ((new-exp (exppoly n))) ;exppoly unrats expon 869 (setq n (car new-exp) ;and rtns deg of expons 870 highest-deg (cdr new-exp))) 871 (cond ((null n) (return nil))) ;nil means expon is not 872 (let ((new-exp (exppoly dn))) ;a rat func. 873 (setq dn (car new-exp) 874 highest-deg (max highest-deg (cdr new-exp)))) 875 (cond ((or (null dn) 876 (= highest-deg 0)) ; prevent infinite recursion 877 (return nil))) 878 (return 879 (do ((answer 1) 880 (degree highest-deg (1- degree)) 881 (numerator n) 882 (denominator dn) 883 (numfactors nil) 884 (denfactors nil)) 885 ((= degree -1) 886 (m* answer 887 (limit (m// numerator denominator) 888 var 889 val 890 'think))) 891 (let ((newnumer-factor (get-newexp&factors 892 numerator 893 degree 894 var))) 895 (setq numerator (car newnumer-factor) 896 numfactors (cdr newnumer-factor))) 897 (let ((newdenom-factor (get-newexp&factors 898 denominator 899 degree 900 var))) 901 (setq denominator (car newdenom-factor) 902 denfactors (cdr newdenom-factor))) 903 (setq answer (simplimit (list '(mexpt) 904 (m* answer 905 (m// numfactors denfactors)) 906 (cond ((> degree 0) var) 907 (t 1))) 908 var 909 val)) 910 (cond ((member answer '($ind $und) :test #'equal) 911 ;; cannot handle limit(exp(x*%i)*x, x, inf); 912 (return nil)) 913 ((member answer '($inf $minf) :test #'equal) 914 ;; 0, zeroa, zerob are passed through to next iteration 915 (return (simplimtimes (list (m// numerator denominator) answer))))))))) 916 917(defun exppoly (exp) ;RETURNS EXPRESSION WITH UNRATTED EXPONENTS 918 (do ((factor nil) 919 (highest-deg 0) 920 (new-exp 1) 921 (exp (cond ((mtimesp exp) 922 (cdr exp)) 923 (t (ncons exp))) 924 (cdr exp))) 925 ((null exp) (cons new-exp highest-deg)) 926 (setq factor (car exp)) 927 (setq new-exp 928 (m* (cond ((or (not (mexptp factor)) 929 (not (ratp (caddr factor) var))) 930 factor) 931 (t (setq highest-deg 932 (max highest-deg 933 (ratdegree (caddr factor)))) 934 (m^ (cadr factor) (unrat (caddr factor))))) 935 new-exp)))) 936 937(defun unrat (exp) ;RETURNS UNRATTED EXPRESION 938 (multiple-value-bind (n d) 939 (numden* exp) 940 (let ((tem ($divide n d))) 941 (m+ (cadr tem) 942 (m// (caddr tem) d))))) 943 944(defun get-newexp&factors (exp degree var) ;RETURNS (CONS NEWEXP FACTORS) 945 (do ((terms (cond ((mtimesp exp)(cdr exp)) ; SUCH THAT 946 (t (ncons exp))) ; NEWEXP*FACTORS^(VAR^DEGREE) 947 (cdr terms)) ; IS EQUAL TO EXP. 948 (factors 1) 949 (newexp 1) 950 (factor nil)) 951 ((null terms) 952 (cons newexp 953 factors)) 954 (setq factor (car terms)) 955 (cond ((not (mexptp factor)) 956 (cond ((= degree 0) 957 (setq factors (m* factor factors))) 958 (t (setq newexp (m* factor newexp))))) 959 ((or (= degree -1) 960 (= (ratdegree (caddr factor)) 961 degree)) 962 (setq factors (m* (m^ (cadr factor) 963 (leading-coef (caddr factor))) 964 factors) 965 newexp (m* (m^ (cadr factor) 966 (m- (caddr factor) 967 (m* (leading-coef (caddr factor)) 968 (m^ var degree)))) 969 newexp))) 970 (t (setq newexp (m* factor newexp)))))) 971 972(defun leading-coef (rat) 973 (ratlim (m// rat (m^ var (ratdegree rat))))) 974 975(defun ratdegree (rat) 976 (multiple-value-bind (n d) 977 (numden* rat) 978 (- (deg n) (deg d)))) 979 980(defun limfact2 (n d var val) 981 (let ((n1 (reflect0 n var val)) 982 (d1 (reflect0 d var val))) 983 (cond ((and (alike1 n n1) 984 (alike1 d d1)) 985 nil) 986 (t (limit (m// n1 d1) var val 'think))))) 987 988;; takes expression and returns operator at front with all flags removed 989;; except array flag. 990;; array flag must match for alike1 to consider two things to be the same. 991;; ((MTIMES SIMP) ... ) => (MTIMES) 992;; ((PSI SIMP ARRAY) 0) => (PSI ARRAY) 993(defun operator-with-array-flag (exp) 994 (cond ((member 'array (car exp) :test #'eq) 995 (list (caar exp) 'array)) 996 (t (list (caar exp))))) 997 998(defun reflect0 (exp var val) 999 (cond ((atom exp) exp) 1000 ((and (eq (caar exp) 'mfactorial) 1001 (let ((argval (limit (cadr exp) var val 'think))) 1002 (or (eq argval '$minf) 1003 (and (numberp argval) 1004 (> 0 argval))))) 1005 (reflect (cadr exp))) 1006 (t (cons (operator-with-array-flag exp) 1007 (mapcar (function 1008 (lambda (term) 1009 (reflect0 term var val))) 1010 (cdr exp)))))) 1011 1012(defun reflect (arg) 1013 (m* -1 1014 '$%pi 1015 (m^ (list (ncons 'mfactorial) 1016 (m+ -1 1017 (m* -1 arg))) 1018 -1) 1019 (m^ (list (ncons '%sin) 1020 (m* '$%pi arg)) 1021 -1))) 1022 1023(defun limfact (n d) 1024 (let ((ans ())) 1025 (setq n (stirling0 n) 1026 d (stirling0 d)) 1027 (setq ans (toplevel-$limit (m// n d) var '$inf)) 1028 (cond ((and (atom ans) 1029 (not (member ans '(und ind ) :test #'eq))) ans) 1030 ((eq (caar ans) '%limit) ()) 1031 (t ans)))) 1032 1033;; substitute asymptotic approximations for gamma, factorial, and 1034;; polylogarithm 1035(defun stirling0 (e) 1036 (cond ((atom e) e) 1037 ((and (setq e (cons (car e) (mapcar 'stirling0 (cdr e)))) 1038 nil)) 1039 ((and (eq (caar e) '%gamma) 1040 (eq (limit (cadr e) var val 'think) '$inf)) 1041 (stirling (cadr e))) 1042 ((and (eq (caar e) 'mfactorial) 1043 (eq (limit (cadr e) var val 'think) '$inf)) 1044 (m* (cadr e) (stirling (cadr e)))) 1045 ((and (eq (caar e) 'mqapply) ;; polylogarithm 1046 (eq (subfunname e) '$li) 1047 (integerp (car (subfunsubs e)))) 1048 (li-asymptotic-expansion (m- (car (subfunsubs e)) 1) 1049 (car (subfunsubs e)) 1050 (car (subfunargs e)))) 1051 (t e))) 1052 1053(defun stirling (x) 1054 (maxima-substitute x '$z 1055 '((mtimes simp) 1056 ((mexpt simp) 2 ((rat simp) 1 2)) 1057 ((mexpt simp) $%pi ((rat simp) 1 2)) 1058 ((mexpt simp) $z ((mplus simp) ((rat simp) -1 2) $z)) 1059 ((mexpt simp) $%e ((mtimes simp) -1 $z))))) 1060 1061(defun no-err-sub (v e &aux ans) 1062 (let ((errorsw t) (*zexptsimp? t) 1063 (errcatch t) 1064 ;; Don't print any error messages 1065 ($errormsg nil)) 1066 (declare (special errcatch)) 1067 ;; Should we just use IGNORE-ERRORS instead HANDLER-CASE here? I 1068 ;; (rtoy) am choosing the latter so that unexpected errors will 1069 ;; actually show up instead of being silently discarded. 1070 (handler-case 1071 (setq ans (catch 'errorsw 1072 (ignore-rat-err 1073 (sratsimp (subin v e))))) 1074 (maxima-$error () 1075 (setq ans nil))) 1076 (cond ((null ans) t) ; Ratfun package returns NIL for failure. 1077 (t ans)))) 1078 1079;; substitute value v for var into expression e. 1080;; if result is defined and e is continuous, we have the limit. 1081(defun simplimsubst (v e) 1082 (let (ans) 1083 (cond ((involve e '(mfactorial)) nil) 1084 1085 ;; functions that are defined at their discontinuities 1086 ((amongl '($atan2 $floor %round $ceiling %signum %integrate 1087 %gamma_incomplete) 1088 e) nil) 1089 1090 ;; substitute value into expression 1091 ((eq (setq ans (no-err-sub (ridofab v) e)) t) 1092 nil) 1093 1094 ((and (member v '($zeroa $zerob) :test #'eq) (=0 ($radcan ans))) 1095 (setq ans (behavior e var v)) 1096 (cond ((equal ans 1) '$zeroa) 1097 ((equal ans -1) '$zerob) 1098 (t nil))) ; behavior can't find direction 1099 (t ans)))) 1100 1101;;;returns (cons numerator denominator) 1102(defun numden* (e) 1103 (let ((e (factor (simplify e))) 1104 (numer ()) 1105 (denom ())) 1106 (cond ((atom e) 1107 (push e numer)) 1108 ((mtimesp e) 1109 (mapc #'forq (cdr e))) 1110 (t 1111 (forq e))) 1112 (cond ((null numer) 1113 (setq numer 1)) 1114 ((null (cdr numer)) 1115 (setq numer (car numer))) 1116 (t 1117 (setq numer (m*l numer)))) 1118 (cond ((null denom) 1119 (setq denom 1)) 1120 ((null (cdr denom)) 1121 (setq denom (car denom))) 1122 (t 1123 (setq denom (m*l denom)))) 1124 (values (factor numer) (factor denom)))) 1125 1126;;;FACTOR OR QUOTIENT 1127;;;Setq's the special vars numer and denom from numden* 1128(defun forq (e) 1129 (cond ((and (mexptp e) 1130 (not (freeof var e)) 1131 (null (pos-neg-p (caddr e)))) 1132 (push (m^ (cadr e) (m* -1. (caddr e))) denom)) 1133 (t (push e numer)))) 1134 1135;;;Predicate to tell whether an expression is pos,zero or neg as var -> val. 1136;;;returns T if pos,zero. () if negative or don't know. 1137(defun pos-neg-p (exp) 1138 (let ((ans (limit exp var val 'think))) 1139 (cond ((and (not (member ans '($und $ind $infinity) :test #'eq)) 1140 (equal ($imagpart ans) 0)) 1141 (let ((sign (getsignl ans))) 1142 (cond ((or (equal sign 1) 1143 (equal sign 0)) 1144 t) 1145 ((equal sign -1) nil)))) 1146 (t 'unknown)))) 1147 1148(declare-top (unspecial n dn)) 1149 1150(defun expp (e) 1151 (cond ((radicalp e var) nil) 1152 ((member (caar e) '(%log %sin %cos %tan %sinh %cosh %tanh mfactorial 1153 %asin %acos %atan %asinh %acosh %atanh) :test #'eq) nil) 1154 ((simplexp e) t) 1155 ((do ((e (cdr e) (cdr e))) 1156 ((null e) nil) 1157 (and (expp (car e)) (return t)))))) 1158 1159(defun simplexp (e) 1160 (and (mexptp e) 1161 (radicalp (cadr e) var) 1162 (among var (caddr e)) 1163 (radicalp (caddr e) var))) 1164 1165 1166(defun gcpower (a b) 1167 ($gcd (getexp a) (getexp b))) 1168 1169(defun getexp (exp) 1170 (cond ((and (mexptp exp) 1171 (free (caddr exp) var) 1172 (eq (ask-integer (caddr exp) '$integer) '$yes)) 1173 (caddr exp)) 1174 ((mtimesp exp) (getexplist (cdr exp))) 1175 (t 1))) 1176 1177(defun getexplist (list) 1178 (cond ((null (cdr list)) 1179 (getexp (car list))) 1180 (t ($gcd (getexp (car list)) 1181 (getexplist (cdr list)))))) 1182 1183(defun limroot (exp power) 1184 (cond ((or (atom exp) (not (member (caar exp) '(mtimes mexpt) :test #'eq))) 1185 (limroot (list '(mexpt) exp 1) power)) ;This is strange-JIM. 1186 ((mexptp exp) (m^ (cadr exp) 1187 (sratsimp (m* (caddr exp) (m^ power -1.))))) 1188 (t (m*l (mapcar #'(lambda (x) 1189 (limroot x power)) 1190 (cdr exp)))))) 1191 1192;;NUMERATOR AND DENOMINATOR HAVE EXPONENTS WITH GCD OF GCP. 1193;;; Used to call simplimit but some of the transformations used here 1194;;; were not stable w.r.t. the simplifier, so try keeping exponent separate 1195;;; from bas. 1196 1197(defun colexpt (n dn gcp) 1198 (let ((bas (m* (limroot n gcp) (limroot dn (m* -1 gcp)))) 1199 (expo gcp) 1200 baslim expolim) 1201 (setq baslim (limit bas var val 'think)) 1202 (setq expolim (limit expo var val 'think)) 1203 (simplimexpt bas expo baslim expolim))) 1204 1205;;; This function will transform an expression such that either all logarithms 1206;;; contain arguments not becoming infinite or are of the form 1207;;; LOG(LOG( ... LOG(VAR))) This reduction takes place only over the operators 1208;;; MPLUS, MTIMES, MEXPT, and %LOG. 1209 1210(defun log-red-contract (facs) 1211 (do ((l facs (cdr l)) 1212 (consts ()) 1213 (log ())) 1214 ((null l) 1215 (if log (cons (cadr log) (m*l consts)) 1216 ())) 1217 (cond ((freeof var (car l)) (push (car l) consts)) 1218 ((mlogp (car l)) 1219 (if (null log) (setq log (car l)) 1220 (return ()))) 1221 (t (return ()))))) 1222 1223(defun log-reduce (x) 1224 (cond ((atom x) x) 1225 ((freeof var x) x) 1226 ((mplusp x) 1227 (do ((l (cdr x) (cdr l)) 1228 (sum ()) 1229 (weak-logs ()) 1230 (strong-logs ()) 1231 (temp)) 1232 ((null l) (m+l `(((%log) ,(m*l strong-logs)) 1233 ((%log) ,(m*l weak-logs)) 1234 ,@sum))) 1235 (setq x (log-reduce (car l))) 1236 (cond ((mlogp x) 1237 (if (infinityp (limit (cadr x) var val 'think)) 1238 (push (cadr x) strong-logs) 1239 (push (cadr x) weak-logs))) 1240 ((and (mtimesp x) (setq temp (log-red-contract (cdr x)))) 1241 (if (infinityp (limit (car temp) var val 'think)) 1242 (push (m^ (car temp) (cdr temp)) strong-logs) 1243 (push (m^ (car temp) (cdr temp)) weak-logs))) 1244 (t (push x sum))))) 1245 ((mtimesp x) 1246 (do ((l (cdr x) (cdr l)) 1247 (ans 1)) 1248 ((null l) ans) 1249 (setq ans ($expand (m* (log-reduce (car l)) ans))))) 1250 ((mexptp x) (m^t (log-reduce (cadr x)) (caddr x))) 1251 ((mlogp x) 1252 (cond ((not (infinityp (limit (cadr x) var val 'think))) x) 1253 (t 1254 (cond ((eq (cadr x) var) x) 1255 ((mplusp (cadr x)) 1256 (let ((strongl (maxi (cdadr x)))) 1257 (m+ (log-reduce `((%log) ,(car strongl))) `((%log) ,(m// (cadr x) (car strongl)))))) 1258 ((mtimesp (cadr x)) 1259 (do ((l (cdadr x) (cdr l)) (ans 0)) ((null l) ans) 1260 (setq ans (m+ (log-reduce (simplify `((%log) ,(log-reduce (car l))))) ans)))) 1261 (t 1262 (let ((red-log (simplify `((%log) ,(log-reduce (cadr x)))))) 1263 (if (alike1 red-log x) x (log-reduce red-log)))))))) 1264 (t x))) 1265 1266;; this function is responsible for the following bug: 1267;; limit(x^2 + %i*x, x, inf) -> inf (should be infinity) 1268(defun ratlim (e) 1269 (cond ((member val '($inf $infinity) :test #'eq) 1270 (setq e (maxima-substitute (m^t 'x -1) var e))) 1271 ((eq val '$minf) 1272 (setq e (maxima-substitute (m^t -1 (m^t 'x -1)) var e))) 1273 ((eq val '$zerob) 1274 (setq e (maxima-substitute (m- 'x) var e))) 1275 ((eq val '$zeroa) 1276 (setq e (maxima-substitute 'x var e))) 1277 ((setq e (maxima-substitute (m+t 'x val) var e)))) 1278 (destructuring-let* ((e (let (($ratfac ())) 1279 ($rat (sratsimp e) 'x))) 1280 ((h n . d) e) 1281 (g (genfind h 'x)) 1282 (nd (lodeg n g)) 1283 (dd (lodeg d g))) 1284 (cond ((and (setq e 1285 (subst var 1286 'x 1287 (sratsimp (m// ($ratdisrep `(,h ,(locoef n g) . 1)) 1288 ($ratdisrep `(,h ,(locoef d g) . 1)))))) 1289 (> nd dd)) 1290 (cond ((not (member val '($zerob $zeroa $inf $minf) :test #'eq)) 1291 0) 1292 ((not (equal ($imagpart e) 0)) 1293 0) 1294 ((null (setq e (getsignl ($realpart e)))) 1295 0) 1296 ((equal e 1) '$zeroa) 1297 ((equal e -1) '$zerob) 1298 (t 0))) 1299 ((equal nd dd) e) 1300 ((not (member val '($zerob $zeroa $infinity $inf $minf) :test #'eq)) 1301 (throw 'limit t)) 1302 ((eq val '$infinity) '$infinity) 1303 ((not (equal ($imagpart e) 0)) '$infinity) 1304 ((null (setq e (getsignl ($realpart e)))) 1305 (throw 'limit t)) 1306 ((equal e 1) '$inf) 1307 ((equal e -1) '$minf) 1308 (t 0)))) 1309 1310(defun lodeg (n x) 1311 (if (or (atom n) (not (eq (car n) x))) 1312 0 1313 (lowdeg (cdr n)))) 1314 1315(defun locoef (n x) 1316 (if (or (atom n) (not (eq (car n) x))) 1317 n 1318 (car (last n)))) 1319 1320(defun behavior (exp var val) ; returns either -1, 0, 1. 1321 (if (= *behavior-count-now* +behavior-count+) 1322 0 1323 (let ((*behavior-count-now* (1+ *behavior-count-now*)) pair sign) 1324 (cond ((real-infinityp val) 1325 (setq val (cond ((eq val '$inf) '$zeroa) 1326 ((eq val '$minf) '$zerob))) 1327 (setq exp (sratsimp (subin (m^ var -1) exp))))) 1328 (cond ((eq val '$infinity) 0) ; Needs more hacking for complex. 1329 ((and (mtimesp exp) 1330 (prog2 (setq pair (partition exp var 1)) 1331 (not (mtimesp (cdr pair))))) 1332 (setq sign (getsignl (car pair))) 1333 (if (not (fixnump sign)) 1334 0 1335 (mul sign (behavior (cdr pair) var val)))) 1336 ((and (=0 (no-err-sub (ridofab val) exp)) 1337 (mexptp exp) 1338 (free (caddr exp) var) 1339 (equal (getsignl (caddr exp)) 1)) 1340 (let ((bas (cadr exp)) (expo (caddr exp))) 1341 (behavior-expt bas expo))) 1342 (t (behavior-by-diff exp var val)))))) 1343 1344(defun behavior-expt (bas expo) 1345 (let ((behavior (behavior bas var val))) 1346 (cond ((= behavior 1) 1) 1347 ((= behavior 0) 0) 1348 ((eq (ask-integer expo '$integer) '$yes) 1349 (cond ((eq (ask-integer expo '$even) '$yes) 1) 1350 (t behavior))) 1351 ((ratnump expo) 1352 (cond ((evenp (cadr expo)) 1) 1353 ((oddp (caddr expo)) behavior) 1354 (t 0))) 1355 (t 0)))) 1356 1357(defun behavior-by-diff (exp var val) 1358 (cond ((not (or (eq val '$zeroa) (eq val '$zerob))) 0) 1359 (t (let ((old-val val) (old-exp exp)) 1360 (setq val (ridofab val)) 1361 (do ((ct 0 (1+ ct)) 1362 (exp (sratsimp (sdiff exp var)) (sratsimp (sdiff exp var))) 1363 (n () (not n)) 1364 (ans ())) ; This do wins by a return. 1365 ((> ct 0) 0) ; This loop used to run up to 5 times, 1366 ;; but the size of some expressions would blow up. 1367 (setq ans (no-err-sub val exp)) ;Why not do an EVENFN and ODDFN 1368 ;test here. 1369 (cond ((eq ans t) 1370 (return (behavior-numden old-exp var old-val))) 1371 ((=0 ans) ()) ;Do it again. 1372 (t (setq ans (getsignl ans)) 1373 (cond (n (return ans)) 1374 ((equal ans 1) 1375 (return (if (eq old-val '$zeroa) 1 -1))) 1376 ((equal ans -1) 1377 (return (if (eq old-val '$zeroa) -1 1))) 1378 (t (return 0)))))))))) 1379 1380(defun behavior-numden (exp var val) 1381 (let ((num ($num exp)) (denom ($denom exp))) 1382 (cond ((equal denom 1) 0) ;Could be hacked more from here. 1383 (t (let ((num-behav (behavior num var val)) 1384 (denom-behav (behavior denom var val))) 1385 (cond ((or (= num-behav 0) (= denom-behav 0)) 0) 1386 ((= num-behav denom-behav) 1) 1387 (t -1))))))) 1388 1389(defun try-lhospital (n d ind) 1390 ;;Make one catch for the whole bunch of lhospital trials. 1391 (let ((ans (lhospital-catch n d ind))) 1392 (cond ((null ans) ()) 1393 ((not (free-infp ans)) (simpinf ans)) 1394 ((not (free-epsilonp ans)) (simpab ans)) 1395 (t ans)))) 1396 1397(defun try-lhospital-quit (n d ind) 1398 (let ((ans (or (lhospital-catch n d ind) 1399 (lhospital-catch (m^ d -1) (m^ n -1) ind)))) 1400 (cond ((null ans) (throw 'limit t)) 1401 ((not (free-infp ans)) (simpinf ans)) 1402 ((not (free-epsilonp ans)) (simpab ans)) 1403 (t ans)))) 1404 1405(defun lhospital-catch (n d ind) 1406 (cond ((> 0 lhcount) 1407 (setq lhcount $lhospitallim) 1408 (throw 'lhospital nil)) 1409 ((equal lhcount $lhospitallim) 1410 (let ((lhcount (m+ lhcount -1))) 1411 (catch 'lhospital (lhospital n d ind)))) 1412 (t (setq lhcount (m+ lhcount -1)) 1413 (prog1 (lhospital n d ind) 1414 (setq lhcount (m+ lhcount 1)))))) 1415;;If this succeeds then raise LHCOUNT. 1416 1417 1418(defun lhospital (n d ind) 1419 (declare (special val lhp?)) 1420 (when (mtimesp n) 1421 (setq n (m*l (mapcar #'(lambda (term) (lhsimp term var val)) (cdr n))))) 1422 (when (mtimesp d) 1423 (setq d (m*l (mapcar #'(lambda (term) (lhsimp term var val)) (cdr d))))) 1424 (multiple-value-bind (n d) 1425 (lhop-numden n d) 1426 (let (const nconst dconst) 1427 (setq lhp? (and (null ind) (cons n d))) 1428 (multiple-value-setq (nconst n) (var-or-const n)) 1429 (multiple-value-setq (dconst d) (var-or-const d)) 1430 1431 (setq n (stirling0 n)) ;; replace factorial and %gamma 1432 (setq d (stirling0 d)) ;; with approximations 1433 1434 (setq n (sdiff n var) ;; take derivatives for l'hospital 1435 d (sdiff d var)) 1436 1437 (if (or (not (free n '%derivative)) (not (free d '%derivative))) 1438 (throw 'lhospital ())) 1439 (setq n (expand-trigs (tansc n) var)) 1440 (setq d (expand-trigs (tansc d) var)) 1441 1442 (multiple-value-setq (const n d) (remove-singularities n d)) 1443 (setq const (m* const (m// nconst dconst))) 1444 (simpinf (let ((ans (if ind 1445 (limit2 n d var val) 1446 (limit-numden n d val)))) 1447 ;; When the limit function returns, it's possible that it will return NIL 1448 ;; (gave up without finding a limit). It's also possible that it will 1449 ;; return something containing UND. We treat that as a failure too. 1450 (when (and ans (freeof '$und ans)) 1451 (m* const ans))))))) 1452 1453;; Try to compute the limit of a quotient NUM/DEN, trying to massage the input 1454;; into a convenient form for LIMIT on the way. 1455(defun limit-numden (n d val) 1456 (let ((expr (cond 1457 ;; For general arguments, the best approach seems to be to use 1458 ;; sratsimp to simplify the quotient as much as we can, then 1459 ;; $multthru, which splits it up into a sum (presumably useful 1460 ;; because limit(a+b) = limit(a) + limit(b) if the limits exist, and 1461 ;; the right hand side might be easier to calculate) 1462 ((not (mplusp n)) 1463 ($multthru (sratsimp (m// n d)))) 1464 1465 ;; If we've already got a sum in the numerator, it seems to be 1466 ;; better not to recombine it. Call LIMIT on the whole lot, though, 1467 ;; because terms with infinite limits might cancel to give a finite 1468 ;; result. 1469 (t 1470 (m+l (mapcar #'(lambda (x) 1471 (sratsimp (m// x d))) 1472 (cdr n))))))) 1473 1474 (limit expr var val 'think))) 1475 1476;; Heuristics for picking the right way to express a LHOSPITAL problem. 1477(defun lhop-numden (num denom) 1478 (declare (special var)) 1479 (cond ((let ((log-num (involve num '(%log)))) 1480 (cond ((null log-num) ()) 1481 ((lessthan (num-of-logs (factor (sratsimp (sdiff (m^ num -1) var)))) 1482 (num-of-logs (factor (sratsimp (sdiff num var))))) 1483 (psetq num (m^ denom -1) denom (m^ num -1)) 1484 t) 1485 (t t)))) 1486 ((let ((log-denom (involve denom '(%log)))) 1487 (cond ((null log-denom) ()) 1488 ((lessthan (num-of-logs (sratsimp (sdiff (m^ denom -1) var))) 1489 (num-of-logs (sratsimp (sdiff denom var)))) 1490 (psetq denom (m^ num -1) num (m^ denom -1)) 1491 t) 1492 (t t)))) 1493 ((let ((exp-num (%einvolve num))) 1494 (cond (exp-num 1495 (cond ((%e-right-placep exp-num) 1496 t) 1497 (t (psetq num (m^ denom -1) 1498 denom (m^ num -1)) t))) 1499 (t ())))) 1500 ((let ((exp-den (%einvolve denom))) 1501 (cond (exp-den 1502 (cond ((%e-right-placep exp-den) 1503 t) 1504 (t (psetq num (m^ denom -1) 1505 denom (m^ num -1)) t))) 1506 (t ())))) 1507 ((let ((scnum (involve num '(%sin)))) 1508 (cond (scnum (cond ((trig-right-placep '%sin scnum) t) 1509 (t (psetq num (m^ denom -1) 1510 denom (m^ num -1)) t))) 1511 (t ())))) 1512 ((let ((scden (involve denom '(%sin)))) 1513 (cond (scden (cond ((trig-right-placep '%sin scden) t) 1514 (t (psetq num (m^ denom -1) 1515 denom (m^ num -1)) t))) 1516 (t ())))) 1517 ((let ((scnum (involve num '(%asin %acos %atan)))) 1518 ;; If the numerator contains an inverse trig and the 1519 ;; denominator or reciprocal of denominator is polynomial, 1520 ;; leave everything as is. If the inverse trig is moved to 1521 ;; the denominator, things get messy, even if the numerator 1522 ;; becomes a polynomial. This is not perfect. 1523 (cond ((and scnum (or (polyinx denom var ()) 1524 (polyinx (m^ denom -1) var ()))) 1525 t) 1526 (t nil)))) 1527 ((or (oscip num) (oscip denom))) 1528 ((frac num) 1529 (psetq num (m^ denom -1) denom (m^ num -1)))) 1530 (values num denom)) 1531 1532;;i don't know what to do here for some cases, may have to be refined. 1533(defun num-of-logs (exp) 1534 (cond ((mapatom exp) 0) 1535 ((equal (caar exp) '%log) 1536 (m+ 1 (num-of-log-l (cdr exp)))) 1537 ((and (mexptp exp) (mnump (caddr exp))) 1538 (m* (simplify `((mabs) ,(caddr exp))) 1539 (num-of-logs (cadr exp)))) 1540 (t (num-of-log-l (cdr exp))))) 1541 1542(defun num-of-log-l (llist) 1543 (do ((temp llist (cdr temp)) (ans 0)) 1544 ((null temp) ans) 1545 (setq ans (m+ ans (num-of-logs (car temp)))))) 1546 1547(defun %e-right-placep (%e-arg) 1548 (let ((%e-arg-diff (sdiff %e-arg var))) 1549 (cond 1550 ((free %e-arg-diff var)) ;simple cases 1551 ((or (and (mexptp denom) 1552 (equal (cadr denom) -1)) 1553 (polyinx (m^ denom -1) var ())) ()) 1554 ((let ((%e-arg-diff-lim (ridofab (limit %e-arg-diff var val 'think))) 1555 (%e-arg-exp-lim (ridofab (limit (m^ '$%e %e-arg) var val 'think)))) 1556 #+nil 1557 (progn 1558 (format t "%e-arg-dif-lim = ~A~%" %e-arg-diff-lim) 1559 (format t "%e-arg-exp-lim = ~A~%" %e-arg-exp-lim)) 1560 (cond ((equal %e-arg-diff-lim %e-arg-exp-lim) 1561 t) 1562 ((and (mnump %e-arg-diff-lim) (mnump %e-arg-exp-lim)) 1563 t) 1564 ((and (mnump %e-arg-diff-lim) (infinityp %e-arg-exp-lim)) 1565 ;; This is meant to make maxima handle bug 1469411 1566 ;; correctly. Undoubtedly, this needs work. 1567 t) 1568 (t ()))))))) 1569 1570(defun trig-right-placep (trig-type arg) 1571 (let ((arglim (ridofab (limit arg var val 'think))) 1572 (triglim (ridofab (limit `((,trig-type) ,arg) var val 'think)))) 1573 (cond ((and (equal arglim 0) (equal triglim 0)) t) 1574 ((and (infinityp arglim) (infinityp triglim)) t) 1575 (t ())))) 1576 1577;;Takes a numerator and a denominator. If they tries all combinations of 1578;;products to try and make a simpler set of subproblems for LHOSPITAL. 1579(defun remove-singularities (numer denom) 1580 (cond ((or (null numer) (null denom) 1581 (atom numer) (atom denom) 1582 (not (mtimesp numer)) ;Leave this here for a while. 1583 (not (mtimesp denom))) 1584 (values 1 numer denom)) 1585 (t 1586 (let ((const 1)) 1587 (multiple-value-bind (num-consts num-vars) 1588 (var-or-const numer) 1589 (multiple-value-bind (denom-consts denom-vars) 1590 (var-or-const denom) 1591 (if (not (mtimesp num-vars)) 1592 (setq num-vars (list num-vars)) 1593 (setq num-vars (cdr num-vars))) 1594 (if (not (mtimesp denom-vars)) 1595 (setq denom-vars (list denom-vars)) 1596 (setq denom-vars (cdr denom-vars))) 1597 (do ((nl num-vars (cdr nl)) 1598 (num-list (copy-list num-vars )) 1599 (den-list denom-vars den-list-temp) 1600 (den-list-temp (copy-list denom-vars))) 1601 ((null nl) (values (m* const (m// num-consts denom-consts)) 1602 (m*l num-list) 1603 (m*l den-list-temp))) 1604 (do ((dl den-list (cdr dl))) 1605 ((null dl) t) 1606 (if (or (%einvolve (car nl)) (%einvolve (car nl))) 1607 t 1608 (let ((lim (catch 'limit (simpinf (simpab (limit (m// (car nl) (car dl)) 1609 var val 'think)))))) 1610 (cond ((or (eq lim t) 1611 (eq lim ()) 1612 (equal (ridofab lim) 0) 1613 (infinityp lim) 1614 (not (free lim '$inf)) 1615 (not (free lim '$minf)) 1616 (not (free lim '$infinity)) 1617 (not (free lim '$ind)) 1618 (not (free lim '$und))) 1619 ()) 1620 (t 1621 (setq const (m* lim const)) 1622 (setq num-list (delete (car nl) num-list :count 1 :test #'equal)) 1623 (setq den-list-temp (delete (car dl) den-list-temp :count 1 :test #'equal)) 1624 (return t))))))))))))) 1625 1626;; separate terms that contain var from constant terms 1627;; returns (const-terms . var-terms) 1628(defun var-or-const (expr) 1629 (setq expr ($factor expr)) 1630 (cond ((atom expr) 1631 (if (eq expr var) 1632 (values 1 expr) 1633 (values expr 1))) 1634 ((free expr var) 1635 (values expr 1)) 1636 ((mtimesp expr) 1637 (do ((l (cdr expr) (cdr l)) 1638 (const 1) 1639 (varl 1)) 1640 ((null l) (values const varl)) 1641 (if (free (car l) var) 1642 (setq const (m* (car l) const)) 1643 (setq varl (m* (car l) varl))))) 1644 (t 1645 (values 1 expr)))) 1646 1647;; if term goes to non-zero constant, replace with constant 1648(defun lhsimp (term var val) 1649 (cond ((atom term) term) 1650 (t 1651 (let ((term-value (ridofab (limit term var val 'think)))) 1652 (cond ((not (member term-value 1653 '($inf $minf $und $ind $infinity 0))) 1654 term-value) 1655 (t term)))))) 1656 1657(defun bylog (expo bas) 1658 (simplimexpt '$%e 1659 (setq bas 1660 (try-lhospital-quit (simplify `((%log) ,(tansc bas))) 1661 (m^ expo -1) 1662 nil)) 1663 '$%e bas)) 1664 1665(defun simplimexpt (bas expo bl el) 1666 (cond ((or (eq bl '$und) (eq el '$und)) '$und) 1667 ((zerop2 bl) 1668 (cond ((eq el '$inf) (if (eq bl '$zeroa) bl 0)) 1669 ((eq el '$minf) (if (eq bl '$zeroa) '$inf '$infinity)) 1670 ((eq el '$ind) '$ind) 1671 ((eq el '$infinity) '$und) 1672 ((zerop2 el) (bylog expo bas)) 1673 (t (cond ((equal (getsignl el) -1) 1674 (cond ((eq bl '$zeroa) '$inf) 1675 ((eq bl '$zerob) 1676 (cond ((even1 el) '$inf) 1677 ((eq (ask-integer el '$integer) '$yes) 1678 (if (eq (ask-integer el '$even) '$yes) 1679 '$inf 1680 '$minf)))) ;Gotta be ODD. 1681 (t (setq bas (behavior bas var val)) 1682 (cond ((equal bas 1) '$inf) 1683 ((equal bas -1) '$minf) 1684 (t (throw 'limit t)))))) 1685 ((and (mnump el) 1686 (member bl '($zeroa $zerob) :test #'eq)) 1687 (cond ((even1 el) '$zeroa) 1688 ((and (eq bl '$zerob) 1689 (ratnump el) 1690 (evenp (caddr el))) 0) 1691 (t bl))) 1692 ((and (equal (getsignl el) 1) 1693 (eq bl '$zeroa)) bl) 1694 ((equal (getsignl el) 0) 1695 1) 1696 (t 0))))) 1697 ((eq bl '$infinity) 1698 (cond ((zerop2 el) (bylog expo bas)) 1699 ((eq el '$minf) 0) 1700 ((eq el '$inf) '$infinity) 1701 ((member el '($infinity $ind) :test #'eq) '$und) 1702 ((equal (setq el (getsignl el)) 1) '$infinity) 1703 ((equal el 0) 1) 1704 ((equal el -1) 0) 1705 (t (throw 'limit t)))) 1706 ((eq bl '$inf) 1707 (cond ((eq el '$inf) '$inf) 1708 ((equal el '$minf) 0) 1709 ((zerop2 el) (bylog expo bas)) 1710 ((member el '($infinity $ind) :test #'eq) '$und) 1711 (t (cond ((eql 0 (getsignl el)) 1) 1712 ((ratgreaterp 0 el) '$zeroa) 1713 ((ratgreaterp el 0) '$inf) 1714 (t (throw 'limit t)))))) 1715 ((eq bl '$minf) 1716 (cond ((zerop2 el) (bylog expo bas)) 1717 ((eq el '$inf) '$und) 1718 ((equal el '$minf) 0) 1719;;;Why not generalize this. We can ask about the number. -Jim 2/23/81 1720 ((mnump el) (cond ((mnegp el) 1721 (if (even1 el) 1722 '$zeroa 1723 (if (eq (ask-integer el '$integer) '$yes) 1724 (if (eq (ask-integer el '$even) '$yes) 1725 '$zeroa 1726 '$zerob) 1727 0))) 1728 (t (cond ((even1 el) '$inf) 1729 ((eq (ask-integer el '$integer) '$yes) 1730 (if (eq (ask-integer el '$even) '$yes) 1731 '$inf 1732 '$minf)) 1733 (t '$infinity))))) 1734 (loginprod? (throw 'lip? 'lip!)) 1735 (t '$und))) 1736 ((equal (simplify (ratdisrep (ridofab bl))) 1) 1737 (if (infinityp el) (bylog expo bas) 1)) 1738 ((and (equal (ridofab bl) -1) 1739 (infinityp el)) '$ind) ;LB 1740 ((eq bl '$ind) (cond ((or (zerop2 el) (infinityp el)) '$und) 1741 ((not (equal (getsignl el) -1)) '$ind) 1742 (t '$und))) 1743 ((eq el '$inf) (cond ((abeq1 bl) 1744 (if (equal (getsignl bl) 1) 1 '$ind)) 1745 ((abless1 bl) 1746 (if (equal (getsignl bl) 1) '$zeroa 0)) 1747 ((equal (getsignl (m1- bl)) 1) '$inf) 1748 ((equal (getsignl (m1- `((mabs) ,bl))) 1) '$infinity) 1749 (t (throw 'limit t)))) 1750 ((eq el '$minf) (cond ((abeq1 bl) 1751 (if (equal (getsignl bl) 1) 1 '$ind)) 1752 ((not (abless1 bl)) 1753 (if (equal (getsignl bl) 1) '$zeroa 0)) 1754 ((ratgreaterp 0 bl) '$infinity) 1755 (t '$inf))) 1756 ((eq el '$infinity) 1757 (if (equal val '$infinity) 1758 '$und ;Not enough info to do anything. 1759 (destructuring-bind (real-el . imag-el) 1760 (trisplit expo) 1761 (setq real-el (limit real-el var origval nil)) 1762 (cond ((eq real-el '$minf) 1763 0) 1764 ((and (eq real-el '$inf) 1765 (not (equal (ridofab (limit imag-el var origval nil)) 0))) 1766 '$infinity) 1767 ((eq real-el '$infinity) 1768 (throw 'limit t)) ;; don't really know real component 1769 (t 1770 '$ind))))) 1771 1772 ((eq el '$ind) '$ind) 1773 ((zerop2 el) 1) 1774 (t (m^ bl el)))) 1775 1776(defun even1 (x) 1777 (cond ((numberp x) (and (integerp x) (evenp x))) 1778 ((and (mnump x) (evenp (cadr x)))))) 1779 1780;; is absolute value less than one? 1781(defun abless1 (bl) 1782 (setq bl (nmr bl)) 1783 (cond ((mnump bl) 1784 (and (ratgreaterp 1. bl) (ratgreaterp bl -1.))) 1785 (t (equal (getsignl (m1- `((mabs) ,bl))) -1.)))) 1786 1787;; is absolute value equal to one? 1788(defun abeq1 (bl) 1789 (setq bl (nmr bl)) 1790 (cond ((mnump bl) 1791 (or (equal 1. bl) (equal bl -1.))) 1792 (t (equal (getsignl (m1- `((mabs) ,bl))) 0)))) 1793 1794(defun simplimit (exp var val &aux op) 1795 (cond 1796 ((eq var exp) val) 1797 ((or (atom exp) (mnump exp)) exp) 1798 ((and (not (infinityp val)) 1799 (not (amongl '(%sin %cos %atanh %cosh %sinh %tanh mfactorial %log) 1800 exp)) 1801 (not (inf-typep exp)) 1802 (simplimsubst val exp))) 1803 ((eq (caar exp) '%limit) (throw 'limit t)) 1804 ((mplusp exp) (simplimplus exp)) 1805 ((mtimesp exp) (simplimtimes (cdr exp))) 1806 ((mexptp exp) (simplimexpt (cadr exp) (caddr exp) 1807 (limit (cadr exp) var val 'think) 1808 (limit (caddr exp) var val 'think))) 1809 ((mlogp exp) (simplimln exp var val)) 1810 ((member (caar exp) '(%sin %cos) :test #'eq) 1811 (simplimsc exp (caar exp) (limit (cadr exp) var val 'think))) 1812 ((eq (caar exp) '%tan) (simplim%tan (cadr exp))) 1813 ((eq (caar exp) '%atan) (simplim%atan (limit (cadr exp) var val 'think))) 1814 ((eq (caar exp) '$atan2) (simplim%atan2 exp)) 1815 ((member (caar exp) '(%sinh %cosh) :test #'eq) 1816 (simplimsch (caar exp) (limit (cadr exp) var val 'think))) 1817 ((eq (caar exp) 'mfactorial) 1818 (simplimfact exp var val)) 1819 ((member (caar exp) '(%erf %tanh) :test #'eq) 1820 (simplim%erf-%tanh (caar exp) (cadr exp))) 1821 ((member (caar exp) '(%acos %asin) :test #'eq) 1822 (simplim%asin-%acos (caar exp) (limit (cadr exp) var val 'think))) 1823 ((eq (caar exp) '%atanh) 1824 (simplim%atanh (limit (cadr exp) var val 'think) val)) 1825 ((eq (caar exp) '%acosh) 1826 (simplim%acosh (limit (cadr exp) var val 'think))) 1827 ((eq (caar exp) '%asinh) 1828 (simplim%asinh (limit (cadr exp) var val 'think))) 1829 ((eq (caar exp) '%inverse_jacobi_ns) 1830 (simplim%inverse_jacobi_ns (limit (cadr exp) var val 'think) (third exp))) 1831 ((eq (caar exp) '%inverse_jacobi_nc) 1832 (simplim%inverse_jacobi_nc (limit (cadr exp) var val 'think) (third exp))) 1833 ((eq (caar exp) '%inverse_jacobi_sc) 1834 (simplim%inverse_jacobi_sc (limit (cadr exp) var val 'think) (third exp))) 1835 ((eq (caar exp) '%inverse_jacobi_cs) 1836 (simplim%inverse_jacobi_cs (limit (cadr exp) var val 'think) (third exp))) 1837 ((eq (caar exp) '%inverse_jacobi_dc) 1838 (simplim%inverse_jacobi_dc (limit (cadr exp) var val 'think) (third exp))) 1839 ((eq (caar exp) '%inverse_jacobi_ds) 1840 (simplim%inverse_jacobi_ds (limit (cadr exp) var val 'think) (third exp))) 1841 ((and (eq (caar exp) 'mqapply) 1842 (eq (subfunname exp) '$li)) 1843 (simplim$li (subfunsubs exp) (subfunargs exp) val)) 1844 ((and (eq (caar exp) 'mqapply) 1845 (eq (subfunname exp) '$psi)) 1846 (simplim$psi (subfunsubs exp) (subfunargs exp) val)) 1847 ((and (eq (caar exp) var) 1848 (member 'array (car exp) :test #'eq) 1849 (every #'(lambda (sub-exp) 1850 (free sub-exp var)) 1851 (cdr exp))) 1852 exp) ;LIMIT(B[I],B,INF); -> B[I] 1853 ((setq op (safe-get (mop exp) 'simplim%function)) 1854 ;; Lookup a simplim%function from the property list 1855 (funcall op exp var val)) 1856 (t (if $limsubst 1857 (simplify (cons (operator-with-array-flag exp) 1858 (mapcar #'(lambda (a) 1859 (limit a var val 'think)) 1860 (cdr exp)))) 1861 (throw 'limit t))))) 1862 1863(defun liminv (e) 1864 (setq e (resimplify (subst (m// 1 var) var e))) 1865 (let ((new-val (cond ((eq val '$zeroa) '$inf) 1866 ((eq val '$zerob) '$minf)))) 1867 (if new-val (let ((preserve-direction t)) 1868 (toplevel-$limit e var new-val)) (throw 'limit t)))) 1869 1870(defun simplimtimes (exp) 1871 ;; The following test 1872 ;; handles (-1)^x * 2^x => (-2)^x => $infinity 1873 ;; wants to avoid (-1)^x * 2^x => $ind * $inf => $und 1874 (let ((try 1875 (and (expfactorp (cons '(mtimes) exp) 1) 1876 (expfactor (cons '(mtimes) exp) 1 var)))) 1877 (when try (return-from simplimtimes try))) 1878 1879 (let ((prod 1) (num 1) (denom 1) 1880 (zf nil) (ind-flag nil) (inf-type nil) 1881 (constant-zero nil) (constant-infty nil)) 1882 (dolist (term exp) 1883 (let* ((loginprod? (involve term '(%log))) 1884 (y (catch 'lip? (limit term var val 'think)))) 1885 (cond 1886 ;; limit failed due to log in product 1887 ((eq y 'lip!) 1888 (return-from simplimtimes (liminv (cons '(mtimes simp) exp)))) 1889 1890 ;; If the limit is infinitesimal or zero 1891 ((zerop2 y) 1892 (setf num (m* num term) 1893 constant-zero (or constant-zero (not (among var term)))) 1894 (case y 1895 ($zeroa 1896 (unless zf (setf zf 1))) 1897 ($zerob 1898 (setf zf (* -1 (or zf 1)))))) 1899 1900 ;; If the limit is not some form of infinity or 1901 ;; undefined/indeterminate. 1902 ((not (member y '($inf $minf $infinity $ind $und) :test #'eq)) 1903 (setq prod (m* prod y))) 1904 1905 ((eq y '$und) (return-from simplimtimes '$und)) 1906 ((eq y '$ind) (setq ind-flag t)) 1907 1908 ;; Some form of infinity 1909 (t 1910 (setf denom (m* denom term) 1911 constant-infty (or constant-infty (not (among var term)))) 1912 (unless (eq inf-type '$infinity) 1913 (cond 1914 ((eq y '$infinity) (setq inf-type '$infinity)) 1915 ((null inf-type) (setf inf-type y)) 1916 ;; minf * minf or inf * inf 1917 ((eq y inf-type) (setf inf-type '$inf)) 1918 ;; minf * inf 1919 (t (setf inf-type '$minf)))))))) 1920 1921 (cond 1922 ;; If there are zeros and infinities among the terms that are free of 1923 ;; VAR, then we have an expression like "inf * zeroa * f(x)" or 1924 ;; similar. That gives an undefined result. Note that we don't 1925 ;; necessarily have something undefined if only the zeros have a term 1926 ;; free of VAR. For example "zeroa * exp(-1/x) * 1/x" as x -> 0. And 1927 ;; similarly for the infinities. 1928 ((and constant-zero constant-infty) '$und) 1929 1930 ;; If num=denom=1, we didn't find any explicit infinities or zeros, so we 1931 ;; either return the simplified product or ind 1932 ((and (eql num 1) (eql denom 1)) 1933 (if ind-flag '$ind prod)) 1934 ;; If denom=1 (and so num != 1), we have some form of zero 1935 ((equal denom 1) 1936 (if (null zf) 1937 0 1938 (let ((sign (getsignl prod))) 1939 (if (or (not sign) (eq sign 'complex)) 1940 0 1941 (ecase (* zf sign) 1942 (0 0) 1943 (1 '$zeroa) 1944 (-1 '$zerob)))))) 1945 ;; If num=1 (and so denom != 1), we have some form of infinity 1946 ((equal num 1) 1947 (let ((sign ($csign prod))) 1948 (cond 1949 (ind-flag '$und) 1950 ((eq sign '$pos) inf-type) 1951 ((eq sign '$neg) (case inf-type 1952 ($inf '$minf) 1953 ($minf '$inf) 1954 (t '$infinity))) 1955 ((member sign '($complex $imaginary)) '$infinity) 1956 ; sign is '$zero, $pnz, $pz, etc 1957 (t (throw 'limit t))))) 1958 ;; Both zeros and infinities 1959 (t 1960 ;; All bets off if there are some infinities or some zeros, but it 1961 ;; needn't be undefined (see above) 1962 (when (or constant-zero constant-infty) (throw 'limit t)) 1963 1964 (let ((ans (limit2 num (m^ denom -1) var val))) 1965 (if ans 1966 (simplimtimes (list prod ans)) 1967 (throw 'limit t))))))) 1968 1969;;;PUT CODE HERE TO ELIMINATE FAKE SINGULARITIES?? 1970 1971(defun simplimplus (exp) 1972 (cond ((memalike exp simplimplus-problems) 1973 (throw 'limit t)) 1974 (t (unwind-protect 1975 (progn (push exp simplimplus-problems) 1976 (let ((ans (catch 'limit (simplimplus1 exp)))) 1977 (cond ((or (eq ans ()) 1978 (eq ans t) 1979 (among '%limit ans)) 1980 (let ((new-exp (sratsimp exp))) 1981 (cond ((not (alike1 exp new-exp)) 1982 (setq ans 1983 (limit new-exp var val 'think)))) 1984 (cond ((or (eq ans ()) 1985 (eq ans t)) 1986 (throw 'limit t)) 1987 (t ans)))) 1988 (t ans)))) 1989 (pop simplimplus-problems))))) 1990 1991(defun simplimplus1 (exp) 1992 (prog (sum y infl infinityl minfl indl) 1993 (setq sum 0.) 1994 (do ((exp (cdr exp) (cdr exp)) (f)) 1995 ((or y (null exp)) nil) 1996 (setq f (limit (car exp) var val 'think)) 1997 (cond ((null f) 1998 (throw 'limit t)) 1999 ((eq f '$und) (setq y t)) 2000 ((not (member f '($inf $minf $infinity $ind) :test #'eq)) 2001 (setq sum (m+ sum f))) 2002 ((eq f '$ind) (push (car exp) indl)) 2003 (infinityl (throw 'limit t)) 2004;;;Don't know what to do with an '$infinity and an $inf or $minf 2005 ((eq f '$inf) (push (car exp) infl)) 2006 ((eq f '$minf) (push (car exp) minfl)) 2007 ((eq f '$infinity) 2008 (cond ((or infl minfl) 2009 (throw 'limit t)) 2010 (t (push (car exp) infinityl)))))) 2011 (cond ((not (or infl minfl indl infinityl)) 2012 (return (cond ((atom sum) sum) 2013 ((or (not (free sum '$zeroa)) 2014 (not (free sum '$zerob))) 2015 (simpab sum)) 2016 (t sum)))) 2017 (t (cond ((null infinityl) 2018 (cond (infl (cond ((null minfl) (return '$inf)) 2019 (t (go oon)))) 2020 (minfl (return '$minf)) 2021 ((> (length indl) 1) 2022 ;; At this point we have a sum of '$ind. We factor 2023 ;; the sum and try again. This way we get the limit 2024 ;; of expressions like (a-b)*ind, where (a-b)--> 0. 2025 (cond ((not (alike1 (setq y ($factorsum exp)) exp)) 2026 (return (limit y var val 'think))) 2027 (t 2028 (return '$ind)))) 2029 (t (return '$ind)))) 2030 (t (setq infl (append infl infinityl)))))) 2031 2032 oon (setq y (m+l (append minfl infl))) 2033 (cond ((alike1 exp (setq y (sratsimp (log-reduce (hyperex y))))) 2034 (cond ((not (infinityp val)) 2035 (setq infl (cnv infl val)) ;THIS IS HORRIBLE!!!! 2036 (setq minfl (cnv minfl val)))) 2037 (let ((val '$inf)) 2038 (cond ((every #'(lambda (j) (radicalp j var)) 2039 (append infl minfl)) 2040 (setq y (rheur infl minfl))) 2041 (t (setq y (sheur infl minfl)))))) 2042 (t (setq y (limit y var val 'think)))) 2043 (cond ((or (eq y ()) 2044 (eq y t)) (return ())) 2045 ((infinityp y) (return y)) 2046 (t (return (m+ sum y)))))) 2047 2048;; Limit n/d, using heuristics on the order of growth. 2049(defun sheur0 (n d) 2050 (let ((orig-n n)) 2051 (cond ((and (free n var) 2052 (free d var)) 2053 (m// n d)) 2054 (t (setq n (cpa n d nil)) 2055 (cond ((equal n 1) 2056 (cond ((oscip orig-n) '$und) 2057 (t '$inf))) 2058 ((equal n -1) '$zeroa) 2059 ((equal n 0) (m// orig-n d))))))) 2060 2061 2062;;;L1 is a list of INF's and L2 is a list of MINF's. Added together 2063;;;it is indeterminate. 2064(defun sheur (l1 l2) 2065 (let ((term (sheur1 l1 l2))) 2066 (cond ((equal term '$inf) '$inf) 2067 ((equal term '$minf) '$minf) 2068 (t (let ((new-num (m+l (mapcar #'(lambda (num-term) 2069 (m// num-term (car l1))) 2070 (append l1 l2))))) 2071 (cond ((limit2 new-num (m// 1 (car l1)) var val)))))))) 2072 2073(defun frac (exp) 2074 (cond ((atom exp) nil) 2075 (t (setq exp (nformat exp)) 2076 (cond ((and (eq (caar exp) 'mquotient) 2077 (among var (caddr exp))) 2078 t))))) 2079 2080(defun zerop2 (z) (=0 (ridofab z))) 2081 2082(defun raise (a) (m+ a '$zeroa)) 2083 2084(defun lower (a) (m+ a '$zerob)) 2085 2086(defun sincoshk (exp1 l sc) 2087 (cond ((equal l 1) (lower l)) 2088 ((equal l -1) (raise l)) 2089 ((among sc l) l) 2090 ((member val '($zeroa $zerob) :test #'eq) (spangside exp1 l)) 2091 (t l))) 2092 2093(defun spangside (e l) 2094 (setq e (behavior e var val)) 2095 (cond ((equal e 1) (raise l)) 2096 ((equal e -1) (lower l)) 2097 (t l))) 2098 2099;; get rid of zeroa and zerob 2100(defun ridofab (e) 2101 (if (among '$zeroa e) (setq e (maxima-substitute 0 '$zeroa e))) 2102 (if (among '$zerob e) (setq e (maxima-substitute 0 '$zerob e))) 2103 e) 2104 2105;; simple radical 2106;; returns true if exp is a polynomial raised to a numeric power 2107(defun simplerd (exp) 2108 (and (mexptp exp) 2109 (mnump (caddr exp)) ;; exponent must be a number - no variables 2110 (polyp (cadr exp)))) 2111 2112(defun branch1 (exp val) 2113 (cond ((polyp exp) nil) 2114 ((simplerd exp) (zerop2 (subin val (cadr exp)))) 2115 (t 2116 (loop for v on (cdr exp) 2117 when (branch1 (car v) val) 2118 do (return v))))) 2119 2120(defun branch (exp val) 2121 (cond ((polyp exp) nil) 2122 ((or (simplerd exp) (mtimesp exp)) 2123 (branch1 exp val)) 2124 ((mplusp exp) 2125 (every #'(lambda (j) (branch j val)) (the list (cdr exp)))))) 2126 2127(defun ser0 (e n d val) 2128 (cond ((and (branch n val) (branch d val)) 2129 (setq nn* nil) 2130 (setq n (ser1 n)) 2131 (setq d (ser1 d)) 2132;;;NN* gets set by POFX, called by SER1, to get a list of exponents. 2133 (setq nn* (ratmin nn*)) 2134 (setq n (sratsimp (m^ n (m^ var nn*)))) 2135 (setq d (sratsimp (m^ d (m^ var nn*)))) 2136 (cond ((member val '($zeroa $zerob) :test #'eq) nil) 2137 (t (setq val 0.))) 2138 (radlim e n d)) 2139 (t (try-lhospital-quit n d nil)))) 2140 2141(defun rheur (l1 l2) 2142 (prog (ans m1 m2) 2143 (setq m1 (mapcar (function asymredu) l1)) 2144 (setq m2 (mapcar (function asymredu) l2)) 2145 (setq ans (m+l (append m1 m2))) 2146 (cond ((rptrouble (m+l (append l1 l2))) 2147 (return (limit (simplify (rdsget (m+l (append l1 l2)))) 2148 var 2149 val 2150 nil))) 2151 ((mplusp ans) (return (sheur m1 m2))) 2152 (t (return (limit ans var val t)))))) 2153 2154(defun rptrouble (rp) 2155 (not (equal (rddeg rp nil) (rddeg (asymredu rp) nil)))) 2156 2157(defun radicalp (exp var) 2158 (cond ((polyinx exp var ())) 2159 ((mexptp exp) (cond ((equal (caddr exp) -1.) 2160 (radicalp (cadr exp) var)) 2161 ((simplerd exp)))) 2162 ((member (caar exp) '(mplus mtimes) :test #'eq) 2163 (every #'(lambda (j) (radicalp j var)) 2164 (cdr exp))))) 2165 2166(defun involve (e nn*) 2167 (declare (special var)) 2168 (cond ((atom e) nil) 2169 ((mnump e) nil) 2170 ((and (member (caar e) nn* :test #'eq) (among var (cdr e))) (cadr e)) 2171 (t (some #'(lambda (j) (involve j nn*)) (cdr e))))) 2172 2173(defun notinvolve (exp nn*) 2174 (cond ((atom exp) t) 2175 ((mnump exp) t) 2176 ((member (caar exp) nn* :test #'eq) (not (among var (cdr exp)))) 2177 ((every #'(lambda (j) (notinvolve j nn*)) 2178 (cdr exp))))) 2179 2180(defun sheur1 (l1 l2) 2181 (prog (ans) 2182 (setq l1 (m+l (maxi l1))) 2183 (setq l2 (m+l (maxi l2))) 2184 (setq ans (cpa l1 l2 t)) 2185 (return (cond ((=0 ans) (m+ l1 l2)) 2186 ((equal ans 1.) '$inf) 2187 (t '$minf))))) 2188 2189(defun zero-lim (cpa-list) 2190 (do ((l cpa-list (cdr l))) 2191 ((null l) ()) 2192 (and (eq (caar l) 'gen) 2193 (zerop2 (limit (cadar l) var val 'think)) 2194 (return t)))) 2195 2196;; Compare order of growth for R1 and R2. The result is 0, -1, +1 2197;; depending on the relative order of growth. 0 is returned if R1 and 2198;; R2 have the same growth; -1 if R1 grows much more slowly than R2; 2199;; +1 if R1 grows much more quickly than R2. 2200(defun cpa (r1 r2 flag) 2201 (let ((t1 r1) 2202 (t2 r2)) 2203 (cond ((alike1 t1 t2) 0.) 2204 ((free t1 var) 2205 (cond ((free t2 var) 0.) 2206 (t (let ((lim-ans (limit1 t2 var val))) 2207 (cond ((not (member lim-ans '($inf $minf $und $ind) :test #'eq)) 0.) 2208 (t -1.)))))) 2209 ((free t2 var) 2210 (let ((lim-ans (limit1 t1 var val))) 2211 (cond ((not (member lim-ans '($inf $minf $und $ind) :test #'eq)) 0.) 2212 (t 1.)))) 2213 (t 2214 ;; Make T1 and T2 be a list of terms that are multiplied 2215 ;; together. 2216 (cond ((mtimesp t1) (setq t1 (cdr t1))) 2217 (t (setq t1 (list t1)))) 2218 (cond ((mtimesp t2) (setq t2 (cdr t2))) 2219 (t (setq t2 (list t2)))) 2220 ;; Find the strengths of each term of T1 and T2 2221 (setq t1 (mapcar (function istrength) t1)) 2222 (setq t2 (mapcar (function istrength) t2)) 2223 ;; Compute the max of the strengths of the terms. 2224 (let ((ans (ismax t1)) 2225 (d (ismax t2))) 2226 (cond ((or (null ans) (null d) 2227 (eq (car ans) 'gen) (eq (car d) 'gen)) 0.)) 2228 (if (eq (car ans) 'var) (setq ans (add-up-deg t1))) 2229 (if (eq (car d) 'var) (setq d (add-up-deg t2))) 2230 ;; Can't just just compare dominating terms if there are 2231 ;; indeterm-inates present; e.g. X-X^2*LOG(1+1/X). So 2232 ;; check for this. 2233 (cond ((or (zero-lim t1) 2234 (zero-lim t2)) 2235 (cpa-indeterm ans d t1 t2 flag)) 2236 ((isgreaterp ans d) 1.) 2237 ((isgreaterp d ans) -1.) 2238 (t 0))))))) 2239 2240(defun cpa-indeterm (ans d t1 t2 flag) 2241 (cond ((not (eq (car ans) 'var)) 2242 (setq ans (gather ans t1) d (gather d t2)))) 2243 (let ((*indicator (and (eq (car ans) 'exp) 2244 flag)) 2245 (test ())) 2246 (setq test (cpa1 ans d)) 2247 (cond ((and (zerop1 test) 2248 (or (equal ($radcan (m// (cadr ans) (cadr d))) 1.) 2249 (and (polyp (cadr ans)) 2250 (polyp (cadr d)) 2251 (equal (limit (m// (cadr ans) (cadr d)) var val 'think) 2252 1.)))) 2253 (let ((new-term1 (m// t1 (cadr ans))) 2254 (new-term2 (m// t2 (cadr d)))) 2255 (cpa new-term1 new-term2 flag))) 2256 (t 0)))) 2257 2258(defun add-up-deg (strengthl) 2259 (do ((stl strengthl (cdr stl)) 2260 (poxl) 2261 (degl)) 2262 ((null stl) (list 'var (m*l poxl) (m+l degl))) 2263 (cond ((eq (caar stl) 'var) 2264 (push (cadar stl) poxl) 2265 (push (caddar stl) degl))))) 2266 2267(defun cpa1 (p1 p2) 2268 (prog (flag s1 s2) 2269 (cond ((eq (car p1) 'gen) (return 0.))) 2270 (setq flag (car p1)) 2271 (setq p1 (cadr p1)) 2272 (setq p2 (cadr p2)) 2273 (cond 2274 ((eq flag 'var) 2275 (setq s1 (istrength p1)) 2276 (setq s2 (istrength p2)) 2277 (return 2278 (cond 2279 ((isgreaterp s1 s2) 1.) 2280 ((isgreaterp s2 s1) -1.) 2281 (*indicator 2282 (setq *indicator nil) 2283 (cond 2284 ((and (poly? p1 var) (poly? p2 var)) 2285 (setq p1 (m- p1 p2)) 2286 (cond ((zerop1 p1) 0.) 2287 (t (getsignl (hot-coef p1))))) 2288 (t 2289 (setq s1 2290 (rheur (list p1) 2291 (list (m*t -1 p2)))) 2292 (cond ((zerop2 s1) 0.) 2293 ((ratgreaterp s1 0.) 1.) 2294 (t -1.))))) 2295 (t 0.)))) 2296 ((eq flag 'exp) 2297 (setq p1 (caddr p1)) 2298 (setq p2 (caddr p2)) 2299 (cond ((and (poly? p1 var) (poly? p2 var)) 2300 (setq p1 (m- p1 p2)) 2301 (return (cond ((or (zerop1 p1) 2302 (not (among var p1))) 2303 0.) 2304 (t (getsignl (hot-coef p1)))))) 2305 ((and (radicalp p1 var) (radicalp p2 var)) 2306 (setq s1 2307 (rheur (list p1) 2308 (list (m*t -1 p2)))) 2309 (return (cond ((eq s1 '$inf) 1.) 2310 ((eq s1 '$minf) -1.) 2311 ((mnump s1) 2312 (cond ((ratgreaterp s1 0.) 1.) 2313 ((ratgreaterp 0. s1) -1.) 2314 (t 0.))) 2315 (t 0.)))) 2316 (t (return (cpa p1 p2 t))))) 2317 ((eq flag 'log) 2318 (setq p1 (try-lhospital (asymredu p1) (asymredu p2) nil)) 2319 (return (cond ((zerop2 p1) -1.) 2320 ((real-infinityp p1) 1.) 2321 (t 0.))))))) 2322 2323;;;EXPRESSIONS TO ISGREATERP ARE OF THE FOLLOWING FORMS 2324;;; ("VAR" POLY DEG) 2325;;; ("EXP" %E^EXP) 2326;;; ("LOG" LOG(EXP)) 2327;;; ("FACT" <A FACTORIAL EXPRESSION>) 2328;;; ("GEN" <ANY OTHER TYPE OF EXPRESSION>) 2329 2330(defun isgreaterp (a b) 2331 (let ((ta (car a)) 2332 (tb (car b))) 2333 (cond ((or (eq ta 'gen) 2334 (eq tb 'gen)) ()) 2335 ((and (eq ta tb) (eq ta 'var)) 2336 (ratgreaterp (caddr a) (caddr b))) 2337 ((and (eq ta tb) (eq ta 'exp)) 2338 ;; Both are exponential order of infinity. Check the 2339 ;; exponents to determine which exponent is bigger. 2340 (eq (limit (m- `((%log) ,(second a)) `((%log) ,(second b))) 2341 var val 'think) 2342 '$inf)) 2343 ((member ta (cdr (member tb '(num log var exp fact gen) :test #'eq)) :test #'eq))))) 2344 2345(defun ismax (l) 2346 ;; Preprocess the list of products. Separate the terms that 2347 ;; exponentials and those that don't. Actually multiply the 2348 ;; exponential terms together to form a single term. Pass this and 2349 ;; the rest to ismax-core to find the max. 2350 (let (exp-terms non-exp-terms) 2351 (dolist (term l) 2352 (if (eq 'exp (car term)) 2353 (push term exp-terms) 2354 (push term non-exp-terms))) 2355 ;; Multiply the exp-terms together 2356 (if exp-terms 2357 (let ((product 1)) 2358 ;;(format t "exp-terms = ~A~%" exp-terms) 2359 (dolist (term exp-terms) 2360 (setf product (simplify (mul product (second term))))) 2361 ;;(format t "product = ~A~%" product) 2362 (setf product `(exp ,($logcontract product))) 2363 ;;(format t "product = ~A~%" product) 2364 (ismax-core (cons product non-exp-terms))) 2365 (ismax-core l)))) 2366 2367(defun ismax-core (l) 2368 (cond ((null l) ()) 2369 ((atom l) ()) 2370 ((= (length l) 1) (car l)) ;If there is only 1 thing give it back. 2371 ((every #'(lambda (x) 2372 (not (eq (car x) 'gen))) l) 2373 2374 (do ((l1 (cdr l) (cdr l1)) 2375 (temp-ans (car l)) 2376 (ans ())) 2377 ((null l1) ans) 2378 (cond ((isgreaterp temp-ans (car l1)) 2379 (setq ans temp-ans)) 2380 ((isgreaterp (car l1) temp-ans) 2381 (setq temp-ans (car l1)) 2382 (setq ans temp-ans)) 2383 (t (setq ans ()))))) 2384 (t ()))) 2385 2386;RETURNS LIST OF HIGH TERMS 2387(defun maxi (all) 2388 (cond ((atom all) nil) 2389 (t (do ((l (cdr all) (cdr l)) 2390 (hi-term (car all)) 2391 (total 1) ; running total constant factor of hi-term 2392 (hi-terms (ncons (car all))) 2393 (compare nil)) 2394 ((null l) (if (zerop2 total) ; if high-order terms cancel each other 2395 all ; keep everything 2396 hi-terms)) ; otherwise return list of high terms 2397 (setq compare (limit (m// (car l) hi-term) var val 'think)) 2398 (cond 2399 ((or (infinityp compare) 2400 (and (eq compare '$und) 2401 (zerop2 (limit (m// hi-term (car l)) var val 'think)))) 2402 (setq total 1) ; have found new high term 2403 (setq hi-terms (ncons (setq hi-term (car l))))) 2404 ((zerop2 compare) nil) 2405 ;; COMPARE IS IND, FINITE-VALUED, or und in both directions 2406 (t ; add to list of high terms 2407 (setq total (m+ total compare)) 2408 (setq hi-terms (append hi-terms (ncons (car l)))))))))) 2409 2410(defun ratmax (l) 2411 (prog (ans) 2412 (cond ((atom l) (return nil))) 2413 l1 (setq ans (car l)) 2414 l2 (setq l (cdr l)) 2415 (cond ((null l) (return ans)) 2416 ((ratgreaterp ans (car l)) (go l2)) 2417 (t (go l1))))) 2418 2419(defun ratmin (l) 2420 (prog (ans) 2421 (cond ((atom l) (return nil))) 2422 l1 (setq ans (car l)) 2423 l2 (setq l (cdr l)) 2424 (cond ((null l) (return ans)) 2425 ((ratgreaterp (car l) ans) (go l2)) 2426 (t (go l1))))) 2427 2428(defun pofx (e) 2429 (cond ((atom e) 2430 (cond ((eq e var) 2431 (push 1 nn*)) 2432 (t ()))) 2433 ((or (mnump e) (not (among var e))) nil) 2434 ((and (mexptp e) (eq (cadr e) var)) 2435 (push (caddr e) nn*)) 2436 ((simplerd e) nil) 2437 (t (mapc #'pofx (cdr e))))) 2438 2439(defun ser1 (e) 2440 (cond ((member val '($zeroa $zerob) :test #'eq) nil) 2441 (t (setq e (subin (m+ var val) e)))) 2442 (setq e (rdfact e)) 2443 (cond ((pofx e) e))) 2444 2445(defun gather (ind l) 2446 (prog (ans) 2447 (setq ind (car ind)) 2448 loop (cond ((null l) 2449 (return (list ind (m*l ans)))) 2450 ((equal (caar l) ind) 2451 (push (cadar l) ans))) 2452 (setq l (cdr l)) 2453 (go loop))) 2454 2455; returns rough class-of-growth of term 2456(defun istrength (term) 2457 (cond ((mnump term) (list 'num term)) 2458 ((atom term) (cond ((eq term var) 2459 (list 'var var 1.)) 2460 (t (list 'num term)))) 2461 ((not (among var term)) (list 'num term)) 2462 ((mplusp term) 2463 (let ((temp (ismax (mapcar #'istrength (cdr term))))) 2464 (cond ((not (null temp)) temp) 2465 (t `(gen ,term))))) 2466 ((mtimesp term) 2467 (let ((temp (mapcar #'istrength (cdr term))) 2468 (temp1 ())) 2469 (setq temp1 (ismax temp)) 2470 (cond ((null temp1) `(gen ,term)) 2471 ((eq (car temp1) 'log) `(log ,temp)) 2472 ((eq (car temp1) 'var) (add-up-deg temp)) 2473 (t `(gen ,temp))))) 2474 ((and (mexptp term) 2475 (real-infinityp (limit term var val t))) 2476 (let ((logterm (logred term))) 2477 (cond ((and (among var (caddr term)) 2478 (member (car (istrength logterm)) 2479 '(var exp fact) :test #'eq) 2480 (real-infinityp (limit logterm var val t))) 2481 (list 'exp (m^ '$%e logterm))) 2482 ((not (among var (caddr term))) 2483 (let ((temp (istrength (cadr term)))) 2484 (cond ((not (alike1 temp term)) 2485 (rplaca (cdr temp) term) 2486 (and (eq (car temp) 'var) 2487 (rplaca (cddr temp) 2488 (m* (caddr temp) (caddr term)))) 2489 temp) 2490 (t `(gen ,term))))) 2491 (t `(gen ,term))))) 2492 ((and (eq (caar term) '%log) 2493 (real-infinityp (limit term var val t))) 2494 (let ((stren (istrength (cadr term)))) 2495 (cond ((member (car stren) '(log var) :test #'eq) 2496 `(log ,term)) 2497 ((and (eq (car stren) 'exp) 2498 (eq (caar (second stren)) 'mexpt)) 2499 (istrength (logred (second stren)))) 2500 (t `(gen ,term))))) 2501 ((eq (caar term) 'mfactorial) 2502 (list 'fact term)) 2503 ((let ((temp (hyperex term))) 2504 (and (not (alike1 term temp)) 2505 (istrength temp)))) 2506 (t (list 'gen term)))) 2507 2508;; log reduce - returns log of s1 2509(defun logred (s1) 2510 (or (and (eq (cadr s1) '$%e) (caddr s1)) 2511 (m* (caddr s1) `((%log) ,(cadr s1))))) 2512 2513(defun asymredu (rd) 2514 (cond ((atom rd) rd) 2515 ((mnump rd) rd) 2516 ((not (among var rd)) rd) 2517 ((polyinx rd var t)) 2518 ((simplerd rd) 2519 (cond ((eq (cadr rd) var) rd) 2520 (t (mabs-subst 2521 (factor ($expand (m^ (polyinx (cadr rd) var t) 2522 (caddr rd)))) 2523 var 2524 val)))) 2525 (t (simplify (cons (list (caar rd)) 2526 (mapcar #'asymredu (cdr rd))))))) 2527 2528(defun rdfact (rd) 2529 (let ((dn** ()) (nn** ())) 2530 (cond ((atom rd) rd) 2531 ((mnump rd) rd) 2532 ((not (among var rd)) rd) 2533 ((polyp rd) 2534 (factor rd)) 2535 ((simplerd rd) 2536 (cond ((eq (cadr rd) var) rd) 2537 (t (setq dn** (caddr rd)) 2538 (setq nn** (factor (cadr rd))) 2539 (cond ((mtimesp nn**) 2540 (m*l (mapcar #'(lambda (j) (m^ j dn**)) 2541 (cdr nn**)))) 2542 (t rd))))) 2543 (t (simplify (cons (ncons (caar rd)) 2544 (mapcar #'rdfact (cdr rd)))))))) 2545 2546(defun cnv (expl val) 2547 (mapcar #'(lambda (e) 2548 (maxima-substitute (cond ((eq val '$zerob) 2549 (m* -1 (m^ var -1))) 2550 ((eq val '$zeroa) 2551 (m^ var -1)) 2552 ((eq val '$minf) 2553 (m* -1 var)) 2554 (t (m^ (m+ var (m* -1 val)) -1.))) 2555 var 2556 e)) 2557 expl)) 2558 2559(defun pwtaylor (exp var l terms) 2560 (prog (coef ans c mc) 2561 (cond ((=0 terms) (return nil)) ((=0 l) (setq mc t))) 2562 (setq c 0.) 2563 (go tag1) 2564 loop (setq c (1+ c)) 2565 (cond ((or (> c 10.) (equal c terms)) 2566 (return (m+l ans))) 2567 (t (setq exp (sdiff exp var)))) 2568 tag1 (setq coef ($radcan (subin l exp))) 2569 (cond ((=0 coef) (setq terms (1+ terms)) (go loop))) 2570 (setq 2571 ans 2572 (append 2573 ans 2574 (list 2575 (m* coef 2576 (m^ `((mfactorial) ,c) -1) 2577 (m^ (if mc var (m+t (m*t -1 l) var)) c))))) 2578 (go loop))) 2579 2580(defun rdsget (e) 2581 (cond ((polyp e) e) 2582 ((simplerd e) (rdtay e)) 2583 (t (cons (list (caar e)) 2584 (mapcar #'rdsget (cdr e)))))) 2585 2586(defun rdtay (rd) 2587 (cond (limit-using-taylor ($ratdisrep ($taylor rd var val 1.))) 2588 (t (lrdtay rd)))) 2589 2590(defun lrdtay (rd) 2591 (prog (varlist p c e d $ratfac) 2592 (setq varlist (ncons var)) 2593 (setq p (ratnumerator (cdr (ratrep* (cadr rd))))) 2594 (cond ((< (length p) 3.) (return rd))) 2595 (setq e (caddr rd)) 2596 (setq d (pdegr p)) 2597 (setq c (m^ var (m* d e))) 2598 (setq d ($ratsimp (varinvert (m* (pdis p) (m^ var (m- d))) 2599 var))) 2600 (setq d (pwtaylor (m^ d e) var 0. 3.)) 2601 (return (m* c (varinvert d var))))) 2602 2603(defun varinvert (e var) (subin (m^t var -1.) e)) 2604 2605(defun deg (p) 2606 (prog ((varlist (list var))) 2607 (return (let (($ratfac nil)) 2608 (newvar p) 2609 (pdegr (cadr (ratrep* p))))))) 2610 2611(defun rat-no-ratfac (e) 2612 (let (($ratfac nil)) 2613 (newvar e) 2614 (ratrep* e))) 2615(setq low* nil) 2616 2617(defun rddeg (rd low*) 2618 (cond ((or (mnump rd) 2619 (not (among var rd))) 2620 0) 2621 ((polyp rd) 2622 (deg rd)) 2623 ((simplerd rd) 2624 (m* (deg (cadr rd)) (caddr rd))) 2625 ((mtimesp rd) 2626 (addn (mapcar #'(lambda (j) 2627 (rddeg j low*)) 2628 (cdr rd)) nil)) 2629 ((and (mplusp rd) 2630 (setq rd (andmapcar #'(lambda (j) (rddeg j low*)) 2631 (cdr rd)))) 2632 (cond (low* (ratmin rd)) 2633 (t (ratmax rd)))))) 2634 2635(defun pdegr (pf) 2636 (cond ((or (atom pf) (not (eq (caadr (ratf var)) (car pf)))) 2637 0) 2638 (low* (cadr (reverse pf))) 2639 (t (cadr pf)))) 2640;;There is some confusion here. We need to be aware of Branch cuts etc.... 2641;;when doing this section of code. It is not very carefully done so there 2642;;are bugs still lurking. Another misfortune is that LIMIT or its inferiors 2643;;sometimes decides to change the limit VAL in midstream. This must be corrected 2644;;since LIMIT's interaction with the data base environment must be maintained. 2645;;I'm not sure that this code can ever be called with VAL other than $INF but 2646;;there is a hook in the first important cond clause to cathc them anyway. 2647 2648(defun asy (n d) 2649 (let ((num-power (rddeg n nil)) 2650 (den-power (rddeg d nil)) 2651 (coef ()) (coef-sign ()) (power ())) 2652 (setq coef (m// ($ratcoef ($expand n) var num-power) 2653 ($ratcoef ($expand d) var den-power))) 2654 (setq coef-sign (getsignl coef)) 2655 (setq power (m// num-power den-power)) 2656 (cond ((eq (ask-integer power '$integer) '$integer) 2657 (cond ((eq (ask-integer power '$even) '$even) '$even) 2658 (t '$odd)))) ;Can be extended from here. 2659 (cond ((or (eq val '$minf) 2660 (eq val '$zerob) 2661 (eq val '$zeroa) 2662 (equal val 0)) ()) ;Can be extended to cover some these. 2663 ((ratgreaterp den-power num-power) 2664 (cond ((equal coef-sign 1.) '$zeroa) 2665 ((equal coef-sign -1) '$zerob) 2666 ((equal coef-sign 0) 0) 2667 (t 0))) 2668 ((ratgreaterp num-power den-power) 2669 (cond ((equal coef-sign 1.) '$inf) 2670 ((equal coef-sign -1) '$minf) 2671 ((equal coef-sign 0) nil) ; should never be zero 2672 ((null coef-sign) '$infinity))) 2673 (t coef)))) 2674 2675(defun radlim (e n d) 2676 (prog (nl dl) 2677 (cond ((eq val '$infinity) (throw 'limit nil)) 2678 ((eq val '$minf) 2679 (setq nl (m* var -1)) 2680 (setq n (subin nl n)) 2681 (setq d (subin nl d)) 2682 (setq val '$inf))) ;This is the Culprit. Doesn't tell the DATABASE. 2683 (cond ((eq val '$inf) 2684 (setq nl (asymredu n)) 2685 (setq dl (asymredu d)) 2686 (cond 2687 ((or (rptrouble n) (rptrouble d)) 2688 (return (limit (m* (rdsget n) (m^ (rdsget d) -1.)) var val t))) 2689 (t (return (asy nl dl)))))) 2690 (setq nl (limit n var val t)) 2691 (setq dl (limit d var val t)) 2692 (cond ((and (zerop2 nl) (zerop2 dl)) 2693 (cond ((or (polyp n) (polyp d)) 2694 (return (try-lhospital-quit n d t))) 2695 (t (return (ser0 e n d val))))) 2696 (t (return ($radcan (ratrad (m// n d) n d nl dl))))))) 2697 2698(defun ratrad (e n d nl dl) 2699 (prog (n1 d1) 2700 (cond ((equal nl 0) (return 0)) 2701 ((zerop2 dl) 2702 (setq n1 nl) 2703 (cond ((equal dl 0) (setq d1 '$infinity)) ;No direction Info. 2704 ((eq dl '$zeroa) 2705 (setq d1 '$inf)) 2706 ((equal (setq d (behavior d var val)) 1) 2707 (setq d1 '$inf)) 2708 ((equal d -1) (setq d1 '$minf)) 2709 (t (throw 'limit nil)))) 2710 ((zerop2 nl) 2711 (setq d1 dl) 2712 (cond ((equal (setq n (behavior n var val)) 1) 2713 (setq n1 '$zeroa)) 2714 ((equal n -1) (setq n1 '$zerob)) 2715 (t (setq n1 0)))) 2716 (t (return ($radcan (ridofab (subin val e)))))) 2717 (return (simplimtimes (list n1 d1))))) 2718 2719;;; Limit of the Logarithm function 2720 2721(defun simplimln (expr var val) 2722 ;; We need to be careful with log because of the branch cut on the 2723 ;; negative real axis. So we look at the imagpart of the argument. If 2724 ;; it's not identically zero, we compute the limit of the real and 2725 ;; imaginary parts and combine them. Otherwise, we can use the 2726 ;; original method for real limits. 2727 (let ((arglim (limit (cadr expr) var val 'think))) 2728 (cond ((eq arglim '$inf) '$inf) 2729 ((member arglim '($minf $infinity) :test #'eq) 2730 '$infinity) 2731 ((member arglim '($ind $und) :test #'eq) '$und) 2732 ((equalp ($imagpart (cadr expr)) 0) 2733 ;; argument is real. 2734 (let* ((real-lim (ridofab arglim))) 2735 (if (equalp real-lim 0) 2736 (cond ((eq arglim '$zeroa) '$minf) 2737 ((eq arglim '$zerob) '$infinity) 2738 (t (let ((dir (behavior (cadr expr) var val))) 2739 (cond ((equal dir 1) '$minf) 2740 ((equal dir -1) '$infinity) 2741 (t (throw 'limit t)))))) 2742 (cond ((equalp arglim 1) 2743 (let ((dir (behavior (cadr expr) var val))) 2744 (if (equal dir 1) '$zeroa 0))) 2745 (t ;; Call simplifier to get value at the limit of the argument. 2746 (simplify `((%log) ,real-lim))))))) 2747 (t ;; argument is complex. 2748 (destructuring-bind (rp . ip) 2749 (trisplit expr) 2750 (if (eq (setq rp (limit rp var val 'think)) '$minf) 2751 ;; Realpart is minf, do not return minf+%i*ip but infinity. 2752 '$infinity 2753 ;; Return a complex limit value. 2754 (add rp (mul '$%i (limit ip var val 'think))))))))) 2755 2756;;; Limit of the Factorial function 2757 2758(defun simplimfact (expr var val) 2759 (let* ((arglim (limit (cadr expr) var val 'think)) ; Limit of the argument. 2760 (arg2 arglim)) 2761 (cond ((eq arglim '$inf) '$inf) 2762 ((member arglim '($minf $infinity $und $ind) :test #'eq) '$und) 2763 ((and (or (maxima-integerp arglim) 2764 (setq arg2 (integer-representation-p arglim))) 2765 (eq ($sign arg2) '$neg)) 2766 ;; A negative integer or float or bigfloat representation. 2767 (let ((dir (limit (add (cadr expr) (mul -1 arg2)) var val 'think)) 2768 (even (mevenp arg2))) 2769 (cond ((or (and even 2770 (eq dir '$zeroa)) 2771 (and (not even) 2772 (eq dir '$zerob))) 2773 '$minf) 2774 ((or (and even 2775 (eq dir '$zerob)) 2776 (and (not even) 2777 (eq dir '$zeroa))) 2778 '$inf) 2779 (t (throw 'limit nil))))) 2780 (t 2781 ;; Call simplifier to get value at the limit of the argument. 2782 (simplify (list '(mfactorial) arglim)))))) 2783 2784(defun simplim%erf-%tanh (fn arg) 2785 (let ((arglim (limit arg var val 'think)) 2786 (ans ()) 2787 (rlim ())) 2788 (cond ((eq arglim '$inf) 1) 2789 ((eq arglim '$minf) -1) 2790 ((eq arglim '$infinity) 2791 (destructuring-bind (rpart . ipart) 2792 (trisplit arg) 2793 (setq rlim (limit rpart var origval 'think)) 2794 (cond ((eq fn '%tanh) 2795 (cond ((equal rlim '$inf) 1) 2796 ((equal rlim '$minf) -1))) 2797 ((eq fn '%erf) 2798 (setq ans (limit (m* rpart (m^t ipart -1)) var origval 'think)) 2799 (setq ans ($asksign (m+ `((mabs) ,ans) -1))) 2800 (cond ((or (eq ans '$pos) (eq ans '$zero)) 2801 (cond ((eq rlim '$inf) 1) 2802 ((eq rlim '$minf) -1) 2803 (t '$und))) 2804 (t '$und)))))) 2805 ((eq arglim '$und) '$und) 2806 ((member arglim '($zeroa $zerob $ind) :test #'eq) arglim) 2807;;;Ignore tanh(%pi/2*%I) and multiples of the argument. 2808 (t 2809 ;; erf (or tanh) of a known value is just erf(arglim). 2810 (simplify (list (ncons fn) arglim)))))) 2811 2812(defun simplim%atan (exp1) 2813 (cond ((zerop2 exp1) exp1) 2814 ((member exp1 '($und $infinity) :test #'eq) 2815 (throw 'limit ())) 2816 ((eq exp1 '$inf) half%pi) 2817 ((eq exp1 '$minf) 2818 (m*t -1. half%pi)) 2819 (t `((%atan) ,exp1)))) 2820 2821;; Most instances of atan2 are simplified to expressions in atan 2822;; by simpatan2 before we get to this point. This routine handles 2823;; tricky cases such as limit(atan2((x^2-2), x^3-2*x), x, sqrt(2), minus). 2824;; Taylor and Gruntz cannot handle the discontinuity at atan(0, -1) 2825(defun simplim%atan2 (exp) 2826 (let* ((exp1 (cadr exp)) 2827 (exp2 (caddr exp)) 2828 (lim1 (limit (cadr exp) var val 'think)) 2829 (lim2 (limit (caddr exp) var val 'think)) 2830 (sign2 ($csign lim2))) 2831 (cond ((and (zerop2 lim1) ;; atan2( 0+, + ) 2832 (eq sign2 '$pos)) 2833 lim1) ;; result is zeroa or zerob 2834 ((and (eq lim1 '$zeroa) 2835 (eq sign2 '$neg)) 2836 '$%pi) 2837 ((and (eq lim1 '$zerob) ;; atan2( 0-, - ) 2838 (eq sign2 '$neg)) 2839 (m- '$%pi)) 2840 ((and (eq lim1 '$zeroa) ;; atan2( 0+, 0 ) 2841 (zerop2 lim2)) 2842 (simplim%atan (limit (m// exp1 exp2) var val 'think))) 2843 ((and (eq lim1 '$zerob) ;; atan2( 0-, 0 ) 2844 (zerop2 lim2)) 2845 (m+ (porm (eq lim2 '$zeroa) '$%pi) 2846 (simplim%atan (limit (m// exp1 exp2) var val 'think)))) 2847 ((member lim1 '($und $infinity) :test #'eq) 2848 (throw 'limit ())) 2849 ((eq lim1 '$inf) half%pi) 2850 ((eq lim1 '$minf) 2851 (m*t -1. half%pi)) 2852 (t (take '($atan2) lim1 lim2))))) 2853 2854(defun simplimsch (sch arg) 2855 (cond ((real-infinityp arg) 2856 (cond ((eq sch '%sinh) arg) (t '$inf))) 2857 ((eq arg '$infinity) '$infinity) 2858 ((eq arg '$ind) '$ind) 2859 ((eq arg '$und) '$und) 2860 (t (let (($exponentialize t)) 2861 (resimplify (list (ncons sch) (ridofab arg))))))) 2862 2863;; simple limit of sin and cos 2864(defun simplimsc (exp fn arg) 2865 (cond ((member arg '($inf $minf $ind) :test #'eq) '$ind) 2866 ((member arg '($und $infinity) :test #'eq) 2867 (throw 'limit ())) 2868 ((member arg '($zeroa $zerob) :test #'eq) 2869 (cond ((eq fn '%sin) arg) 2870 (t (m+ 1 '$zerob)))) 2871 ((sincoshk exp 2872 (simplify (list (ncons fn) (ridofab arg))) 2873 fn)))) 2874 2875(defun simplim%tan (arg) 2876 (let ((arg1 (ridofab (limit arg var val 'think)))) 2877 (cond 2878 ((member arg1 '($inf $minf $infinity $ind $und) :test #'eq) '$und) 2879 ((pip arg1) 2880 (let ((c (trigred (pip arg1)))) 2881 (cond ((not (equal ($imagpart arg1) 0)) '$infinity) 2882 ((and (eq (caar c) 'rat) 2883 (equal (caddr c) 2) 2884 (> (cadr c) 0)) 2885 (setq arg1 (behavior arg var val)) 2886 (cond ((= arg1 1) '$inf) 2887 ((= arg1 -1) '$minf) 2888 (t '$und))) 2889 ((and (eq (caar c) 'rat) 2890 (equal (caddr c) 2) 2891 (< (cadr c) 0)) 2892 (setq arg1 (behavior arg var val)) 2893 (cond ((= arg1 1) '$minf) 2894 ((= arg1 -1) '$inf) 2895 (t '$und))) 2896 (t (throw 'limit ()))))) 2897 ((equal arg1 0) 2898 (setq arg1 (behavior arg var val)) 2899 (cond ((equal arg1 1) '$zeroa) 2900 ((equal arg1 -1) '$zerob) 2901 (t 0))) 2902 (t (simp-%tan (list '(%tan) arg1) 1 nil))))) 2903 2904(defun simplim%asinh (arg) 2905 (cond ((member arg '($inf $minf $zeroa $zerob $ind $und) :test #'eq) 2906 arg) 2907 ((eq arg '$infinity) '$und) 2908 (t (simplify (list '(%asinh) (ridofab arg)))))) 2909 2910(defun simplim%acosh (arg) 2911 (cond ((equal (ridofab arg) 1) '$zeroa) 2912 ((eq arg '$inf) arg) 2913 ((eq arg '$minf) '$infinity) 2914 ((member arg '($und $ind $infinity) :test #'eq) '$und) 2915 (t (simplify (list '(%acosh) (ridofab arg)))))) 2916 2917(defun simplim%atanh (arg dir) 2918 ;; Compute limit(atanh(x),x,arg). If ARG is +/-1, we need to take 2919 ;; into account which direction we're approaching ARG. 2920 (cond ((zerop2 arg) arg) 2921 ((member arg '($ind $und $infinity $minf $inf) :test #'eq) 2922 '$und) 2923 ((equal (setq arg (ridofab arg)) 1.) 2924 ;; The limit at 1 should be complex infinity because atanh(x) 2925 ;; is complex for x > 1, but inf if we're approaching 1 from 2926 ;; below. 2927 (if (eq dir '$zerob) 2928 '$inf 2929 '$infinity)) 2930 ((equal arg -1.) 2931 ;; Same as above, except for the limit is at -1. 2932 (if (eq dir '$zeroa) 2933 '$minf 2934 '$infinity)) 2935 (t (simplify (list '(%atanh) arg))))) 2936 2937(defun simplim%asin-%acos (fn arg) 2938 (cond ((member arg '($und $ind $inf $minf $infinity) :test #'eq) 2939 '$und) 2940 ((and (eq fn '%asin) 2941 (member arg '($zeroa $zerob) :test #'eq)) 2942 arg) 2943 (t (simplify (list (ncons fn) (ridofab arg)))))) 2944 2945(defun simplim$li (order arg val) 2946 (if (and (not (equal (length order) 1)) 2947 (not (equal (length arg) 1))) 2948 (throw 'limit ()) 2949 (setq order (car order) 2950 arg (car arg))) 2951 (if (not (equal order 2)) 2952 (throw 'limit ()) 2953 (destructuring-bind (rpart . ipart) 2954 (trisplit arg) 2955 (cond ((not (equal ipart 0)) 2956 (throw 'limit ())) 2957 (t 2958 (setq rpart (limit rpart var val 'think)) 2959 (cond ((eq rpart '$zeroa) '$zeroa) 2960 ((eq rpart '$zerob) '$zerob) 2961 ((eq rpart '$minf) '$minf) 2962 ((eq rpart '$inf) '$infinity) 2963 (t (simplify (subfunmake '$li (list order) 2964 (list rpart)))))))))) 2965 2966(defun simplim$psi (order arg val) 2967 (if (and (not (equal (length order) 1)) 2968 (not (equal (length arg) 1))) 2969 (throw 'limit ()) 2970 (setq order (car order) 2971 arg (car arg))) 2972 (cond ((equal order 0) 2973 (destructuring-bind (rpart . ipart) 2974 (trisplit arg) 2975 (cond ((not (equal ipart 0)) (throw 'limit ())) 2976 (t (setq rpart (limit rpart var val 'think)) 2977 (cond ((eq rpart '$zeroa) '$minf) 2978 ((eq rpart '$zerob) '$inf) 2979 ((eq rpart '$inf) '$inf) 2980 ((eq rpart '$minf) '$und) 2981 ((equal (getsignl rpart) -1) (throw 'limit ())) 2982 (t (simplify (subfunmake '$psi (list order) 2983 (list rpart))))))))) 2984 ((and (integerp order) (> order 0) 2985 (equal (limit arg var val 'think) '$inf)) 2986 (cond ((mevenp order) '$zerob) 2987 ((moddp order) '$zeroa) 2988 (t (throw 'limit ())))) 2989 (t (throw 'limit ())))) 2990 2991(defun simplim%inverse_jacobi_ns (arg m) 2992 (if (or (eq arg '$inf) (eq arg '$minf)) 2993 0 2994 `((%inverse_jacobi_ns) ,arg ,m))) 2995 2996(defun simplim%inverse_jacobi_nc (arg m) 2997 (if (or (eq arg '$inf) (eq arg '$minf)) 2998 `((%elliptic_kc) ,m) 2999 `((%inverse_jacobi_nc) ,arg ,m))) 3000 3001(defun simplim%inverse_jacobi_sc (arg m) 3002 (if (or (eq arg '$inf) (eq arg '$minf)) 3003 `((%elliptic_kc) ,m) 3004 `((%inverse_jacobi_sc) ,arg ,m))) 3005 3006(defun simplim%inverse_jacobi_dc (arg m) 3007 (if (or (eq arg '$inf) (eq arg '$minf)) 3008 `((%elliptic_kc) ,m) 3009 `((%inverse_jacobi_dc) ,arg ,m))) 3010 3011(defun simplim%inverse_jacobi_cs (arg m) 3012 (if (or (eq arg '$inf) (eq arg '$minf)) 3013 0 3014 `((%inverse_jacobi_cs) ,arg ,m))) 3015 3016(defun simplim%inverse_jacobi_ds (arg m) 3017 (if (or (eq arg '$inf) (eq arg '$minf)) 3018 0 3019 `((%inverse_jacobi_ds) ,arg ,m))) 3020 3021(setf (get '%signum 'simplim%function) 'simplim%signum) 3022 3023(defun simplim%signum (e x pt) 3024 (let* ((e (limit (cadr e) x pt 'think)) (sgn (mnqp e 0))) 3025 (cond ((eq t sgn) (take '(%signum) e)) ;; limit of argument of signum is not zero 3026 ((eq nil sgn) '$und) ;; limit of argument of signum is zero (noncontinuous) 3027 (t (throw 'limit nil))))) ;; don't know 3028 3029;; more functions for limit to handle 3030 3031(defun lfibtophi (e) 3032 (cond ((not (involve e '($fib))) e) 3033 ((eq (caar e) '$fib) 3034 (let ((lnorecurse t)) 3035 ($fibtophi (list '($fib) (lfibtophi (cadr e))) lnorecurse))) 3036 (t (cons (car e) 3037 (mapcar #'lfibtophi (cdr e)))))) 3038 3039;;; FOLLOWING CODE MAKES $LDEFINT WORK 3040 3041(defmfun $ldefint (exp var ll ul &aux $logabs ans a1 a2) 3042 (setq $logabs t ans (sinint exp var) 3043 a1 (toplevel-$limit ans var ul '$minus) 3044 a2 (toplevel-$limit ans var ll '$plus)) 3045 (and (member a1 '($inf $minf $infinity $und $ind) :test #'eq) 3046 (setq a1 (nounlimit ans var ul))) 3047 (and (member a2 '($inf $minf $infinity $und $ind) :test #'eq) 3048 (setq a2 (nounlimit ans var ll))) 3049 ($expand (m- a1 a2))) 3050 3051(defun nounlimit (exp var val) 3052 (setq exp (restorelim exp)) 3053 (nconc (list '(%limit) exp var (ridofab val)) 3054 (cond ((eq val '$zeroa) '($plus)) 3055 ((eq val '$zerob) '($minus))))) 3056 3057;; replace noun form of %derivative and indefinite %integrate with gensym. 3058;; prevents substitution x -> x+1 for limit('diff(x+2,x), x, 1) 3059;; 3060;; however, this doesn't work for limit('diff(x+2,x)/x, x, inf) 3061;; because the rest of the limit code thinks the gensym is const wrt x. 3062(defun hide (exp) 3063 (cond ((atom exp) exp) 3064 ((or (eq '%derivative (caar exp)) 3065 (and (eq '%integrate (caar exp)) ; indefinite integral 3066 (null (cdddr exp)))) 3067 (hidelim exp (caar exp))) 3068 (t (cons (car exp) (mapcar 'hide (cdr exp)))))) 3069 3070(defun hidelim (exp func) 3071 (setq func (gensym)) 3072 (putprop func 3073 (hidelima exp) 3074 'limitsub) 3075 func) 3076 3077(defun hidelima (e) 3078 (if (among var e) 3079 (nounlimit e var val) 3080 e)) 3081 3082;;;Used by Defint also. 3083(defun oscip (e) 3084 (or (involve e '(%sin %cos %tan)) 3085 (among '$%i (%einvolve e)))) 3086 3087(defun %einvolve (e) 3088 (when (among '$%e e) (%einvolve01 e))) 3089 3090(defun %einvolve01 (e) 3091 (cond ((atom e) nil) 3092 ((mnump e) nil) 3093 ((and (mexptp e) 3094 (eq (cadr e) '$%e) 3095 (among var (caddr e))) 3096 (caddr e)) 3097 (t (some #'%einvolve (cdr e))))) 3098 3099(declare-top (unspecial *indicator nn* dn* exp var val origval taylored 3100 $tlimswitch logcombed lhp? lhcount $ratfac)) 3101 3102 3103;; GRUNTZ ALGORITHM 3104 3105;; Dominik Gruntz 3106;; "On Computing Limits in a Symbolic Manipulation System" 3107;; PhD Dissertation ETH Zurich 1996 3108 3109;; The algorithm identifies the most rapidly varying (MRV) subexpression, 3110;; replaces it with a new variable w, rewrites the expression in terms 3111;; of the new variable, and then repeats. 3112 3113;; The algorithm doesn't handle oscillating functions, so it can't do things like 3114;; limit(sin(x)/x, x, inf). 3115 3116;; To handle limits involving functions like gamma(x) and erf(x), the 3117;; gruntz algorithm requires them to be written in terms of asymptotic 3118;; expansions, which maxima cannot currently do. 3119 3120;; The algorithm assumes that everything is real, so it can't 3121;; currently handle limit((-2)^x, x, inf). 3122 3123;; This is one of the methods used by maxima's $limit. 3124;; It is also directly available to the user as $gruntz. 3125 3126 3127;; most rapidly varying subexpression of expression exp with respect to limit variable var. 3128;; returns a list of subexpressions which are in the same MRV equivalence class. 3129(defun mrv (exp var) 3130 (cond ((freeof var exp) 3131 nil) 3132 ((eq var exp) 3133 (list var)) 3134 ((mtimesp exp) 3135 (mrv-max (mrv (cadr exp) var) 3136 (mrv (m*l (cddr exp)) var) 3137 var)) 3138 ((mplusp exp) 3139 (mrv-max (mrv (cadr exp) var) 3140 (mrv (m+l (cddr exp)) var) 3141 var)) 3142 ((mexptp exp) 3143 (cond ((freeof var (caddr exp)) 3144 (mrv (cadr exp) var)) 3145 ((member (limitinf (logred exp) var) '($inf $minf) :test #'eq) 3146 (mrv-max (list exp) (mrv (caddr exp) var) var)) 3147 (t (mrv-max (mrv (cadr exp) var) (mrv (caddr exp) var) var)))) 3148 ((mlogp exp) 3149 (mrv (cadr exp) var)) 3150 ((equal (length (cdr exp)) 1) 3151 (mrv (cadr exp) var)) 3152 ((equal (length (cdr exp)) 2) 3153 (mrv-max (mrv (cadr exp) var) 3154 (mrv (caddr exp) var) 3155 var)) 3156 (t (tay-error "mrv not implemented" exp)))) 3157 3158;; takes two lists of expresions, f and g, and limit variable var. 3159;; members in each list are assumed to be in same MRV equivalence 3160;; class. returns MRV set of the union of the inputs - either f or g 3161;; or the union of f and g. 3162(defun mrv-max (f g var) 3163 (prog () 3164 (cond ((not f) 3165 (return g)) 3166 ((not g) 3167 (return f)) 3168 ((intersection f g) 3169 (return (union f g)))) 3170 (let ((c (mrv-compare (car f) (car g) var))) 3171 (cond ((eq c '>) 3172 (return f)) 3173 ((eq c '<) 3174 (return g)) 3175 ((eq c '=) 3176 (return (union f g))) 3177 (t (merror "MRV-MAX: expected '>' '<' or '='; found: ~M" c)))))) 3178 3179(defun mrv-compare (a b var) 3180 (let ((c (limitinf (m// `((%log) ,a) `((%log) ,b)) var))) 3181 (cond ((equal c 0) 3182 '<) 3183 ((member c '($inf $minf) :test #'eq) 3184 '>) 3185 (t '=)))) 3186 3187;; rewrite expression exp by replacing members of MRV set omega with 3188;; expressions in terms of new variable wsym. return cons pair of new 3189;; version of exp and the log of the new variable wsym. 3190(defun mrv-rewrite (exp omega var wsym) 3191 (setq omega (sort omega (lambda (x y) (> (length (mrv x var)) 3192 (length (mrv y var)))))) 3193 (let* ((g (car (last omega))) 3194 (logg (logred g)) 3195 (sig (equal (mrv-sign logg var) 1)) 3196 (w (if sig (m// 1 wsym) wsym)) 3197 (logw (if sig (m* -1 logg) logg))) 3198 (mapcar (lambda (x y) 3199 ;;(mtell "y:~M x:~M exp:~M~%" y x exp) 3200 (setq exp (syntactic-substitute y x exp))) 3201 omega 3202 (mapcar (lambda (f) ;; rewrite each element of omega 3203 (let* ((logf (logred f)) 3204 (c (mrv-leadterm (m// logf logg) var nil))) 3205 (cond ((not (equal (cadr c) 0)) 3206 (merror "MRV-REWRITE: expected leading term to be constant in ~M" c))) 3207 ;;(mtell "logg: ~M logf: ~M~%" logg logf) 3208 (m* (m^ w (car c)) 3209 (m^ '$%e (m- logf 3210 (m* (car c) logg)))))) 3211 omega)) 3212 (cons exp logw))) 3213 3214;;; if log w(x) = h(x), rewrite all subexpressions of the form 3215;;; log(f(x)) as log(w^-c f(x)) + c h(x) with c the unique constant 3216;;; such that w^-c f(x) is strictly less rapidly varying than w. 3217(defun mrv-rewrite-logs (exp wsym logw) 3218 (cond ((atom exp) exp) 3219 ((and (mlogp exp) 3220 (not (freeof wsym exp))) 3221 (let* ((f (cadr exp)) 3222 (c ($lopow (calculate-series f wsym) 3223 wsym))) 3224 (m+ (list (car exp) 3225 (m* (m^ wsym (m- c)) 3226 (mrv-rewrite-logs f wsym logw))) 3227 (m* c logw)))) 3228 (t 3229 (cons (car exp) 3230 (mapcar (lambda (e) 3231 (mrv-rewrite-logs e wsym logw)) 3232 (cdr exp)))))) 3233 3234;; returns list of two elements: coeff and exponent of leading term of exp, 3235;; after rewriting exp in term of its MRV set omega. 3236(defun mrv-leadterm (exp var omega) 3237 (prog ((new-omega ())) 3238 (cond ((freeof var exp) 3239 (return (list exp 0)))) 3240 (dolist (term omega) 3241 (cond ((subexp exp term) 3242 (push term new-omega)))) 3243 (setq omega new-omega) 3244 (cond ((not omega) 3245 (setq omega (mrv exp var)))) 3246 (cond ((member var omega :test #'eq) 3247 (let* ((omega-up (mrv-moveup omega var)) 3248 (e-up (car (mrv-moveup (list exp) var))) 3249 (mrv-leadterm-up (mrv-leadterm e-up var omega-up))) 3250 (return (mrv-movedown mrv-leadterm-up var))))) 3251 (destructuring-let* ((wsym (gensym "w")) 3252 lo 3253 coef 3254 ((f . logw) (mrv-rewrite exp omega var wsym)) 3255 (series (calculate-series (mrv-rewrite-logs f wsym logw) 3256 wsym))) 3257 (setq series (maxima-substitute logw `((%log) ,wsym) series)) 3258 (setq lo ($lopow series wsym)) 3259 (when (or (not ($constantp lo)) 3260 (not (free series '%derivative))) 3261 ;; (mtell "series: ~M lo: ~M~%" series lo) 3262 (tay-error "error in series expansion" f)) 3263 (setq coef ($coeff series wsym lo)) 3264 ;;(mtell "exp: ~M f: ~M~%" exp f) 3265 ;;(mtell "series: ~M~%coeff: ~M~%pow: ~M~%" series coef lo) 3266 (return (list coef lo))))) 3267 3268(defun mrv-moveup (l var) 3269 (mapcar (lambda (exp) 3270 (simplify-log-of-exp 3271 (syntactic-substitute `((mexpt) $%e ,var) var exp))) 3272 l)) 3273 3274(defun mrv-movedown (l var) 3275 (mapcar (lambda (exp) (syntactic-substitute `((%log simp) ,var) var exp)) 3276 l)) 3277 3278;; test whether sub is a subexpression of exp 3279(defun subexp (exp sub) 3280 (not (equal (maxima-substitute 'dummy 3281 sub 3282 exp) 3283 exp))) 3284 3285;; Generate $lhospitallim terms of taylor expansion. 3286;; Ideally we would use a lazy series representation that generates 3287;; more terms as higher order terms cancel. 3288(defun calculate-series (exp var) 3289 (assume `((mgreaterp) ,var 0)) 3290 (putprop var t 'internal);; keep var from appearing in questions to user 3291 (let ((series ($taylor exp var 0 $lhospitallim))) 3292 (forget `((mgreaterp) ,var 0)) 3293 series)) 3294 3295(defun mrv-sign (exp var) 3296 (cond ((freeof var exp) 3297 (let ((sign ($sign ($radcan exp)))) 3298 (cond ((eq sign '$zero) 3299 0) 3300 ((eq sign '$pos) 3301 1) 3302 ((eq sign '$neg) 3303 -1) 3304 (t (tay-error " cannot determine mrv-sign" exp))))) 3305 ((eq exp var) 3306 1) 3307 ((mtimesp exp) 3308 (* (mrv-sign (cadr exp) var) 3309 (mrv-sign (m*l (cddr exp)) var))) 3310 ((and (mexptp exp) 3311 (equal (mrv-sign (cadr exp) var) 1)) 3312 1) 3313 ((mlogp exp) 3314 (cond ((equal (mrv-sign (cadr exp) var) -1) 3315 (tay-error " complex expression in gruntz limit" exp))) 3316 (mrv-sign (m+ -1 (cadr exp)) var)) 3317 ((mplusp exp) 3318 (mrv-sign (limitinf exp var) var)) 3319 (t (tay-error " cannot determine mrv-sign" exp)))) 3320 3321;; gruntz algorithm for limit of exp as var goes to positive infinity 3322(defun limitinf (exp var) 3323 (prog (($exptsubst nil)) 3324 (cond ((freeof var exp) 3325 (return exp))) 3326 (destructuring-let* ((c0-e0 (mrv-leadterm exp var nil)) 3327 (c0 (car c0-e0)) 3328 (e0 (cadr c0-e0)) 3329 (sig (mrv-sign e0 var))) 3330 (cond ((equal sig 1) 3331 (return 0)) 3332 ((equal sig -1) 3333 (cond ((equal (mrv-sign c0 var) 1) 3334 (return '$inf)) 3335 ((equal (mrv-sign c0 var) -1) 3336 (return '$minf)))) 3337 ((equal sig 0) 3338 (return (limitinf c0 var))))))) 3339 3340;; user-level function equivalent to $limit. 3341;; direction must be specified if limit point is not infinite 3342;; The arguments are checked and a failure of taylor is catched. 3343 3344(defmfun $gruntz (expr var val &rest rest) 3345 (let (ans dir) 3346 (when (> (length rest) 1) 3347 (merror 3348 (intl:gettext "gruntz: too many arguments; expected just 3 or 4"))) 3349 (setq dir (car rest)) 3350 (when (and (not (member val '($inf $minf $zeroa $zerob))) 3351 (not (member dir '($plus $minus)))) 3352 (merror 3353 (intl:gettext "gruntz: direction must be 'plus' or 'minus'"))) 3354 (setq ans 3355 (catch 'taylor-catch 3356 (let ((silent-taylor-flag t)) 3357 (declare (special silent-taylor-flag)) 3358 (gruntz1 expr var val dir)))) 3359 (if (or (null ans) (eq ans t)) 3360 (if dir 3361 `(($gruntz simp) ,expr ,var, val ,dir) 3362 `(($gruntz simp) ,expr ,var ,val)) 3363 ans))) 3364 3365;; This function is for internal use in $limit. 3366(defun gruntz1 (exp var val &rest rest) 3367 (cond ((> (length rest) 1) 3368 (merror (intl:gettext "gruntz: too many arguments; expected just 3 or 4")))) 3369 (let (($logexpand t) ; gruntz needs $logexpand T 3370 (newvar (gensym "w")) 3371 (dir (car rest))) 3372 (cond ((eq val '$inf) 3373 (setq newvar var)) 3374 ((eq val '$minf) 3375 (setq exp (maxima-substitute (m* -1 newvar) var exp))) 3376 ((eq val '$zeroa) 3377 (setq exp (maxima-substitute (m// 1 newvar) var exp))) 3378 ((eq val '$zerob) 3379 (setq exp (maxima-substitute (m// -1 newvar) var exp))) 3380 ((eq dir '$plus) 3381 (setq exp (maxima-substitute (m+ val (m// 1 newvar)) var exp))) 3382 ((eq dir '$minus) 3383 (setq exp (maxima-substitute (m+ val (m// -1 newvar)) var exp))) 3384 (t (merror (intl:gettext "gruntz: direction must be 'plus' or 'minus'; found: ~M") dir))) 3385 (limitinf exp newvar))) 3386 3387;; substitute y for x in exp 3388;; similar to maxima-substitute but does not simplify result 3389(defun syntactic-substitute (y x exp) 3390 (cond ((alike1 x exp) y) 3391 ((atom exp) exp) 3392 (t (cons (car exp) 3393 (mapcar (lambda (exp) 3394 (syntactic-substitute y x exp)) 3395 (cdr exp)))))) 3396 3397;; log(exp(subexpr)) -> subexpr 3398;; without simplifying entire exp 3399(defun simplify-log-of-exp (exp) 3400 (cond ((atom exp) exp) 3401 ((and (mlogp exp) 3402 (mexptp (cadr exp)) 3403 (eq '$%e (cadadr exp))) 3404 (caddr (cadr exp))) 3405 (t (cons (car exp) 3406 (mapcar #'simplify-log-of-exp 3407 (cdr exp)))))) 3408