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;;; ************************************************************** 9;;; ***** HAYAT ******* Finite Power Series Routines ************* 10;;; ************************************************************** 11;;; ** (c) Copyright 1982 Massachusetts Institute of Technology ** 12;;; ****** This is a read-only file! (All writes reserved) ******* 13;;; ************************************************************** 14 15(in-package :maxima) 16 17;;; TOP LEVEL STRUCTURE 18 19;;; Power series have the following format when seen outside the power 20;;; series package: 21;;; 22;;; ((MRAT SIMP <varlist> <genvar> <tlist> trunc) <poly-form>) 23;;; 24;;; This is the form of the output of the expressions, to 25;;; be displayed they are RATDISREPed and passed to DISPLA. 26 27;;; The <poly-forms> consist of a header and list of exponent-coefficient 28;;; pairs as shown below. The PS is used to distinguish power series 29;;; from their coefficients which have a similar representation. 30;;; 31;;; (PS (<var> . <ord-num>) (<trunc-lvl>) 32;;; (<exponent> . <coeff>) (<exponent> . <coeff>) . . .) 33;;; 34;;; The <var> component of the power series is a gensym which represents the 35;;; kernel of the power series. If the package is called with the arguments: 36;;; Taylor(<expr>, x, a, n) then the kernel will be (x - a). 37;;; The <ord-num> is a relative ordering for the various kernels in a 38;;; multivariate expansion. 39;;; <trunc-lvl> is the highest degree of the variable <var> which is retained 40;;; in the current power series. 41;;; The terms in the list of exponent-coefficient pairs are ordered by 42;;; increasing degree. 43 44;;; Problem: fix expansion of logs so that taylor(log(1+exp(-1/x)),x,0,3) 45;;; works. Done. 46;;; 47;;; Problem: taylor(log(1+exp(-1/x)),x,0,5) loses because, while 48;;; taylor_simplify_recurse'ing exp(-3/x) get trunc level = -3. FIxed. 49;;; 50;;; Problem: Need to fix things so that asymptotic kernels aren't put onto 51;;; tvars via tlist merge etc. in taylor1. Done. 52;;; 53;;; Problem: get-series returns 0 for taylor(log(1+exp(-1/x)),x,0,5) and 54;;; need to make log(exp(1/x)) -> 1/x. Fixed. 55;;; 56;;; Problem: Fix psexpt-fn so that it doesn't lose via the invert-var 57;;; scheme, e.g. try taylor(exp(exp(-1/x)+x),x,0,5). Note that if it did 58;;; just that scheme. Done. 59;;; 60;;; Problem: fix adjoin-tvar so that the new tvars are ordered correctly 61;;; according to their strength. This is necessary in order to read the limit 62;;; directly from the leading term. E.g. see the misordered: 63;;; taylor(subst(1/x,x,part(screw2,1)),x,0,2) from ALJABR;SCREW2 LIMIT. 64;;; Note that the answer given for this appear to be incorrect when the 65;;; truncation on x is < 4. Is this due to the misordering? 66;;; Also taylor(screwa,x,0,4)+taylor(screwb,x,0,8) doesn't agree with 67;;; taylor(screw,x,0,8) where it should (here screwa = part(screw,1), 68;;; screwb = part(screw, 2); is this a truncation problem on the 69;;; gvar exp(1/x)?). 70;;; 71;;; Problem: new gvars have to be intro'd for logs just as for exp's instead 72;;; of treating them like constants as currently done. For example, 73;;; taylor(log(1+1/log(x)),x,0,2) currently doesn't expand. Done. 74;;; 75;;; Problem: The display routines need pieces of the taylor environment 76;;; (tvar-limits, tvars, tlist, etc.) to figure out how to order terms. 77;;; This means we'll probably have to store some of this on the local tlist. 78;;; When this is done the commented out code in psdisrep and psdisrep2 can 79;;; be re-inserted. Psdisrep2expand will also need to be modified. 80;;; I just fixed srdisrep to get the local env it needs; psdisrep2expand 81;;; still needs to be updated. The display order problem is still around 82;;; however: try taylor(exp(exp(-1/x)+x),x,0,3). After more investigation, 83;;; it seems that the term reversal always occurs for ps's that are coeff's 84;;; of terms whose expt is < 0. Probably the psdisrep routines should reverse 85;;; these terms to account for this (the bug is somewhere in the DISPLA 86;;; routines, possible DIM-MPLUS). 87;;; 88;;; Problem: Since gvar's like exp(-1/x) can't be put on genvar, they have 89;;; to be saved somewhere locally and restored by everyone who needs to setup 90;;; disrep info etc. Done for re-taylor. 91;;; 92;;; Problem: All new code needs to be checked to ensure it does the correct 93;;; thing when the expansion point is infinite (e.g. see the code in 94;;; TSEXPT-RED which handles this case). 95;;; 96;;; Perhaps the code for exp's and log's which pushes trunc degrees 97;;; can be done by first computing exp(c0) or log(c0) first and see 98;;; how much to push by looking at this series. Done for exp in tsexpt-red. 99;;; 100;;; Problems: taylor(part(screwa,2)-2/x,x,0,1) shouldn't be exact. 101;;; taylor(screwa,x,0,-2) misses the degree -2 term. This part is now fixed. 102;;; 103;;; Tvar-limits should be stored locally so that psdisrep need not recompute 104;;; each gvar limit when disrepping. 105 106(macsyma-module hayat) 107 108(defmvar tlist nil) 109 110(defvar *within-srf?* nil) 111 112(load-macsyma-macros mhayat rzmac ratmac) 113 114;;; Subtitle Special Stuff for Compiling 115 116(declare-top 117 (special vlist 118 varlist ;List of all the variables occurring in a power 119 ;series, the power series variables at the end 120 genvar ;The list of gensyms corresponding to varlist 121 modulus ; 122 *a* ;Temporary special 123 silent-taylor-flag ;If true indicates that errors will be 124 ;returned via a throw to TAY-ERR 125 tlist ;An association list which contains the 126 ;relevant information for the expansion which 127 ;is passed in at toplevel invocation. 128 $float ;Indicates whether to convert rational numbers 129 ;to floating point numbers. 130 $keepfloat ;When true retains floatin point numbers 131 ;internal to Taylor. 132 $radexpand ; 133 log-1 ;What log(-1) should be log(-1) or pi*i. 134 log%i ;Similarly for log(i) 135 exact-poly ;Inicates whether polynomials are to be 136 ;considered exact or not. True within SRF, 137 ;false within TAYLOR. 138 tvars ; 139 half%pi ;Has pi/2 to save space. 140 const-exp-funs ; 141 tay-const-expand ;For rediculousness like csch(log(x)) 142 $exponentialize ;which we do by exponentiation. 143 tay-pole-expand ; 144 trigdisp ; 145 last-exp ;last-expression through taylor2 146 $taylordepth ; 147 $ratexpand ; 148 genpairs ;List of dotted pairs 149 ps-bmt-disrep ; 150 ivars ;Pairlist if gensym and disreped version 151 key-vars ;Pairlist of gensym and key var (for searching 152 ;TLIST) 153 $algebraic ; 154 *psacirc ; 155 *pscirc ; 156 full-log ; 157 $logarc ; 158 trunclist ; 159 *within-srf?* ;flag for in srf 160 mainvar-datum ; 161 least_term? ; If non-null then the addition routines 162 ; are adding or subtracting coeff's of the 163 ; least term of a sum so they should do 164 ; zero checking on it if it is desired. 165 taylor_simplifier ; This is set by taylor1 to the object 166 ; which will be funcalled whenever 167 ; coefficient simplification is desired. 168 zerolist ; A list of constant expressions which have 169 ; been verified to be zero by a call to 170 ; $TAYLOR_SIMPLIFIER in taylor2. It is used to 171 ; suppress the message that TAYLOR is assumming 172 ; an expression to be zero. 173 ; 0p-funord lexp-non0 ; referenced only in commented-out code, so comment out here too 174 $zerobern $simp) 175 ) ;Don't want to see closed compilation notes. 176 177(defmvar $psexpand () 178 "When TRUE extended rational function expressions will be displayed fully 179 expanded. (RATEXPAND will also cause this.) If FALSE, multivariate 180 expressions will be displayed just as in the rational function package. 181 If PSEXPAND:MULTI, then terms with the same total degree in the variables 182 are grouped together.") 183 184(defmvar $maxtayorder t 185 "When true TAYLOR retains as many terms as are certain to be correct 186 during power series arithmetic. Otherwise, truncation is controlled 187 by the arguments specified to TAYLOR.") 188 189(defmvar $taylor_truncate_polynomials t 190 "When FALSE polynomials input to TAYLOR are considered to have infinite 191 precison; otherwise (the default) they are truncated based upon the input 192 truncation levels.") 193 194(defmvar $taylor_logexpand t 195 "Unless FALSE log's of products will be expanded fully in TAYLOR (the default) 196 to avoid identically-zero constant terms which involve log's. When FALSE, 197 only expansions necessary to produce a formal series will be executed.") 198 199(defmvar $taylor_simplifier 'simplify 200 "A function of one argument which TAYLOR uses to simplify coefficients 201 of power series.") 202 203(defvar taylor_simplifier nil) 204 205;;; Subtitle General Macsyma Free Predicates 206 207(defun zfree (e x) 208 (cond ((equal e x) () ) 209 ((atom e) 't) 210 ((eq (caar e) 'mrat) 211 (null (member x (cdr ($listofvars e)) :test #'equal))) 212 ('t (do ((l (cdr e) (cdr l))) ((null l) 't) 213 (or (zfree (car l) x) (return () )))))) 214 215(defun mfree (exp varl) 216 (declare (special dummy-variable-operators)) 217 (cond ((atom exp) (not (member exp varl :test #'eq))) 218 ((eq (caar exp) 'mrat) 219 (do ((l (mrat-varlist exp) (cdr l))) 220 ((null l) 't) 221 (unless (mfree (car l) varl) (return () )))) 222 ((or (member (caar exp) dummy-variable-operators :test #'eq) 223 (member 'array (cdar exp) :test #'eq)) 224 (do ((vars varl (cdr vars))) 225 ((null vars) 't) 226 (unless (freeof (car vars) exp) (return () )))) 227 ('t (and (mfree (caar exp) varl) (mfreel (cdr exp) varl))))) 228 229(defun mfreel (l varl) 230 (or (null l) (and (mfree (car l) varl) (mfreel (cdr l) varl)))) 231 232;;; Subtitle Coefficient Arithmetic 233 234(defun rcexpt (x y) 235 (cond ((equal x (rcone)) (rcone)) 236 ((rczerop y) (rcone)) 237 ((and (equal (cdr y) 1) (fixnump (car y))) 238 (ratexpt x (car y))) 239 ((and $radexpand (numberp (car y)) (numberp (cdr y))) 240 (if (floatp (car y)) 241 (setq y (maxima-rationalize (quot (car y) (cdr y))))) 242 (ratexpt (rcquo (rcexpt1 (car x) (cdr y)) 243 (rcexpt1 (cdr x) (cdr y))) 244 (car y))) 245 (t (let ($keepfloat) 246 (prep1 (m^ (rcdisrep x) (rcdisrep y))))))) 247 248(defun rcexpt1 (p n) 249 (cond ((equal p 1) (rcone)) 250 ((pcoefp p) (prep1 (m^ (pdis p) (*red 1 n)))) 251 ;; psfr does a square-free decom on p yielding (p1 e1 p2 e2 ... pn en) 252 ;; where p = p1^e1 p2^e2 ... pn^en, the pi being square-free 253 (t (do ((l (psqfr p) (cddr l)) 254 (ans (rcone))) 255 ((null l) ans) 256 (if (not (equal (rem (cadr l) n) 0)) 257 (setq ans (rctimes ans (prep1 (m^ (pdis (car l)) 258 (*red (cadr l) n))))) 259 ;; If pi<0, n=2m and n|ei then ei=2e and 260 ;; (pi^ei)^(1/(2m)) = (-pi)^(e/m) 261 (progn 262 (when (and (evenp n) (eq ($sign (pdis (car l))) '$neg)) 263 (rplaca l (pminus (car l)))) 264 (setq ans (rctimes ans (ratexpt (cons (car l) 1) 265 (truncate (cadr l) n)))))))))) 266 267(defun rccoefp (e) ;a sure check, but expensive 268 (and (null (atom e)) 269 (or (atom (car e)) 270 (member (caar e) genvar :test #'eq)) 271 (or (atom (cdr e)) 272 (member (cadr e) genvar :test #'eq)))) 273 274;;; Subtitle Exponent arithmetic 275 276(defun ezerop (x) 277 (and (not (infp x)) (signp e (car x)))) 278 279(defun e+ (x y) 280 (cond ((or (infp x) (infp y)) (inf)) 281 ((and (equal (cdr x) 1) (equal (cdr y) 1)) 282 (cons (+ (car x) (car y)) 1)) 283 (t (ereduce (+ (* (car x) (cdr y)) (* (cdr x) (car y))) 284 (* (cdr x) (cdr y)))))) 285 286(defun ediff (x y) 287 (cond ((infp x) (inf)) 288 ((and (equal (cdr x) 1) (equal (cdr y) 1)) 289 (cons (- (car x) (car y)) 1)) 290 (t (ereduce (- (* (car x) (cdr y)) (* (cdr x) (car y))) 291 (* (cdr x) (cdr y)))))) 292 293(defun emin (x y) 294 (cond ((infp x) y) 295 ((infp y) x) 296 ((equal (cdr x) (cdr y)) (cons (min (car x) (car y)) (cdr x))) 297 ((< (* (car x) (cdr y)) (* (cdr x) (car y))) x) 298 (t y))) 299 300(defun emax (x y) 301 (cond ((or (infp x) (infp y)) (inf)) 302 ((equal (cdr x) (cdr y)) (cons (max (car x) (car y)) (cdr x))) 303 ((> (* (car x) (cdr y)) (* (cdr x) (car y))) x) 304 (t y))) 305 306(defun e* (x y) 307 (cond ((or (infp x) (infp y)) (inf)) 308 ((and (equal (cdr x) 1) (equal (cdr y) 1)) 309 (cons (* (car x) (car y)) 1)) 310 (t (ereduce (* (car x) (car y)) (* (cdr x) (cdr y)))))) 311 312(defun erecip (e) 313 (if (minusp (car e)) 314 (cons (- (cdr e)) (- (car e))) 315 (cons (cdr e) (car e)))) 316 317(defun equo (x y) 318 (cond ((infp x) (inf)) 319 ((infp y) (rczero)) 320 (t (ereduce (* (car x) (cdr y)) 321 (* (cdr x) (car y)))))) 322 323(defun e1+ (x) 324 (cond ((infp x) (inf)) 325 ((= (cdr x) 1) (cons (1+ (car x)) 1)) 326 (t (cons (+ (cdr x) (car x)) (cdr x))))) 327 328(defun e1- (x) 329 (cond ((infp x) (inf)) 330 ((equal (cdr x) 1) (cons (1- (car x)) 1)) 331 (t (cons (- (car x) (cdr x)) (cdr x))))) 332 333(defun e> (x y) 334 (cond ((infp x) t) 335 ((infp y) ()) 336 ((equal (cdr x) (cdr y)) (> (car x) (car y))) 337 (t (> (* (car x) (cdr y)) (* (car y) (cdr x)))))) 338 339(defun e= (e1 e2) 340 (cond ((eq e1 e2) t) 341 ((or (null e1) (null e2)) ()) 342 (t (and (equal (car e1) (car e2)) 343 (equal (cdr e1) (cdr e2)))))) 344 345(defun ereduce (n d) 346 (if (signp l d) (setq d (- d) n (- n))) 347 (if (zerop n) (rczero) 348 (let ((gcd (gcd n d))) 349 (cons (/ n gcd) (/ d gcd))))) 350 351(defun egcd (x y) 352 (let ((xn (abs (car x))) (xd (cdr x)) 353 (yn (abs (car y))) (yd (cdr y))) 354 (cons (gcd xn yn) (* xd (/ yd (gcd xd yd)))))) 355 356;;; Subtitle polynomial arithmetic 357 358(declare-top (special vars)) 359 360(defun ord-vector (p) 361 (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum))) tlist))) 362 (declare (special vars)) 363 (cond ((not (cdr vars)) (ncons (ps-le* p))) 364 (t (ord-vect1 p) (mapcar #'(lambda (x) (or (cdr x) (rczero))) vars))))) 365 366(defun ord-vect1 (p) 367 (declare (special vars)) 368 (unless (pscoefp p) 369 (let ((data (assoc (gvar p) vars :test #'eq)) 370 (le (ps-le p))) 371 (rplacd data (cond ((not (cdr data)) le) 372 (t (emin (cdr data) le)))) 373 (mapl #'(lambda (l) (ord-vect1 (lc l))) (terms p))))) 374 375(defun trunc-vector (p min?) 376 (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum))) tlist))) 377 (declare (special vars)) 378 (if (null (cdr vars)) (ncons (if (psp p) (trunc-lvl p) () )) 379 (progn 380 (trunc-vect1 p min?) 381 (mapcar 'cdr vars))))) 382 383(defun trunc-vect1 (p min?) 384 (declare (special vars)) 385 (unless (pscoefp p) 386 (let ((data (assoc (gvar p) vars :test #'eq)) 387 (trunc (trunc-lvl p))) 388 (when trunc 389 (rplacd data (if (null (cdr data)) trunc 390 (if min? (emin (cdr data) trunc) 391 (emax (cdr data) trunc)))))) 392 (dolist (term (terms p)) 393 (trunc-vect1 (c term) min?)))) 394 395(declare-top (unspecial vars)) 396 397(defun psplus (x y) 398 (cond ((pscoefp x) 399 (cond ((pscoefp y) (rcplus x y)) 400 ((rczerop x) y) 401 (t (pscplus x y)))) 402 ((pscoefp y) (if (rczerop y) x (pscplus y x))) 403 ((eqgvar (gvar-o x) (gvar-o y)) (psplus1 x y)) 404 ((pointerp (gvar-o x) (gvar-o y)) (pscplus y x)) 405 (t (pscplus x y)))) 406 407(defun rcplus! (x y) 408 (if (not (and least_term? taylor_simplifier)) (rcplus x y) 409 (prep1 (funcall taylor_simplifier (m+ (rcdisrep x) (rcdisrep y)))))) 410 411(defun psdiff (x y) 412 (cond ((pscoefp x) (cond ((pscoefp y) (rcdiff x y)) 413 ((rczerop x) (pstimes (rcmone) y)) 414 (t (pscdiff x y () )))) 415 ((pscoefp y) (if (rczerop y) x (pscdiff y x t))) 416 ((eqgvar (gvar-o x) (gvar-o y)) (psdiff1 x y)) 417 ((pointerp (gvar-o x) (gvar-o y)) (pscdiff y x t)) 418 (t (pscdiff x y () )))) 419 420(defun rcdiff! (x y) 421 (if (not (and least_term? taylor_simplifier)) (rcdiff x y) 422 (prep1 (funcall taylor_simplifier (m- (rcdisrep x) (rcdisrep y)))))) 423 424(defun psplus1 (x y) 425 (let ((ans (cons () () ))) 426 (psplus2 (gvar-o x) (emin (trunc-lvl x) (trunc-lvl y)) 427 (cons 0 (terms x)) (cons 0 (terms y)) ans ans))) 428 429(defun pscplus (c p) 430 (if (e> (rczero) (trunc-lvl p)) p 431 (pscheck (gvar-o p) (poly-data p) (pscplus1 c (terms p))))) 432 433(defun pscdiff (c p fl) 434 (if (e> (rczero) (trunc-lvl p)) 435 (if fl p (psminus p)) 436 (pscheck (gvar-o p) (poly-data p) 437 (cond ((not fl) (pscplus1 c (psminus-terms (terms p)))) 438 (t (pscplus1 (psminus c) (terms p))))))) 439 440(defun strip-zeroes (terms ps?) 441 (cond ((or (null terms) (null taylor_simplifier)) terms) 442 ((null ps?) 443 (do ((terms terms (n-term terms))) 444 ((null terms) () ) 445 (change-lc terms (strip-zeroes (lc terms) 't)) 446 (unless (rczerop (lc terms)) (return terms)))) 447 ((pscoefp terms) 448 (if (null taylor_simplifier) terms 449 (let ((exp (rcdisrep terms))) 450 ;; If a pscoeff is not free of tvars then the ps is a 451 ;; multivar series and we can't handle a recursive 452 ;; call to taylor (as opposed to a call to prep1, as below) 453 ;; because this would be circuler (e.g. try 454 ;; taylor(x/ (x^2+1),[x],%i,-1) ). Besides, in this case 455 ;; the pscoeff contains a tvar hence should not be 0. 456 (if (not (mfree exp tvars)) terms 457 (prep1 (funcall taylor_simplifier exp)))))) 458 (t (pscheck (gvar-o terms) (poly-data terms) 459 (strip-zeroes (terms terms) () ))))) 460 461(defun pscplus1 (c l) 462 (cond ((null l) (list (term (rczero) c))) 463 ((rczerop (le l)) (setq c (psplus c (lc l))) 464 (if (rczerop c) (strip-zeroes (n-term l) () ) 465 (cons (term (rczero) c) (n-term l)))) 466 ((e> (le l) (rczero)) (cons (term (rczero) c) l)) 467 (t (cons (lt l) (let ((least_term?)) (pscplus1 c (n-term l))))))) 468 469;;; Both here and in psdiff2 xx and yy point one before where one 470;;; might think they should point so that extensions will be retained. 471 472(defun psplus2 (varh trunc xx yy ans a) 473 (prog (c) 474 a (cond ((mono-term? xx) 475 (if (mono-term? yy) (go end) (go null))) 476 ((mono-term? yy) (setq yy xx) (go null))) 477 (cond ((equal (le (n-term xx)) (le (n-term yy))) 478 (setq xx (n-term xx) yy (n-term yy)) 479 (setq c (let ((least_term? (null (n-term ans)))) 480 (psplus (lc xx) (lc yy)))) 481 (if (rczerop c) (go a) (add-term a (le xx) c))) 482 ((e> (le (n-term xx)) (le (n-term yy))) 483 (setq yy (n-term yy)) 484 (add-term a (lt yy))) 485 (t (setq xx (n-term xx)) 486 (add-term a (lt xx)))) 487 (setq a (n-term a)) 488 (go a) 489 null (if (or (mono-term? yy) (e> (le (n-term yy)) trunc)) 490 (go end) 491 (progn 492 (setq yy (n-term yy)) 493 (add-term-&-pop a (lt yy)) 494 (go null))) 495 end (return (pscheck varh (list trunc) (cdr ans))))) 496 497(defun psdiff1 (x y) 498 (let ((ans (cons () () ))) 499 (psdiff2 (gvar-o x) (emin (trunc-lvl x) (trunc-lvl y)) 500 (cons 0 (terms x)) (cons 0 (terms y)) ans ans))) 501 502(defun psdiff2 (varh trunc xx yy ans a) 503 (prog (c) 504 a (cond ((mono-term? xx) 505 (if (mono-term? yy) 506 (go end) 507 (progn 508 (setq yy 509 (cons 0 (mapcar #'(lambda (q) 510 (term (e q) (psminus (c q)))) 511 (cdr yy)))) 512 (go null)))) 513 ((mono-term? yy) 514 (setq yy xx) (go null))) 515 (cond ((equal (le (n-term xx)) (le (n-term yy))) 516 (setq xx (n-term xx) yy (n-term yy)) 517 (setq c (let ((least_term? (null (n-term ans)))) 518 (psdiff (lc xx) (lc yy)))) 519 (if (rczerop c) (go a) 520 (add-term a (le xx) c))) 521 ((e> (le (n-term xx)) (le (n-term yy))) 522 (setq yy (n-term yy)) 523 (add-term a (le yy) (psminus (lc yy)))) 524 (t (setq xx (n-term xx)) 525 (add-term a (lt xx)))) 526 (setq a (n-term a)) 527 (go a) 528 null (if (or (mono-term? yy) (e> (le (n-term yy)) trunc)) 529 (go end) 530 (progn 531 (setq yy (n-term yy)) 532 (add-term-&-pop a (le yy) (lc yy)) 533 (go null))) 534 end (return (pscheck varh (list trunc) (cdr ans))))) 535 536(defun psminus (x) 537 (if (psp x) (make-ps x (psminus-terms (terms x))) 538 (rcminus x))) 539 540(defun psminus-terms (terms) 541 (let ((ans (cons () () ))) 542 (do ((q terms (n-term q)) 543 (a ans (cdr a))) 544 ((null q) (cdr ans)) 545 (add-term a (le q) (psminus (lc q)))))) 546 547(defun pscheck (a b terms) 548 (cond ((null terms) (rczero)) 549 ((and (mono-term? terms) (rczerop (le terms))) 550 (lc terms)) 551 (t (make-ps a b terms)))) 552 553(defun pstrim-terms (terms e) 554 (do () (()) 555 (cond ((null terms) (return () )) 556 ((null (e> e (le terms))) (return terms)) 557 (t (setq terms (n-term terms)))))) 558 559(defun psterm (terms e) 560 (psterm1 (pstrim-terms terms e) e)) 561 562(defun psterm1 (l e) 563 (cond ((null l) (rczero)) 564 ((e= (le l) e) (lc l)) 565 (t (rczero)))) 566 567(defun pscoeff1 (a b c) ;a is an mrat!!! 568 (let ((tlist (mrat-tlist a))) 569 (cons (nconc (list 'mrat 'simp (mrat-varlist a) (mrat-genvar a)) 570 (do ((l (mrat-tlist a) (cdr l)) 571 (ans () (cons (car l) ans))) 572 ((null l) ans) 573 (when (alike1 (caar l) b) 574 (return 575 (and (or ans (cdr l)) 576 (list (nreconc ans (cdr l)) 'trunc)))))) 577 (pscoef (mrat-ps a) (int-gvar (get-datum b)) (prep1 c))))) 578 579(defun pscoef (a b c) 580 (cond ((pscoefp a) (if (rczerop c) a (rczero))) 581 ((eq b (gvar a)) (psterm (terms a) c)) 582 (t (do ((gvar-o (gvar-o a)) 583 (poly-data (poly-data a)) 584 (ans (rczero)) 585 (terms (terms a) (n-term terms)) 586 (temp)) 587 ((null terms) ans) 588 (unless (rczerop (setq temp (pscoef (lc terms) b c))) 589 (setq ans (psplus ans 590 (make-ps gvar-o poly-data 591 (ncons (term (le terms) 592 temp)))))))))) 593 594(defun psdisextend (p) 595 (cond ((not (psp p)) p) 596 (t (make-ps p (mapcar #'(lambda (q) (cons (car q) (psdisextend (cdr q)))) 597 (terms p)))))) 598 599(defun psfloat (p) 600 (if (psp p) (psfloat1 p (trunc-lvl p) (terms p) (ncons 0)) 601 (rctimes (rcfone) p))) 602 603(defun psfloat1 (p trunc l ans) 604 (do (($float 't) 605 (a (last ans) (n-term a))) 606 ((or (null l) (e> (le l) trunc)) 607 (pscheck (gvar-o p) (poly-data p) (cdr ans))) 608 (add-term a (le l) (psfloat (lc l))) 609 (setq l (n-term l)))) 610 611(defun pstrunc (p) 612 (pstrunc1 p (mapcar #'(lambda (q) (cons (int-gvar q) (current-trunc q))) 613 tlist))) 614 615(defun pstrunc1 (p trlist) 616 (cond ((not (psp p)) 617 p) 618 (t 619 (let ((trnc (cdr (assoc (gvar p) trlist :test #'eq))) (trunc-ps) (a nil)) 620 (do ((l (terms p) (n-term l))) 621 ((null l) (pscheck (gvar-o p) (ncons (trunc-lvl p)) (nreverse a))) 622 (when (e> (le l) trnc) 623 (return (pscheck (gvar-o p) (ncons trnc) (nreverse a)))) 624 (unless (rczerop (setq trunc-ps (pstrunc1 (lc l) trlist))) 625 (push (term (le l) trunc-ps) a))))))) 626 627(defun pstimes (x y) 628 (cond ((or (rczerop x) (rczerop y)) (rczero)) 629 ((pscoefp x) (cond ((pscoefp y) (rctimes x y)) 630 ((equal x (rcone)) y) 631 (t (psctimes* x y)))) 632 ((pscoefp y) (if (equal y (rcone)) x (psctimes* y x))) 633 ((eqgvar (gvar-o x) (gvar-o y)) (pstimes*1 x y)) 634 ((pointerp (gvar-o x) (gvar-o y)) (psctimes* y x)) 635 (t (psctimes* x y)))) 636 637(defun psctimes* (c p) 638 (make-ps p (maplist #'(lambda (l) 639 (term (le l) (pstimes c (lc l)))) 640 (terms p)))) 641 642(defun pstimes*1 (xa ya) 643 (let ((ans (cons () () )) 644 (trunc (let ((lex (ps-le xa)) (ley (ps-le ya))) 645 (e+ (emin (e- (trunc-lvl xa) lex) (e- (trunc-lvl ya) ley)) 646 (e+ lex ley))))) 647 (unless $maxtayorder 648 (setq trunc (emin trunc (t-o-var (gvar xa))))) 649 (pstimes*2 xa ya trunc ans))) 650 651(defun pstimes*2 (xa ya trunc ans) 652 (prog (a c e x y yy) 653 (setq x (terms xa) y (setq yy (terms ya)) a ans) 654 a (cond ((or (null y) (e> (setq e (e+ (le x) (le y))) trunc)) 655 (go b)) 656 ((not (rczerop (setq c (pstimes (lc x) (lc y))))) 657 (add-term-&-pop a e c))) 658 (setq y (n-term y)) 659 (go a) 660 b (unless (setq x (n-term x)) 661 (return (pscheck (gvar-o xa) (list trunc) (cdr ans)))) 662 (setq y yy a ans) 663 c (when (or (null y) (e> (setq e (e+ (le x) (le y))) trunc)) 664 (go b)) 665 (setq c (pstimes (lc x) (lc y))) 666 d (cond ((or (mono-term? a) (e> (le (n-term a)) e)) 667 (add-term-&-pop a e c)) 668 ((e> e (le (n-term a))) 669 (setq a (n-term a)) 670 (go d)) 671 (t (setq c (psplus c (lc (n-term a)))) 672 (if (rczerop c) 673 (rplacd a (n-term (n-term a))) 674 (progn 675 (change-lc (n-term a) c) 676 (setq a (n-term a)))))) 677 (setq y (n-term y)) 678 (go c))) 679 680(defun pscsubst (c v p) 681 (cond ((pscoefp p) p) 682 ((eq v (gvar p)) (pscsubst1 c p)) 683 ((pointerp v (gvar p)) p) 684 (t (make-ps p (maplist 685 #'(lambda (q) (term (le q) 686 (pscsubst c v (lc q)))) 687 (terms p)))))) 688 689(defun pscsubst1 (v u) 690 (do ((a (rczero)) 691 (ul (terms u) (n-term ul))) 692 ((null ul) a) 693 (setq a (psplus a (pstimes (lc ul) (psexpt v (le ul))))))) 694 695(defun get-series (func trunc var e c) 696 (let ((pw (e// trunc e))) 697 (setq e (if (and (equal e (rcone)) (equal c (rcone))) 698 (getexp-fun func var pw) 699 (psmonsubst (getexp-fun func var pw) trunc e c))) 700 (if (and $float $keepfloat) (psfloat e) e))) 701 702(defun psmonsubst (p trunc e c) 703 (if (psp p) 704 (psmonsubst1 p trunc e c 705 `(() . ,(terms p)) (ncons () ) (rcone) (rczero)) 706 p)) 707 708 709(defun psmonsubst1 (p trunc e c l ans cc el) 710 ;; We set $MAXTAYORDER to () here so that the calls to psexpt below 711 ;; won't do needless extra work, e.g. see rwg's complaint of 9/7/82. 712 (prog (a ee varh $maxtayorder) 713 (setq a ans varh (gvar-o p)) 714 a (cond ((or (mono-term? l) 715 (e> (setq ee (e* e (le (n-term l)))) trunc)) 716 (go end)) 717 ((rczerop (setq cc 718 (pstimes cc 719 (psexpt c (e- (le (setq l (n-term l))) 720 el)))))) 721 ((mono-term? a) 722 (add-term a ee (pstimes cc (lc l))))) 723 (setq a (n-term a) el (le l)) 724 (go a) 725 end (return (pscheck varh (list trunc) (cdr ans))))) 726 727(defun psexpon-gcd (terms) 728 (do ((gcd (le terms) (egcd (le l) gcd)) 729 (l (n-term terms) (n-term l))) 730 ((null l) gcd))) 731 732(defun psfind-s (p) 733 (if (psp p) (psfind-s (psterm (terms p) (rczero))) 734 (psfind-s1 p))) 735 736(defun psfind-s1 (r) 737 (cond ((null (atom (cdr r))) (rczero)) 738 ((atom (car r)) r) 739 (t (do ((p (ptterm (cdar r) 0) (ptterm (cdr p) 0))) 740 ((atom p) (cons p (cdr r))))))) 741 742(defun psexpt (p n) 743 (cond ((rczerop n) ;; p^0 744 (if (rczerop p) ;; 0^0 745 (merror (intl:gettext "taylor: 0^0 is undefined.")) 746 (rcone))) ;; Otherwise can let p^0 = 1 747 ((or (equal n (rcone)) (equal n (rcfone))) p) ;; p^1 cases 748 ((pscoefp p) (rcexpt p n)) 749 ((mono-term? (terms p)) ;; A monomial to a power 750 (let ((s (psfind-s n)) (n-s) (x) (l (terms p))) 751 ;; s is the numeric part of the exponent 752 (if (floatp (car s)) ;; Perhaps we souldn't 753 ;; rationalize if $keepfloat is true? 754 (setq s (maxima-rationalize (quot (car s) (cdr s))))) 755 (setq n-s (psdiff n s) ;; the non-numeric part of exponent 756 x (e* s (le l))) ;; the degree of the lowest term 757 (setq x (if (and (null $maxtayorder) ;; if not getting all terms 758 (e> x (t-o-var (gvar p)))) 759 ;; and result is of high order 760 (rczero) ;; then zero is enough 761 (pscheck (gvar-o p) ;; otherwise 762 (ncons (e+ (trunc-lvl p) ;; new trunc-level 763 (e- x (le l)))) ;; kick exponent 764 (ncons (term x (psexpt (lc l) n)))))) 765 ;; x is now p^s 766 (if (or (rczerop n-s) (rczerop x)) ;; is that good enough? 767 x ;; yes! The rest is bletcherous. 768 (pstimes x (psexpt (prep1 (m^ (get-inverse (gvar p)) 769 (rcdisrep n-s))) 770 (ps-le p)))))) 771 (t (prog (l lc le inc trunc s lt mr lim lcinv ans) 772 (setq lc (lc (setq l (terms p))) 773 le (le l) lt (lt l) trunc (trunc-lvl p) 774 inc (psexpon-gcd l) s (psfind-s n)) 775 (when (floatp (car s)) 776 (setq s (maxima-rationalize (quot (car s) (cdr s))))) 777 (setq ans (psexpt (setq lt (pscheck (gvar-o p) (list trunc) 778 (list lt))) n) 779 lcinv (psexpt lc (rcmone)) 780 mr (e+ inc (e* s le)) 781 lim (if (and (infp trunc) (not (e> s (rczero)))) 782 (t-o-var (gvar p)) 783 ;; See the comment in PSEXPT1 below which tells 784 ;; why we don't allow inf. trunc's here. 785 (e+ (if (and (infp trunc) (not (rcintegerp s))) 786 (if (infp (setq lim (t-o-var (gvar p)))) 787 (infin-ord-err) 788 lim) 789 trunc) 790 (e* (e1- s) le))) 791 ans 792 (if (or (pscoefp ans) (null (eq (gvar p) (gvar ans)))) 793 (list 0 (term (rczero) ans)) 794 (cons 0 (terms ans)))) 795 (and (null $maxtayorder) 796 (or (not (infp lim)) 797 (not (rcintegerp s)) 798 (e> (e* s (le (last l))) (t-o-var (gvar p)))) 799 (setq lim (emin lim (t-o-var (gvar p))))) 800 ;;(and (infp lim) (n-term l) (e> (rczero) n) 801 ;; (infin-ord-err)) 802 (return (psexpt1 (gvar-o p) 803 lim l n s inc 1 mr ans le lcinv)))))) 804 805(defun psexpt1 (varh trunc l n s inc m mr ans r linv) 806 ;; n is the power we are raising the series to 807 ;; inc is the exponent increment 808 ;; mr is the current exponent 809 ;; tr is the truncation level desired 810 ;; m is the term index 811 (declare (fixnum m)) 812 (prog (a (k 0) ak cm-k c ma0 sum kr tr) 813 (declare (fixnum k)) 814 ;; truly unfortunate that we need so many variables in this hack 815 (setq a (last ans) tr trunc) 816 ;; I don't see what's wrong with truncating exact series when 817 ;; raising them to fractional powers so we'll allow it for now. 818 ;; This is accomplished above in PSEXPT (see the comment). Thus, 819 ;; presumably, this check should never be needed anymore. 820 ;; Bugs catching this clause were sqrt(1-x)*taylor(f1,x,0,0) 821 ;; and sqrt(taylor(x+x^2,x,0,2)),taylor_truncate_polynomials=false. 822 (when (infp tr) 823 (if (rcintegerp s) 824 (setq tr (e* s (le (last l)))) 825 (merror (intl:gettext "taylor: expected an integer, instead found: ~:M") s))) 826 (when (infp tr) (setq tr (t-o-var (car varh)))) 827 b (and (e> mr tr) (go end)) 828 (setq kr inc ak l ma0 (pstimes (cons 1 m) linv) 829 k 1 sum (rczero)) 830 a (if (or (> k m) (null (setq cm-k (psterm (cdr ans) (e- mr kr))))) 831 (go add-term)) 832 (setq ak (or (pstrim-terms ak (e+ kr r)) (go add-term)) 833 c (pstimes (psdiff (pstimes (cons k 1) n) 834 (cons (- m k) 1)) 835 (pstimes (if (e= (e+ kr r) (le ak)) 836 (lc ak) 837 (rczero)) 838 cm-k))) 839 (setq sum (psplus sum c) 840 k (1+ k) kr (e+ kr inc)) 841 (go a) 842 add-term 843 (and (null (rczerop sum)) 844 (add-term-&-pop a mr (pstimes ma0 sum))) 845 (setq m (1+ m) mr (e+ mr inc)) 846 (go b) 847 end (return (pscheck varh (list trunc) (cdr ans))))) 848 849(defun psderivative (p v) 850 (cond ((pscoefp p) (rcderiv p v)) 851 ((eq v (gvar p)) 852 (if (prog1 (rczerop (ps-le p)) 853 (setq p (psderiv1 (gvar-o p) 854 (trunc-lvl p) (cons 0 (terms p)) (list 0)))) 855 (strip-zeroes p 't) p)) 856 (t (psderiv2 (gvar-o p) 857 (trunc-lvl p) v (cons 0 (terms p)) (list 0))))) 858 859(defun psderiv1 (varh trunc l ans) 860 (do ((a (last ans))) 861 ((or (mono-term? l) (e> (le (n-term l)) trunc)) 862 (pscheck varh (list (e1- trunc)) (cdr ans))) 863 (setq l (n-term l)) 864 (when (not (rczerop (le l))) 865 (add-term-&-pop a (e1- (le l)) (pstimes (le l) (lc l)))))) 866 867(defun psderiv2 (varh trunc v l ans) 868 (do ((a (last ans) (n-term a)) (c)) 869 ((or (mono-term? l) (e> (le (n-term l)) trunc)) 870 (pscheck varh (list trunc) (cdr ans))) 871 (setq l (n-term l)) 872 (or (rczerop (setq c (psderivative (lc l) v))) 873 (add-term a (le l) c)))) 874 875(defun psdp (p) 876 (let (temp temp2) 877 (cond ((pscoefp p) (rcderivx p)) 878 ((or (rczerop (setq temp (getdiff (gvar-o p)))) 879 (eq (car temp) 'multi)) 880 (setq temp2 (psdp2 (gvar-o p) (trunc-lvl p) 881 (cons 0 (terms p)) (list 0))) 882 (if (eq (car temp) 'multi) 883 (pstimes temp2 884 (make-ps (gvar-o p) (ncons (inf)) 885 (list (term (cdr temp) (rcone))))) 886 temp2)) 887 (t (psdp1 (gvar-o p) 888 (trunc-lvl p) (cons 0 (terms p)) 889 (list 0) temp))))) 890 891(defun psdp1 (varh trunc l ans dx) 892 (do ((a (last ans)) (c (rczero))) 893 ((or (mono-term? l) (e> (le (n-term l)) trunc)) 894 (psplus c (pscheck varh (list (e1- trunc)) (cdr ans)))) 895 (setq l (n-term l)) 896 (if (rczerop (le l)) (setq c (psdp (lc l))) 897 (add-term-&-pop 898 a (e1- (le l)) (pstimes (le l) (pstimes dx (lc l))))))) 899 900(defun psdp2 (varh trunc l ans) 901 (do ((a (last ans)) (c)) 902 ((or (mono-term? l) (e> (le (n-term l)) trunc)) 903 (pscheck varh (list trunc) (cdr ans))) 904 (setq l (n-term l)) 905 (when (null (rczerop (setq c (psdp (lc l))))) 906 (add-term-&-pop a (le l) c)))) 907 908;;; Currently unused 909;;; 910;;; (defun psintegrate (p v) 911;;; (cond ((rczerop p) (rczero)) 912;;; ((pscoefp p) 913;;; (pstimes p (taylor2 (get-inverse (car v))))) 914;;; ((eqgvar v (gvar-o p)) 915;;; (psinteg1 (gvar-o p) 916;;; (trunc-lvl p) (cons 0 (terms p)) (list 0))) 917;;; (t (psinteg2 (gvar-o p) 918;;; (trunc-lvl p) v (cons 0 (terms p)) (list 0))))) 919;;; 920;;; (defun psinteg1 (varh trunc l ans) 921;;; (prog (a) 922;;; (setq a (last ans)) 923;;; a (if (or (null (n-term l)) (e> (le (n-term l)) trunc)) 924;;; (go end) 925;;; (add-term a (e1+ (le (setq l (n-term l)))) 926;;; (pstimes (le l) 927;;; (if (e= (le l) (rcmone)) 928;;; (prep1 (list '(%LOG) 929;;; (get-inverse 930;;; (car varh)))) 931;;; (lc l)))) 932;;; (setq a (n-term a))) 933;;; (go a) 934;;; end (return (pscheck varh (list (e1+ trunc)) (cdr ans))))) 935 936;;; (defun psinteg2 (varh trunc v l ans) 937;;; (prog (a) 938;;; (setq a (last ans)) 939;;; a (if (or (null (n-term l)) (e> (le (n-term l)) trunc)) 940;;; (go end) 941;;; (add-term a (le l) 942;;; (psintegrate (lc (setq l (n-term l))) v)) 943;;; (setq a (n-term a))) 944;;; (go a) 945;;; end (return (pscheck varh (list trunc) (cdr ans))))) 946 947(defun psexpt-log-ord (p) 948 (cond ((null $maxtayorder) (emin (trunc-lvl p) (t-o-var (gvar p)))) 949 ((infp (trunc-lvl p)) (t-o-var (gvar p))) 950 (t (trunc-lvl p)))) 951 952;(defun ps-infp (p) 953; (if (pscoefp p) () 954; (get-) "...")) 955 956(defun psexpt-fn (p) 957 (let (ans ord<0?) 958 (cond ((pscoefp p) (psexpt-fn2 (rcdisrep p))) 959 ((ps-lim-infp p) (psexpt-fn-sing p)) 960 ((prog2 (setq ord<0? (e> (rczero) (ps-le p))) 961 (null (n-term (terms p)))) 962 (setq ans (get-series '%exp (psexpt-log-ord p) (gvar-o p) 963 (if ord<0? (e- (ps-le p)) (ps-le p)) 964 (ps-lc p))) 965 (if ord<0? (ps-invert-var ans) ans)) 966 ((if ord<0? 967 (when (e= (rczero) (e (setq ans (ps-gt p)))) 968 (pstimes (psexpt-fn (pscheck (gvar-o p) (list (trunc-lvl p)) 969 (delete ans (terms p) :test #'eq))) 970 (psexpt-fn2 (srdis (c ans))))) 971 (when (e= (rczero) (ps-le p)) 972 (pstimes (psexpt-fn2 (srdis (lc (terms p)))) 973 (psexpt-fn (pscheck (gvar-o p) (list (trunc-lvl p)) 974 (n-term (terms p))))))) ) 975 (t (prog (l inc trunc ea0 ans) 976 (setq l (terms p)) 977 (when ord<0? 978 ;(return (ps-invert-var (psexpt-fn (ps-invert-var p)))) 979 (setq l (invert-terms l))) 980 (setq trunc (trunc-lvl p) 981 inc (psexpon-gcd l) ea0 (rcone)) 982 (unless (e> (le l) (rczero)) 983 ;; MEANING OF FOLLOWING MESSAGE IS OBSCURE 984 (merror "PSEXPT-FN: unreachable point.")) 985 (setq ans 986 (if (or (pscoefp ea0) (null (eq (gvar p) (gvar ea0)))) 987 (list 0 (term (rczero) ea0)) 988 (cons 0 (terms ea0)))) 989 (unless $maxtayorder 990 (setq trunc (emin trunc (t-o-var (gvar p))))) 991 (when (infp trunc) (setq trunc (t-o-var (gvar p)))) 992 (setq ans (psexpt-fn1 (gvar-o p) trunc l inc 1 inc ans)) 993 (return (if ord<0? (ps-invert-var ans) ans))))))) 994 995(defun psexpt-fn-sing (p) 996 (let ((inf-var? (member (gvar-lim (gvar p)) '($inf $minf) :test #'eq)) 997 (c*logs (c*logs (lt-poly p))) c strongest-term) 998 ;; Must pull out out logs here: exp(ci*log(ui)+x) -> ui^ci*exp(x) 999 ;; since its much harder for adjoin-tvar to do this transformation 1000 ;; below after things have been disrepped. 1001 (setq c (exp-c*logs c*logs) p (psdiff p (sum-c*logs c*logs))) 1002 (if (not (ps-lim-infp p)) 1003 ;; Here we just subtracted the only infinite term, e.g. 1004 ;; p = 1/2*log(x)+1/log(x)+... 1005 (pstimes c (psexpt-fn p)) 1006 (progn 1007 (setq strongest-term (if inf-var? (ps-gt p) (ps-lt p))) 1008 ;; If the strongest term has degree 0 in the mainvar then the singular 1009 ;; terms occur in some other weaker var. There may be terms in this 1010 ;; coef which arent singular (e.g. 1 in (1/x+1+...)+exp(-1/x)+...) so 1011 ;; we must recursively psexpt-fn this term to get only what we need. 1012 (if (rczerop (e strongest-term)) 1013 (setq c (pstimes c (psexpt-fn (c strongest-term)))) 1014 (dolist (exp (expand-and-disrep strongest-term p)) 1015 (setq c (pstimes c (adjoin-tvar (m^ '$%e exp)))))) 1016 (pstimes c (psexpt-fn (pscheck (gvar-o p) (list (trunc-lvl p)) 1017 (if inf-var? 1018 (delete strongest-term (terms p) :test #'eq) 1019 (n-term (terms p)))))))))) 1020 1021(defun gvar-logp (gvar) 1022 (let ((var (get-inverse gvar))) 1023 (and (consp var) (eq (caar var) 'mexpt) (equal (caddr var) -1) 1024 (consp (setq var (cadr var))) (eq (caar var) '%log) 1025 var))) 1026 1027(defun c*logs (p) 1028 (if (pscoefp p) () 1029 (let ((log (gvar-logp (gvar p))) c) 1030 (if (not log) 1031 () 1032 (progn 1033 (setq c (psconst (psterm (terms p) (rcmone)))) 1034 ;; We don't want e.g. exp(x^a*log(x)) -> x^x^a 1035 (if (not (mfree (rcdisrep c) tvars)) () 1036 (cons (cons c (cons log p)) 1037 (c*logs (psterm (terms p) (rczero)))))))))) 1038 1039(defun psconst (p) 1040 (if (pscoefp p) p (psconst (psterm (terms p) (rczero))))) 1041 1042(defun exp-c*logs (c*logs) 1043 (if (null c*logs) (rcone) 1044 (pstimes (taylor2 `((mexpt) ,(cadr (cadr (car c*logs))) 1045 ,(rcdisrep (caar c*logs)))) 1046 (exp-c*logs (cdr c*logs))))) 1047 1048(defun sum-c*logs (c*logs) 1049 (if (null c*logs) (rczero) 1050 (let ((ps (cddr (car c*logs)))) 1051 (psplus (make-ps ps (ncons (term (ps-le ps) (caar c*logs)))) 1052 (sum-c*logs (cdr c*logs)))))) 1053 1054;; Calculatest the limit of a series at the expansion point. Returns one of 1055;; {$zeroa, $zerob, $pos, $neg, $inf, $minf}. 1056 1057(defvar tvar-limits () 1058 "A list of the form ((gvar . limit(gvar)) ...)") 1059 1060(defun ps-lim-infp (ps) 1061 (if (pscoefp ps) () 1062 ;; Assume taylor vars at 0+ for now. Should handle the cases when 1063 ;; the expansion point is INF, MINF,etc. 1064 (let* ((lim (gvar-lim (gvar ps))) 1065 (strongest-term 1066 (if (member lim '($inf $minf) :test #'eq) (ps-gt ps) (ps-lt ps)))) 1067 (if (ezerop (e strongest-term)) 1068 (ps-lim-infp (c strongest-term)) 1069 (progn 1070 (setq lim (lim-power lim (e strongest-term))) 1071 (and (lim-infp lim) (not (eq lim '$infinity)))))))) 1072 1073(defun lim-zerop (lim) 1074 (member lim '($zeroa $zerob $zeroim) :test #'eq)) 1075 1076(defun lim-plusp (lim) 1077 (member lim '($zeroa $pos $inf $finite) :test #'eq)) 1078 1079(defun lim-finitep (lim) 1080 (member lim '($pos $neg $im $finite) :test #'eq)) 1081 1082(defun lim-infp (lim) 1083 (member lim '($inf $minf $infinity) :test #'eq)) 1084 1085(defun lim-imagp (lim) 1086 (member lim '($im $infinity) :test #'eq)) 1087 1088(defun lim-minus (lim) 1089 (cdr (assoc lim '(($zeroa . $zerob) ($zerob . $zeroa) ($pos . $neg) ($zero . $zero) 1090 ($neg . $pos) ($inf . $minf) ($minf . $inf) 1091 ($im . $im) ($infinity . $infinity) ($finite . $finite)) :test #'eq))) 1092(defun lim-abs (lim) 1093 (or (cdr (assoc lim '(($zerob . $zeroa) ($neg . $pos) ($minf . $inf)) :test #'eq)) 1094 lim)) 1095 1096(defun lim-times (lim1 lim2) 1097 (let (lim) 1098 (cond ((or (eq lim1 '$zero) (eq lim2 '$zero)) (setq lim '$zero)) 1099 ((and (lim-infp lim1) (lim-infp lim2)) (setq lim '$inf)) 1100 ((and (lim-zerop lim1) (lim-zerop lim2)) (setq lim '$pos)) 1101 ((or (when (lim-finitep lim2) (rotatef lim1 lim2) 't) 1102 (lim-finitep lim1)) 1103 (when (and (eq lim1 '$finite) (lim-infp lim1)) 1104 (break "Undefined finite*inf in lim-times")) 1105 (setq lim (lim-abs lim2))) 1106 (t (break "Undefined limit product ~A * ~A in lim-times" lim1 lim2))) 1107 (if (or (lim-imagp lim1) (lim-imagp lim2)) 1108 (if (lim-infp lim) '$infinity '$im) 1109 (if (and (lim-plusp lim1) (lim-plusp lim2)) lim (lim-minus lim))))) 1110 1111(defun lim-power (lim power) 1112 (cond ((ezerop power) '$pos) 1113 ((e> (rczero) power) (lim-recip (lim-power lim (e- power)))) 1114 ((not (oddp (car power))) 1115 (if (lim-plusp lim) lim (lim-minus lim))) 1116 (t lim))) 1117 1118(defun lim-recip (lim) 1119 (or (cdr (assoc lim '(($zeroa . $inf) ($zerob . $minf) 1120 ($inf . $zeroa) ($minf . $zerob)) :test #'eq)) 1121 (if (eq lim '$finite) (break "inverting $finite?") 1122 lim))) 1123 1124(defun lim-exp (lim) 1125 (case lim 1126 (($zeroa $zerob $zero $pos $neg $minf) '$zeroa) 1127 (($inf $finite) lim) 1128 ($infinity '$infinity) ; actually only if Re lim = inf 1129 (t (break "Unhandled limit in lim-exp")))) 1130 1131(defun lim-log (lim) 1132 (case lim 1133 ($zeroa '$minf) 1134 ($inf '$inf) 1135 ($minf '$infinity) 1136 ($zerob '$infinity) 1137 (t (break "Unhandled limit in lim-log")))) 1138 1139(defun expand-and-disrep (term p) 1140 (let ((x^n (list '(mexpt) (get-inverse (gvar p)) (edisrep (e term)))) 1141 (a (c term))) 1142 (if (pscoefp a) (ncons (m* (srdis a) x^n)) 1143 (mapcar #'(lambda (subterm) 1144 (m* (cons '(mtimes) (expand-and-disrep subterm a)) x^n)) 1145 (terms a))))) 1146 1147(defun adjoin-sing-datum (d) 1148 (let ((r (prep1 (datum-var d))) (g (gensym)) (kernel (datum-var d)) 1149 (no (1+ (cdr (int-var (car (last tlist))))))) 1150 (unless (and (equal (car r) 1) (equal (cddr r) '(1 1))) 1151 (break "bad singular datum")) 1152 (putprop g kernel 'disrep) 1153 (rplacd (cdddr d) (cons g no)) 1154 (adjoin-datum d) 1155 (push (cons (cadr r) kernel) key-vars) 1156 (push (cons g kernel) key-vars) 1157 (push (car key-vars) ivars) 1158 ;(push (cons kernel (cons (pget g) 1)) genpairs) 1159 (push (cons g (exp-pt d)) tvar-limits))) 1160 1161(defun adjoin-tvar (exp) (rat->ps (prep1 exp))) 1162 1163(defun rat->ps (rat) 1164 (pstimes (poly->ps (car rat)) 1165 (psexpt (poly->ps (cdr rat)) (rcmone)))) 1166 1167(defun poly->ps (poly) 1168 (if (or (pcoefp poly) (mfree (pdis poly) tvars)) (prep1 poly) 1169 (let ((g (p-var poly)) datum (pow (rcone))) 1170 (if (setq datum (key-var-pow g)) (desetq (g . pow) datum) 1171 (desetq (g . pow) (adjoin-pvar g))) 1172 (if (and (not (atom g)) (psp g)) 1173 g 1174 (progn 1175 (setq datum (gvar-data g)) 1176 (do ((po-terms (p-terms poly) (p-red po-terms)) 1177 (ps-terms () 1178 (push (term (e* pow (prep1 (pt-le po-terms))) 1179 (poly->ps (pt-lc po-terms))) 1180 ps-terms))) 1181 ((null po-terms) 1182 ;; This must be exact so that when we invert in rat-ps above 1183 ;; we don't lose any terms. E.g. try 1184 ;; taylor(log(1+exp(-1/x)),x,0,5). When taylor2'ing exp(-1/x), 1185 ;; if you used current trunc here this would return exp(1/x)...5 1186 ;; which would then be trunc'd to degree 3 by psexpt. 1187 (make-ps (int-var datum) 1188 (ncons (current-trunc datum)) 1189 (if (eq g (data-gvar datum)) ps-terms 1190 (invert-terms ps-terms)))))))))) 1191 1192(defun key-var-pow (g) 1193 (let ((var (get-key-var g)) datum) 1194 (when var 1195 (setq datum (get-datum var)) 1196 (if (eq g (setq g (data-gvar datum))) (cons g (rcone)) 1197 (cons g (rcmone)))))) 1198 1199(defun adjoin-pvar (g) 1200 (let ((kernel (get g 'disrep)) g* lim datum ans 1201 (no (1+ (cdr (int-var (car (last tlist)))))) (pow (rcone)) expt) 1202 (when (assol kernel tlist) (break "bad1")) 1203 (if (and (eq (caar kernel) 'mexpt) (eq (cadr kernel) '$%e) 1204 (not (atom (setq expt (caddr kernel)))) 1205 (eq (caar expt) 'mtimes) (not (mfree expt (ncons '$%i)))) 1206 (destructuring-let (((rpart . ipart) (trisplit expt))) 1207 (cons (pstimes (prep1 (m^ '$%e rpart)) 1208 (psplus (adjoin-tvar `((%cos) ,ipart)) 1209 (pstimes (prep1 '$%i) 1210 (adjoin-tvar `((%sin) ,ipart))))) 1211 pow)) 1212 (progn 1213 (when (eq (caar kernel) 'mexpt) 1214 (when (and (not (atom (setq expt (caddr kernel)))) 1215 (eq (caar expt) 'mtimes) 1216 ($ratnump (cadr expt))) 1217 (setq pow (cadr expt) kernel (m^ kernel (m// pow)) 1218 g (prep1 kernel) pow (prep1 pow)) 1219 (unless (and (equal (cdr g) 1) (equal (cdar g) '(1 1))) 1220 (break "Illegal kernel in `adjoin-pvar'")) 1221 (setq g (caar g) kernel (get g 'disrep)))) 1222 (if (setq ans (key-var-pow g)) 1223 (cons (car ans) (e* pow (cdr ans))) 1224 (progn 1225 (when (lim-infp (or lim (setq lim (tvar-lim kernel)))) 1226 (setq g* g g (gensym) kernel (m// kernel) 1227 lim (lim-recip lim) pow (e* (rcmone) pow)) 1228 (putprop g kernel 'disrep) 1229 ;(push g genvar) (push kernel varlist) 1230 (push (cons g* kernel) key-vars)) 1231 (when (assol kernel tlist) (break "bad2")) 1232 (setq datum (list* kernel 1233 ;(mapcar #'(lambda (e) (emax e (rczero))) 1234 ; (trunc-stack (car tlist))) 1235 (copy-list (trunc-stack (car tlist))) 1236 lim () g no)) 1237 ;(setq tlist (nconc tlist (ncons datum))) 1238 (adjoin-datum datum) 1239 (push (cons g kernel) key-vars) 1240 (push (car key-vars) ivars) 1241 ;(push (cons kernel (cons (pget g) 1)) genpairs) 1242 (push (cons g lim) tvar-limits) 1243 (cons g pow))))))) 1244 1245(defun adjoin-datum (datum) 1246 (do ((tlist* tlist (cdr tlist*)) 1247 (tlist** () tlist*)) 1248 ((null tlist*) (setq tlist (nconc tlist (ncons datum)))) 1249 (when (stronger-var? (datum-var (car tlist*)) (datum-var datum)) 1250 (return (if (null tlist**) 1251 (progn 1252 (push datum tlist) 1253 (renumber-tlist tlist)) 1254 (progn 1255 (rplacd tlist** (cons datum tlist*)) 1256 (renumber-tlist (cdr tlist**)))))))) 1257 1258;; Maybe this should just permute the numbering in case it isn't sequential? 1259 1260(defun renumber-tlist (tlist) 1261 (rplacd (data-gvar-o (car tlist)) (cdr (data-gvar-o (cadr tlist)))) 1262 (do ((tlist* (cdr tlist) (cdr tlist*))) 1263 ((null tlist*)) 1264 (rplacd (data-gvar-o (car tlist*)) 1265 (1+ (cdr (data-gvar-o (car tlist*))))))) 1266 1267(defun tvar? (var) (or (atom var) (member 'array (cdar var) :test #'eq))) 1268 1269;; Needs to be extended to handle singular tvars in > 1 var's. 1270 1271(defun stronger-var? (v1 v2) 1272 (let ((e1 (rcone)) (e2 (rcone)) reverse? ans) 1273 (when (alike1 v1 v2) 1274 (tay-err (intl:gettext "taylor: stronger-var? called on equal vars."))) 1275 (when (and (mexptp v1) ($ratnump (caddr v1))) 1276 (setq e1 (prep1 (caddr v1)) v1 (cadr v1))) 1277 (when (and (mexptp v2) ($ratnump (caddr v2))) 1278 (setq e2 (prep1 (caddr v2)) v2 (cadr v2))) 1279 (if (alike1 v1 v2) 1280 (if (equal e1 e2) 1281 (tay-err 1282 (intl:gettext "taylor: stronger-var? called on equal vars.")) 1283 (e> e1 e2)) 1284 (progn 1285 (when (eq (tvar-lim v2) '$finite) 1286 (rotatef v1 v2) 1287 (rotatef e1 e2) 1288 (setq reverse? (not reverse?))) 1289 (if (eq (tvar-lim v1) '$finite) 1290 (if (eq (tvar-lim v2) '$finite) 1291 (great v1 v2) reverse?) 1292 (progn 1293 (when (mtimesp v2) 1294 (rotatef v1 v2) 1295 (rotatef e1 e2) 1296 (setq reverse? (not reverse?))) 1297 (setq ans 1298 (if (mtimesp v1) 1299 (stronger-vars? (order-vars-by-strength (cdr v1)) 1300 (order-vars-by-strength (if (mtimesp v2) (cdr v2) 1301 (ncons (m^ v2 (edisrep e2)))))) 1302 (progn 1303 (when (tvar? v2) 1304 (rotatef v1 v2) 1305 (rotatef e1 e2) 1306 (setq reverse? (not reverse?))) 1307 (if (tvar? v1) 1308 (cond ((tvar? v2) 1309 (let ((n1 (cdr (data-gvar-o (get-datum v1 t)))) 1310 (n2 (cdr (data-gvar-o (get-datum v2 t))))) 1311 (> n1 n2))) 1312 ((mfree v2 (ncons v1)) 1313 (tay-err 1314 (intl:gettext "taylor: Unhandled multivar datum comparison."))) 1315 ((eq (caar v2) '%log) 't) 1316 ((and (eq (caar v2) 'mexpt) (eq (cadr v2) '$%e)) 1317 (stronger-var? `((%log) ,v1) (caddr v2))) 1318 (t (tay-err (intl:gettext "taylor: Unhandled var in stronger-var?.")))) 1319 (progn 1320 (when (eq (caar v2) '%log) 1321 (rotatef v1 v2) 1322 (rotatef e1 e2) 1323 (setq reverse? (not reverse?))) 1324 (if (eq (caar v1) '%log) 1325 (cond ((eq (caar v2) '%log) 1326 (stronger-var? (cadr v1) (cadr v2))) 1327 ((and (eq (caar v2) 'mexpt) (eq (cadr v2) '$%e)) 1328 (stronger-var? `((%log) ,v1) (caddr v2))) 1329 (t (tay-err (intl:gettext "taylor: Unhandled var in stronger-var?")))) 1330 (if (and (eq (caar v1) 'mexpt) (eq (cadr v1) '$%e) 1331 (eq (caar v2) 'mexpt) (eq (cadr v2) '$%e)) 1332 (stronger-var? (caddr v1) (caddr v2)) 1333 (tay-err (intl:gettext "taylor: Unhandled var in stronger-var?"))))))))) 1334 (if reverse? (not ans) ans))))))) 1335 1336(defun neg-monom? (exp) 1337 (and (mtimesp exp) (equal (cadr exp) -1) (null (cdddr exp)) 1338 (caddr exp))) 1339 1340(defun order-vars-by-strength (vars) 1341 (do ((vars* vars (cdr vars*)) (ordvars () )) 1342 ((null vars*) ordvars) 1343 (unless (mfree (car vars*) tvars) ; ignore constants 1344 (do ((ordvars* ordvars (cdr ordvars*))) 1345 ((null ordvars*) 1346 (if (null ordvars) (setq ordvars (ncons (car vars*))) 1347 (rplacd (last ordvars) (ncons (car vars*))))) 1348 (when (stronger-var? (car vars*) (car ordvars*)) 1349 (rplacd ordvars* (cons (car ordvars*) (cdr ordvars*))) 1350 (rplaca ordvars* (car vars*)) 1351 (return () )))))) 1352 1353(defun stronger-vars? (vars1 vars2) 1354 (do ((vars1* vars1 (cdr vars1*)) 1355 (vars2* vars2 (cdr vars2*))) 1356 (()) 1357 (cond ((null vars1*) 1358 (if (null vars2*) 1359 ;; two equal vars generated 1360 (return 't) 1361 (let ((lim (tvar-lim (car vars2*)))) 1362 (return 1363 (cond ((lim-infp lim) ()) 1364 ((lim-zerop lim) 't) 1365 (t (break "var with non-zero finite lim?"))))))) 1366 ((null vars2*) 1367 (let ((lim (tvar-lim (car vars1*)))) 1368 (return 1369 (cond ((lim-infp lim) 't) 1370 ((lim-zerop lim) ()) 1371 (t (break "var with non-zero finite lim?")))))) 1372 ((alike1 (car vars1*) (car vars2*)) ) 1373 ((return (stronger-var? (car vars1*) (car vars2*))))))) 1374 1375(defun stronger-datum? (d1 d2) 1376 (setq d1 (datum-var d1) d2 (datum-var d2)) 1377 (do ((end-flag) (answer)) 1378 (end-flag (member answer '($yes $y) :test #'eq)) 1379 (setq answer (retrieve `((mtext) |Is | ,d1 | stronger than | ,d2 |?|) 1380 nil)) 1381 (if (member answer '($yes $y $no $n) :test #'eq) (setq end-flag 't) 1382 (mtell "~%Acceptable answers are: yes, y, no, n~%")))) 1383 1384(defun datum-lim (datum) 1385 (if (not (tvar? (datum-var datum))) 1386 (exp-pt datum) 1387 (let ((pt (exp-pt datum))) 1388 (if (member pt '($inf $minf) :test #'eq) pt '$zeroa)))) 1389 1390(defun tvar-lim (kernel) 1391 (if (mfree kernel tvars) (coef-sign kernel) 1392 (let ((datum (get-datum kernel t)) lim) 1393 (or (and datum (datum-lim datum)) 1394 (and (setq datum (get-datum (m// kernel) t)) 1395 (setq lim (datum-lim datum)) 1396 (lim-recip lim)) 1397 (progn 1398 (setq lim 1399 (cond ((eq (caar kernel) 'mexpt) 1400 (cond ((and (setq datum (get-datum (cadr kernel) t)) 1401 ($ratnump (caddr kernel))) 1402 (lim-power (datum-lim datum) 1403 (prep1 (caddr kernel)))) 1404 (($ratnump (caddr kernel)) 1405 (lim-power (tvar-lim (cadr kernel)) 1406 (prep1 (caddr kernel)))) 1407 ((eq (cadr kernel) '$%e) 1408 (lim-exp (tvar-lim (caddr kernel)))) 1409 (t (tay-error "Unhandled case in tvar-lim" kernel)))) 1410 ((eq (caar kernel) 'mtimes) 1411 (do ((ans (tvar-lim (cadr kernel)) 1412 (lim-times ans (tvar-lim (car facs)))) 1413 (facs (cddr kernel) (cdr facs))) 1414 ((null facs) ans))) 1415 ((eq (caar kernel) '%log) 1416 (lim-log (datum-lim (get-datum (cadr kernel) t)))) 1417 ((member (caar kernel) '(%sin %cos) :test #'eq) 1418 (unless (lim-infp (tvar-lim (cadr kernel))) 1419 (tay-error "Invalid trig kernel in tvar-lim" kernel)) 1420 '$finite) 1421 (t (tay-error "Unhandled kernel in tvar-lim" kernel)))) 1422 lim))))) 1423 1424(defun coef-sign (coef) 1425 (if (not ($freeof '$%i ($rectform coef))) 1426 '$im 1427 ($asksign coef))) 1428 1429(defun gvar-lim (gvar) 1430 (or (cdr (assoc gvar tvar-limits :test #'eq)) 1431 (if (member (gvar->var gvar) tvars :test #'eq) '$zeroa ; user tvars assumed 0+ now 1432 (break "Invalid gvar")))) 1433 1434(defun psexpt-fn1 (varh trunc l inc m mr ans) 1435 (declare (fixnum m )) 1436 (prog (a (k 0) ak cm-k c sum kr lim) 1437 (declare (fixnum k )) 1438 ;; truly unfortunate that we need so many variables in this hack 1439 (setq a (last ans)) 1440 b (and (e> mr trunc) (go end)) 1441 (setq kr inc ak l k 1 sum (rczero) lim m) 1442 a (cond ((or (> k lim) 1443 (null (setq cm-k (psterm (cdr ans) (e- mr kr))))) 1444 (go add-term))) 1445 (setq ak (or (pstrim-terms ak kr) 1446 (go add-term)) 1447 c (pstimes (ereduce k m) 1448 (pstimes (psterm1 ak kr) cm-k)) 1449 sum (psplus sum c)) 1450 (setq k (1+ k) kr (e+ kr inc)) 1451 (go a) 1452 add-term 1453 (unless (rczerop sum) (add-term-&-pop a mr sum)) 1454 (setq m (1+ m) mr (e+ mr inc)) 1455 (go b) 1456 end 1457 (return (pscheck varh (list trunc) (cdr ans))))) 1458 1459;;; PSEXPT-FN2 and RED-MONO-LOG are needed to reduce exponentials of logs. 1460 1461(defun psexpt-fn2 (p) 1462 (cond ((atom p) (if (get-datum p) 1463 (psexpt-fn (taylor2 p)) 1464 (prep1 `((mexpt) $%e ,p)))) 1465 ((eq (caar p) '%log) 1466 (if (get-datum (cadr p)) (taylor2 (cadr p)) (prep1 (cadr p)))) 1467 ((or (eq (caar p) 'mplus) (eq (caar p) 'mtimes)) 1468 (let ((e ($ratexpand p)) temp) 1469 (cond ((not (and (consp e) (member (caar e) '(mplus mtimes) :test #'eq))) 1470 (psexpt-fn2 e)) 1471 (t 1472 (if (eq (caar e) 'mplus) 1473 (do ((sumnds (cdr e) (cdr sumnds)) (log-facs) (l)) 1474 ((null sumnds) 1475 (cond ((not log-facs) (tsexpt '$%e p)) 1476 (t (tstimes (cons (m^t '$%e (m+l l)) log-facs))))) 1477 (if (setq temp (red-mono-log (car sumnds))) 1478 (push temp log-facs) 1479 (push (car sumnds) l))) 1480 (progn 1481 (setq temp (red-mono-log e)) 1482 (if temp 1483 (taylor2 temp) 1484 (prep1 (power '$%e p))))))))) 1485 (t (prep1 (power '$%e p))))) 1486 1487(defun red-mono-log (e) 1488 (cond ((atom e) ()) 1489 ((eq (caar e) '%log) (cadr e)) 1490 ((mtimesp e) 1491 (do ((facs (cdr e) (cdr facs)) (log-term)) 1492 ((null facs) 1493 (when log-term 1494 (m^t (cadr log-term) (m*l (remove log-term (cdr e) :test #'eq))))) 1495 (if (and (null (atom (car facs))) (eq (caaar facs) '%log)) 1496 (if log-term (return ()) (setq log-term (car facs))) 1497 (unless (mfree (car facs) tvars) (return nil))))) 1498 (t nil ))) 1499 1500(defun pslog (p) 1501 (if (pscoefp p) (pslog2 (rcdisrep p)) 1502 (let ((terms (terms p))) 1503 (cond ((mono-term? terms) ; log(c x^n) = log(c) + n log(x) 1504 ;; do this always for now 1505 (if 't ;$TAYLOR_LOGEXPAND 1506 ;(psplus (pslog (lc terms)) 1507 ; (pstimes (le terms) (pslog-of-gvar (gvar p)))) 1508 (pslog-monom p) 1509 ;(prep1 `((%LOG) ,(term-disrep (lt terms) p))) 1510 )) 1511 ;; expand log(1+ax^n) directly by series substitution 1512 ((not (or (n-term (setq terms (terms (psplus p (rcmone))))) 1513 ;(e> (rczero) (le terms)) 1514 (ps-lim-infp p))) 1515 (setq p (get-series '%log (psexpt-log-ord p) (gvar-o p) 1516 (if (e> (rczero) (le terms)) (e- (le terms)) 1517 (le terms)) 1518 (lc terms))) 1519 (if (e> (rczero) (le terms)) (ps-invert-var p) p)) 1520 (t (prog (l inc trunc lt ans lterm $maxtayorder gvar-lim gt) 1521 ;; log(lt+y) = log(lt) + log(1 + y/lt) = lterm + p 1522 (setq trunc (trunc-lvl p)) 1523 (if (not (member (setq gvar-lim (gvar-lim (gvar p))) 1524 '($zeroa $zerob $inf $minf) :test #'eq)) 1525 (tay-error "bad gvar lim" gvar-lim) 1526 (if (member gvar-lim '($inf $minf) :test #'eq) 1527 (setq lt (ps-gt p) gt lt) 1528 (setq lt (ps-lt p) gt () ))) 1529 (setq lterm (pslog 1530 (setq lt (pscheck (gvar-o p) 1531 (ncons trunc) 1532 (ncons lt)))) 1533 p (pstimes p (let (($maxtayorder 't)) 1534 (psexpt lt (rcmone))))) 1535 (when (and (member gvar-lim '($inf $minf) :test #'eq) 1536 (e> (le terms) (rczero))) 1537 (return (psplus lterm (pslog p)))) 1538 (when (pscoefp p) 1539 (unless (equal p (rcone)) 1540 (merror "PSLOG: internal error.")) 1541 (return lterm)) 1542 (setq l (terms p) inc (psexpon-gcd l)) 1543 (if gt (setq l (delete (last l) l :test #'equal)) 1544 (setq l (n-term l))) 1545 (setq ans (ncons 0)) 1546 (unless $maxtayorder 1547 (setq trunc (emin trunc (t-o-var (gvar p))))) 1548 ;; When we've divided by the greatest term, all terms 1549 ;; have non-positive exponents and we must perform the 1550 ;; transformation x -> 1/x befor calling pslog1 and then 1551 ;; perform the inverse afterwards. 1552 (when gt (setq l (invert-terms l))) 1553 (when (e> (rczero) inc) (setq inc (e- inc))) 1554 (setq ans (psplus lterm 1555 (pslog1 (gvar-o p) trunc l inc 1 inc ans))) 1556 (return 1557 (if (and gt (psp ans) (eq (gvar ans) (gvar p))) 1558 (ps-invert-var ans) 1559 ans)))))))) 1560 1561(defun invert-terms (terms) 1562 (nreverse (mapc #'(lambda (x) (rplaca x (e- (e x)))) terms))) 1563 1564(defun ps-invert-var (ps) 1565 (when (psp ps) (rplacd (cddr ps) (invert-terms (terms ps)))) 1566 ps) 1567 1568(defun ps-gt (ps) 1569 (if (pscoefp ps) (term (rczero) ps) 1570 (lt (last (terms ps))))) 1571 1572(defun pslog1 (varh trunc l inc m mr ans) 1573 (declare (fixnum m )) 1574 (prog (a (k 0) ak cm-k c sum kr m-kr) 1575 (declare (fixnum k )) 1576 ;; truly unfortunate that we need so many variables in this hack 1577 ;; 1578 (setq a (last ans)) 1579 b (and (e> mr trunc) (go end)) 1580 (setq kr inc ak l k 1 sum (rczero)) 1581 a (cond ((or (= k m) 1582 (null (setq cm-k (psterm (cdr ans) 1583 (setq m-kr (e- mr kr)))))) 1584 (go add-term))) 1585 (setq ak (or (pstrim-terms ak kr) 1586 (go add-term)) 1587 c (pstimes m-kr (pstimes (psterm1 ak kr) cm-k)) 1588 sum (psplus sum c) 1589 k (1+ k) kr (e+ kr inc)) 1590 (go a) 1591 add-term 1592 (cond ((setq c (pstrim-terms ak mr)) 1593 (setq c (psterm1 c mr))) 1594 ((setq c (rczero)))) 1595 (setq sum (psdiff c (pstimes sum (e// mr)))) 1596 (unless (rczerop sum) (add-term-&-pop a mr sum)) 1597 (setq m (1+ m) mr (e+ mr inc)) 1598 (go b) 1599 end 1600 (return (pscheck varh (list trunc) (cdr ans))))) 1601 1602;; Computes log(monom), where monom = c x^n. Is extra careful trying to keep 1603;; singular logs going to INF and not generating log(-1)'s unless it is 1604;; necessary to transform a log at MINF to INF. 1605 1606(defun pslog-monom (monom) 1607 (let* ((gvar (gvar monom)) 1608 (datum (gvar-data gvar)) var pt logvar c) 1609 (if (switch 'multivar datum) 1610 (pslog (ps-lc monom)) 1611 (progn 1612 (setq var (datum-var datum)) 1613 (if (tvar? var) 1614 (if (not (member (setq pt (exp-pt datum)) '($inf $minf) :test #'eq)) 1615 (setq logvar (adjoin-tvar `((%log) ,(m- var pt)))) 1616 (progn 1617 ;; At x = inf: log(c (1/x)^n) -> log(c) - n log(x) 1618 ;; At x = minf: log(c (-1/x)^n) -> log(c (-1)^n) - n log(x) 1619 (setq logvar (psminus (adjoin-tvar `((%log) ,var)))) 1620 (when (eq pt '$minf) 1621 (setq c (rcexpt (rcmone) (ps-le monom)))))) 1622 (if (eq (caar var) 'mexpt) 1623 (if (equal (caddr var) -1);; var must be 1/log(y) 1624 ;; Try to keep inf's real. Here we want 1625 ;; log(c (1/log(x))^n) -> log(c (-1)^n) - n log(-log(x)) 1626 (if (equal (tvar-lim (cadr var)) '$minf) 1627 (setq c (rcexpt (rcmone) (ps-le monom)) 1628 logvar 1629 (psminus (adjoin-tvar 1630 `((%log) ,(m- (cadr var)))))) 1631 (setq logvar (psminus 1632 (adjoin-tvar `((%log) ,(cadr var)))))) 1633 (if (equal (cadr var) '$%e) 1634 (setq logvar (taylor2 (caddr var))) 1635 (break "Unhandled gvar in `pslog-of-gvar'"))))) 1636 (psplus (pslog (if c (pstimes c (ps-lc monom)) (ps-lc monom))) 1637 (pstimes (ps-le monom) logvar)))))) 1638 1639;; Computes log(p), where p is an rcdisrep'd pscoef. 1640 1641(defun pslog2 (p) (let ($logarc) (pslog3 p))) 1642 1643(defun pslog3 (p) 1644 (cond ((atom p) 1645 (prep1 (cond ((equal p 1) 0) 1646 ((equal p -1) log-1) 1647 ((eq p '$%i) log%i) 1648 ((eq p '$%e) 1) 1649 ((equal p 0) 1650 (merror (intl:gettext "taylor: log(0) encountered while processing ~:M") last-exp)) 1651 (t `((%log) ,p))))) 1652 ((eq (caar p) 'rat) 1653 (prep1 (cond ((not $taylor_logexpand) `((%log) ,p)) 1654 (t (m- `((%log) ,(cadr p)) `((%log) ,(caddr p))))))) 1655 ((and full-log (not (free p '$%i))) 1656 (let ((full-log () )) (pslog3 ($polarform p)))) 1657 ((eq (caar p) 'mexpt) 1658 ;; Must handle things like x^a, %e^(a*x), etc. which are pscoef's. 1659 (pstimes (taylor2 (caddr p)) (pslog (taylor2 (cadr p))))) 1660 ((and (eq (caar p) 'mtimes) $taylor_logexpand) 1661 (do ((l (cddr p) (cdr l)) 1662 (ans (pslog3 (cadr p)) (psplus ans (pslog3 (car l))))) 1663 ((null l) ans))) 1664 (t (prep1 `((%log) ,p))))) 1665 1666;;; Subtitle Extending Routines 1667 1668(defun getfun-lt (fun) 1669 (let ((exp-datum (get (oper-name fun) 'exp-form))) 1670 (cond (exp-datum 1671 ;; Info not needed yet. 1672 ;; (or (atom (car exp-datum)) 1673 ;; (setq 0p-funord (copy-tree (cdar exp-datum)))) 1674 (exp-datum-lt fun exp-datum)) 1675 ((setq exp-datum (get (oper-name fun) 'sp2)) 1676 (setq exp-datum (get-lexp (subst (dummy-var) 'sp2var exp-datum) 1677 (rcone) ())) 1678 ;; Info not needed yet; need to bind lexp-non0 to T when 1679 ;; this is used though so n-term will be there. 1680 ;; (and (rczerop (le exp-datum)) 1681 ;; (setq 0p-funord (le (n-term exp-datum)))) 1682 (if (psp exp-datum) (ps-lt exp-datum) 1683 (term (rczero) exp-datum))) 1684 (t (merror "GETFUN-LT: unknown function ~A" fun))))) 1685 1686(declare-top (special var)) 1687 1688(defun getexp-fun (fun var pw) 1689 (declare (special var)) 1690 (let ((exp-datum (copy-tree (get (oper-name fun) 'exp-form)))) 1691 (cond ((infp pw) (infin-ord-err)) 1692 ((null exp-datum) 1693 (if (null (setq exp-datum 1694 (get-ps-form (if (atom fun) fun (caar fun))))) 1695 (merror (intl:gettext "taylor: power series unavailable for function ~A") fun) 1696 (progn 1697 (unless (atom fun) 1698 (do ((subvals (cdr fun) (cdr subvals)) 1699 (subs (safe-get (caar fun) 'sp2subs) (cdr subs))) 1700 ((or (null subvals) (null subs)) 1701 (when (or subvals subs) 1702 (merror (intl:gettext "taylor: incorrect number of subscripts to the deftaylor'd function ~A") (caar fun)))) 1703 (setq exp-datum (maxima-substitute (car subvals) (car subs) 1704 exp-datum)))) 1705 (ts-formula exp-datum var pw)))) 1706 ((e> (exp-datum-le fun exp-datum) pw) (pszero var pw)) 1707 ((setq exp-datum 1708 (apply (exp-fun exp-datum) 1709 (if (atom fun) (cons pw (cdr exp-datum)) 1710 (cons pw (cons (cdr fun) (cdr exp-datum)))))) 1711 (cond ((null exp-datum) (pszero var pw)) 1712 ((psp exp-datum) exp-datum) 1713 (t (make-ps var (ncons pw) exp-datum))))))) 1714 1715(declare-top (unspecial var)) 1716 1717(defun expexp-funs (pw l sign chng inc) 1718 (prog (e lt-l) 1719 (setq e (e l) lt-l (setq l (ncons l))) 1720 a (cond ((e> (setq e (e+ e inc)) pw) (return l)) 1721 (t (add-term-&-pop 1722 lt-l 1723 e 1724 (rctimes (e// sign 1725 (cond ((e= inc (rcone)) e) 1726 ((e* e (e1- e))))) 1727 (cons 1 (cdr (lc lt-l))))) 1728 (setq sign (e* sign chng)))) 1729 (go a))) 1730 1731;; returns series expansion of %expintegral_si 1732;; from term x^i to term x^max 1733;; sign is sign of term x^i (1 or -1) 1734;; ifac is i! (to avoid repeated computation) 1735(defun expsi_series (i sign ifac max) 1736 (if (> i max) 1737 nil ; we have all powers up to max 1738 (cons (cons (cons i 1) ; this is i'th term 1739 (cons sign (* i ifac))); i'th coefficient 1740 (expsi_series (+ 2 i) ; advance to next even/odd power 1741 (* -1 sign) ; flip sign 1742 (* (+ 2 i) (+ 1 i) ifac) ; update factorial 1743 max)))) ; pass through the max 1744 1745;; returns series expansion of %expintegral_si up to term with x^pw 1746(defun exp_%expintegral_si (pw l) 1747 (expsi_series 1 1748 1 1749 1 1750 (/ (float (car pw)) (float (cdr pw))))) 1751 1752(defun explog-funs (pw l sign chng inc) 1753 (prog (e lt-l) 1754 (setq e (e l) lt-l (setq l (ncons l))) 1755 a (cond ((e> (setq e (e+ e inc)) pw) (return l)) 1756 (t (add-term lt-l e (e// sign e)) 1757 (setq lt-l (n-term lt-l) 1758 sign (e* sign chng)))) 1759 (go a))) 1760 1761(defun exptan-funs (pw l chng) 1762 (prog (e lt-l sign fact pow) 1763 (setq e (e l) lt-l (setq l (ncons l)) 1764 sign (rcone) fact '(1 . 2) pow '(4 . 1)) 1765 a (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l)) 1766 (t (setq fact (e// fact (e* e (e1+ e))) 1767 pow (e* '(4 . 1) pow) 1768 sign (e* chng sign)) 1769 (add-term lt-l e (e* (e* sign fact) 1770 (e* (prep1 1771 ($bern (rcdisrep (e1+ e)))) 1772 (e* pow (e1- pow))))) 1773 (setq lt-l (n-term lt-l)))) 1774 (go a))) 1775 1776(defun expcot-funs (pw l sign chng plus) 1777 (prog (e lt-l fact pow) 1778 (setq e (e l) lt-l (setq l (ncons l)) 1779 fact (rcone) pow (rcone)) 1780 a (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l)) 1781 (t (setq fact (e// fact (e* e (e1+ e))) 1782 pow (e* '(4 . 1) pow) 1783 sign (e* chng sign)) 1784 (add-term lt-l e (e* (e* sign fact) 1785 (e* (prep1 1786 ($bern (rcdisrep (e1+ e)))) 1787 (e+ pow plus)))) 1788 (setq lt-l (n-term lt-l)))) 1789 (go a))) 1790 1791(defun expsec-funs (pw l chng) 1792 (prog (e lt-l sign fact) 1793 (setq e (e l) lt-l (setq l (ncons l)) 1794 sign (rcone) fact (rcone)) 1795 a (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l)) 1796 (t (setq fact (e// fact (e* e (e1- e))) 1797 sign (e* chng sign)) 1798 (add-term lt-l e (e* (e* sign fact) 1799 (prep1 ($euler (rcdisrep e))))) 1800 (setq lt-l (n-term lt-l)))) 1801 (go a))) 1802 1803(defun expasin-funs (pw l chng) 1804 (prog (e lt-l sign n d) 1805 (setq e (e l) lt-l (setq l (ncons l)) sign 1 n 1 d 1) 1806 a (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l)) 1807 (t (setq n (* n (car (e- e (rctwo)))) 1808 d (* d (car (e1- e))) 1809 sign (* sign chng)) 1810 (add-term lt-l e ; need to reduce here ? - check this. 1811 (let ((x (*red (* n sign) 1812 (* d (car e))))) 1813 (if (atom x) x 1814 (cons (cadr x) (caddr x))))) 1815 (setq lt-l (n-term lt-l)))) 1816 (go a))) 1817 1818;;; This is the table of expansion data for known functions. 1819;;; The format of the EXP-FORM property is as follows: 1820;;; (<name of the expanding routine for the function or 1821;;; (name . le of n-term) if expansion is of order 0> 1822;;; <first term in the expansion or the name of a routine which 1823;;; computes the order when it may depend on parameters (e.g subsripts)> 1824;;; <data for the expanding routine>) 1825 1826 1827(loop for (fun exp) on 1828 '(%exp ((expexp-funs 1 . 1) ((0 . 1) 1 . 1) (1 . 1) (1 . 1) (1 . 1)) 1829 %sin (expexp-funs ((1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (2 . 1)) 1830 %cos ((expexp-funs 2 . 1) ((0 . 1) 1 . 1) (-1 . 1) (-1 . 1) (2 . 1)) 1831 %sinh (expexp-funs ((1 . 1) 1 . 1) (1 . 1) (1 . 1) (2 . 1)) 1832 %cosh ((expexp-funs 2 . 1) ((0 . 1) 1 . 1) (1 . 1) (1 . 1) (2 . 1)) 1833 %log (explog-funs ((1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (1 . 1)) 1834 %atan (explog-funs ((1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (2 . 1)) 1835 %atanh (explog-funs ((1 . 1) 1 . 1) (1 . 1) (1 . 1) (2 . 1)) 1836 %cot (expcot-funs ((-1 . 1) 1 . 1) (1 . 1) (-1 . 1) (0 . 1)) 1837 %csc (expcot-funs ((-1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (-2 . 1)) 1838 %csch (expcot-funs ((-1 . 1) 1 . 1) (-1 . 1) (1 . 1) (-2 . 1)) 1839 %coth (expcot-funs ((-1 . 1) 1 . 1) (1 . 1) (1 . 1) (0 . 1)) 1840 %tan (exptan-funs ((1 . 1) 1 . 1) (-1 . 1)) 1841 %tanh (exptan-funs ((1 . 1) 1 . 1) (1 . 1)) 1842 %sec ((expsec-funs 2 . 1) ((0 . 1) 1 . 1) (-1 . 1)) 1843 %sech ((expsec-funs 2 . 1) ((0 . 1) 1 . 1) (1 . 1)) 1844 %asin (expasin-funs ((1 . 1) 1 . 1) 1) 1845 %asinh (expasin-funs ((1 . 1) 1 . 1) -1) 1846 %gamma (expgam-fun ((-1 . 1) 1 . 1)) 1847 $li (exp$li-fun li-ord) 1848 %expintegral_si (exp_%expintegral_si ((1 . 1) 1 . 1)) 1849 $psi (expplygam-funs plygam-ord)) 1850 by #'cddr 1851 do (putprop fun exp 'exp-form)) 1852 1853 1854(defun known-ps (fun) 1855 (getl fun '(exp-form sp2))) 1856 1857;;; Autoload Properties 1858 1859;;; Taylor series expansion routines 1860 1861;;; SRF is only called externally; by RATF and SIMPEXPT. 1862 1863(defun srf (x) 1864 (let ((exact-poly t) (tlist) (*within-srf?* 't)) 1865 (setq x (taylor1 x ()) tlist (mrat-tlist x)) 1866 ;; Set trunc levels in the local tlist to correspond to the maximum 1867 ;; level occurring in any series. 1868 (do ((data tlist (cdr data)) 1869 (truncs (trunc-vector (mrat-ps x) () ))) 1870 ((null data)) 1871 (when (and (car truncs) (e> (car truncs) (current-trunc (car data)))) 1872 (setf (current-trunc (car data)) (car truncs)))) 1873 x)) 1874 1875;;; [var, pt, order, asymp] 1876 1877(defmfun $taylor (e &rest args) 1878 (when (not ($ratp e)) 1879 ;; Not a mrat expression. Remove the special representation. 1880 (setq e (specrepcheck e))) 1881 (taylor* e args)) 1882 1883(defun taylor* (arg l) 1884 ;; We must bind $MAXTAYORDER to () below because of the problem of constants 1885 ;; not retaining their truncation level. This means that when we add a 1886 ;; series which has more terms than the user-specified truncation to a 1887 ;; constant we must truncate the series with more terms down to the user 1888 ;; specified level because, in the worst case, the constant could be a 1889 ;; series no better than to the user-specified level. Hence $MAXTAYORDER 1890 ;; is essentially useless until the constant problem is fixed. If one 1891 ;; decides to not bind $MAXTAYORDER below then the sum routines must 1892 ;; be updated to truncate series with more terms than the user-specified 1893 ;; level down to that level---taylor(sin(x)^2-cos(x)^2-1,x,0,1) would 1894 ;; give x^2+... in this case if the sum routines weren't updated. 1895 ;; Also, batch(mquery,160,aljabr) for another truncation bug which crops 1896 ;; up when $maxtayorder isn't bound here. Similarly, loadfile(taybad,rl, 1897 ;; aljabr) and see tomh's bug note of 4/15/81. 1898 (let ((tlist () ) ($maxtayorder () ) (*within-srf?* () ) 1899 (exact-poly (if l (not $taylor_truncate_polynomials) 'user-specified))) 1900 (declare (special *within-srf?*)) 1901 1902 (parse-tay-args l) 1903 (taylor1 arg (ncons tlist)))) 1904 1905(defun tay-order (n) 1906 (let (($float) (modulus)) 1907 (cond ((eq n '$inf) (ncons (inf))) 1908 ((null n) (wna-err '$taylor)) 1909 ((null (mnump n)) 1910 (merror (intl:gettext "taylor: expansion order must be a number; found: ~:M") n)) 1911 (t (ncons (prep1 n)))))) 1912 1913(defun re-erat (head exp) 1914 (taylor1 exp (list (cadddr (cdr head))))) 1915 1916(defun parse-tay-args (l) 1917 (cond ((null l) ) 1918 ((numberp (car l)) 1919 (merror (intl:gettext "taylor: variable of expansion cannot be a number: ~M") (car l))) 1920 ((or (symbolp (car l)) (not (eq (caaar l) 'mlist))) 1921 (parse-tay-args1 (list (car l) ($ratdisrep (cadr l)) (caddr l))) 1922 (parse-tay-args (cdddr l))) 1923 ((do ((l (cddar l) (cdr l))) 1924 ((null l) () ) 1925 (and (or (mnump (car l)) (eq (car l) '$inf)) 1926 (return 't))) 1927 (parse-tay-args1 (cdar l)) 1928 (parse-tay-args (cdr l))) 1929 (t (parse-tay-args2 (list (car l) (cadr l) (caddr l))) 1930 (parse-tay-args (cdddr l))))) 1931 1932(defun parse-tay-args1 (l) 1933 (if ($listp (car l)) (parse-tay-args2 l) 1934 (let ((v (car l)) 1935 (pt ($ratdisrep (cadr l))) 1936 (ord (tay-order (caddr l))) 1937 (switches (make-switch-list (cdddr l)))) 1938 (push (list v ord pt switches) tlist)))) 1939 1940(defun parse-tay-args2 (l) 1941 (let ((label (gensym)) 1942 (vs (cdar l)) 1943 (pts (make-long-list (if ($listp (cadr l)) 1944 (copy-list (cdadr l)) 1945 (ncons (ratdisrep (cadr l)))))) 1946 (ord (caddr l)) 1947 (switches (make-switch-list (cdddr l))) 1948 (lcm 1) 1949 (max 0)) 1950 (if (atom ord) 1951 (setq lcm ord max ord ord (make-long-list (ncons ord))) 1952 (do ((a vs (cdr a)) 1953 (l (cdr ord) (cdr l))) 1954 ((null a) (setq ord (cdr ord))) 1955 (cond ((not l) (merror "PARSE-TAY-ARGS2: ran out of truncation levels.")) 1956 (t (setq lcm (lcm lcm (car l)) max (max max (car l))))))) 1957 (push (list label (tay-order max) 0 1958 (ncons (list 'multivar lcm vs))) 1959 tlist) 1960 (do ((vl vs (cdr vl)) 1961 (ordl ord (cdr ordl)) 1962 (ptl pts (cdr ptl))) 1963 ((null vl) ) 1964 (cond ((not ptl) (merror "PARSE-TAY-ARGS2: ran out of matching points of expansion.")) 1965 (t 1966 (push 1967 (list (car vl) (tay-order (car ordl)) (car ptl) 1968 (cons (list 'multi label (timesk lcm (expta (car ordl) -1))) switches)) 1969 tlist)))))) 1970 1971(defun make-switch-list (l) 1972 (mapcar #'(lambda (q) (cons q t)) l)) 1973 1974(defun make-long-list (q) 1975 (nconc q q)) 1976 1977;;; This checks to ensure that there isn't more than one set of multi- 1978;;; dependent variables with different orders of expansion, e.g. 1979;;; taylor(u+v+x+y,[u,v],0,[2,3],[x,y],0,[5,7]) is one. 1980 1981(defun ratwtsetup (l) 1982 (do ((l l (cdr l)) (a) (sw)) 1983 ((null l) ) 1984 (when (setq a (switch 'multivar (car l))) 1985 (do ((ll (cadr a) (cdr ll))) 1986 ((null ll) ) 1987 (cond ((equal (cadr (switch 'multi (get-datum (car ll)))) 1) ) 1988 (sw (merror (intl:gettext "taylor: multiple dependent variables must all have the same order of expansion."))) 1989 ('t (setq sw 't) (return 't))))))) 1990 1991(defmvar $taylor_order_coefficients t 1992 "When `true', coefficients of taylor series will be ordered canonically.") 1993 1994(defun taylor1 (e tlist) 1995 (declare (special *within-srf?* )) 1996 (setq tlist (tlist-merge (nconc (find-tlists e) tlist))) 1997 (prog ($zerobern $simp $algebraic genpairs varlist tvars sing-tvars 1998 log-1 log%i ivars key-vars ans full-log last-exp 1999 mainvar-datum zerolist taylor_simplifier least_term? tvar-limits 2000 genvar) 2001 (setq tlist (mapcan #'(lambda (d) 2002 (if (tvar? (datum-var d)) 2003 (ncons d) 2004 (progn 2005 (push d sing-tvars) 2006 () ))) 2007 tlist)) 2008 (setq $zerobern t $simp t $algebraic t last-exp e least_term? 't 2009 log-1 '((%log simp) -1) log%i '((%log simp) $%i) 2010 tvars (mapcar 'car tlist) varlist (copy-list tvars)) 2011 (when $taylor_simplifier 2012 ; This symbolp/fboundp check is presumably for efficiency (so it 2013 ; can be directly funcalled). 2014 (setq taylor_simplifier 2015 (if (and (symbolp $taylor_simplifier) 2016 (fboundp $taylor_simplifier)) 2017 $taylor_simplifier 2018 'taylor_simplifier_caller))) 2019 ;; Ensure that the expansion points don't depend on the expansion vars. 2020 ;; This could cause an infinite loop, e.g. taylor(x,x,x,1). 2021 (do ((tl tlist (cdr tl))) 2022 ((null tl) ) 2023 (unless (mfree (exp-pt (car tl)) tvars) 2024 (merror (intl:gettext "taylor: attempt to expand ~M at a point depending on ~M.") e (caar tl)))) 2025 ;; This drastic initialization ensures that ALGEBRAIC, TELLRAT, DISREP, 2026 ;; etc. prop's are removed from our gensyms. RATSETUP does not appear 2027 ;; to do this correctly, e.g. see ASB's bug of 1/10/83 (MQUERY 17). 2028 (mapc #'(lambda (g) (setf (symbol-plist g) nil)) genvar) 2029 (ratsetup varlist genvar) 2030 (when (and $taylor_order_coefficients (not *within-srf?*)) (newvar e)) 2031 (orderpointer varlist) 2032 (maplist #'(lambda (q g) 2033 (push (cons (car g) (car q)) key-vars) 2034 (let ((data (get-datum (car q)))) 2035 (rplaca q (transform-tvar (car q) data)) 2036 (push (cons (car g) (car q)) ivars) 2037 ;(setf (data-gvar-o data) 2038 ; (cons (car g) (valget (car g)))) 2039 (rplacd (cdddr data) 2040 (cons (car g) (valget (car g)))))) 2041 (do ((v varlist (cdr v))) 2042 ((eq (car v) (car tvars)) v)) 2043 (do ((v varlist (cdr v)) 2044 (g genvar (cdr g))) 2045 ((eq (car v) (car tvars)) g))) 2046 (setq genpairs (mapcar #'(lambda (y z) 2047 (putprop z y 'disrep) 2048 (cons y (cons (pget z) 1))) 2049 varlist genvar)) 2050 (ratwtsetup tlist) 2051 (setup-multivar-disrep () ) 2052 (setq mainvar-datum (car (last tlist))) 2053 (mapc #'(lambda (d) (adjoin-sing-datum d)) sing-tvars) 2054 (setq ans (catch 'tay-err (taylor3 e))) 2055 (return 2056 (if (atom (car ans)) (tay-error (car ans) (cadr ans)) ans)))) 2057 2058(defun transform-tvar (var data) 2059 (if (not (tvar? var)) var 2060 (cond ((and (signp e (exp-pt data)) (null (switch '$asymp data))) 2061 var) ;Simple case 2062 ((member (exp-pt data) '($inf infinity) :test #'eq) 2063 (m^ var -1)) 2064 ((eq (exp-pt data) '$minf) 2065 (m- (m^ var -1))) 2066 ((let ((temp (m- var (exp-pt data)))) 2067 (if (switch '$asymp data) (m^ temp -1) temp)))))) 2068 2069(defun taylor_simplifier_caller (exp) 2070 (mcall $taylor_simplifier exp)) 2071 2072(defun taylor_simplify_recurse (ps) 2073 (if (pscoefp ps) (taylor2 (funcall taylor_simplifier (rcdisrep ps))) 2074 (let ((datum (ps-data ps)) (var () )) 2075 ;; We must treat multivars like 1, since they'll reappear again 2076 ;; when we call taylor2 on their disrep'd coeff's. 2077 (if (switch 'multivar datum) 2078 (setq datum '()) 2079 (progn 2080 (setq var (getdisrep (gvar-o ps))) 2081 ;; Don't push pw's < 0, else constants will be truncated 2082 (push-pw datum (emax (trunc-lvl ps) (rczero))))) 2083 (do ((terms (terms ps) (n-term terms)) 2084 (ans (rczero) (psplus (if (null datum) 2085 (taylor_simplify_recurse (lc terms)) 2086 (pstimes (taylor_simplify_recurse 2087 (lc terms)) 2088 ;; Don't do 2089 ;; (taylor2 (funcall taylor_simplifier 2090 ;; (m^ var (edisrep (le terms))))) 2091 ;; causes terms to be lost when inverting. E.g. 2092 ;; taylor(log(1+exp(-1/x)),x,0,5) calls psexpt(<exp(1/x)^3>...3,-1) 2093 ;; which must return a series good to 3+3(-1-1)=-3 which, when added 2094 ;; to other terms will truncate them to degree -3 also. 2095 (if (ezerop (le terms)) (rcone) 2096 (make-ps ps 2097 (ncons 2098 (term (le terms) (rcone))))))) 2099 2100 ans))) 2101 ((null terms) 2102 (when datum (pop-pw datum)) 2103 ans))))) 2104 2105(defun push-pw (datum pw) 2106 (push pw (trunc-stack datum)) 2107 ;; When changing the truncation on an internal multivar we must also 2108 ;; propagate the change to all var's which depend upon it. See WGD's 2109 ;; bug report of 9/15/82 which exhibits the necessity of this. 2110 (when (setq datum (switch 'multivar datum)) 2111 (do ((vars (cadr datum) (cdr vars))) 2112 ((null vars) ) 2113 (push pw (trunc-stack (get-datum (car vars))))))) 2114 2115(defun pop-pw (datum) 2116 (pop (trunc-stack datum)) 2117 ;; See the comment above in push-pw; here we must undo the propagation. 2118 (when (setq datum (switch 'multivar datum)) 2119 (do ((vars (cadr datum) (cdr vars))) 2120 ((null vars) ) 2121 (pop (trunc-stack (get-datum (car vars))))))) 2122 2123(defun setup-multivar-disrep (mrat?) 2124 (let ((varlist varlist) (genvar genvar) (multivars () )) 2125 (when mrat? 2126 (setq varlist (mrat-varlist mrat?) genvar (mrat-genvar mrat?))) 2127 (mapc #'(lambda (datum) 2128 (when (switch 'multivar datum) 2129 (push (car datum) multivars) 2130 ;; All genvar's corresponding to internally generated 2131 ;; multivars must "disappear" when disrep'd. If this 2132 ;; were not done then disrep'ing gx*gt would give x*t 2133 ;; which, upon, re-tayloring would give (gx*gt)*gt, 2134 ;; where t is the internal multivar for x, and gt, gx 2135 ;; are their genvars. An example where this would occur is 2136 ;; taylor(sin(x+y),[x],0,f1,[y],0,1). 2137 (putprop (int-gvar datum) 1 'disrep))) 2138 (if mrat? (mrat-tlist mrat?) tlist)) 2139 ;; Here we must substitute 1 for any genvars which depend on multivars. 2140 ;; For example, taylor(x^n,[x],0,0) generates a multivar^n. 2141 (when multivars 2142 (do ((expl varlist (cdr expl)) 2143 (gvarl genvar (cdr gvarl))) 2144 ((null expl) ) 2145 (unless (mfree (car expl) multivars) 2146 (putprop (car gvarl) 1 'disrep)))))) 2147 2148;; An example showing why this flag is need is given by 2149;; taylor(-exp(exp(-1/x)+2/x),x,0,-1). Without it, tstimes and 2150;; taylor_simplify_recurse end up trunc'ing the -1. 2151 2152(defvar trunc-constants? 't) 2153 2154(defun taylor3 (e) 2155 (cond ((mbagp e) (cons (car e) (mapcar #'taylor3 (cdr e)))) 2156 ((and (null tlist) (not (eq exact-poly 'user-specified))) 2157 (xcons (prep1 e) 2158 (list 'mrat 'simp varlist genvar))) 2159 (t (xcons (if (null taylor_simplifier) 2160 (taylor2 e) 2161 (progn 2162 (setq e (taylor2 e)) 2163 (let ((exact-poly () ) (trunc-constants? () )) 2164 (taylor_simplify_recurse e)))) 2165 (list 'mrat 'simp varlist genvar tlist 'trunc))))) 2166 2167(defun find-tlists (e) (let (*a*) (findtl1 e) *a*)) 2168 2169(defun findtl1 (e) 2170 (cond ((or (atom e) (mnump e)) ) 2171 ((eq (caar e) 'mrat) 2172 (when (member 'trunc (car e) :test #'eq) 2173 (push (mapcar #'copy-tree (mrat-tlist e)) *a*))) 2174 (t (mapc #'findtl1 (cdr e))))) 2175 2176(defun tlist-merge (tlists) 2177 (do ((tlists tlists (cdr tlists)) (tlist () )) 2178 ((null tlists) (nreverse tlist)) 2179 (do ((a_tlist (car tlists) (cdr a_tlist)) (temp nil)) 2180 ((null a_tlist) ) 2181 (if (null (setq temp (get-datum (datum-var (car a_tlist)) t))) 2182 (if (prog2 (setq temp (car a_tlist)) 2183 (or (tvar? (datum-var temp)) 2184 (member (caar (datum-var temp)) '(mexpt %log) :test #'eq))) 2185 (push (list (datum-var temp) (trunc-stack temp) 2186 (exp-pt temp) (switches temp)) 2187 tlist) 2188 (merror (intl:gettext "taylor: ~M cannot be a variable.") (datum-var temp))) 2189 (progn 2190 (if $maxtayorder 2191 ;; We must take the max truncation level when $maxtayorder 2192 ;; is T, cf. JPG's bug of 9/15/82. 2193 (when (e> (current-trunc (car a_tlist)) (current-trunc temp)) 2194 (setf (current-trunc temp) (current-trunc (car a_tlist)))) 2195 (unless (e> (current-trunc (car a_tlist)) (current-trunc temp)) 2196 (setf (current-trunc temp) (current-trunc (car a_tlist))))) 2197 (unless (alike1 (exp-pt temp) (exp-pt (car a_tlist))) 2198 (merror (intl:gettext "taylor: cannot combine expressions expanded at different points."))) 2199 (setf (switches temp) 2200 (union* (switches temp) (switches (car a_tlist))))))))) 2201 2202(defun compattlist (list) 2203 (do ((l list (cdr l))) 2204 ((null l) t) 2205 (or (alike1 (exp-pt (get-datum (datum-var (car l)))) (exp-pt (car l))) 2206 (return () )))) 2207 2208(defun taylor2 (e) 2209 (let ((last-exp e)) ;; lexp-non0 should be bound here when needed 2210 (cond ((assolike e tlist) (var-expand e 1 () )) 2211 ((or (mnump e) (atom e) (mfree e tvars)) 2212 (if (or (e> (rczero) (current-trunc mainvar-datum)) 2213 (lim-zerop e)) 2214 (pszero (data-gvar-o mainvar-datum) 2215 (current-trunc mainvar-datum)) 2216 (if (and taylor_simplifier (not (atom e))) 2217 (let ((e-simp (prep1 (funcall taylor_simplifier e)))) 2218 (when (and (rczerop e-simp) (not (member e-simp zerolist :test #'eq))) 2219 (push e zerolist)) 2220 e-simp) 2221 (prep1 e)))) 2222 ((null (atom (caar e))) (merror "TAYLOR2: internal error.")) 2223 (($taylorp e) 2224 (if (and (compatvarlist varlist (mrat-varlist e) 2225 genvar (mrat-genvar e)) 2226 (compattlist (mrat-tlist e))) 2227 (pstrunc (cdr e)) 2228 (let ((exact-poly () )) (re-taylor e)))) 2229 ((eq (caar e) 'mplus) (tsplus (cdr e))) 2230 ((eq (caar e) 'mtimes) (tstimes (cdr e))) 2231 ((eq (caar e) 'mexpt) (tsexpt (cadr e) (caddr e))) 2232 ((eq (caar e) '%log) (tslog (cadr e))) 2233 ((and (or (known-ps (caar e)) (get (caar e) 'tay-trans)) 2234 (not (member 'array (cdar e) :test #'eq)) 2235 (try-expansion (if (cddr e) (cdr e) (cadr e)) 2236 (caar e))) ) 2237 ((and (mqapplyp e) 2238 (cond ((get (subfunname e) 'spec-trans) 2239 (funcall (get (subfunname e) 'spec-trans) e)) 2240 ((known-ps (subfunname e)) 2241 (try-expansion (caddr e) (cadr e))))) ) 2242 ((and (member (caar e) '(%sum %product) :test #'eq) 2243 (mfreel (cddr e) tvars)) 2244 (tsprsum (cadr e) (cddr e) (caar e))) 2245 ((eq (caar e) '%derivative) (tsdiff (cadr e) (cddr e) e)) 2246 ((or (eq (caar e) '%at) 2247 (do ((l (mapcar 'car tlist) (cdr l))) 2248 ((null l) t) 2249 (or (free e (car l)) (return ())))) 2250 (newsym e)) 2251 (t (let ((exact-poly () )) ; Taylor series aren't exact 2252 (taylor2 (diff-expand e tlist))))))) 2253 2254(defun compatvarlist (a b c d) 2255 (cond ((null a) t) 2256 ((or (null b) (null c) (null d)) () ) 2257 ((alike1 (car a) (car b)) 2258 (if (not (eq (car c) (car d))) () 2259 (compatvarlist (cdr a) (cdr b) (cdr c) (cdr d)))) 2260 (t (compatvarlist a (cdr b) c (cdr d))))) 2261 2262 2263(defun re-taylor (mrat) 2264 (let ((old-tlist (mrat-tlist mrat)) (old-varlist (mrat-varlist mrat)) 2265 (old-genvar (mrat-genvar mrat)) old-ivars) 2266 (declare (special old-tlist old-ivars)) 2267 ;; Put back the old disrpes so rcdisrep's will work correctly. 2268 (mapc #'(lambda (g v) (putprop g v 'disrep)) old-genvar old-varlist) 2269 (setup-multivar-disrep mrat) 2270 (setq old-ivars (mapcar #'(lambda (g v) (cons g v)) 2271 old-genvar old-varlist)) 2272 (prog1 (re-taylor-recurse (mrat-ps mrat)) 2273 ;; Restore the correct disreps. 2274 (mapc #'(lambda (g v) (putprop g v 'disrep)) genvar varlist) 2275 (setup-multivar-disrep () )))) 2276 2277(defun re-taylor-recurse (ps) 2278 (declare (special old-tlist old-ivars)) 2279 (if (not (psp ps)) (taylor2 (rcdisrep ps)) 2280 (let (var (datum () )) 2281 (setq var (cdr (assoc (gvar ps) old-ivars :test #'eq))) 2282 ;; We must treat multivars like 1, since they'll reappear again 2283 ;; when we call taylor2 or var-expand below. 2284 (if (switch 'multivar (assoc var old-tlist :test #'equal)) 2285 (setq var () ) 2286 (when (setq datum (var-data var)) 2287 (push-pw datum (trunc-lvl ps)))) 2288 (prog1 2289 (do ((terms (terms ps) (n-term terms)) 2290 (ans (rczero) 2291 (psplus (if (null var) (re-taylor-recurse (lc terms)) 2292 (pstimes (re-taylor-recurse (lc terms)) 2293 (if datum 2294 (var-expand (car datum) 2295 (edisrep (le terms)) 2296 () ) 2297 (taylor2 2298 (m^t var (edisrep (le terms))))))) 2299 ans))) 2300 ((null terms) ans)) 2301 (when datum (pop-pw datum)))))) 2302 2303(defun var-expand (var exp dont-truncate?) 2304 (let (($keepfloat) ($float) (modulus)) 2305 (setq exp (prep1 exp))) ;; exp must be a rational integer 2306 (let ((temp (get-datum var 't))) 2307 (cond ((null temp) (merror "VAR-EXPAND: invalid call.")) 2308 ((member (exp-pt temp) '($inf $minf $infinity) :test #'eq) 2309 (cond ((switch '$asymp temp) 2310 (merror (intl:gettext "taylor: cannot create an asymptotic expansion at infinity."))) 2311 ((e> (setq exp (rcminus exp)) (current-trunc temp)) 2312 (rczero)) 2313 (t (make-ps (int-var temp) 2314 (ncons (if exact-poly (inf) (current-trunc temp))) 2315 (ncons (term exp 2316 (if (eq (exp-pt temp) '$minf) 2317 (rcmone) 2318 (rcone)))))))) 2319 ;; multivar expansion does not work at infinity, so 2320 ;; expansion at infinity is handled by above clause even if doing multivar. 2321 ((switch 'multi temp) ;; multivar expansion 2322 (psexpt (psplus 2323 ;; The reason we call var-expand below instead of taylor2 2324 ;; is that we must be sure the call is not truncated to 2325 ;; 0 which would cause an error in psexpt if exp < 0. 2326 ;; For example, this occurred in TAYLOR(X^2/Y,[X,Y],0,2). 2327 (pstimes 2328 ;; Must ensure that we get back a series truncated 2329 ;; to at least what is specified by tlist. This means 2330 ;; we'll have to push-pw unless exp>0 since psexpt'n 2331 ;; kills (exp-1) terms. The bug that discovered this 2332 ;; is taylor(li[2](x+1/2)/x,[x],0,0) missing 2*log(2). 2333 (if (not (e> exp (rczero))) 2334 (let-pw (get-datum (car (switch 'multi temp))) 2335 (e+ (current-trunc temp) (e- (e1- exp))) 2336 (var-expand (car (switch 'multi temp)) 1 't)) 2337 (var-expand (car (switch 'multi temp)) 1 't)) 2338 (cons (list (int-gvar temp) 1 1) 1)) 2339 (taylor2 (exp-pt temp))) 2340 exp)) 2341 ((signp e (exp-pt temp)) 2342 (let ((exp>trunc? () )) 2343 (if (and (e> exp (current-trunc temp)) (setq exp>trunc? 't) 2344 (not dont-truncate?)) 2345 (rczero) 2346 (make-ps (int-var temp) 2347 (ncons (if exact-poly (inf) 2348 (if exp>trunc? exp (current-trunc temp)))) 2349 (ncons (term (if (switch '$asymp temp) (rcminus exp) 2350 exp) 2351 (rcone))))))) 2352 (t (psexpt (psplus 2353 (make-ps (int-var temp) 2354 (ncons (if exact-poly (inf) (current-trunc temp))) 2355 (ncons (term (if (switch '$asymp temp) 2356 (rcmone) 2357 (rcone)) 2358 (rcone)))) 2359 (taylor2 (exp-pt temp))) 2360 exp))))) 2361 2362(defun expand (arg func) 2363 (or (try-expansion arg func) (exp-pt-err))) 2364 2365(defun try-expansion (arg func) 2366 (prog (funame funord fun-lc argord psarg arg-trunc temp exact-poly) 2367 ;; We bind exact-poly to () since we don't want psexpt retaining 2368 ;; higher order terms when subst'ing into series (which aren't exact). 2369 ;; Try diff-expanding unknown subsripted functions. 2370 (unless (or (atom func) (known-ps (caar func))) 2371 (taylor2 (diff-expand `((mqapply) ,func ,arg) tlist))) 2372 (when (setq temp (get (setq funame (oper-name func)) 'tay-trans)) 2373 (return (funcall temp arg func))) 2374 (let ((lterm (getfun-lt func))) 2375 (setq funord (e lterm) fun-lc (c lterm))) 2376 begin-expansion 2377 (when (rczerop (or psarg (setq psarg (get-lexp arg (rcone) () )))) 2378 (if (e> (rczero) funord) 2379 (if (rczerop (setq psarg (get-lexp arg (rcone) 't))) 2380 (tay-depth-err) 2381 (go begin-expansion)) 2382 (return (cond ((setq temp (assoc funame tay-pole-expand :test #'eq)) 2383 (funcall (cdr temp) arg psarg func)) 2384 ((rczerop funord) fun-lc) 2385 (t (rczero)))))) 2386 (when (pscoefp psarg) (setq psarg (taylor2 arg))) 2387 (when (pscoefp psarg) 2388 (return 2389 (cond ((null (mfree (rdis psarg) tvars)) 2390 (symbolic-expand arg psarg func)) 2391 ((setq temp (assoc funame tay-pole-expand :test #'eq)) 2392 (funcall (cdr temp) arg psarg func)) 2393 (t (prep1 (simplify 2394 (if (atom func) `((,func) ,(rcdisrep psarg)) 2395 `((mqapply) ,func ,(rcdisrep psarg))))))))) 2396 (when (e> (rczero) (setq argord (ps-le psarg))) 2397 (cond ((not (member funame '(%atan %asin %asinh %atanh) :test #'eq)) 2398 (if (e> (rczero) (ps-le* (setq psarg (get-lexp arg (rcone) 't)))) 2399 (essen-sing-err) 2400 (go begin-expansion))) 2401 (t 2402 (if (and (eq funame '%atan) 2403 (eq (coef-sign arg) '$neg)) 2404 (return (psplus (atrigh arg func) (taylor2 (m- '$%pi)))) 2405 (return (atrigh arg func)))))) 2406 (setq temp (t-o-var (gvar psarg))) 2407 (when (e> (e* funord argord) temp) (return (rczero))) 2408 ;; the following form need not be executed if psarg is really exact. 2409 ;; The constant problem does not allow one to determine this now, 2410 ;; so we always have to execute this currently. 2411 ;; This really should be 2412 ;; (unless (infp (trunc-lvl psarg)) ... ) 2413 ;; Likewise, the infp checks shouldn't be there; have to assume 2414 ;; nothing is exact until constant problem is fixed. 2415 (setq arg-trunc (if (and (not (infp (trunc-lvl psarg))) 2416 (e= funord (rcone))) 2417 temp 2418 (e- temp (e* (e1- funord) argord))) 2419 psarg (let-pw (get-datum (get-key-var (gvar psarg))) 2420 arg-trunc 2421 (if (or (infp (trunc-lvl psarg)) 2422 (e> arg-trunc (trunc-lvl psarg))) 2423 (taylor2 arg) 2424 (pstrunc psarg))) 2425 ;; We must recalculate argord since pstrunc may have "picked" 2426 ;; a coeff out of a constant monomial; e.g. this occurs in 2427 ;; taylor(sin(x+y),x,0,0,y,0,1) where psarg is (Y+...)*X^0+... 2428 ;; which truncates to Y+... of order 1. 2429 argord (ps-le* psarg)) 2430 (if (rczerop argord) 2431 (cond ((member funame '(%atan %asin %asinh %atanh) :test #'eq) 2432 (return (atrigh arg func))) 2433 ((setq temp (assoc funame const-exp-funs :test #'eq)) 2434 (return (funcall (cdr temp) arg psarg func))) 2435 ((rczerop (ps-le* (setq psarg (get-lexp arg (rcone) 't)))) 2436 (return () )) ; Don't know an addition formula 2437 (t (go begin-expansion))) 2438 (return 2439 (if (mono-term? (terms psarg)) 2440 (get-series func (current-trunc 2441 (get-datum (get-key-var (gvar psarg)))) 2442 (gvar-o psarg) (ps-le psarg) (ps-lc psarg)) 2443 (progn 2444 (setq temp (get-series func 2445 (e// temp argord) (gvar-o psarg) 2446 (rcone) (rcone))) 2447 (cond ((not (psp temp)) temp) 2448 (t (pscsubst1 psarg temp))))))))) 2449 2450(defun symbolic-expand (ign psarg func) ; should be much stronger 2451 (declare (ignore ign)) 2452 (prep1 (simplifya (if (atom func) 2453 `((,func) ,(rcdisrep psarg)) 2454 `((mqapply) ,func ,(rcdisrep psarg))) 2455 () ))) 2456 2457(defun expand-sing-trig? (arg func) 2458 (cond ((member func *pscirc :test #'eq) (tay-exponentialize arg func)) 2459 ((member func *psacirc :test #'eq) (atrigh arg func)) 2460 (t (essen-sing-err)))) 2461 2462(defun trig-const (a arg func) 2463 (let ((const (ps-lc* arg)) (temp (cdr (assoc func trigdisp :test #'eq)))) 2464 (cond ((and (pscoefp const) 2465 (member func '(%tan %cot) :test #'eq) 2466 (multiple-%pi a (srdis const) func))) 2467 (temp (funcall temp (setq const (psdisrep const)) 2468 (m- a const))) 2469 (t (tsexpt `((,(get func 'recip)) ,(srdis arg)) -1))))) 2470 2471(defun multiple-%pi (a const func) 2472 (let (coef) 2473 (and (equal ($hipow const '$%pi) 1) 2474 ($ratnump (setq coef ($ratcoef const '$%pi 1))) 2475 (cond ((numberp coef) (expand (m- a const) func)) 2476 ((equal (caddr coef) 2) 2477 (psminus (expand (m- a const) 2478 (cond ((eq func '%tan) '%cot) 2479 ((eq func '%cot) '%tan) 2480 (t (merror "MULTIPLE-%PI: internal error in Taylor expansion.")))))))))) 2481 2482(setq *pscirc '(%cot %tan %csc %sin %sec %cos %coth 2483 %tanh %csch %sinh %sech %cosh) 2484 2485 *psacirc '(%acot %atan %acsc %asin %asec %acos %acoth 2486 %atanh %acsch %asinh %asech %acosh)) 2487 2488(setq const-exp-funs 2489 `((%gamma . gam-const) ($psi . plygam-const) 2490 . ,(mapcar #'(lambda (q) (cons q 'trig-const)) *pscirc)) 2491 2492 trigdisp '((%sin . psina+b) (%cos . pscosa+b) (%tan . pstana+b) 2493 (%sinh . psinha+b) (%cosh . pscosha+b) (%tanh . pstanha+b)) 2494 2495 tay-pole-expand '((%gamma . plygam-pole) ($psi . plygam-pole)) 2496 2497 tay-const-expand ; !these should be handled by symbolic-expand 2498 ; be sure to change this with tay-exponentialize! 2499 (append (mapcar #'(lambda (q) (cons q 'tay-exponentialize)) *pscirc) 2500 (mapcar #'(lambda (q) (cons q 'tay-exponentialize)) *psacirc))) 2501 2502(mapc #'(lambda (q) (putprop q 'atrig-trans 'tay-trans)) 2503 '(%acos %acot %asec %acsc %acosh %acoth %asech %acsch)) 2504 2505(defprop mfactorial factorial-trans tay-trans) 2506 2507(defun factorial-trans (arg func) 2508 (declare (ignore func)) 2509 (taylor2 `((%gamma) ,(m1+ arg)))) 2510 2511 2512(defprop %gamma_incomplete gamma-upper-trans tay-trans) 2513(defprop $gamma_incomplete gamma-upper-trans tay-trans) 2514(defprop %gamma_incomplete_lower gamma-lower-trans tay-trans) 2515(defprop $gamma_incomplete_lower gamma-lower-trans tay-trans) 2516 2517;; for gamma_incomplete(s,z) 2518;; translate into gamma_incomplete_lower if s>0 and z=0 2519(defun gamma-upper-trans (arg func) 2520 (let ((s (car arg)) 2521 (z (cadr arg))) 2522 (if (and 2523 (eq ($sign s) '$pos) 2524 (zerop1 ($limit z (caar tlist) (exp-pt (car tlist))))) 2525 (taylor2 `((mplus) ((%gamma) ,s) 2526 ((mtimes) -1 ((%gamma_incomplete_lower) ,s ,z)))) 2527 (taylor2 (diff-expand `((,func) . ,arg) 2528 tlist))))) 2529 2530;; for gamma_incomplete_lower(s,z) 2531;; if z=0, use A&S 6.5.29 2532;;; inf 2533;;; === 2534;;; \ (-z)^k 2535;;; gamma_incomplete_lower(s,z) = z^s * > ------------ 2536;;; / (s+k) k! 2537;;; === 2538;;; k=0 2539(defun gamma-lower-trans (arg func) 2540 (let ((s (car arg)) 2541 (z (cadr arg))) 2542 (if (zerop1 ($limit z (caar tlist) (exp-pt (car tlist)))) 2543 (taylor2 `((mtimes) 2544 ((mexpt) ,z ,s) 2545 ((%sum) 2546 ((mtimes) 2547 ((mexpt) ((mtimes) -1 ,z) k) 2548 ((mexpt) ((mtimes) ((mfactorial) k) 2549 ((mplus) ,s k)) 2550 -1)) 2551 k 2552 0 2553 $inf))) 2554 (taylor2 (diff-expand `((,func) . ,arg) 2555 tlist))))) 2556 2557 2558;;; Not done properly yet 2559;;; 2560;;; (defprop $BETA BETA-TRANS TAY-TRANS) 2561 2562(defun psina+b (a b) 2563 (psplus (pstimes (expand a '%sin) (expand b '%cos)) 2564 (pstimes (expand a '%cos) (expand b '%sin)))) 2565 2566(defun pscosa+b (a b) 2567 (psdiff (pstimes (expand a '%cos) (expand b '%cos)) 2568 (pstimes (expand a '%sin) (expand b '%sin)))) 2569 2570(defun pstana+b (a b) 2571 (setq a (expand a '%tan) b (expand b '%tan)) 2572 (pstimes (psplus a b) 2573 (psexpt (psdiff (rcone) (pstimes a b)) 2574 (rcmone)))) 2575 2576(defun psinha+b (a b) 2577 (psplus (pstimes (expand a '%sinh) (expand b '%cosh)) 2578 (pstimes (expand a '%cosh) (expand b '%sinh)))) 2579 2580(defun pscosha+b (a b) 2581 (psplus (pstimes (expand a '%cosh) (expand b '%cosh)) 2582 (pstimes (expand a '%sinh) (expand b '%sinh)))) 2583 2584(defun pstanha+b (a b) 2585 (setq a (expand a '%tanh) b (expand b '%tanh)) 2586 (pstimes (psplus a b) 2587 (psexpt (psplus (rcone) (pstimes a b)) 2588 (rcmone)))) 2589 2590(defun atrig-trans (arg func) 2591 (taylor2 2592 (cond ((eq func '%acos) 2593 `((mplus) ,half%pi ((mtimes) -1 ((%asin) ,arg)))) 2594 2595 ((eq func '%acosh) 2596 `((mtimes) -1 $%i ((mplus) ,half%pi ((mtimes) -1 ((%asin) ,arg))))) 2597 2598 (t 2599 `((,(cdr (assoc func '((%acsc . %asin) (%asec . %acos) 2600 (%acot . %atan) (%acsch . %asinh) 2601 (%asech . %acosh) (%acoth . %atanh)) :test #'eq))) 2602 ,(m^ arg -1)))))) 2603 2604(defun atrigh (arg func) 2605 (let ((full-log t) ($logarc t) (log-1 '((mtimes) $%i $%pi)) 2606 (log%i '((mtimes) ((rat) 1 2) $%i $%pi))) 2607 (taylor2 (simplify `((,func) ,arg))))) 2608 2609(defun tay-exponentialize (arg fun) ; !this should be in symbolic-expand! 2610 (let (($exponentialize t) ($logarc t)) 2611 (setq arg (meval `((,fun) ,arg)))) 2612 (taylor2 arg)) 2613 2614(defun tsplus (l) 2615 (do ((l (cdr l) (cdr l)) 2616 (ans (taylor2 (car l)) 2617 (psplus ans (taylor2 (car l))))) 2618 ((null l) ans))) 2619 2620(defun ts-formula (form var pw) 2621 (let ((datum (get-datum (get-key-var (car var))))) 2622 (let-pw datum pw 2623 (taylor2 (subst (get-inverse (car var)) 'sp2var form))))) 2624 2625(defmacro next-series (l) `(cdadr ,l)) 2626 2627(defun tstimes-get-pw (l pw) 2628 (do ((l l (cdr l)) (vect)) 2629 ((null l) pw) 2630 (setq pw (mapcar #'(lambda (pwq ple) (e+ pwq ple)) 2631 pw (setq vect (ord-vector (cdar l))))) 2632 (rplacd (car l) (cons (cdar l) vect)))) 2633 2634(defun tstimes-l-mult (a) 2635 (do ((l (cdr a) (cdr l)) ($maxtayorder t) 2636 (a (car a) (pstimes a (car l)))) 2637 ((null l) a))) 2638 2639(defun mzfree (e l) 2640 (do ((l l (cdr l))) 2641 ((null l) 't) 2642 (or (zfree e (car l)) (return () )))) 2643 2644;;; The lists posl, negl and zerl have the following format: 2645;;; 2646;;; ( (<expression> <expansion> <order vector>) . . . ) 2647 2648(defun tstimes (l) 2649 (*bind* ((funl) (expl) (negl) (zerl) (posl) 2650 (pw) (negfl) (temp) (fixl (rcone))) 2651 (dolist (fun l) ;; find the exponentials 2652 (if (mexptp fun) 2653 (push (if (free (caddr fun) (car tvars)) fun 2654 `((mexpt) $%e ,(m* (caddr fun) 2655 `((%log) ,(cadr fun))))) 2656 expl) 2657 (push fun funl))) 2658 (when expl 2659 (setq expl (tsexp-comb expl)) ;; simplify exps 2660 (setq expl (tsbase-comb expl))) ;; and again 2661 (setq l (nconc expl funl)) ;; now try expanding 2662 (let ((trunc-constants? () )) 2663 (setq expl (cons 0 (mapcar #'(lambda (exp) 2664 (cons exp (taylor2 exp))) 2665 l))) ) 2666 ;; EXPL is now of the form (0 ( <form> . <taylor2(form)> ) ...) 2667 ;; L points behind the cons considered for destructive updating. 2668 (do ((l expl) (tem)) 2669 ((null (cdr l)) ) 2670 (cond ((rczerop (cdadr l)) 2671 ;; Consider taylor((a+1/x)*1/x,x,0,-2). Each factor will be on 2672 ;; zerl. Each factor will also appear to have le = 0 since its 2673 ;; series is 0, which would fool the get-pw routines below if 2674 ;; they tried to handle this case. The easiest fix for now 2675 ;; appears to be to always call get-lexp here, killing this: 2676 (cond ;((null $maxtayorder) 2677 ; (setq zerl (cons (cadr l) zerl)) 2678 ; (rplacd l (cddr l))) 2679 ((rczerop (setq tem (get-lexp (caadr l) (rcone) ()))) 2680 (return (setq zerl 0))) 2681 ('t (setq posl (cons (cons (caadr l) tem) posl)) 2682 (rplacd l (cddr l))))) 2683 ((pscoefp (cdadr l)) 2684 (cond ((mzfree (caadr l) tvars) ;must be zfree to permit ratfuns 2685 (setq fixl (pstimes (cdadr l) fixl)) 2686 (rplacd l (cddr l))) 2687 ((setq zerl (cons (cadr l) zerl)) 2688 (rplacd l (cddr l))))) 2689 ((rczerop (ps-le (cdadr l))) 2690 (setq zerl (cons (cadr l) zerl)) 2691 (rplacd l (cddr l))) 2692 ((e> (ps-le (cdadr l)) (rczero)) 2693 (setq posl (cons (cadr l) posl)) 2694 (rplacd l (cddr l))) 2695 ('t (setq l (cdr l))))) 2696 (when (equal zerl 0) (return (rczero))) 2697 (setq negl (cdr expl) temp (ord-vector fixl)) 2698 (mapcar #'(lambda (x) (and (e> (rczero) x) (setq negfl t))) temp) 2699 (tstimes-get-pw zerl temp) 2700 (setq pw (tstimes-get-pw posl (tstimes-get-pw negl temp))) 2701 (if (or negl negfl) 2702 (setq posl 2703 (mapcar #'(lambda (x) 2704 (prog2 (mapcar #'(lambda (datum lel pwl) 2705 (push-pw datum 2706 (e+ (current-trunc datum) 2707 (e- lel pwl)))) 2708 tlist (cddr x) pw) 2709 (taylor2 (car x)) 2710 (mapcar #'(lambda (datum) (pop-pw datum)) 2711 tlist))) 2712 (nconc posl zerl negl))) 2713 (setq posl (nconc (mapcar 'cadr posl) (mapcar 'cadr zerl) 2714 (mapcar 'cadr negl)))) 2715 (setq posl (tstimes-l-mult posl)) 2716 (let ((ans (cond ((null fixl) posl) 2717 ((null posl) fixl) 2718 ('t (pstimes fixl posl))))) 2719 (if $maxtayorder ans (pstrunc ans))))) 2720 2721;;; This routine transforms a list of exponentials as follows: 2722;;; 2723;;; a^c*b^(n*c) ===> (a*b^n)^c, where n is free of var. 2724;;; 2725;;; This transformation is only applicable when c is not free of var. 2726 2727(defun tsexp-comb (l) ;; ***** clobbers l ***** 2728 (setq l (cons '* l)) 2729 (do ((a l) (e)) ;; updated by a rplacd or cdr. 2730 ((null (cddr a)) (cdr l)) ;; get rid of the * 2731 (rplaca (cddadr a) (setq e ($ratsimp (caddr (cadr a))))) 2732 ;; Must delete e^0 lest we divide by the 0 below. RWG's byzero bug 2733 ;; of 3/1/78 used to cause this. 2734 (if (equal e 0) (rplacd a (cddr a)) 2735 (if (mfree (caddr (cadr a)) tvars) (pop a) 2736 (do ((b (cddr a) (cdr b)) (n)) 2737 ((null b) (setq a (cdr a))) 2738 (when (mfree (setq n ($ratsimp (m// (caddar b) 2739 (caddr (cadr a))))) 2740 tvars) 2741 (rplaca b (list '(mexpt simp) 2742 (m* (cadadr a) 2743 (m^ (cadar b) n)) ;; b^n 2744 (caddr (cadr a)))) 2745 (rplacd a (cddr a)) ;; delete a^c 2746 (return () ))))))) 2747 2748;;; This routine transforms a list of exponentials as follows: 2749;;; 2750;;; a^b*a^c ===> a^(b+c), 2751;;; 2752;;; this is only necessary when b and c depend on "var." 2753 2754(defun tsbase-comb (l) ;;; *******clobbers l******** 2755 (setq l (cons '* l)) 2756 (do ((a l)) ;;; updated by a rplacd or cdr 2757 ((null (cddr a)) (cdr l)) 2758 (do ((b (cddr a) (cdr b))) 2759 ((null b) (pop a)) ;;; did not return early so pop. 2760 (when (alike1 (cadar b) (cadadr a)) 2761 (rplaca b (m^ (cadar b) (m+ (caddar b) (caddr (cadr a))))) 2762 (rplacd a (cddr a)) 2763 (return 't))))) 2764 2765(defun tsexpt (b e) 2766 (cond ((and (atom b) (mnump e) 2767 (get-datum b) 2768 (not (eq (exp-pt (get-datum b)) '$minf))) 2769 ;; one could remove this clause and let this case be handled by tsexpt1 2770 (var-expand b e () )) 2771 ((mfree e tvars) (tsexpt1 b e)) 2772 ((eq '$%e b) (tsexpt-red (list e))) 2773 (t (tsexpt-red (list (list '(%log) b) e))))) 2774 2775(defun tsexpt-red (l) 2776 (*bind* ((free) (nfree) (full-log) ($logarc t) (expt) (ps) (e) 2777 (log-1 '((mtimes) $%i $%pi)) 2778 (log%i '((mtimes) ((rat) 1 2) $%i $%pi))) 2779 (declare (special e)) 2780 a (do ((l l (cdr l))) 2781 ((null l) ) 2782 (cond ((mtimesp (car l)) (setq l (append l (cdar l)))) 2783 ((mfree (car l) tvars) (push (car l) free)) 2784 (t (push (car l) nfree)))) 2785 (cond ((or (cdr nfree) (atom (car nfree))) ) 2786 ((eq (caaar nfree) '%log) 2787 (return (tsexpt1 (cadar nfree) (m*l free)))) 2788 ((member (caaar nfree) *psacirc :test #'eq) 2789 (setq l (ncons (simplifya ;; simplify after removing simp flag 2790 (cons (ncons (caaar nfree)) (cdar nfree)) 2791 () )) 2792 nfree (cdr nfree)) 2793 (go a))) 2794 ;; Must have truncs > 0 so that logs in the expt aren't trunc'd. 2795 ;; E.g, consider taylor(x^(x-1),x,0,-1). 2796 (tlist-mapc d (push-pw d (emax (current-trunc d) (rcone)))) 2797 (setq ps (taylor2 (setq expt (m*l (append nfree free))))) 2798 (tlist-mapc d (pop-pw d)) 2799 ;; Here we must account for the truncation gain or lossage that 2800 ;; is encountered in exp(c*log(x)+y) -> x^c*exp(y). 2801 (let ((c0 (if (pscoefp ps) ps (psterm (terms ps) (rczero)))) 2802 e^c0 ord-e^c0) 2803 (unless (rczerop c0) 2804 (setq ord-e^c0 (ord-vector (setq e^c0 (psexpt-fn c0)))) 2805 ;; Must emax with 0 so that new singular kernals won't be trunc'd 2806 ;; e.g exp(1/x+...) to degree -2 should be exp(-1/x)+... 2807 ;; Also try taylor(screwa,x,0,-2). 2808 (mapc #'(lambda (d o) (push-pw d (emax (e- (current-trunc d) o) 2809 (rczero)))) 2810 tlist ord-e^c0) 2811 (setq ps (psdiff (taylor2 expt) c0))) 2812 (setq ps (psexpt-fn ps)) 2813 (when e^c0 2814 (tlist-mapc d (pop-pw d)) 2815 (setq ps (pstimes e^c0 ps))) 2816 (pstrunc ps)))) 2817 2818;; Taylor's b^e, where e is independent of tvars. 2819 2820(defun tsexpt1 (b e) 2821 (prog (s le pw tb) 2822 (setq e (let ((modulus () )) ; Don't mod exponents! See WGM's bug 2823 (prog2 (mapcar ; of 3/6/83 for an example. 2824 #'(lambda (datum) 2825 (push-pw datum 2826 (emax (current-trunc datum) (rczero)))) 2827 tlist) 2828 (taylor2 e) 2829 (mapcar #'(lambda (datum) (pop-pw datum)) tlist))) 2830 s (psfind-s e) 2831 tb (taylor2 b) 2832 pw (if (psp tb) (current-trunc (get-datum 2833 (get-key-var (gvar tb)))) 2834 ;; Constant problem kludge. 2835 (if (rczerop tb) (current-trunc (car tlist)) (rczero)))) 2836 (if (floatp (car s)) 2837 (setq s (maxima-rationalize (quot (car s) (cdr s))))) 2838 ;; We must ensure that the lc is non-zero since it will be inverted in 2839 ;; psexpt. 2840 (setq tb (strip-zeroes tb 't)) 2841 (cond ((rczerop tb) 2842 (when (or ;; When 1 > s we need more terms since -le*(s-1) > 0. 2843 (e> (rcone) s) 2844 (and (e> (rczero) pw) (e> s (rcone)))) 2845 (setq tb (get-lexp b () 't))) 2846 (setq le (ps-le* tb))) 2847 ((psp tb) (setq le (ps-le tb))) 2848 (t (return (rcexpt tb e)))) 2849 (and (e> (e* s le) pw) (null $maxtayorder) (return (rczero))) 2850 (setq s (e- pw (e* le (e1- s)))) 2851 ;(setq le (increment-truncs tb)) 2852 (return 2853 (psexpt 2854 (if (e> pw s) 2855 (if $maxtayorder tb 2856 (pstrunc1 tb (list (cons (gvar tb) s)))) 2857 ;; because of constants not retaining info, have to 2858 ;; just keep the constant here 2859 (cond ((not (psp tb)) tb) 2860 (t (let-pw (get-datum (get-key-var (gvar tb))) s (strip-zeroes (taylor2 b) 't))))) 2861 e)))) 2862 2863;;; the method of calculating truncation levels below is incorrect. 2864;;; (i.e. increment-truncs & decrement-truncs, also used above) 2865;;; Examples which exhibit this incorrectness are: 2866;;; taylor(log(sin(y)+x),x,0,2,y,0,1) is missing a y/6*x and -1/6*x^2 2867;;; taylor(log(sin(z)+sin(y)+x),x,0,f1,y,0,3,z,0,5) misses a z^5*y^3 term. 2868 2869;;; TSLOG must find the lowest degree term in the expansion of the 2870;;; log arg, then expand with the orders of all var's in this low term 2871;;; incremented by their order in this low term. Note that this is 2872;;; only necessary for var's with ord > 0, since otherwise we have 2873;;; already expanded to a higher ord than required. Also we must 2874;;; not do this for var's with trunc < 0, since this may incorrectly 2875;;; truncate terms which should end up as logs. 2876 2877(defun increment-truncs (ps) 2878 (do ((ps ps (ps-lc ps)) (trunc (t-o-var (gvar ps))) (data () )) 2879 ((pscoefp ps) data) 2880 (when (e> (ps-le ps) (rczero)) 2881 (push (assoc (get-key-var (gvar ps)) tlist :test #'eq) data) 2882 (push-pw (car data) (e+ (e* (e+ trunc (rctwo)) (ps-le ps)) 2883 (current-trunc (car data)))) 2884 (setq trunc (e+ trunc (current-trunc (car data)))) 2885 ))) 2886 2887(defun decrement-truncs (data) 2888 (mapc #'(lambda (data) (pop-pw data)) data)) 2889 2890(defun tslog (arg) 2891 (let ((psarg (taylor2 arg)) datum) 2892 (when (rczerop psarg) (setq psarg (get-lexp arg () 't))) 2893 ;; We must ensure that the lc is non-zero since it will be inverted in pslog 2894 (setq psarg (strip-zeroes psarg 't)) 2895 (do ((ps psarg (ps-lc ps)) (shift (rcone) (e* shift (rctwo)))) 2896 ((pscoefp ps) 2897 (when datum 2898 (when (rczerop (setq psarg (taylor2 arg))) 2899 (setq psarg (get-lexp arg () 't))) 2900 (mapc #'(lambda (data) (pop-pw data)) datum)) 2901 (pslog psarg)) 2902 (push (get-datum (get-key-var (gvar ps))) datum) 2903 (if (and (e> (ps-le ps) (rczero)) 2904 (e> (current-trunc (car datum)) (rczero))) 2905 (push-pw (car datum) (e+ (e* shift (ps-le ps)) 2906 (current-trunc (car datum)))) 2907 (pop datum))))) 2908 2909;; When e-start is non-null we start expanding at order e-start, ... , 2^m, 2910;; then 2^m*pow, instead of the normal sequence pow, ... , 2^m*pow 2911;; (where m = $taylordepth, pow = ord of var). This is done because it is 2912;; usually much more efficient for large, non-trivial expansions when we only 2913;; want the lowest order term. 2914 2915(defun get-lexp (exp e-start zerocheck?) 2916 (if (equal exp 0) 2917 (if zerocheck? 2918 (tay-depth-err) 2919 (rczero)) 2920 (progn 2921 (tlist-mapc d (push-pw d (or e-start (emax (orig-trunc d) (rcone))))) 2922 (do ((psexp) (i (1+ $taylordepth) (1- i))) 2923 ((signp e i) 2924 (tlist-mapc d (pop-pw d)) 2925 (if zerocheck? 2926 (tay-depth-err) 2927 (progn 2928 (unless silent-taylor-flag (zero-warn exp)) 2929 (rczero)))) 2930 (declare (fixnum i)) 2931 (cond ((and (rczerop (setq psexp (if zerocheck? 2932 (strip-zeroes (taylor2 exp) 't) 2933 (taylor2 exp)))) 2934 (not (member exp zerolist :test #'eq))) ) 2935 ;; Info not needed yet. 2936 ;; ((and lexp-non0 (rczerop (le (terms psexp))) 2937 ;; (mono-term? (terms psexp)))) 2938 (t (tlist-mapc d (pop-pw d)) 2939 (return psexp))) 2940 (cond ((and (= i 1) e-start) 2941 (setq e-start () i 2) 2942 (tlist-mapc d (push-pw d (prog1 (e* (orig-trunc d) (current-trunc d)) 2943 (pop-pw d))))) 2944 (t (tlist-mapc d (push-pw d (prog1 (e* (rctwo) (current-trunc d)) 2945 (pop-pw d)))))))))) 2946 2947(defun 1p (x) 2948 (or (equal x 1) (equal x 1.0))) 2949 2950(defun [max-trunc] () 2951 (do ((l tlist (cdr l)) (emax (rczero))) 2952 ((null l) (1+ (truncate (car emax) (cdr emax)))) 2953 (when (e> (current-trunc (car l)) emax) 2954 (setq emax (orig-trunc (car l)))))) 2955 2956(defun tsprsum (f l type) 2957 (if (mfree f tvars) (newsym f) 2958 (let ((li (ncons (car l))) (hi (caddr l)) (lv (ncons (cadr l))) a aa 2959 ($maxtayorder () ));; needed to determine when terms are 0 2960 (if (and (numberp (car lv)) (numberp hi) (> (car lv) hi)) 2961 (if (eq type '%sum) (taylor2 0) (taylor2 1)) 2962 (progn 2963 (if (eq type '%sum) (setq type '())) 2964 (do ((m (* ([max-trunc]) (ash 1 $taylordepth))) 2965 (k 0 (1+ k)) 2966 (ans (taylor2 (maxima-substitute (car lv) (car li) f)))) 2967 ((equal hi (car lv)) ans) 2968 (rplaca lv (m1+ (car lv))) 2969 ;; A cheap heuristic to catch infinite recursion when 2970 ;; possible, should be improved in the future 2971 (if (> k m) (exp-pt-err) 2972 (setq a ;(mlet li lv (taylor2 (setq aa (meval f)))) 2973 (taylor2 (maxima-substitute (car lv) (car li) f)))) 2974 (if type 2975 (if (and (1p (car a)) (1p (cdr a)) (not (1p aa))) 2976 (return ans) 2977 (setq ans (pstimes a ans))) 2978 (if (and (rczerop a) (not (signp e aa))) 2979 (return ans) 2980 (setq ans (psplus ans a)))))))))) 2981 2982(defun tsdiff (e l check) 2983 (*bind* ((n) (v) (u)) 2984 (do ((l l (cddr l))) 2985 ((null l)) 2986 (if (and (atom (car l)) (numberp (cadr l)) 2987 (assoc (car l) tlist :test #'eq)) 2988 (setq n (cons (cadr l) n) v (cons (car l) v)) 2989 (setq u (cons (car l) (cons (cadr l) u))))) 2990 (or n (return (prep1 check))) 2991 (if u (setq e (meval (cons '($diff) (cons e l))))) 2992 (setq l (mapcar #'(lambda (x) (get-datum x)) v)) 2993 (mapcar #'(lambda (datum pw) 2994 (push-pw datum (e+ (current-trunc datum) (prep1 pw)))) 2995 l n) 2996 (setq e (taylor2 e)) 2997 (mapc #'(lambda (datum) (pop-pw datum)) l) 2998 (do ((vl v (cdr vl)) 2999 (nl n (cdr nl))) 3000 ((null vl ) e) 3001 (do ((i 1 (1+ i))) 3002 ((> i (car nl)) ) 3003 (mapc #'(lambda (a b) 3004 (putprop a (prep1 (sdiff b (car v))) 3005 'diff)) 3006 genvar varlist) 3007 (setq e (psdp e)))))) 3008 3009 3010(defun no-sing-err (x) ;; try to catch all singularities 3011 (let ((errorsw t)) 3012 (declare (special errorsw)) 3013 (let ((ans (catch 'errorsw (eval x)))) 3014 (if (eq ans t) (unfam-sing-err) ans)))) 3015 3016;; evaluate deriv at location var=pt 3017(defun eval-deriv (deriv var pt) 3018 (let ((errorsw t)) 3019 (declare (special errorsw)) 3020 (let ((ans (no-sing-err `(meval '(($at) ,deriv ((mequal) ,var ,pt)))))) 3021 ans))) 3022 3023(defun check-inf-sing (pt-list) ; don't know behavior of random fun's @ inf 3024 (and (or (member '$inf pt-list :test #'eq) (member '$minf pt-list :test #'eq)) 3025 (unfam-sing-err))) 3026 3027(defun diff-expand (exp l) ;l is tlist 3028 (check-inf-sing (mapcar (function caddr) l)) 3029 (cond ((not l) exp) 3030 (t 3031 (setq exp (diff-expand exp (cdr l))) 3032 (do ((deriv (sdiff exp (caar l)) (sdiff deriv var)) 3033 (var (caar l)) 3034 (coef 1 (* coef (1+ cnt))) 3035 (cnt 1 (1+ cnt)) 3036 (pt (exp-pt (car l))) 3037 (lim (rcdisrep (current-trunc (car l)))) 3038 (ans (list (no-sing-err `(meval '(($at) ,exp ((mequal) ,(caar l) ,(exp-pt (car l))))))) 3039 (cons `((mtimes) ((rat simp) 1 ,coef) 3040 ,(eval-deriv deriv var pt) 3041 ((mexpt) ,(sub* var pt) ,cnt)) 3042 ans))) 3043 ((or (great cnt lim) (equal deriv 0)) (cons '(mplus) ans)))))) 3044 3045;;; subtitle disreping routines 3046 3047(defun edisrep (e) 3048 (if (= (cdr e) 1) (car e) (list '(rat) (car e) (cdr e)))) 3049 3050(defun striptimes (a) 3051 (if (mtimesp a) (cdr a) (ncons a))) 3052 3053(defun srdis (x) 3054 (let (($psexpand () )) ; Called only internally, no need to expand. 3055 ($ratdisrep 3056 (cons (list 'mrat 'simp varlist genvar tlist 'trunc) 3057 x)))) 3058 3059(defun srdisrep (r) 3060 (let ((varlist (mrat-varlist r)) (genvar (mrat-genvar r))) 3061 (mapc #'(lambda (exp genv) (putprop genv exp 'disrep)) 3062 varlist genvar) 3063 (setup-multivar-disrep r) 3064 ;; This used to return 0 if psdisrep returned () but this is wrong 3065 ;; since taylor(false,x,0,0) would lose. If psdisrep really wants to 3066 ;; return () for 0 then we will probably find out soon. 3067 (if (eq $psexpand '$multi) (psdisexpand (cdr r)) 3068 (psdisrep (cdr r))))) 3069 3070(defun psdisrep (p) 3071 (if (psp p) 3072 (psdisrep+ (psdisrep2 (terms p) (getdisrep (gvar-o p)) (trunc-lvl p)) 3073 (if (or $psexpand (trunc-lvl p)) '(mplus trunc) 3074 '(mplus exact))) 3075 (rcdisrep p))) 3076 3077(defun psdisrep^ (n var) 3078 ;; If var = () then it is an internal var generated in a multivariate 3079 ;; expansion so it shouldn't be displayed. If var = 1 then it probably 3080 ;; resulted from the substitution in srdisrep, so it depends on an 3081 ;; internal var and likewise shouldn't be displayed. 3082 (cond ((or (rczerop n) (null var) (equal var 1)) 1) 3083 ((equal n (rcone)) var) 3084 ((and ps-bmt-disrep (mexptp var) (equal (caddr var) -1)) 3085 (psdisrep^ (e- n) (cadr var))) 3086 ('t `((mexpt ratsimp) ,var ,(edisrep n))))) 3087 3088;;; There used to be a hack below that would print a series consisting 3089;;; of merely one term as exact polynomial (i.e. no trailing "..."'s). 3090;;; This is, of course, wrong but the problem with the fix is that 3091;;; now exact things like taylor(y*x,x,0,f1,y,0,1) will display like 3092;;; (y+...) x+... because of the problem with $MAXTAYORDER being internally 3093;;; bound to ()---which causes exact things to look inexact, such as 3094;;; x and y above. See the comment above taylor* for the $MAXTAYORDER problem. 3095 3096(defun psdisrep+ (p plush &aux lowest-degree-term) 3097 (if;; An exact sum of one arg is just that arg. 3098 (and (null (cdr p)) (eq (cadr plush) 'exact)) 3099 (car p) 3100 (progn 3101 ;; Since the DISPLAY package prints trunc'd sum's arguments 3102 ;; from right to left we must put the terms of any constant term 3103 ;; in decreasing order. Note that only a constant (wrt to the 3104 ;; mainvar) term can be a term which is a sum. 3105 (when (mplusp (setq lowest-degree-term (car (last p)))) 3106 (rplacd lowest-degree-term (nreverse (cdr lowest-degree-term)))) 3107 (cons plush p)))) 3108 3109(defun psdisrep* (a b) 3110 (cond ((equal a 1) b) 3111 ((equal b 1) a) 3112 (t (cons '(mtimes ratsimp) 3113 (nconc (striptimes a) (striptimes b)))))) 3114 3115(defun psdisrep2 (p var trunc) 3116 (if (or $ratexpand $psexpand) (psdisrep2expand p var) 3117 (do ((a () (cons (psdisrep* (psdisrep (lc p)) (psdisrep^ (le p) var)) 3118 a)) 3119 (p p (cdr p))) 3120 ((or (null p) (e> (le p) trunc)) a)))) 3121 3122(defun psdisrep2expand (p var) 3123 (do ((p p (cdr p)) 3124 (l () (nconc (psdisrep*expand (psdisrep (lc p)) (psdisrep^ (le p) var)) 3125 l))) 3126 ((null p) l))) 3127 3128(defun psdisrep*expand (a b) 3129 (cond ((equal a 1) (list b)) 3130 ((equal b 1) (list a)) 3131 ((null (mplusp a)) 3132 (list (cons '(mtimes ratimes) (nconc (striptimes a) (striptimes b))))) 3133 ('t (mapcar #'(lambda (z) (psdisrep* z b)) 3134 (cdr a))))) 3135 3136 3137(defun psdisexpand (p) 3138 (let ((ans (ncons ()))) 3139 (declare (special ans)) ;used in pans-add 3140 (psdisexcnt p () (rczero)) 3141 (setq ans 3142 (nreverse 3143 (mapcar #'(lambda (x) (cond ((not (cddr x)) (cadr x)) 3144 (t (cons '(mplus trunc) (cdr x))))) 3145 (cdr ans)))) 3146 (cond ((not (cdr ans)) (car ans)) 3147 (t (cons '(mplus trunc) ans))))) 3148 3149(defun psdisexcnt (p l n) 3150 (if (psp p) 3151 (do ((var (getdisrep (gvar-o p))) (ll (terms p) (n-term ll))) 3152 ((null ll) ()) 3153 (if (rczerop (le ll)) (psdisexcnt (lc ll) l n) 3154 (psdisexcnt (lc ll) 3155 (cons (psdisrep^ (le ll) var) l) 3156 (e+ (le ll) n)))) 3157 (psans-add (cond ((not l) (rcdisrep p)) 3158 (t (psdisrep* (rcdisrep p) 3159 (cond ((not (cdr l)) (car l)) 3160 (t (cons '(mtimes trunc) l)))))) 3161 n))) 3162 3163(defun psans-add (exp n) 3164 (declare (special ans)) ;bound in psdisexpand 3165 (do ((l ans (cdr l))) 3166 ((cond ((null (cdr l)) (rplacd l (ncons (list n exp)))) 3167 ((e= (caadr l) n) (rplacd (cadr l) (cons exp (cdadr l)))) 3168 ((e> (caadr l) n) (rplacd l (cons (list n exp) (cdr l)))))))) 3169 3170(defun srconvert (r) 3171 (cond ((not (atom (caadr (cdddar r)))) 3172 (cons (car r) (psdisextend (cdr r)))) 3173 (t 3174 (*bind* ((trunclist (cadr (cdddar r))) 3175 (tlist) 3176 (gps) 3177 (temp) 3178 (vs (caddar r)) 3179 (gens (cadddr (car r)))) 3180 (setq gps (mapcar #'cons gens vs)) 3181 (do ((tl (cdr trunclist) (cddr tl))) 3182 ((null tl) (cons (list 'mrat 'simp vs gens tlist 'trunc) (srconvert1 (cdr r)))) 3183 (setq temp (cdr (assoc (car tl) gps :test #'eq))) 3184 (cond ((null (member (car tl) (cdr trunclist) :test #'eq))) 3185 ((mplusp temp) (merror "SRCONVERT: internal error.")) 3186 (t 3187 (setq tlist 3188 (cons (list* temp (tay-order (cadr tl)) 0 nil 3189 (cons (car tl) (symbol-value (car tl)))) 3190 tlist))))))))) 3191 3192(defun srconvert1 (p) 3193 (cond ((not (member (car p) genvar :test #'eq)) p) 3194 (t 3195 (do ((l (cdr p) (cddr l)) 3196 (a nil (cons (term (prep1 (car l)) (srconvert1 (cadr l))) a))) 3197 ((null l) 3198 (make-ps (cons (car p) (symbol-value (car p))) 3199 (tay-order (zl-get trunclist (car p))) a)))))) 3200 3201;;; subtitle error handling 3202 3203(defun tay-error (msg exp) 3204 (if silent-taylor-flag (throw 'taylor-catch ()) 3205 (if exp 3206 (merror "taylor: ~A~%~M" msg exp) 3207 (merror "taylor: ~A" msg)))) 3208 3209(defun exp-pt-err () 3210 (tay-err (intl:gettext "unable to expand at a point specified in:"))) 3211 3212(defun essen-sing-err () 3213 (tay-err (intl:gettext "encountered an essential singularity in:"))) 3214 3215(defun unfam-sing-err () 3216 (tay-err (intl:gettext "encountered an unfamiliar singularity in:"))) 3217 3218(defun infin-ord-err () 3219 (tay-err (intl:gettext "expansion to infinite order?"))) 3220 3221(defun tay-depth-err () 3222 (tay-err (intl:gettext "'taylordepth' exceeded while expanding:"))) 3223 3224;;; Subtitle TAYLORINFO 3225 3226(defun taylor-trunc (q) 3227 (setq q (current-trunc q)) 3228 (cond ((null q) '$inf) 3229 ((equal (cdr q) 1) (car q)) 3230 (t `((rat) ,(car q) ,(cdr q))))) 3231 3232(defun taylor-info (q) 3233 (let ((acc-var nil) (acc-pt nil) (acc-ord nil) (qk) (acc)) 3234 (cond ((null q) nil) 3235 (t 3236 (setq qk (pop q)) 3237 (cond ((and (fourth qk) (consp (fourth qk)) (eq (caar (fourth qk)) 'multivar)) nil) 3238 ((and (fourth qk) (consp (fourth qk)) (eq (caar (fourth qk)) 'multi)) 3239 (while (and (fourth qk) (consp (fourth qk)) (eq (caar (fourth qk)) 'multi)) 3240 (setq acc nil) 3241 (push (taylor-trunc qk) acc-ord) 3242 (push (exp-pt qk) acc-pt) 3243 (push (datum-var qk) acc-var) 3244 (setq qk (pop q))) 3245 (push '(mlist) acc-ord) 3246 (push '(mlist) acc-pt) 3247 (push '(mlist) acc-var) 3248 (setq q (taylor-info q)) 3249 (if (null q) (list acc-var acc-pt acc-ord) (append q (list acc-var acc-pt acc-ord)))) 3250 3251 (t 3252 (setq acc (if (and (fourth qk) (consp (fourth qk)) (eq '$asympt (caar (fourth qk)))) 3253 (list '$asympt) nil)) 3254 (push (taylor-trunc qk) acc) 3255 (push (exp-pt qk) acc) 3256 (push (datum-var qk) acc) 3257 (push '(mlist) acc) 3258 (setq q (taylor-info q)) 3259 (if (null q) (list acc) (append q (list acc))))))))) 3260 3261(defmfun $taylorinfo (x) 3262 (if (and (consp x) (member 'trunc (first x) :test #'eq)) 3263 (cons '(mlist) (taylor-info (mrat-tlist x))) 3264 nil)) 3265 3266 3267;;; Local Modes: 3268;;; Lisp let-pw Indent:2 3269;;; End: 3270