1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; **************************************************************
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 ;;; **************************************************************
17 ;;; TOP LEVEL STRUCTURE
19 ;;; Power series have the following format when seen outside the power
22 ;;; ((MRAT SIMP <varlist> <genvar> <tlist> trunc) <poly-form>)
24 ;;; This is the form of the output of the expressions, to
25 ;;; be displayed they are RATDISREPed and passed to DISPLA.
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.
31 ;;; (PS (<var> . <ord-num>) (<trunc-lvl>)
32 ;;; (<exponent> . <coeff>) (<exponent> . <coeff>) . . .)
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.
44 ;;; Problem: fix expansion of logs so that taylor(log(1+exp(-1/x)),x,0,3)
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.
50 ;;; Problem: Need to fix things so that asymptotic kernels aren't put onto
51 ;;; tvars via tlist merge etc. in taylor1. Done.
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.
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.
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
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.
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).
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.
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).
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.
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.
103 ;;; Tvar-limits should be stored locally so that psdisrep need not recompute
104 ;;; each gvar limit when disrepping.
106 (macsyma-module hayat
)
110 (defvar *within-srf?
* nil
)
112 (load-macsyma-macros mhayat rzmac ratmac
)
114 ;;; Subtitle Special Stuff for Compiling
118 *a
* ;Temporary special
119 tlist
;An association list which contains the
120 ;relevant information for the expansion which
121 ;is passed in at toplevel invocation.
122 log-1
;What log(-1) should be log(-1) or pi*i.
123 log%i
;Similarly for log(i)
124 exact-poly
;Inicates whether polynomials are to be
125 ;considered exact or not. True within SRF,
126 ;false within TAYLOR.
129 tay-const-expand
;For rediculousness like csch(log(x))
132 last-exp
;last-expression through taylor2
133 ivars
;Pairlist if gensym and disreped version
134 key-vars
;Pairlist of gensym and key var (for searching
140 *within-srf?
* ;flag for in srf
142 least_term?
; If non-null then the addition routines
143 ; are adding or subtracting coeff's of the
144 ; least term of a sum so they should do
145 ; zero checking on it if it is desired.
146 taylor_simplifier
; This is set by taylor1 to the object
147 ; which will be funcalled whenever
148 ; coefficient simplification is desired.
149 zerolist
; A list of constant expressions which have
150 ; been verified to be zero by a call to
151 ; $TAYLOR_SIMPLIFIER in taylor2. It is used to
152 ; suppress the message that TAYLOR is assuming
153 ; an expression to be zero.
154 ; 0p-funord lexp-non0 ; referenced only in commented-out code, so comment out here too
156 ) ;Don't want to see closed compilation notes.
158 (defmvar $psexpand
()
159 "When TRUE extended rational function expressions will be displayed fully
160 expanded. (RATEXPAND will also cause this.) If FALSE, multivariate
161 expressions will be displayed just as in the rational function package.
162 If PSEXPAND:MULTI, then terms with the same total degree in the variables
163 are grouped together.")
165 (defmvar $maxtayorder t
166 "When true TAYLOR retains as many terms as are certain to be correct
167 during power series arithmetic. Otherwise, truncation is controlled
168 by the arguments specified to TAYLOR.")
170 (defmvar $taylor_truncate_polynomials t
171 "When FALSE polynomials input to TAYLOR are considered to have infinite
172 precision; otherwise (the default) they are truncated based upon the input
175 (defmvar $taylor_logexpand t
176 "Unless FALSE log's of products will be expanded fully in TAYLOR (the default)
177 to avoid identically-zero constant terms which involve log's. When FALSE,
178 only expansions necessary to produce a formal series will be executed.")
180 (defmvar $taylor_simplifier
'simplify
181 "A function of one argument which TAYLOR uses to simplify coefficients
184 (defvar taylor_simplifier nil
)
186 ;;; Subtitle General Macsyma Free Predicates
189 (cond ((equal e x
) () )
192 (null (member x
(cdr ($listofvars e
)) :test
#'equal
)))
193 ('t
(do ((l (cdr e
) (cdr l
))) ((null l
) 't
)
194 (or (zfree (car l
) x
) (return () ))))))
196 (defun mfree (exp varl
)
197 (cond ((atom exp
) (not (member exp varl
:test
#'eq
)))
198 ((eq (caar exp
) 'mrat
)
199 (do ((l (mrat-varlist exp
) (cdr l
)))
201 (unless (mfree (car l
) varl
) (return () ))))
202 ((or (member (caar exp
) dummy-variable-operators
:test
#'eq
)
203 (member 'array
(cdar exp
) :test
#'eq
))
204 (do ((vars varl
(cdr vars
)))
206 (unless (freeof (car vars
) exp
) (return () ))))
207 ('t
(and (mfree (caar exp
) varl
) (mfreel (cdr exp
) varl
)))))
209 (defun mfreel (l varl
)
210 (or (null l
) (and (mfree (car l
) varl
) (mfreel (cdr l
) varl
))))
212 ;;; Subtitle Coefficient Arithmetic
215 (cond ((equal x
(rcone)) (rcone))
216 ((rczerop y
) (rcone))
217 ((and (equal (cdr y
) 1) (fixnump (car y
)))
219 ((and $radexpand
(numberp (car y
)) (numberp (cdr y
)))
221 (setq y
(maxima-rationalize (quot (car y
) (cdr y
)))))
222 (ratexpt (rcquo (rcexpt1 (car x
) (cdr y
))
223 (rcexpt1 (cdr x
) (cdr y
)))
226 (prep1 (m^
(rcdisrep x
) (rcdisrep y
)))))))
229 (cond ((equal p
1) (rcone))
230 ((pcoefp p
) (prep1 (m^
(pdis p
) (*red
1 n
))))
231 ;; psfr does a square-free decom on p yielding (p1 e1 p2 e2 ... pn en)
232 ;; where p = p1^e1 p2^e2 ... pn^en, the pi being square-free
233 (t (do ((l (psqfr p
) (cddr l
))
236 (if (not (equal (rem (cadr l
) n
) 0))
237 (setq ans
(rctimes ans
(prep1 (m^
(pdis (car l
))
238 (*red
(cadr l
) n
)))))
239 ;; If pi<0, n=2m and n|ei then ei=2e and
240 ;; (pi^ei)^(1/(2m)) = (-pi)^(e/m)
242 (when (and (evenp n
) (eq ($sign
(pdis (car l
))) '$neg
))
243 (rplaca l
(pminus (car l
))))
244 (setq ans
(rctimes ans
(ratexpt (cons (car l
) 1)
245 (truncate (cadr l
) n
))))))))))
247 (defun rccoefp (e) ;a sure check, but expensive
250 (member (caar e
) genvar
:test
#'eq
))
252 (member (cadr e
) genvar
:test
#'eq
))))
254 ;;; Subtitle Exponent arithmetic
257 (and (not (infp x
)) (signp e
(car x
))))
260 (cond ((or (infp x
) (infp y
)) (inf))
261 ((and (equal (cdr x
) 1) (equal (cdr y
) 1))
262 (cons (+ (car x
) (car y
)) 1))
263 (t (ereduce (+ (* (car x
) (cdr y
)) (* (cdr x
) (car y
)))
264 (* (cdr x
) (cdr y
))))))
267 (cond ((infp x
) (inf))
268 ((and (equal (cdr x
) 1) (equal (cdr y
) 1))
269 (cons (- (car x
) (car y
)) 1))
270 (t (ereduce (- (* (car x
) (cdr y
)) (* (cdr x
) (car y
)))
271 (* (cdr x
) (cdr y
))))))
276 ((equal (cdr x
) (cdr y
)) (cons (min (car x
) (car y
)) (cdr x
)))
277 ((< (* (car x
) (cdr y
)) (* (cdr x
) (car y
))) x
)
281 (cond ((or (infp x
) (infp y
)) (inf))
282 ((equal (cdr x
) (cdr y
)) (cons (max (car x
) (car y
)) (cdr x
)))
283 ((> (* (car x
) (cdr y
)) (* (cdr x
) (car y
))) x
)
287 (cond ((or (infp x
) (infp y
)) (inf))
288 ((and (equal (cdr x
) 1) (equal (cdr y
) 1))
289 (cons (* (car x
) (car y
)) 1))
290 (t (ereduce (* (car x
) (car y
)) (* (cdr x
) (cdr y
))))))
294 (cons (- (cdr e
)) (- (car e
)))
295 (cons (cdr e
) (car e
))))
298 (cond ((infp x
) (inf))
300 (t (ereduce (* (car x
) (cdr y
))
301 (* (cdr x
) (car y
))))))
304 (cond ((infp x
) (inf))
305 ((= (cdr x
) 1) (cons (1+ (car x
)) 1))
306 (t (cons (+ (cdr x
) (car x
)) (cdr x
)))))
309 (cond ((infp x
) (inf))
310 ((equal (cdr x
) 1) (cons (1- (car x
)) 1))
311 (t (cons (- (car x
) (cdr x
)) (cdr x
)))))
316 ((equal (cdr x
) (cdr y
)) (> (car x
) (car y
)))
317 (t (> (* (car x
) (cdr y
)) (* (car y
) (cdr x
))))))
321 ((or (null e1
) (null e2
)) ())
322 (t (and (equal (car e1
) (car e2
))
323 (equal (cdr e1
) (cdr e2
))))))
326 (if (signp l d
) (setq d
(- d
) n
(- n
)))
327 (if (zerop n
) (rczero)
328 (let ((gcd (gcd n d
)))
329 (cons (/ n gcd
) (/ d gcd
)))))
332 (let ((xn (abs (car x
))) (xd (cdr x
))
333 (yn (abs (car y
))) (yd (cdr y
)))
334 (cons (gcd xn yn
) (* xd
(/ yd
(gcd xd yd
))))))
336 ;;; Subtitle polynomial arithmetic
338 (declare-top (special vars
))
340 (defun ord-vector (p)
341 (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum
))) tlist
)))
342 (declare (special vars
))
343 (cond ((not (cdr vars
)) (ncons (ps-le* p
)))
344 (t (ord-vect1 p
) (mapcar #'(lambda (x) (or (cdr x
) (rczero))) vars
)))))
347 (declare (special vars
))
349 (let ((data (assoc (gvar p
) vars
:test
#'eq
))
351 (rplacd data
(cond ((not (cdr data
)) le
)
352 (t (emin (cdr data
) le
))))
353 (mapl #'(lambda (l) (ord-vect1 (lc l
))) (terms p
)))))
355 (defun trunc-vector (p min?
)
356 (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum
))) tlist
)))
357 (declare (special vars
))
358 (if (null (cdr vars
)) (ncons (if (psp p
) (trunc-lvl p
) () ))
361 (mapcar 'cdr vars
)))))
363 (defun trunc-vect1 (p min?
)
364 (declare (special vars
))
366 (let ((data (assoc (gvar p
) vars
:test
#'eq
))
367 (trunc (trunc-lvl p
)))
369 (rplacd data
(if (null (cdr data
)) trunc
370 (if min?
(emin (cdr data
) trunc
)
371 (emax (cdr data
) trunc
))))))
372 (dolist (term (terms p
))
373 (trunc-vect1 (c term
) min?
))))
375 (declare-top (unspecial vars
))
379 (cond ((pscoefp y
) (rcplus x y
))
382 ((pscoefp y
) (if (rczerop y
) x
(pscplus y x
)))
383 ((eqgvar (gvar-o x
) (gvar-o y
)) (psplus1 x y
))
384 ((pointerp (gvar-o x
) (gvar-o y
)) (pscplus y x
))
388 (if (not (and least_term? taylor_simplifier
)) (rcplus x y
)
389 (prep1 (funcall taylor_simplifier
(m+ (rcdisrep x
) (rcdisrep y
))))))
392 (cond ((pscoefp x
) (cond ((pscoefp y
) (rcdiff x y
))
393 ((rczerop x
) (pstimes (rcmone) y
))
394 (t (pscdiff x y
() ))))
395 ((pscoefp y
) (if (rczerop y
) x
(pscdiff y x t
)))
396 ((eqgvar (gvar-o x
) (gvar-o y
)) (psdiff1 x y
))
397 ((pointerp (gvar-o x
) (gvar-o y
)) (pscdiff y x t
))
398 (t (pscdiff x y
() ))))
401 (if (not (and least_term? taylor_simplifier
)) (rcdiff x y
)
402 (prep1 (funcall taylor_simplifier
(m- (rcdisrep x
) (rcdisrep y
))))))
405 (let ((ans (cons () () )))
406 (psplus2 (gvar-o x
) (emin (trunc-lvl x
) (trunc-lvl y
))
407 (cons 0 (terms x
)) (cons 0 (terms y
)) ans ans
)))
410 (if (e> (rczero) (trunc-lvl p
)) p
411 (pscheck (gvar-o p
) (poly-data p
) (pscplus1 c
(terms p
)))))
413 (defun pscdiff (c p fl
)
414 (if (e> (rczero) (trunc-lvl p
))
415 (if fl p
(psminus p
))
416 (pscheck (gvar-o p
) (poly-data p
)
417 (cond ((not fl
) (pscplus1 c
(psminus-terms (terms p
))))
418 (t (pscplus1 (psminus c
) (terms p
)))))))
420 (defun strip-zeroes (terms ps?
)
421 (cond ((or (null terms
) (null taylor_simplifier
)) terms
)
423 (do ((terms terms
(n-term terms
)))
425 (change-lc terms
(strip-zeroes (lc terms
) 't
))
426 (unless (rczerop (lc terms
)) (return terms
))))
428 (if (null taylor_simplifier
) terms
429 (let ((exp (rcdisrep terms
)))
430 ;; If a pscoeff is not free of tvars then the ps is a
431 ;; multivar series and we can't handle a recursive
432 ;; call to taylor (as opposed to a call to prep1, as below)
433 ;; because this would be circuler (e.g. try
434 ;; taylor(x/ (x^2+1),[x],%i,-1) ). Besides, in this case
435 ;; the pscoeff contains a tvar hence should not be 0.
436 (if (not (mfree exp tvars
)) terms
437 (prep1 (funcall taylor_simplifier exp
))))))
438 (t (pscheck (gvar-o terms
) (poly-data terms
)
439 (strip-zeroes (terms terms
) () )))))
441 (defun pscplus1 (c l
)
442 (cond ((null l
) (list (term (rczero) c
)))
443 ((rczerop (le l
)) (setq c
(psplus c
(lc l
)))
444 (if (rczerop c
) (strip-zeroes (n-term l
) () )
445 (cons (term (rczero) c
) (n-term l
))))
446 ((e> (le l
) (rczero)) (cons (term (rczero) c
) l
))
447 (t (cons (lt l
) (let ((least_term?
)) (pscplus1 c
(n-term l
)))))))
449 ;;; Both here and in psdiff2 xx and yy point one before where one
450 ;;; might think they should point so that extensions will be retained.
452 (defun psplus2 (varh trunc xx yy ans a
)
454 a
(cond ((mono-term? xx
)
455 (if (mono-term? yy
) (go end
) (go null
)))
456 ((mono-term? yy
) (setq yy xx
) (go null
)))
457 (cond ((equal (le (n-term xx
)) (le (n-term yy
)))
458 (setq xx
(n-term xx
) yy
(n-term yy
))
459 (setq c
(let ((least_term?
(null (n-term ans
))))
460 (psplus (lc xx
) (lc yy
))))
461 (if (rczerop c
) (go a
) (add-term a
(le xx
) c
)))
462 ((e> (le (n-term xx
)) (le (n-term yy
)))
463 (setq yy
(n-term yy
))
464 (add-term a
(lt yy
)))
465 (t (setq xx
(n-term xx
))
466 (add-term a
(lt xx
))))
469 null
(if (or (mono-term? yy
) (e> (le (n-term yy
)) trunc
))
472 (setq yy
(n-term yy
))
473 (add-term-&-pop a
(lt yy
))
475 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
478 (let ((ans (cons () () )))
479 (psdiff2 (gvar-o x
) (emin (trunc-lvl x
) (trunc-lvl y
))
480 (cons 0 (terms x
)) (cons 0 (terms y
)) ans ans
)))
482 (defun psdiff2 (varh trunc xx yy ans a
)
484 a
(cond ((mono-term? xx
)
489 (cons 0 (mapcar #'(lambda (q)
490 (term (e q
) (psminus (c q
))))
494 (setq yy xx
) (go null
)))
495 (cond ((equal (le (n-term xx
)) (le (n-term yy
)))
496 (setq xx
(n-term xx
) yy
(n-term yy
))
497 (setq c
(let ((least_term?
(null (n-term ans
))))
498 (psdiff (lc xx
) (lc yy
))))
499 (if (rczerop c
) (go a
)
500 (add-term a
(le xx
) c
)))
501 ((e> (le (n-term xx
)) (le (n-term yy
)))
502 (setq yy
(n-term yy
))
503 (add-term a
(le yy
) (psminus (lc yy
))))
504 (t (setq xx
(n-term xx
))
505 (add-term a
(lt xx
))))
508 null
(if (or (mono-term? yy
) (e> (le (n-term yy
)) trunc
))
511 (setq yy
(n-term yy
))
512 (add-term-&-pop a
(le yy
) (lc yy
))
514 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
517 (if (psp x
) (make-ps x
(psminus-terms (terms x
)))
520 (defun psminus-terms (terms)
521 (let ((ans (cons () () )))
522 (do ((q terms
(n-term q
))
525 (add-term a
(le q
) (psminus (lc q
))))))
527 (defun pscheck (a b terms
)
528 (cond ((null terms
) (rczero))
529 ((and (mono-term? terms
) (rczerop (le terms
)))
531 (t (make-ps a b terms
))))
533 (defun pstrim-terms (terms e
)
535 (cond ((null terms
) (return () ))
536 ((null (e> e
(le terms
))) (return terms
))
537 (t (setq terms
(n-term terms
))))))
539 (defun psterm (terms e
)
540 (psterm1 (pstrim-terms terms e
) e
))
543 (cond ((null l
) (rczero))
544 ((e= (le l
) e
) (lc l
))
547 (defun pscoeff1 (a b c
) ;a is an mrat!!!
548 (let ((tlist (mrat-tlist a
)))
549 (cons (nconc (list 'mrat
'simp
(mrat-varlist a
) (mrat-genvar a
))
550 (do ((l (mrat-tlist a
) (cdr l
))
551 (ans () (cons (car l
) ans
)))
553 (when (alike1 (caar l
) b
)
555 (and (or ans
(cdr l
))
556 (list (nreconc ans
(cdr l
)) 'trunc
))))))
557 (pscoef (mrat-ps a
) (int-gvar (get-datum b
)) (prep1 c
)))))
559 (defun pscoef (a b c
)
560 (cond ((pscoefp a
) (if (rczerop c
) a
(rczero)))
561 ((eq b
(gvar a
)) (psterm (terms a
) c
))
562 (t (do ((gvar-o (gvar-o a
))
563 (poly-data (poly-data a
))
565 (terms (terms a
) (n-term terms
))
568 (unless (rczerop (setq temp
(pscoef (lc terms
) b c
)))
569 (setq ans
(psplus ans
570 (make-ps gvar-o poly-data
571 (ncons (term (le terms
)
574 (defun psdisextend (p)
575 (cond ((not (psp p
)) p
)
576 (t (make-ps p
(mapcar #'(lambda (q) (cons (car q
) (psdisextend (cdr q
))))
580 (if (psp p
) (psfloat1 p
(trunc-lvl p
) (terms p
) (ncons 0))
581 (rctimes (rcfone) p
)))
583 (defun psfloat1 (p trunc l ans
)
585 (a (last ans
) (n-term a
)))
586 ((or (null l
) (e> (le l
) trunc
))
587 (pscheck (gvar-o p
) (poly-data p
) (cdr ans
)))
588 (add-term a
(le l
) (psfloat (lc l
)))
589 (setq l
(n-term l
))))
592 (pstrunc1 p
(mapcar #'(lambda (q) (cons (int-gvar q
) (current-trunc q
)))
595 (defun pstrunc1 (p trlist
)
599 (let ((trnc (cdr (assoc (gvar p
) trlist
:test
#'eq
))) (trunc-ps) (a nil
))
600 (do ((l (terms p
) (n-term l
)))
601 ((null l
) (pscheck (gvar-o p
) (ncons (trunc-lvl p
)) (nreverse a
)))
602 (when (e> (le l
) trnc
)
603 (return (pscheck (gvar-o p
) (ncons trnc
) (nreverse a
))))
604 (unless (rczerop (setq trunc-ps
(pstrunc1 (lc l
) trlist
)))
605 (push (term (le l
) trunc-ps
) a
)))))))
608 (cond ((or (rczerop x
) (rczerop y
)) (rczero))
609 ((pscoefp x
) (cond ((pscoefp y
) (rctimes x y
))
610 ((equal x
(rcone)) y
)
611 (t (psctimes* x y
))))
612 ((pscoefp y
) (if (equal y
(rcone)) x
(psctimes* y x
)))
613 ((eqgvar (gvar-o x
) (gvar-o y
)) (pstimes*1 x y
))
614 ((pointerp (gvar-o x
) (gvar-o y
)) (psctimes* y x
))
615 (t (psctimes* x y
))))
617 (defun psctimes* (c p
)
618 (make-ps p
(maplist #'(lambda (l)
619 (term (le l
) (pstimes c
(lc l
))))
622 (defun pstimes*1 (xa ya
)
623 (let ((ans (cons () () ))
624 (trunc (let ((lex (ps-le xa
)) (ley (ps-le ya
)))
625 (e+ (emin (e- (trunc-lvl xa
) lex
) (e- (trunc-lvl ya
) ley
))
628 (setq trunc
(emin trunc
(t-o-var (gvar xa
)))))
629 (pstimes*2 xa ya trunc ans
)))
631 (defun pstimes*2 (xa ya trunc ans
)
633 (setq x
(terms xa
) y
(setq yy
(terms ya
)) a ans
)
634 a
(cond ((or (null y
) (e> (setq e
(e+ (le x
) (le y
))) trunc
))
636 ((not (rczerop (setq c
(pstimes (lc x
) (lc y
)))))
637 (add-term-&-pop a e c
)))
640 b
(unless (setq x
(n-term x
))
641 (return (pscheck (gvar-o xa
) (list trunc
) (cdr ans
))))
643 c
(when (or (null y
) (e> (setq e
(e+ (le x
) (le y
))) trunc
))
645 (setq c
(pstimes (lc x
) (lc y
)))
646 d
(cond ((or (mono-term? a
) (e> (le (n-term a
)) e
))
647 (add-term-&-pop a e c
))
648 ((e> e
(le (n-term a
)))
651 (t (setq c
(psplus c
(lc (n-term a
))))
653 (rplacd a
(n-term (n-term a
)))
655 (change-lc (n-term a
) c
)
656 (setq a
(n-term a
))))))
660 (defun pscsubst (c v p
)
661 (cond ((pscoefp p
) p
)
662 ((eq v
(gvar p
)) (pscsubst1 c p
))
663 ((pointerp v
(gvar p
)) p
)
664 (t (make-ps p
(maplist
665 #'(lambda (q) (term (le q
)
666 (pscsubst c v
(lc q
))))
669 (defun pscsubst1 (v u
)
671 (ul (terms u
) (n-term ul
)))
673 (setq a
(psplus a
(pstimes (lc ul
) (psexpt v
(le ul
)))))))
675 (defun get-series (func trunc var e c
)
676 (let ((pw (e// trunc e
)))
677 (setq e
(if (and (equal e
(rcone)) (equal c
(rcone)))
678 (getexp-fun func var pw
)
679 (psmonsubst (getexp-fun func var pw
) trunc e c
)))
680 (if (and $float $keepfloat
) (psfloat e
) e
)))
682 (defun psmonsubst (p trunc e c
)
684 (psmonsubst1 p trunc e c
685 `(() .
,(terms p
)) (ncons () ) (rcone) (rczero))
689 (defun psmonsubst1 (p trunc e c l ans cc el
)
690 ;; We set $MAXTAYORDER to () here so that the calls to psexpt below
691 ;; won't do needless extra work, e.g. see rwg's complaint of 9/7/82.
692 (prog (a ee varh $maxtayorder
)
693 (setq a ans varh
(gvar-o p
))
694 a
(cond ((or (mono-term? l
)
695 (e> (setq ee
(e* e
(le (n-term l
)))) trunc
))
699 (psexpt c
(e- (le (setq l
(n-term l
)))
702 (add-term a ee
(pstimes cc
(lc l
)))))
703 (setq a
(n-term a
) el
(le l
))
705 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
707 (defun psexpon-gcd (terms)
708 (do ((gcd (le terms
) (egcd (le l
) gcd
))
709 (l (n-term terms
) (n-term l
)))
713 (if (psp p
) (psfind-s (psterm (terms p
) (rczero)))
717 (cond ((null (atom (cdr r
))) (rczero))
719 (t (do ((p (ptterm (cdar r
) 0) (ptterm (cdr p
) 0)))
720 ((atom p
) (cons p
(cdr r
)))))))
723 (cond ((rczerop n
) ;; p^0
724 (if (rczerop p
) ;; 0^0
725 (merror (intl:gettext
"taylor: 0^0 is undefined."))
726 (rcone))) ;; Otherwise can let p^0 = 1
727 ((or (equal n
(rcone)) (equal n
(rcfone))) p
) ;; p^1 cases
728 ((pscoefp p
) (rcexpt p n
))
729 ((mono-term?
(terms p
)) ;; A monomial to a power
730 (let ((s (psfind-s n
)) (n-s) (x) (l (terms p
)))
731 ;; s is the numeric part of the exponent
732 (if (floatp (car s
)) ;; Perhaps we shouldn't
733 ;; rationalize if $keepfloat is true?
734 (setq s
(maxima-rationalize (quot (car s
) (cdr s
)))))
735 (setq n-s
(psdiff n s
) ;; the non-numeric part of exponent
736 x
(e* s
(le l
))) ;; the degree of the lowest term
737 (setq x
(if (and (null $maxtayorder
) ;; if not getting all terms
738 (e> x
(t-o-var (gvar p
))))
739 ;; and result is of high order
740 (rczero) ;; then zero is enough
741 (pscheck (gvar-o p
) ;; otherwise
742 (ncons (e+ (trunc-lvl p
) ;; new trunc-level
743 (e- x
(le l
)))) ;; kick exponent
744 (ncons (term x
(psexpt (lc l
) n
))))))
746 (if (or (rczerop n-s
) (rczerop x
)) ;; is that good enough?
747 x
;; yes! The rest is bletcherous.
748 (pstimes x
(psexpt (prep1 (m^
(get-inverse (gvar p
))
751 (t (prog (l lc le inc trunc s lt mr lim lcinv ans
)
752 (setq lc
(lc (setq l
(terms p
)))
753 le
(le l
) lt
(lt l
) trunc
(trunc-lvl p
)
754 inc
(psexpon-gcd l
) s
(psfind-s n
))
755 (when (floatp (car s
))
756 (setq s
(maxima-rationalize (quot (car s
) (cdr s
)))))
757 (setq ans
(psexpt (setq lt
(pscheck (gvar-o p
) (list trunc
)
759 lcinv
(psexpt lc
(rcmone))
760 mr
(e+ inc
(e* s le
))
761 lim
(if (and (infp trunc
) (not (e> s
(rczero))))
763 ;; See the comment in PSEXPT1 below which tells
764 ;; why we don't allow inf. trunc's here.
765 (e+ (if (and (infp trunc
) (not (rcintegerp s
)))
766 (if (infp (setq lim
(t-o-var (gvar p
))))
772 (if (or (pscoefp ans
) (null (eq (gvar p
) (gvar ans
))))
773 (list 0 (term (rczero) ans
))
774 (cons 0 (terms ans
))))
775 (and (null $maxtayorder
)
778 (e> (e* s
(le (last l
))) (t-o-var (gvar p
))))
779 (setq lim
(emin lim
(t-o-var (gvar p
)))))
780 ;;(and (infp lim) (n-term l) (e> (rczero) n)
782 (return (psexpt1 (gvar-o p
)
783 lim l n s inc
1 mr ans le lcinv
))))))
785 (defun psexpt1 (varh trunc l n s inc m mr ans r linv
)
786 ;; n is the power we are raising the series to
787 ;; inc is the exponent increment
788 ;; mr is the current exponent
789 ;; tr is the truncation level desired
790 ;; m is the term index
792 (prog (a (k 0) ak cm-k c ma0 sum kr tr
)
794 ;; truly unfortunate that we need so many variables in this hack
795 (setq a
(last ans
) tr trunc
)
796 ;; I don't see what's wrong with truncating exact series when
797 ;; raising them to fractional powers so we'll allow it for now.
798 ;; This is accomplished above in PSEXPT (see the comment). Thus,
799 ;; presumably, this check should never be needed anymore.
800 ;; Bugs catching this clause were sqrt(1-x)*taylor(f1,x,0,0)
801 ;; and sqrt(taylor(x+x^2,x,0,2)),taylor_truncate_polynomials=false.
804 (setq tr
(e* s
(le (last l
))))
805 (merror (intl:gettext
"taylor: expected an integer, instead found: ~:M") s
)))
806 (when (infp tr
) (setq tr
(t-o-var (car varh
))))
807 b
(and (e> mr tr
) (go end
))
808 (setq kr inc ak l ma0
(pstimes (cons 1 m
) linv
)
810 a
(if (or (> k m
) (null (setq cm-k
(psterm (cdr ans
) (e- mr kr
)))))
812 (setq ak
(or (pstrim-terms ak
(e+ kr r
)) (go add-term
))
813 c
(pstimes (psdiff (pstimes (cons k
1) n
)
815 (pstimes (if (e= (e+ kr r
) (le ak
))
819 (setq sum
(psplus sum c
)
820 k
(1+ k
) kr
(e+ kr inc
))
823 (and (null (rczerop sum
))
824 (add-term-&-pop a mr
(pstimes ma0 sum
)))
825 (setq m
(1+ m
) mr
(e+ mr inc
))
827 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
829 (defun psderivative (p v
)
830 (cond ((pscoefp p
) (rcderiv p v
))
832 (if (prog1 (rczerop (ps-le p
))
833 (setq p
(psderiv1 (gvar-o p
)
834 (trunc-lvl p
) (cons 0 (terms p
)) (list 0))))
835 (strip-zeroes p
't
) p
))
836 (t (psderiv2 (gvar-o p
)
837 (trunc-lvl p
) v
(cons 0 (terms p
)) (list 0)))))
839 (defun psderiv1 (varh trunc l ans
)
841 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
842 (pscheck varh
(list (e1- trunc
)) (cdr ans
)))
844 (when (not (rczerop (le l
)))
845 (add-term-&-pop a
(e1- (le l
)) (pstimes (le l
) (lc l
))))))
847 (defun psderiv2 (varh trunc v l ans
)
848 (do ((a (last ans
) (n-term a
)) (c))
849 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
850 (pscheck varh
(list trunc
) (cdr ans
)))
852 (or (rczerop (setq c
(psderivative (lc l
) v
)))
853 (add-term a
(le l
) c
))))
857 (cond ((pscoefp p
) (rcderivx p
))
858 ((or (rczerop (setq temp
(getdiff (gvar-o p
))))
859 (eq (car temp
) 'multi
))
860 (setq temp2
(psdp2 (gvar-o p
) (trunc-lvl p
)
861 (cons 0 (terms p
)) (list 0)))
862 (if (eq (car temp
) 'multi
)
864 (make-ps (gvar-o p
) (ncons (inf))
865 (list (term (cdr temp
) (rcone)))))
868 (trunc-lvl p
) (cons 0 (terms p
))
871 (defun psdp1 (varh trunc l ans dx
)
872 (do ((a (last ans
)) (c (rczero)))
873 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
874 (psplus c
(pscheck varh
(list (e1- trunc
)) (cdr ans
))))
876 (if (rczerop (le l
)) (setq c
(psdp (lc l
)))
878 a
(e1- (le l
)) (pstimes (le l
) (pstimes dx
(lc l
)))))))
880 (defun psdp2 (varh trunc l ans
)
881 (do ((a (last ans
)) (c))
882 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
883 (pscheck varh
(list trunc
) (cdr ans
)))
885 (when (null (rczerop (setq c
(psdp (lc l
)))))
886 (add-term-&-pop a
(le l
) c
))))
890 ;;; (defun psintegrate (p v)
891 ;;; (cond ((rczerop p) (rczero))
893 ;;; (pstimes p (taylor2 (get-inverse (car v)))))
894 ;;; ((eqgvar v (gvar-o p))
895 ;;; (psinteg1 (gvar-o p)
896 ;;; (trunc-lvl p) (cons 0 (terms p)) (list 0)))
897 ;;; (t (psinteg2 (gvar-o p)
898 ;;; (trunc-lvl p) v (cons 0 (terms p)) (list 0)))))
900 ;;; (defun psinteg1 (varh trunc l ans)
902 ;;; (setq a (last ans))
903 ;;; a (if (or (null (n-term l)) (e> (le (n-term l)) trunc))
905 ;;; (add-term a (e1+ (le (setq l (n-term l))))
907 ;;; (if (e= (le l) (rcmone))
908 ;;; (prep1 (list '(%LOG)
912 ;;; (setq a (n-term a)))
914 ;;; end (return (pscheck varh (list (e1+ trunc)) (cdr ans)))))
916 ;;; (defun psinteg2 (varh trunc v l ans)
918 ;;; (setq a (last ans))
919 ;;; a (if (or (null (n-term l)) (e> (le (n-term l)) trunc))
921 ;;; (add-term a (le l)
922 ;;; (psintegrate (lc (setq l (n-term l))) v))
923 ;;; (setq a (n-term a)))
925 ;;; end (return (pscheck varh (list trunc) (cdr ans)))))
927 (defun psexpt-log-ord (p)
928 (cond ((null $maxtayorder
) (emin (trunc-lvl p
) (t-o-var (gvar p
))))
929 ((infp (trunc-lvl p
)) (t-o-var (gvar p
)))
938 (cond ((pscoefp p
) (psexpt-fn2 (rcdisrep p
)))
939 ((ps-lim-infp p
) (psexpt-fn-sing p
))
940 ((prog2 (setq ord
<0?
(e> (rczero) (ps-le p
)))
941 (null (n-term (terms p
))))
942 (setq ans
(get-series '%exp
(psexpt-log-ord p
) (gvar-o p
)
943 (if ord
<0?
(e- (ps-le p
)) (ps-le p
))
945 (if ord
<0?
(ps-invert-var ans
) ans
))
947 (when (e= (rczero) (e (setq ans
(ps-gt p
))))
948 (pstimes (psexpt-fn (pscheck (gvar-o p
) (list (trunc-lvl p
))
949 (delete ans
(terms p
) :test
#'eq
)))
950 (psexpt-fn2 (srdis (c ans
)))))
951 (when (e= (rczero) (ps-le p
))
952 (pstimes (psexpt-fn2 (srdis (lc (terms p
))))
953 (psexpt-fn (pscheck (gvar-o p
) (list (trunc-lvl p
))
954 (n-term (terms p
))))))) )
955 (t (prog (l inc trunc ea0 ans
)
958 ;(return (ps-invert-var (psexpt-fn (ps-invert-var p))))
959 (setq l
(invert-terms l
)))
960 (setq trunc
(trunc-lvl p
)
961 inc
(psexpon-gcd l
) ea0
(rcone))
962 (unless (e> (le l
) (rczero))
963 ;; MEANING OF FOLLOWING MESSAGE IS OBSCURE
964 (merror "PSEXPT-FN: unreachable point."))
966 (if (or (pscoefp ea0
) (null (eq (gvar p
) (gvar ea0
))))
967 (list 0 (term (rczero) ea0
))
968 (cons 0 (terms ea0
))))
970 (setq trunc
(emin trunc
(t-o-var (gvar p
)))))
971 (when (infp trunc
) (setq trunc
(t-o-var (gvar p
))))
972 (setq ans
(psexpt-fn1 (gvar-o p
) trunc l inc
1 inc ans
))
973 (return (if ord
<0?
(ps-invert-var ans
) ans
)))))))
975 (defun psexpt-fn-sing (p)
976 (let ((inf-var?
(member (gvar-lim (gvar p
)) '($inf $minf
) :test
#'eq
))
977 (c*logs
(c*logs
(lt-poly p
))) c strongest-term
)
978 ;; Must pull out out logs here: exp(ci*log(ui)+x) -> ui^ci*exp(x)
979 ;; since its much harder for adjoin-tvar to do this transformation
980 ;; below after things have been disrepped.
981 (setq c
(exp-c*logs c
*logs
) p
(psdiff p
(sum-c*logs c
*logs
)))
982 (if (not (ps-lim-infp p
))
983 ;; Here we just subtracted the only infinite term, e.g.
984 ;; p = 1/2*log(x)+1/log(x)+...
985 (pstimes c
(psexpt-fn p
))
987 (setq strongest-term
(if inf-var?
(ps-gt p
) (ps-lt p
)))
988 ;; If the strongest term has degree 0 in the mainvar then the singular
989 ;; terms occur in some other weaker var. There may be terms in this
990 ;; coef which are not singular (e.g. 1 in (1/x+1+...)+exp(-1/x)+...) so
991 ;; we must recursively psexpt-fn this term to get only what we need.
992 (if (rczerop (e strongest-term
))
993 (setq c
(pstimes c
(psexpt-fn (c strongest-term
))))
994 (dolist (exp (expand-and-disrep strongest-term p
))
995 (setq c
(pstimes c
(adjoin-tvar (m^
'$%e exp
))))))
996 (pstimes c
(psexpt-fn (pscheck (gvar-o p
) (list (trunc-lvl p
))
998 (delete strongest-term
(terms p
) :test
#'eq
)
999 (n-term (terms p
))))))))))
1001 (defun gvar-logp (gvar)
1002 (let ((var (get-inverse gvar
)))
1003 (and (consp var
) (eq (caar var
) 'mexpt
) (equal (caddr var
) -
1)
1004 (consp (setq var
(cadr var
))) (eq (caar var
) '%log
)
1009 (let ((log (gvar-logp (gvar p
))) c
)
1013 (setq c
(psconst (psterm (terms p
) (rcmone))))
1014 ;; We don't want e.g. exp(x^a*log(x)) -> x^x^a
1015 (if (not (mfree (rcdisrep c
) tvars
)) ()
1016 (cons (cons c
(cons log p
))
1017 (c*logs
(psterm (terms p
) (rczero))))))))))
1020 (if (pscoefp p
) p
(psconst (psterm (terms p
) (rczero)))))
1022 (defun exp-c*logs
(c*logs
)
1023 (if (null c
*logs
) (rcone)
1024 (pstimes (taylor2 `((mexpt) ,(cadr (cadr (car c
*logs
)))
1025 ,(rcdisrep (caar c
*logs
))))
1026 (exp-c*logs
(cdr c
*logs
)))))
1028 (defun sum-c*logs
(c*logs
)
1029 (if (null c
*logs
) (rczero)
1030 (let ((ps (cddr (car c
*logs
))))
1031 (psplus (make-ps ps
(ncons (term (ps-le ps
) (caar c
*logs
))))
1032 (sum-c*logs
(cdr c
*logs
))))))
1034 ;; Calculatest the limit of a series at the expansion point. Returns one of
1035 ;; {$zeroa, $zerob, $pos, $neg, $inf, $minf}.
1037 (defvar tvar-limits
()
1038 "A list of the form ((gvar . limit(gvar)) ...)")
1040 (defun ps-lim-infp (ps)
1042 ;; Assume taylor vars at 0+ for now. Should handle the cases when
1043 ;; the expansion point is INF, MINF,etc.
1044 (let* ((lim (gvar-lim (gvar ps
)))
1046 (if (member lim
'($inf $minf
) :test
#'eq
) (ps-gt ps
) (ps-lt ps
))))
1047 (if (ezerop (e strongest-term
))
1048 (ps-lim-infp (c strongest-term
))
1050 (setq lim
(lim-power lim
(e strongest-term
)))
1051 (and (lim-infp lim
) (not (eq lim
'$infinity
))))))))
1053 (defun lim-zerop (lim)
1054 (member lim
'($zeroa $zerob $zeroim
) :test
#'eq
))
1056 (defun lim-plusp (lim)
1057 (member lim
'($zeroa $pos $inf $finite
) :test
#'eq
))
1059 (defun lim-finitep (lim)
1060 (member lim
'($pos $neg $im $finite
) :test
#'eq
))
1062 (defun lim-infp (lim)
1063 (member lim
'($inf $minf $infinity
) :test
#'eq
))
1065 (defun lim-imagp (lim)
1066 (member lim
'($im $infinity
) :test
#'eq
))
1068 (defun lim-minus (lim)
1069 (cdr (assoc lim
'(($zeroa . $zerob
) ($zerob . $zeroa
) ($pos . $neg
) ($zero . $zero
)
1070 ($neg . $pos
) ($inf . $minf
) ($minf . $inf
)
1071 ($im . $im
) ($infinity . $infinity
) ($finite . $finite
)) :test
#'eq
)))
1072 (defun lim-abs (lim)
1073 (or (cdr (assoc lim
'(($zerob . $zeroa
) ($neg . $pos
) ($minf . $inf
)) :test
#'eq
))
1076 (defun lim-times (lim1 lim2
)
1078 (cond ((or (eq lim1
'$zero
) (eq lim2
'$zero
)) (setq lim
'$zero
))
1079 ((and (lim-infp lim1
) (lim-infp lim2
)) (setq lim
'$inf
))
1080 ((and (lim-zerop lim1
) (lim-zerop lim2
)) (setq lim
'$pos
))
1081 ((or (when (lim-finitep lim2
) (rotatef lim1 lim2
) 't
)
1083 (when (and (eq lim1
'$finite
) (lim-infp lim1
))
1084 (tay-error "Undefined finite*inf in lim-times" lim2
))
1085 (setq lim
(lim-abs lim2
)))
1086 (t (tay-error "Undefined limit product in lim-times" (list (list 'mtimes
) lim1 lim2
))))
1087 (if (or (lim-imagp lim1
) (lim-imagp lim2
))
1088 (if (lim-infp lim
) '$infinity
'$im
)
1089 (if (and (lim-plusp lim1
) (lim-plusp lim2
)) lim
(lim-minus lim
)))))
1090 (defun lim-power (lim power
)
1091 (cond ((ezerop power
) '$pos
)
1092 ((e> (rczero) power
) (lim-recip (lim-power lim
(e- power
))))
1093 ((not (oddp (car power
)))
1094 (if (lim-plusp lim
) lim
(lim-minus lim
)))
1097 (defun lim-recip (lim)
1098 (or (cdr (assoc lim
'(($zeroa . $inf
) ($zerob . $minf
)
1099 ($inf . $zeroa
) ($minf . $zerob
)) :test
#'eq
))
1100 (if (eq lim
'$finite
) (break "inverting $finite?")
1103 (defun lim-exp (lim)
1105 (($zeroa $zerob $zero $pos $neg $minf
) '$zeroa
)
1106 (($inf $finite
) lim
)
1107 ($infinity
'$infinity
) ; actually only if Re lim = inf
1108 (t (break "Unhandled limit in lim-exp"))))
1110 (defun lim-log (lim)
1116 (t (break "Unhandled limit in lim-log"))))
1118 (defun expand-and-disrep (term p
)
1119 (let ((x^n
(list '(mexpt) (get-inverse (gvar p
)) (edisrep (e term
))))
1121 (if (pscoefp a
) (ncons (m* (srdis a
) x^n
))
1122 (mapcar #'(lambda (subterm)
1123 (m* (cons '(mtimes) (expand-and-disrep subterm a
)) x^n
))
1126 (defun adjoin-sing-datum (d)
1127 (let ((r (prep1 (datum-var d
))) (g (gensym)) (kernel (datum-var d
))
1128 (no (1+ (cdr (int-var (car (last tlist
)))))))
1129 (unless (and (equal (car r
) 1) (equal (cddr r
) '(1 1)))
1130 (break "bad singular datum"))
1131 (putprop g kernel
'disrep
)
1132 (rplacd (cdddr d
) (cons g no
))
1134 (push (cons (cadr r
) kernel
) key-vars
)
1135 (push (cons g kernel
) key-vars
)
1136 (push (car key-vars
) ivars
)
1137 ;(push (cons kernel (cons (pget g) 1)) genpairs)
1138 (push (cons g
(exp-pt d
)) tvar-limits
)))
1140 (defun adjoin-tvar (exp) (rat->ps
(prep1 exp
)))
1142 (defun rat->ps
(rat)
1143 (pstimes (poly->ps
(car rat
))
1144 (psexpt (poly->ps
(cdr rat
)) (rcmone))))
1146 (defun poly->ps
(poly)
1147 (if (or (pcoefp poly
) (mfree (pdis poly
) tvars
)) (prep1 poly
)
1148 (let ((g (p-var poly
)) datum
(pow (rcone)))
1149 (if (setq datum
(key-var-pow g
)) (desetq (g . pow
) datum
)
1150 (desetq (g . pow
) (adjoin-pvar g
)))
1151 (if (and (not (atom g
)) (psp g
))
1154 (setq datum
(gvar-data g
))
1155 (do ((po-terms (p-terms poly
) (p-red po-terms
))
1157 (push (term (e* pow
(prep1 (pt-le po-terms
)))
1158 (poly->ps
(pt-lc po-terms
)))
1161 ;; This must be exact so that when we invert in rat-ps above
1162 ;; we don't lose any terms. E.g. try
1163 ;; taylor(log(1+exp(-1/x)),x,0,5). When taylor2'ing exp(-1/x),
1164 ;; if you used current trunc here this would return exp(1/x)...5
1165 ;; which would then be trunc'd to degree 3 by psexpt.
1166 (make-ps (int-var datum
)
1167 (ncons (current-trunc datum
))
1168 (if (eq g
(data-gvar datum
)) ps-terms
1169 (invert-terms ps-terms
))))))))))
1171 (defun key-var-pow (g)
1172 (let ((var (get-key-var g
)) datum
)
1174 (setq datum
(get-datum var
))
1175 (if (eq g
(setq g
(data-gvar datum
))) (cons g
(rcone))
1176 (cons g
(rcmone))))))
1178 (defun adjoin-pvar (g)
1179 (let ((kernel (get g
'disrep
)) g
* lim datum ans
1180 (no (1+ (cdr (int-var (car (last tlist
)))))) (pow (rcone)) expt
)
1181 (when (assol kernel tlist
) (break "bad1"))
1182 (if (and (eq (caar kernel
) 'mexpt
) (eq (cadr kernel
) '$%e
)
1183 (not (atom (setq expt
(caddr kernel
))))
1184 (eq (caar expt
) 'mtimes
) (not (mfree expt
(ncons '$%i
))))
1185 (destructuring-let (((rpart . ipart
) (trisplit expt
)))
1186 (cons (pstimes (prep1 (m^
'$%e rpart
))
1187 (psplus (adjoin-tvar `((%cos
) ,ipart
))
1188 (pstimes (prep1 '$%i
)
1189 (adjoin-tvar `((%sin
) ,ipart
)))))
1192 (when (eq (caar kernel
) 'mexpt
)
1193 (when (and (not (atom (setq expt
(caddr kernel
))))
1194 (eq (caar expt
) 'mtimes
)
1195 ($ratnump
(cadr expt
)))
1196 (setq pow
(cadr expt
) kernel
(m^ kernel
(m// pow
))
1197 g
(prep1 kernel
) pow
(prep1 pow
))
1198 (unless (and (equal (cdr g
) 1) (equal (cdar g
) '(1 1)))
1199 (break "Illegal kernel in `adjoin-pvar'"))
1200 (setq g
(caar g
) kernel
(get g
'disrep
))))
1201 (if (setq ans
(key-var-pow g
))
1202 (cons (car ans
) (e* pow
(cdr ans
)))
1204 (when (lim-infp (or lim
(setq lim
(tvar-lim kernel
))))
1205 (setq g
* g g
(gensym) kernel
(m// kernel
)
1206 lim
(lim-recip lim
) pow
(e* (rcmone) pow
))
1207 (putprop g kernel
'disrep
)
1208 ;(push g genvar) (push kernel varlist)
1209 (push (cons g
* kernel
) key-vars
))
1210 (when (assol kernel tlist
) (break "bad2"))
1211 (setq datum
(list* kernel
1212 ;(mapcar #'(lambda (e) (emax e (rczero)))
1213 ; (trunc-stack (car tlist)))
1214 (copy-list (trunc-stack (car tlist
)))
1216 ;(setq tlist (nconc tlist (ncons datum)))
1217 (adjoin-datum datum
)
1218 (push (cons g kernel
) key-vars
)
1219 (push (car key-vars
) ivars
)
1220 ;(push (cons kernel (cons (pget g) 1)) genpairs)
1221 (push (cons g lim
) tvar-limits
)
1224 (defun adjoin-datum (datum)
1225 (do ((tlist* tlist
(cdr tlist
*))
1226 (tlist** () tlist
*))
1227 ((null tlist
*) (setq tlist
(nconc tlist
(ncons datum
))))
1228 (when (stronger-var?
(datum-var (car tlist
*)) (datum-var datum
))
1229 (return (if (null tlist
**)
1232 (renumber-tlist tlist
))
1234 (rplacd tlist
** (cons datum tlist
*))
1235 (renumber-tlist (cdr tlist
**))))))))
1237 ;; Maybe this should just permute the numbering in case it isn't sequential?
1239 (defun renumber-tlist (tlist)
1240 (rplacd (data-gvar-o (car tlist
)) (cdr (data-gvar-o (cadr tlist
))))
1241 (do ((tlist* (cdr tlist
) (cdr tlist
*)))
1243 (rplacd (data-gvar-o (car tlist
*))
1244 (1+ (cdr (data-gvar-o (car tlist
*)))))))
1246 (defun tvar?
(var) (or (atom var
) (member 'array
(cdar var
) :test
#'eq
)))
1248 ;; Needs to be extended to handle singular tvars in > 1 var's.
1250 (defun stronger-var?
(v1 v2
)
1251 (let ((e1 (rcone)) (e2 (rcone)) reverse? ans
)
1252 (when (alike1 v1 v2
)
1253 (tay-err (intl:gettext
"taylor: stronger-var? called on equal vars.")))
1254 (when (and (mexptp v1
) ($ratnump
(caddr v1
)))
1255 (setq e1
(prep1 (caddr v1
)) v1
(cadr v1
)))
1256 (when (and (mexptp v2
) ($ratnump
(caddr v2
)))
1257 (setq e2
(prep1 (caddr v2
)) v2
(cadr v2
)))
1261 (intl:gettext
"taylor: stronger-var? called on equal vars."))
1264 (when (eq (tvar-lim v2
) '$finite
)
1267 (setq reverse?
(not reverse?
)))
1268 (if (eq (tvar-lim v1
) '$finite
)
1269 (if (eq (tvar-lim v2
) '$finite
)
1270 (great v1 v2
) reverse?
)
1275 (setq reverse?
(not reverse?
)))
1278 (stronger-vars?
(order-vars-by-strength (cdr v1
))
1279 (order-vars-by-strength (if (mtimesp v2
) (cdr v2
)
1280 (ncons (m^ v2
(edisrep e2
))))))
1285 (setq reverse?
(not reverse?
)))
1288 (let ((n1 (cdr (data-gvar-o (get-datum v1 t
))))
1289 (n2 (cdr (data-gvar-o (get-datum v2 t
)))))
1291 ((mfree v2
(ncons v1
))
1293 (intl:gettext
"taylor: Unhandled multivar datum comparison.")))
1294 ((eq (caar v2
) '%log
) 't
)
1295 ((and (eq (caar v2
) 'mexpt
) (eq (cadr v2
) '$%e
))
1296 (stronger-var?
`((%log
) ,v1
) (caddr v2
)))
1297 (t (tay-err (intl:gettext
"taylor: Unhandled var in stronger-var?."))))
1299 (when (eq (caar v2
) '%log
)
1302 (setq reverse?
(not reverse?
)))
1303 (if (eq (caar v1
) '%log
)
1304 (cond ((eq (caar v2
) '%log
)
1305 (stronger-var?
(cadr v1
) (cadr v2
)))
1306 ((and (eq (caar v2
) 'mexpt
) (eq (cadr v2
) '$%e
))
1307 (stronger-var?
`((%log
) ,v1
) (caddr v2
)))
1308 (t (tay-err (intl:gettext
"taylor: Unhandled var in stronger-var?"))))
1309 (if (and (eq (caar v1
) 'mexpt
) (eq (cadr v1
) '$%e
)
1310 (eq (caar v2
) 'mexpt
) (eq (cadr v2
) '$%e
))
1311 (stronger-var?
(caddr v1
) (caddr v2
))
1312 (tay-err (intl:gettext
"taylor: Unhandled var in stronger-var?")))))))))
1313 (if reverse?
(not ans
) ans
)))))))
1315 (defun neg-monom?
(exp)
1316 (and (mtimesp exp
) (equal (cadr exp
) -
1) (null (cdddr exp
))
1319 (defun order-vars-by-strength (vars)
1320 (do ((vars* vars
(cdr vars
*)) (ordvars () ))
1321 ((null vars
*) ordvars
)
1322 (unless (mfree (car vars
*) tvars
) ; ignore constants
1323 (do ((ordvars* ordvars
(cdr ordvars
*)))
1325 (if (null ordvars
) (setq ordvars
(ncons (car vars
*)))
1326 (rplacd (last ordvars
) (ncons (car vars
*)))))
1327 (when (stronger-var?
(car vars
*) (car ordvars
*))
1328 (rplacd ordvars
* (cons (car ordvars
*) (cdr ordvars
*)))
1329 (rplaca ordvars
* (car vars
*))
1332 (defun stronger-vars?
(vars1 vars2
)
1333 (do ((vars1* vars1
(cdr vars1
*))
1334 (vars2* vars2
(cdr vars2
*)))
1336 (cond ((null vars1
*)
1338 ;; two equal vars generated
1340 (let ((lim (tvar-lim (car vars2
*))))
1342 (cond ((lim-infp lim
) ())
1343 ((lim-zerop lim
) 't
)
1344 (t (break "var with non-zero finite lim?")))))))
1346 (let ((lim (tvar-lim (car vars1
*))))
1348 (cond ((lim-infp lim
) 't
)
1349 ((lim-zerop lim
) ())
1350 (t (break "var with non-zero finite lim?"))))))
1351 ((alike1 (car vars1
*) (car vars2
*)) )
1352 ((return (stronger-var?
(car vars1
*) (car vars2
*)))))))
1354 (defun stronger-datum?
(d1 d2
)
1355 (setq d1
(datum-var d1
) d2
(datum-var d2
))
1356 (do ((end-flag) (answer))
1357 (end-flag (member answer
'($yes $y
) :test
#'eq
))
1358 (setq answer
(retrieve `((mtext) |Is |
,d1 | stronger than |
,d2 |?|
)
1360 (if (member answer
'($yes $y $no $n
) :test
#'eq
) (setq end-flag
't
)
1361 (mtell "~%Acceptable answers are: yes, y, no, n~%"))))
1363 (defun datum-lim (datum)
1364 (if (not (tvar?
(datum-var datum
)))
1366 (let ((pt (exp-pt datum
)))
1367 (if (member pt
'($inf $minf
) :test
#'eq
) pt
'$zeroa
))))
1369 (defun tvar-lim (kernel)
1370 (if (mfree kernel tvars
) (coef-sign kernel
)
1371 (let ((datum (get-datum kernel t
)) lim
)
1372 (or (and datum
(datum-lim datum
))
1373 (and (setq datum
(get-datum (m// kernel
) t
))
1374 (setq lim
(datum-lim datum
))
1378 (cond ((eq (caar kernel
) 'mexpt
)
1379 (cond ((and (setq datum
(get-datum (cadr kernel
) t
))
1380 ($ratnump
(caddr kernel
)))
1381 (lim-power (datum-lim datum
)
1382 (prep1 (caddr kernel
))))
1383 (($ratnump
(caddr kernel
))
1384 (lim-power (tvar-lim (cadr kernel
))
1385 (prep1 (caddr kernel
))))
1386 ((eq (cadr kernel
) '$%e
)
1387 (lim-exp (tvar-lim (caddr kernel
))))
1388 (t (tay-error "Unhandled case in tvar-lim" kernel
))))
1389 ((eq (caar kernel
) 'mtimes
)
1390 (do ((ans (tvar-lim (cadr kernel
))
1391 (lim-times ans
(tvar-lim (car facs
))))
1392 (facs (cddr kernel
) (cdr facs
)))
1394 ((eq (caar kernel
) '%log
)
1395 (lim-log (datum-lim (get-datum (cadr kernel
) t
))))
1396 ((member (caar kernel
) '(%sin %cos
) :test
#'eq
)
1397 (unless (lim-infp (tvar-lim (cadr kernel
)))
1398 (tay-error "Invalid trig kernel in tvar-lim" kernel
))
1400 (t (tay-error "Unhandled kernel in tvar-lim" kernel
))))
1403 (defun coef-sign (coef)
1404 (if (not ($freeof
'$%i
($rectform coef
)))
1408 (defun gvar-lim (gvar)
1409 (or (cdr (assoc gvar tvar-limits
:test
#'eq
))
1410 (if (member (gvar->var gvar
) tvars
:test
#'eq
) '$zeroa
; user tvars assumed 0+ now
1411 (break "Invalid gvar"))))
1413 (defun psexpt-fn1 (varh trunc l inc m mr ans
)
1414 (declare (fixnum m
))
1415 (prog (a (k 0) ak cm-k c sum kr lim
)
1416 (declare (fixnum k
))
1417 ;; truly unfortunate that we need so many variables in this hack
1419 b
(and (e> mr trunc
) (go end
))
1420 (setq kr inc ak l k
1 sum
(rczero) lim m
)
1421 a
(cond ((or (> k lim
)
1422 (null (setq cm-k
(psterm (cdr ans
) (e- mr kr
)))))
1424 (setq ak
(or (pstrim-terms ak kr
)
1426 c
(pstimes (ereduce k m
)
1427 (pstimes (psterm1 ak kr
) cm-k
))
1429 (setq k
(1+ k
) kr
(e+ kr inc
))
1432 (unless (rczerop sum
) (add-term-&-pop a mr sum
))
1433 (setq m
(1+ m
) mr
(e+ mr inc
))
1436 (return (pscheck varh
(list trunc
) (cdr ans
)))))
1438 ;;; PSEXPT-FN2 and RED-MONO-LOG are needed to reduce exponentials of logs.
1440 (defun psexpt-fn2 (p)
1441 (cond ((atom p
) (if (get-datum p
)
1442 (psexpt-fn (taylor2 p
))
1443 (prep1 `((mexpt) $%e
,p
))))
1444 ((eq (caar p
) '%log
)
1445 (if (get-datum (cadr p
)) (taylor2 (cadr p
)) (prep1 (cadr p
))))
1446 ((or (eq (caar p
) 'mplus
) (eq (caar p
) 'mtimes
))
1447 (let ((e ($ratexpand p
)) temp
)
1448 (cond ((not (and (consp e
) (member (caar e
) '(mplus mtimes
) :test
#'eq
)))
1451 (if (eq (caar e
) 'mplus
)
1452 (do ((sumnds (cdr e
) (cdr sumnds
)) (log-facs) (l))
1454 (cond ((not log-facs
) (tsexpt '$%e p
))
1455 (t (tstimes (cons (m^t
'$%e
(m+l l
)) log-facs
)))))
1456 (if (setq temp
(red-mono-log (car sumnds
)))
1457 (push temp log-facs
)
1458 (push (car sumnds
) l
)))
1460 (setq temp
(red-mono-log e
))
1463 (prep1 (power '$%e p
)))))))))
1464 (t (prep1 (power '$%e p
)))))
1466 (defun red-mono-log (e)
1468 ((eq (caar e
) '%log
) (cadr e
))
1470 (do ((facs (cdr e
) (cdr facs
)) (log-term))
1473 (m^t
(cadr log-term
) (m*l
(remove log-term
(cdr e
) :test
#'eq
)))))
1474 (if (and (null (atom (car facs
))) (eq (caaar facs
) '%log
))
1475 (if log-term
(return ()) (setq log-term
(car facs
)))
1476 (unless (mfree (car facs
) tvars
) (return nil
)))))
1480 (if (pscoefp p
) (pslog2 (rcdisrep p
))
1481 (let ((terms (terms p
)))
1482 (cond ((mono-term? terms
) ; log(c x^n) = log(c) + n log(x)
1483 ;; do this always for now
1484 (if 't
;$TAYLOR_LOGEXPAND
1485 ;(psplus (pslog (lc terms))
1486 ; (pstimes (le terms) (pslog-of-gvar (gvar p))))
1488 ;(prep1 `((%LOG) ,(term-disrep (lt terms) p)))
1490 ;; expand log(1+ax^n) directly by series substitution
1491 ((not (or (n-term (setq terms
(terms (psplus p
(rcmone)))))
1492 ;(e> (rczero) (le terms))
1494 (setq p
(get-series '%log
(psexpt-log-ord p
) (gvar-o p
)
1495 (if (e> (rczero) (le terms
)) (e- (le terms
))
1498 (if (e> (rczero) (le terms
)) (ps-invert-var p
) p
))
1499 (t (prog (l inc trunc lt ans lterm $maxtayorder gvar-lim gt
)
1500 ;; log(lt+y) = log(lt) + log(1 + y/lt) = lterm + p
1501 (setq trunc
(trunc-lvl p
))
1502 (if (not (member (setq gvar-lim
(gvar-lim (gvar p
)))
1503 '($zeroa $zerob $inf $minf
) :test
#'eq
))
1504 (tay-error "bad gvar lim" gvar-lim
)
1505 (if (member gvar-lim
'($inf $minf
) :test
#'eq
)
1506 (setq lt
(ps-gt p
) gt lt
)
1507 (setq lt
(ps-lt p
) gt
() )))
1509 (setq lt
(pscheck (gvar-o p
)
1512 p
(pstimes p
(let (($maxtayorder
't
))
1513 (psexpt lt
(rcmone)))))
1514 (when (and (member gvar-lim
'($inf $minf
) :test
#'eq
)
1515 (e> (le terms
) (rczero)))
1516 (return (psplus lterm
(pslog p
))))
1518 (unless (equal p
(rcone))
1519 (merror "PSLOG: internal error."))
1521 (setq l
(terms p
) inc
(psexpon-gcd l
))
1522 (if gt
(setq l
(delete (last l
) l
:test
#'equal
))
1523 (setq l
(n-term l
)))
1524 (setq ans
(ncons 0))
1525 (unless $maxtayorder
1526 (setq trunc
(emin trunc
(t-o-var (gvar p
)))))
1527 ;; When we've divided by the greatest term, all terms
1528 ;; have non-positive exponents and we must perform the
1529 ;; transformation x -> 1/x before calling pslog1 and then
1530 ;; perform the inverse afterwards.
1531 (when gt
(setq l
(invert-terms l
)))
1532 (when (e> (rczero) inc
) (setq inc
(e- inc
)))
1533 (setq ans
(psplus lterm
1534 (pslog1 (gvar-o p
) trunc l inc
1 inc ans
)))
1536 (if (and gt
(psp ans
) (eq (gvar ans
) (gvar p
)))
1540 (defun invert-terms (terms)
1541 (nreverse (mapc #'(lambda (x) (rplaca x
(e- (e x
)))) terms
)))
1543 (defun ps-invert-var (ps)
1544 (when (psp ps
) (rplacd (cddr ps
) (invert-terms (terms ps
))))
1548 (if (pscoefp ps
) (term (rczero) ps
)
1549 (lt (last (terms ps
)))))
1551 (defun pslog1 (varh trunc l inc m mr ans
)
1552 (declare (fixnum m
))
1553 (prog (a (k 0) ak cm-k c sum kr m-kr
)
1554 (declare (fixnum k
))
1555 ;; truly unfortunate that we need so many variables in this hack
1558 b
(and (e> mr trunc
) (go end
))
1559 (setq kr inc ak l k
1 sum
(rczero))
1560 a
(cond ((or (= k m
)
1561 (null (setq cm-k
(psterm (cdr ans
)
1562 (setq m-kr
(e- mr kr
))))))
1564 (setq ak
(or (pstrim-terms ak kr
)
1566 c
(pstimes m-kr
(pstimes (psterm1 ak kr
) cm-k
))
1568 k
(1+ k
) kr
(e+ kr inc
))
1571 (cond ((setq c
(pstrim-terms ak mr
))
1572 (setq c
(psterm1 c mr
)))
1573 ((setq c
(rczero))))
1574 (setq sum
(psdiff c
(pstimes sum
(e// mr
))))
1575 (unless (rczerop sum
) (add-term-&-pop a mr sum
))
1576 (setq m
(1+ m
) mr
(e+ mr inc
))
1579 (return (pscheck varh
(list trunc
) (cdr ans
)))))
1581 ;; Computes log(monom), where monom = c x^n. Is extra careful trying to keep
1582 ;; singular logs going to INF and not generating log(-1)'s unless it is
1583 ;; necessary to transform a log at MINF to INF.
1585 (defun pslog-monom (monom)
1586 (let* ((gvar (gvar monom
))
1587 (datum (gvar-data gvar
)) var pt logvar c
)
1588 (if (switch 'multivar datum
)
1589 (pslog (ps-lc monom
))
1591 (setq var
(datum-var datum
))
1593 (if (not (member (setq pt
(exp-pt datum
)) '($inf $minf
) :test
#'eq
))
1594 (setq logvar
(adjoin-tvar `((%log
) ,(m- var pt
))))
1596 ;; At x = inf: log(c (1/x)^n) -> log(c) - n log(x)
1597 ;; At x = minf: log(c (-1/x)^n) -> log(c (-1)^n) - n log(x)
1598 (setq logvar
(psminus (adjoin-tvar `((%log
) ,var
))))
1599 (when (eq pt
'$minf
)
1600 (setq c
(rcexpt (rcmone) (ps-le monom
))))))
1601 (if (eq (caar var
) 'mexpt
)
1602 (if (equal (caddr var
) -
1);; var must be 1/log(y)
1603 ;; Try to keep inf's real. Here we want
1604 ;; log(c (1/log(x))^n) -> log(c (-1)^n) - n log(-log(x))
1605 (if (equal (tvar-lim (cadr var
)) '$minf
)
1606 (setq c
(rcexpt (rcmone) (ps-le monom
))
1608 (psminus (adjoin-tvar
1609 `((%log
) ,(m- (cadr var
))))))
1610 (setq logvar
(psminus
1611 (adjoin-tvar `((%log
) ,(cadr var
))))))
1612 (if (equal (cadr var
) '$%e
)
1613 (setq logvar
(taylor2 (caddr var
)))
1614 (break "Unhandled gvar in `pslog-of-gvar'")))))
1615 (psplus (pslog (if c
(pstimes c
(ps-lc monom
)) (ps-lc monom
)))
1616 (pstimes (ps-le monom
) logvar
))))))
1618 ;; Computes log(p), where p is an rcdisrep'd pscoef.
1620 (defun pslog2 (p) (let ($logarc
) (pslog3 p
)))
1624 (prep1 (cond ((equal p
1) 0)
1625 ((equal p -
1) log-1
)
1629 (merror (intl:gettext
"taylor: log(0) encountered while processing ~:M") last-exp
))
1632 (prep1 (cond ((not $taylor_logexpand
) `((%log
) ,p
))
1633 (t (m- `((%log
) ,(cadr p
)) `((%log
) ,(caddr p
)))))))
1634 ((and full-log
(not (free p
'$%i
)))
1635 (let ((full-log () )) (pslog3 ($polarform p
))))
1636 ((eq (caar p
) 'mexpt
)
1637 ;; Must handle things like x^a, %e^(a*x), etc. which are pscoef's.
1638 (pstimes (taylor2 (caddr p
)) (pslog (taylor2 (cadr p
)))))
1639 ((and (eq (caar p
) 'mtimes
) $taylor_logexpand
)
1640 (do ((l (cddr p
) (cdr l
))
1641 (ans (pslog3 (cadr p
)) (psplus ans
(pslog3 (car l
)))))
1643 (t (prep1 `((%log
) ,p
)))))
1645 ;;; Subtitle Extending Routines
1647 (defun getfun-lt (fun)
1648 (let ((exp-datum (get (oper-name fun
) 'exp-form
)))
1650 ;; Info not needed yet.
1651 ;; (or (atom (car exp-datum))
1652 ;; (setq 0p-funord (copy-tree (cdar exp-datum))))
1653 (exp-datum-lt fun exp-datum
))
1654 ((setq exp-datum
(get (oper-name fun
) 'sp2
))
1655 (setq exp-datum
(get-lexp (subst (dummy-var) 'sp2var exp-datum
)
1657 ;; Info not needed yet; need to bind lexp-non0 to T when
1658 ;; this is used though so n-term will be there.
1659 ;; (and (rczerop (le exp-datum))
1660 ;; (setq 0p-funord (le (n-term exp-datum))))
1661 (if (psp exp-datum
) (ps-lt exp-datum
)
1662 (term (rczero) exp-datum
)))
1663 (t (merror "GETFUN-LT: unknown function ~A" fun
)))))
1665 (declare-top (special var
))
1667 (defun getexp-fun (fun var pw
)
1668 (declare (special var
))
1669 (let ((exp-datum (copy-tree (get (oper-name fun
) 'exp-form
))))
1670 (cond ((infp pw
) (infin-ord-err))
1672 (if (null (setq exp-datum
1673 (get-ps-form (if (atom fun
) fun
(caar fun
)))))
1674 (merror (intl:gettext
"taylor: power series unavailable for function ~A") fun
)
1677 (do ((subvals (cdr fun
) (cdr subvals
))
1678 (subs (safe-get (caar fun
) 'sp2subs
) (cdr subs
)))
1679 ((or (null subvals
) (null subs
))
1680 (when (or subvals subs
)
1681 (merror (intl:gettext
"taylor: incorrect number of subscripts to the deftaylor'd function ~A") (caar fun
))))
1682 (setq exp-datum
(maxima-substitute (car subvals
) (car subs
)
1684 (ts-formula exp-datum var pw
))))
1685 ((e> (exp-datum-le fun exp-datum
) pw
) (pszero var pw
))
1687 (apply (exp-fun exp-datum
)
1688 (if (atom fun
) (cons pw
(cdr exp-datum
))
1689 (cons pw
(cons (cdr fun
) (cdr exp-datum
))))))
1690 (cond ((null exp-datum
) (pszero var pw
))
1691 ((psp exp-datum
) exp-datum
)
1692 (t (make-ps var
(ncons pw
) exp-datum
)))))))
1694 (declare-top (unspecial var
))
1696 (defun expexp-funs (pw l sign chng inc
)
1698 (setq e
(e l
) lt-l
(setq l
(ncons l
)))
1699 a
(cond ((e> (setq e
(e+ e inc
)) pw
) (return l
))
1704 (cond ((e= inc
(rcone)) e
)
1706 (cons 1 (cdr (lc lt-l
)))))
1707 (setq sign
(e* sign chng
))))
1710 ;; returns series expansion of %expintegral_si
1711 ;; from term x^i to term x^max
1712 ;; sign is sign of term x^i (1 or -1)
1713 ;; ifac is i! (to avoid repeated computation)
1714 (defun expsi_series (i sign ifac max
)
1716 nil
; we have all powers up to max
1717 (cons (cons (cons i
1) ; this is i'th term
1718 (cons sign
(* i ifac
))); i'th coefficient
1719 (expsi_series (+ 2 i
) ; advance to next even/odd power
1720 (* -
1 sign
) ; flip sign
1721 (* (+ 2 i
) (+ 1 i
) ifac
) ; update factorial
1722 max
)))) ; pass through the max
1724 ;; returns series expansion of %expintegral_si up to term with x^pw
1725 (defun exp_%expintegral_si
(pw l
)
1729 (/ (float (car pw
)) (float (cdr pw
)))))
1731 (defun explog-funs (pw l sign chng inc
)
1733 (setq e
(e l
) lt-l
(setq l
(ncons l
)))
1734 a
(cond ((e> (setq e
(e+ e inc
)) pw
) (return l
))
1735 (t (add-term lt-l e
(e// sign e
))
1736 (setq lt-l
(n-term lt-l
)
1737 sign
(e* sign chng
))))
1740 (defun exptan-funs (pw l chng
)
1741 (prog (e lt-l sign fact pow
)
1742 (setq e
(e l
) lt-l
(setq l
(ncons l
))
1743 sign
(rcone) fact
'(1 .
2) pow
'(4 .
1))
1744 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1745 (t (setq fact
(e// fact
(e* e
(e1+ e
)))
1746 pow
(e* '(4 .
1) pow
)
1747 sign
(e* chng sign
))
1748 (add-term lt-l e
(e* (e* sign fact
)
1750 (ftake '%bern
(rcdisrep (e1+ e
))))
1751 (e* pow
(e1- pow
)))))
1752 (setq lt-l
(n-term lt-l
))))
1755 (defun expcot-funs (pw l sign chng plus
)
1756 (prog (e lt-l fact pow
)
1757 (setq e
(e l
) lt-l
(setq l
(ncons l
))
1758 fact
(rcone) pow
(rcone))
1759 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1760 (t (setq fact
(e// fact
(e* e
(e1+ e
)))
1761 pow
(e* '(4 .
1) pow
)
1762 sign
(e* chng sign
))
1763 (add-term lt-l e
(e* (e* sign fact
)
1765 (ftake '%bern
(rcdisrep (e1+ e
))))
1767 (setq lt-l
(n-term lt-l
))))
1770 (defun expsec-funs (pw l chng
)
1771 (prog (e lt-l sign fact
)
1772 (setq e
(e l
) lt-l
(setq l
(ncons l
))
1773 sign
(rcone) fact
(rcone))
1774 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1775 (t (setq fact
(e// fact
(e* e
(e1- e
)))
1776 sign
(e* chng sign
))
1777 (add-term lt-l e
(e* (e* sign fact
)
1778 (prep1 (ftake '%euler
(rcdisrep e
)))))
1779 (setq lt-l
(n-term lt-l
))))
1782 (defun expasin-funs (pw l chng
)
1783 (prog (e lt-l sign n d
)
1784 (setq e
(e l
) lt-l
(setq l
(ncons l
)) sign
1 n
1 d
1)
1785 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1786 (t (setq n
(* n
(car (e- e
(rctwo))))
1787 d
(* d
(car (e1- e
)))
1789 (add-term lt-l e
; need to reduce here ? - check this.
1790 (let ((x (*red
(* n sign
)
1793 (cons (cadr x
) (caddr x
)))))
1794 (setq lt-l
(n-term lt-l
))))
1797 ;;; This is the table of expansion data for known functions.
1798 ;;; The format of the EXP-FORM property is as follows:
1799 ;;; (<name of the expanding routine for the function or
1800 ;;; (name . le of n-term) if expansion is of order 0>
1801 ;;; <first term in the expansion or the name of a routine which
1802 ;;; computes the order when it may depend on parameters (e.g subscripts)>
1803 ;;; <data for the expanding routine>)
1806 (loop for
(fun exp
) on
1807 '(%exp
((expexp-funs 1 .
1) ((0 .
1) 1 .
1) (1 .
1) (1 .
1) (1 .
1))
1808 %sin
(expexp-funs ((1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (2 .
1))
1809 %cos
((expexp-funs 2 .
1) ((0 .
1) 1 .
1) (-1 .
1) (-1 .
1) (2 .
1))
1810 %sinh
(expexp-funs ((1 .
1) 1 .
1) (1 .
1) (1 .
1) (2 .
1))
1811 %cosh
((expexp-funs 2 .
1) ((0 .
1) 1 .
1) (1 .
1) (1 .
1) (2 .
1))
1812 %log
(explog-funs ((1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (1 .
1))
1813 %atan
(explog-funs ((1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (2 .
1))
1814 %atanh
(explog-funs ((1 .
1) 1 .
1) (1 .
1) (1 .
1) (2 .
1))
1815 %cot
(expcot-funs ((-1 .
1) 1 .
1) (1 .
1) (-1 .
1) (0 .
1))
1816 %csc
(expcot-funs ((-1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (-2 .
1))
1817 %csch
(expcot-funs ((-1 .
1) 1 .
1) (-1 .
1) (1 .
1) (-2 .
1))
1818 %coth
(expcot-funs ((-1 .
1) 1 .
1) (1 .
1) (1 .
1) (0 .
1))
1819 %tan
(exptan-funs ((1 .
1) 1 .
1) (-1 .
1))
1820 %tanh
(exptan-funs ((1 .
1) 1 .
1) (1 .
1))
1821 %sec
((expsec-funs 2 .
1) ((0 .
1) 1 .
1) (-1 .
1))
1822 %sech
((expsec-funs 2 .
1) ((0 .
1) 1 .
1) (1 .
1))
1823 %asin
(expasin-funs ((1 .
1) 1 .
1) 1)
1824 %asinh
(expasin-funs ((1 .
1) 1 .
1) -
1)
1825 %gamma
(expgam-fun ((-1 .
1) 1 .
1))
1826 $li
(exp$li-fun li-ord
)
1827 %expintegral_si
(exp_%expintegral_si
((1 .
1) 1 .
1))
1828 $psi
(expplygam-funs plygam-ord
))
1830 do
(putprop fun exp
'exp-form
))
1833 (defun known-ps (fun)
1834 (getl fun
'(exp-form sp2
)))
1836 ;;; Autoload Properties
1838 ;;; Taylor series expansion routines
1840 ;;; SRF is only called externally; by RATF and SIMPEXPT.
1843 (let ((exact-poly t
) (tlist) (*within-srf?
* 't
))
1844 (setq x
(taylor1 x
()) tlist
(mrat-tlist x
))
1845 ;; Set trunc levels in the local tlist to correspond to the maximum
1846 ;; level occurring in any series.
1847 (do ((data tlist
(cdr data
))
1848 (truncs (trunc-vector (mrat-ps x
) () )))
1850 (when (and (car truncs
) (e> (car truncs
) (current-trunc (car data
))))
1851 (setf (current-trunc (car data
)) (car truncs
))))
1854 ;;; [var, pt, order, asymp]
1856 (defmfun $taylor
(e &rest args
)
1857 (when (not ($taylorp e
))
1858 ;; Not a taylor expression. Remove the special representation.
1859 (setq e
(specrepcheck e
)))
1862 (defun taylor* (arg l
)
1863 ;; We must bind $MAXTAYORDER to () below because of the problem of constants
1864 ;; not retaining their truncation level. This means that when we add a
1865 ;; series which has more terms than the user-specified truncation to a
1866 ;; constant we must truncate the series with more terms down to the user
1867 ;; specified level because, in the worst case, the constant could be a
1868 ;; series no better than to the user-specified level. Hence $MAXTAYORDER
1869 ;; is essentially useless until the constant problem is fixed. If one
1870 ;; decides to not bind $MAXTAYORDER below then the sum routines must
1871 ;; be updated to truncate series with more terms than the user-specified
1872 ;; level down to that level---taylor(sin(x)^2-cos(x)^2-1,x,0,1) would
1873 ;; give x^2+... in this case if the sum routines weren't updated.
1874 ;; Also, batch(mquery,160,aljabr) for another truncation bug which crops
1875 ;; up when $maxtayorder isn't bound here. Similarly, loadfile(taybad,rl,
1876 ;; aljabr) and see tomh's bug note of 4/15/81.
1877 (let ((tlist () ) ($maxtayorder
() ) (*within-srf?
* () )
1878 (exact-poly (if l
(not $taylor_truncate_polynomials
) 'user-specified
)))
1879 (declare (special *within-srf?
*))
1882 (taylor1 arg
(ncons tlist
))))
1884 (defun tay-order (n)
1885 (let (($float
) (modulus))
1886 (cond ((eq n
'$inf
) (ncons (inf)))
1887 ((null n
) (wna-err '$taylor
))
1889 (merror (intl:gettext
"taylor: expansion order must be a number; found: ~:M") n
))
1890 (t (ncons (prep1 n
))))))
1892 (defun re-erat (head exp
)
1893 (taylor1 exp
(list (cadddr (cdr head
)))))
1895 (defun parse-tay-args (l)
1898 (merror (intl:gettext
"taylor: variable of expansion cannot be a number: ~M") (car l
)))
1899 ((or (symbolp (car l
)) (not (eq (caaar l
) 'mlist
)))
1900 (parse-tay-args1 (list (car l
) ($ratdisrep
(cadr l
)) (caddr l
)))
1901 (parse-tay-args (cdddr l
)))
1902 ((do ((l (cddar l
) (cdr l
)))
1904 (and (or (mnump (car l
)) (eq (car l
) '$inf
))
1906 (parse-tay-args1 (cdar l
))
1907 (parse-tay-args (cdr l
)))
1908 (t (parse-tay-args2 (list (car l
) (cadr l
) (caddr l
)))
1909 (parse-tay-args (cdddr l
)))))
1911 (defun parse-tay-args1 (l)
1912 (if ($listp
(car l
)) (parse-tay-args2 l
)
1914 (pt ($ratdisrep
(cadr l
)))
1915 (ord (tay-order (caddr l
)))
1916 (switches (make-switch-list (cdddr l
))))
1917 (push (list v ord pt switches
) tlist
))))
1919 (defun parse-tay-args2 (l)
1920 (let ((label (gensym))
1922 (pts (make-long-list (if ($listp
(cadr l
))
1923 (copy-list (cdadr l
))
1924 (ncons (ratdisrep (cadr l
))))))
1926 (switches (make-switch-list (cdddr l
)))
1930 (setq lcm ord max ord ord
(make-long-list (ncons ord
)))
1932 (l (cdr ord
) (cdr l
)))
1933 ((null a
) (setq ord
(cdr ord
)))
1934 (cond ((not l
) (merror "PARSE-TAY-ARGS2: ran out of truncation levels."))
1935 (t (setq lcm
(lcm lcm
(car l
)) max
(max max
(car l
)))))))
1936 (push (list label
(tay-order max
) 0
1937 (ncons (list 'multivar lcm vs
)))
1939 (do ((vl vs
(cdr vl
))
1940 (ordl ord
(cdr ordl
))
1941 (ptl pts
(cdr ptl
)))
1943 (cond ((not ptl
) (merror "PARSE-TAY-ARGS2: ran out of matching points of expansion."))
1946 (list (car vl
) (tay-order (car ordl
)) (car ptl
)
1947 (cons (list 'multi label
(timesk lcm
(expta (car ordl
) -
1))) switches
))
1950 (defun make-switch-list (l)
1951 (mapcar #'(lambda (q) (cons q t
)) l
))
1953 (defun make-long-list (q)
1956 ;;; This checks to ensure that there isn't more than one set of multi-
1957 ;;; dependent variables with different orders of expansion, e.g.
1958 ;;; taylor(u+v+x+y,[u,v],0,[2,3],[x,y],0,[5,7]) is one.
1960 (defun ratwtsetup (l)
1961 (do ((l l
(cdr l
)) (a) (sw))
1963 (when (setq a
(switch 'multivar
(car l
)))
1964 (do ((ll (cadr a
) (cdr ll
)))
1966 (cond ((equal (cadr (switch 'multi
(get-datum (car ll
)))) 1) )
1967 (sw (merror (intl:gettext
"taylor: multiple dependent variables must all have the same order of expansion.")))
1968 ('t
(setq sw
't
) (return 't
)))))))
1970 (defmvar $taylor_order_coefficients t
1971 "When `true', coefficients of taylor series will be ordered canonically.")
1973 (defun taylor1 (e tlist
)
1974 (declare (special *within-srf?
* ))
1975 (setq tlist
(tlist-merge (nconc (find-tlists e
) tlist
)))
1976 (prog ($zerobern $simp $algebraic genpairs varlist tvars sing-tvars
1977 log-1 log%i ivars key-vars ans full-log last-exp
1978 mainvar-datum zerolist taylor_simplifier least_term? tvar-limits
1980 (setq tlist
(mapcan #'(lambda (d)
1981 (if (tvar?
(datum-var d
))
1987 (setq $zerobern t $simp t $algebraic t last-exp e least_term?
't
1988 log-1
'((%log simp
) -
1) log%i
'((%log simp
) $%i
)
1989 tvars
(mapcar 'car tlist
) varlist
(copy-list tvars
))
1990 (when $taylor_simplifier
1991 ; This symbolp/fboundp check is presumably for efficiency (so it
1992 ; can be directly funcalled).
1993 (setq taylor_simplifier
1994 (if (and (symbolp $taylor_simplifier
)
1995 (fboundp $taylor_simplifier
))
1997 'taylor_simplifier_caller
)))
1998 ;; Ensure that the expansion points don't depend on the expansion vars.
1999 ;; This could cause an infinite loop, e.g. taylor(x,x,x,1).
2000 (do ((tl tlist
(cdr tl
)))
2002 (unless (mfree (exp-pt (car tl
)) tvars
)
2003 (merror (intl:gettext
"taylor: attempt to expand ~M at a point depending on ~M.") e
(caar tl
))))
2004 ;; This drastic initialization ensures that ALGEBRAIC, TELLRAT, DISREP,
2005 ;; etc. prop's are removed from our gensyms. RATSETUP does not appear
2006 ;; to do this correctly, e.g. see ASB's bug of 1/10/83 (MQUERY 17).
2007 (mapc #'(lambda (g) (setf (symbol-plist g
) nil
)) genvar
)
2008 (ratsetup varlist genvar
)
2009 (when (and $taylor_order_coefficients
(not *within-srf?
*)) (newvar e
))
2010 (orderpointer varlist
)
2011 (maplist #'(lambda (q g
)
2012 (push (cons (car g
) (car q
)) key-vars
)
2013 (let ((data (get-datum (car q
))))
2014 (rplaca q
(transform-tvar (car q
) data
))
2015 (push (cons (car g
) (car q
)) ivars
)
2016 ;(setf (data-gvar-o data)
2017 ; (cons (car g) (valget (car g))))
2018 (rplacd (cdddr data
)
2019 (cons (car g
) (valget (car g
))))))
2020 (do ((v varlist
(cdr v
)))
2021 ((eq (car v
) (car tvars
)) v
))
2022 (do ((v varlist
(cdr v
))
2024 ((eq (car v
) (car tvars
)) g
)))
2025 (setq genpairs
(mapcar #'(lambda (y z
)
2026 (putprop z y
'disrep
)
2027 (cons y
(cons (pget z
) 1)))
2030 (setup-multivar-disrep () )
2031 (setq mainvar-datum
(car (last tlist
)))
2032 (mapc #'(lambda (d) (adjoin-sing-datum d
)) sing-tvars
)
2033 (setq ans
(catch 'tay-err
(taylor3 e
)))
2035 (if (atom (car ans
)) (tay-error (car ans
) (cadr ans
)) ans
))))
2037 (defun transform-tvar (var data
)
2038 (if (not (tvar? var
)) var
2039 (cond ((and (signp e
(exp-pt data
)) (null (switch '$asymp data
)))
2041 ((member (exp-pt data
) '($inf infinity
) :test
#'eq
)
2043 ((eq (exp-pt data
) '$minf
)
2045 ((let ((temp (m- var
(exp-pt data
))))
2046 (if (switch '$asymp data
) (m^ temp -
1) temp
))))))
2048 (defun taylor_simplifier_caller (exp)
2049 (mcall $taylor_simplifier exp
))
2051 (defun taylor_simplify_recurse (ps)
2052 (if (pscoefp ps
) (taylor2 (funcall taylor_simplifier
(rcdisrep ps
)))
2053 (let ((datum (ps-data ps
)) (var () ))
2054 ;; We must treat multivars like 1, since they'll reappear again
2055 ;; when we call taylor2 on their disrep'd coeff's.
2056 (if (switch 'multivar datum
)
2059 (setq var
(getdisrep (gvar-o ps
)))
2060 ;; Don't push pw's < 0, else constants will be truncated
2061 (push-pw datum
(emax (trunc-lvl ps
) (rczero)))))
2062 (do ((terms (terms ps
) (n-term terms
))
2063 (ans (rczero) (psplus (if (null datum
)
2064 (taylor_simplify_recurse (lc terms
))
2065 (pstimes (taylor_simplify_recurse
2068 ;; (taylor2 (funcall taylor_simplifier
2069 ;; (m^ var (edisrep (le terms)))))
2070 ;; causes terms to be lost when inverting. E.g.
2071 ;; taylor(log(1+exp(-1/x)),x,0,5) calls psexpt(<exp(1/x)^3>...3,-1)
2072 ;; which must return a series good to 3+3(-1-1)=-3 which, when added
2073 ;; to other terms will truncate them to degree -3 also.
2074 (if (ezerop (le terms
)) (rcone)
2077 (term (le terms
) (rcone)))))))
2081 (when datum
(pop-pw datum
))
2084 (defun push-pw (datum pw
)
2085 (push pw
(trunc-stack datum
))
2086 ;; When changing the truncation on an internal multivar we must also
2087 ;; propagate the change to all var's which depend upon it. See WGD's
2088 ;; bug report of 9/15/82 which exhibits the necessity of this.
2089 (when (setq datum
(switch 'multivar datum
))
2090 (do ((vars (cadr datum
) (cdr vars
)))
2092 (push pw
(trunc-stack (get-datum (car vars
)))))))
2094 (defun pop-pw (datum)
2095 (pop (trunc-stack datum
))
2096 ;; See the comment above in push-pw; here we must undo the propagation.
2097 (when (setq datum
(switch 'multivar datum
))
2098 (do ((vars (cadr datum
) (cdr vars
)))
2100 (pop (trunc-stack (get-datum (car vars
)))))))
2102 (defun setup-multivar-disrep (mrat?
)
2103 (let ((varlist varlist
) (genvar genvar
) (multivars () ))
2105 (setq varlist
(mrat-varlist mrat?
) genvar
(mrat-genvar mrat?
)))
2106 (mapc #'(lambda (datum)
2107 (when (switch 'multivar datum
)
2108 (push (car datum
) multivars
)
2109 ;; All genvar's corresponding to internally generated
2110 ;; multivars must "disappear" when disrep'd. If this
2111 ;; were not done then disrep'ing gx*gt would give x*t
2112 ;; which, upon, re-tayloring would give (gx*gt)*gt,
2113 ;; where t is the internal multivar for x, and gt, gx
2114 ;; are their genvars. An example where this would occur is
2115 ;; taylor(sin(x+y),[x],0,f1,[y],0,1).
2116 (putprop (int-gvar datum
) 1 'disrep
)))
2117 (if mrat?
(mrat-tlist mrat?
) tlist
))
2118 ;; Here we must substitute 1 for any genvars which depend on multivars.
2119 ;; For example, taylor(x^n,[x],0,0) generates a multivar^n.
2121 (do ((expl varlist
(cdr expl
))
2122 (gvarl genvar
(cdr gvarl
)))
2124 (unless (mfree (car expl
) multivars
)
2125 (putprop (car gvarl
) 1 'disrep
))))))
2127 ;; An example showing why this flag is need is given by
2128 ;; taylor(-exp(exp(-1/x)+2/x),x,0,-1). Without it, tstimes and
2129 ;; taylor_simplify_recurse end up trunc'ing the -1.
2131 (defvar trunc-constants?
't
)
2134 (cond ((mbagp e
) (cons (car e
) (mapcar #'taylor3
(cdr e
))))
2135 ((and (null tlist
) (not (eq exact-poly
'user-specified
)))
2137 (list 'mrat
'simp varlist genvar
)))
2138 (t (xcons (if (null taylor_simplifier
)
2141 (setq e
(taylor2 e
))
2142 (let ((exact-poly () ) (trunc-constants?
() ))
2143 (taylor_simplify_recurse e
))))
2144 (list 'mrat
'simp varlist genvar tlist
'trunc
)))))
2146 (defun find-tlists (e) (let (*a
*) (findtl1 e
) *a
*))
2149 (cond ((or (atom e
) (mnump e
)) )
2150 ((eq (caar e
) 'mrat
)
2151 (when (member 'trunc
(car e
) :test
#'eq
)
2152 (push (mapcar #'copy-tree
(mrat-tlist e
)) *a
*)))
2153 (t (mapc #'findtl1
(cdr e
)))))
2155 (defun tlist-merge (tlists)
2156 (do ((tlists tlists
(cdr tlists
)) (tlist () ))
2157 ((null tlists
) (nreverse tlist
))
2158 (do ((a_tlist (car tlists
) (cdr a_tlist
)) (temp nil
))
2160 (if (null (setq temp
(get-datum (datum-var (car a_tlist
)) t
)))
2161 (if (prog2 (setq temp
(car a_tlist
))
2162 (or (tvar?
(datum-var temp
))
2163 (member (caar (datum-var temp
)) '(mexpt %log
) :test
#'eq
)))
2164 (push (list (datum-var temp
) (trunc-stack temp
)
2165 (exp-pt temp
) (switches temp
))
2167 (merror (intl:gettext
"taylor: ~M cannot be a variable.") (datum-var temp
)))
2170 ;; We must take the max truncation level when $maxtayorder
2171 ;; is T, cf. JPG's bug of 9/15/82.
2172 (when (e> (current-trunc (car a_tlist
)) (current-trunc temp
))
2173 (setf (current-trunc temp
) (current-trunc (car a_tlist
))))
2174 (unless (e> (current-trunc (car a_tlist
)) (current-trunc temp
))
2175 (setf (current-trunc temp
) (current-trunc (car a_tlist
)))))
2176 (unless (alike1 (exp-pt temp
) (exp-pt (car a_tlist
)))
2177 (merror (intl:gettext
"taylor: cannot combine expressions expanded at different points.")))
2178 (setf (switches temp
)
2179 (union* (switches temp
) (switches (car a_tlist
)))))))))
2181 (defun compattlist (list)
2182 (do ((l list
(cdr l
)))
2184 (or (alike1 (exp-pt (get-datum (datum-var (car l
)))) (exp-pt (car l
)))
2188 (let ((last-exp e
)) ;; lexp-non0 should be bound here when needed
2189 (cond ((assolike e tlist
) (var-expand e
1 () ))
2190 ((or (mnump e
) (atom e
) (mfree e tvars
))
2191 (if (or (e> (rczero) (current-trunc mainvar-datum
))
2193 (pszero (data-gvar-o mainvar-datum
)
2194 (current-trunc mainvar-datum
))
2195 (if (and taylor_simplifier
(not (atom e
)))
2196 (let ((e-simp (prep1 (funcall taylor_simplifier e
))))
2197 (when (and (rczerop e-simp
) (not (member e-simp zerolist
:test
#'eq
)))
2201 ((null (atom (caar e
))) (merror "TAYLOR2: internal error."))
2203 (if (and (compatvarlist varlist
(mrat-varlist e
)
2204 genvar
(mrat-genvar e
))
2205 (compattlist (mrat-tlist e
)))
2207 (let ((exact-poly () )) (re-taylor e
))))
2208 ((eq (caar e
) 'mplus
) (tsplus (cdr e
)))
2209 ((eq (caar e
) 'mtimes
) (tstimes (cdr e
)))
2210 ((eq (caar e
) 'mexpt
) (tsexpt (cadr e
) (caddr e
)))
2211 ((eq (caar e
) '%log
) (tslog (cadr e
)))
2212 ((and (or (known-ps (caar e
)) (get (caar e
) 'tay-trans
))
2213 (not (member 'array
(cdar e
) :test
#'eq
))
2214 (try-expansion (if (cddr e
) (cdr e
) (cadr e
))
2217 (cond ((get (subfunname e
) 'spec-trans
)
2218 (funcall (get (subfunname e
) 'spec-trans
) e
))
2219 ((known-ps (subfunname e
))
2220 (try-expansion (caddr e
) (cadr e
))))) )
2221 ((and (member (caar e
) '(%sum %product
) :test
#'eq
)
2222 (mfreel (cddr e
) tvars
))
2223 (tsprsum (cadr e
) (cddr e
) (caar e
)))
2224 ;; If TSDIFF returns NIL, it means that the derivative is wrt to some variable
2225 ;; other than one of the taylor variables.
2226 ;; If so, keep going and handle E as if it were a general expression.
2227 ((and (eq (caar e
) '%derivative
) (tsdiff (cadr e
) (cddr e
) e
)))
2228 ((or (eq (caar e
) '%at
)
2229 (do ((l (mapcar 'car tlist
) (cdr l
)))
2231 (or (free e
(car l
)) (return ()))))
2233 (t (let ((exact-poly () )) ; Taylor series aren't exact
2234 (let ((ee (diff-expand e tlist
)))
2235 (cond ((equal (sratsimp e
) (sratsimp ee
)) (return-from taylor2
(prep1 ee
)))
2236 (t (taylor2 ee
)))))))))
2238 (defun compatvarlist (a b c d
)
2240 ((or (null b
) (null c
) (null d
)) () )
2241 ((alike1 (car a
) (car b
))
2242 (if (not (eq (car c
) (car d
))) ()
2243 (compatvarlist (cdr a
) (cdr b
) (cdr c
) (cdr d
))))
2244 (t (compatvarlist a
(cdr b
) c
(cdr d
)))))
2247 (defun re-taylor (mrat)
2248 (let ((old-tlist (mrat-tlist mrat
)) (old-varlist (mrat-varlist mrat
))
2249 (old-genvar (mrat-genvar mrat
)) old-ivars
)
2250 (declare (special old-tlist old-ivars
))
2251 ;; Put back the old disrpes so rcdisrep's will work correctly.
2252 (mapc #'(lambda (g v
) (putprop g v
'disrep
)) old-genvar old-varlist
)
2253 (setup-multivar-disrep mrat
)
2254 (setq old-ivars
(mapcar #'(lambda (g v
) (cons g v
))
2255 old-genvar old-varlist
))
2256 (prog1 (re-taylor-recurse (mrat-ps mrat
))
2257 ;; Restore the correct disreps.
2258 (mapc #'(lambda (g v
) (putprop g v
'disrep
)) genvar varlist
)
2259 (setup-multivar-disrep () ))))
2261 (defun re-taylor-recurse (ps)
2262 (declare (special old-tlist old-ivars
))
2263 (if (not (psp ps
)) (taylor2 (rcdisrep ps
))
2264 (let (var (datum () ))
2265 (setq var
(cdr (assoc (gvar ps
) old-ivars
:test
#'eq
)))
2266 ;; We must treat multivars like 1, since they'll reappear again
2267 ;; when we call taylor2 or var-expand below.
2268 (if (switch 'multivar
(assoc var old-tlist
:test
#'equal
))
2270 (when (setq datum
(var-data var
))
2271 (push-pw datum
(trunc-lvl ps
))))
2273 (do ((terms (terms ps
) (n-term terms
))
2275 (psplus (if (null var
) (re-taylor-recurse (lc terms
))
2276 (pstimes (re-taylor-recurse (lc terms
))
2278 (var-expand (car datum
)
2279 (edisrep (le terms
))
2282 (m^t var
(edisrep (le terms
)))))))
2285 (when datum
(pop-pw datum
))))))
2287 (defun var-expand (var exp dont-truncate?
)
2288 (let (($keepfloat
) ($float
) (modulus))
2289 (setq exp
(prep1 exp
))) ;; exp must be a rational integer
2290 (let ((temp (get-datum var
't
)))
2291 (cond ((null temp
) (merror "VAR-EXPAND: invalid call."))
2292 ((member (exp-pt temp
) '($inf $minf $infinity
) :test
#'eq
)
2293 (cond ((switch '$asymp temp
)
2294 (merror (intl:gettext
"taylor: cannot create an asymptotic expansion at infinity.")))
2295 ((e> (setq exp
(rcminus exp
)) (current-trunc temp
))
2297 (t (make-ps (int-var temp
)
2298 (ncons (if exact-poly
(inf) (current-trunc temp
)))
2300 (if (eq (exp-pt temp
) '$minf
)
2303 ;; multivar expansion does not work at infinity, so
2304 ;; expansion at infinity is handled by above clause even if doing multivar.
2305 ((switch 'multi temp
) ;; multivar expansion
2307 ;; The reason we call var-expand below instead of taylor2
2308 ;; is that we must be sure the call is not truncated to
2309 ;; 0 which would cause an error in psexpt if exp < 0.
2310 ;; For example, this occurred in TAYLOR(X^2/Y,[X,Y],0,2).
2312 ;; Must ensure that we get back a series truncated
2313 ;; to at least what is specified by tlist. This means
2314 ;; we'll have to push-pw unless exp>0 since psexpt'n
2315 ;; kills (exp-1) terms. The bug that discovered this
2316 ;; is taylor(li[2](x+1/2)/x,[x],0,0) missing 2*log(2).
2317 (if (not (e> exp
(rczero)))
2318 (let-pw (get-datum (car (switch 'multi temp
)))
2319 (e+ (current-trunc temp
) (e- (e1- exp
)))
2320 (var-expand (car (switch 'multi temp
)) 1 't
))
2321 (var-expand (car (switch 'multi temp
)) 1 't
))
2322 (cons (list (int-gvar temp
) 1 1) 1))
2323 (taylor2 (exp-pt temp
)))
2325 ((signp e
(exp-pt temp
))
2326 (let ((exp>trunc?
() ))
2327 (if (and (e> exp
(current-trunc temp
)) (setq exp
>trunc?
't
)
2328 (not dont-truncate?
))
2330 (make-ps (int-var temp
)
2331 (ncons (if exact-poly
(inf)
2332 (if exp
>trunc? exp
(current-trunc temp
))))
2333 (ncons (term (if (switch '$asymp temp
) (rcminus exp
)
2337 (make-ps (int-var temp
)
2338 (ncons (if exact-poly
(inf) (current-trunc temp
)))
2339 (ncons (term (if (switch '$asymp temp
)
2343 (taylor2 (exp-pt temp
)))
2346 (defun expand (arg func
)
2347 (or (try-expansion arg func
) (exp-pt-err)))
2349 (defun try-expansion (arg func
)
2350 (prog (funame funord fun-lc argord psarg arg-trunc temp exact-poly
)
2351 ;; We bind exact-poly to () since we don't want psexpt retaining
2352 ;; higher order terms when subst'ing into series (which aren't exact).
2353 ;; Try diff-expanding unknown subsripted functions.
2354 (unless (or (atom func
) (known-ps (caar func
)))
2355 (taylor2 (diff-expand `((mqapply) ,func
,arg
) tlist
)))
2356 (when (setq temp
(get (setq funame
(oper-name func
)) 'tay-trans
))
2357 (return (funcall temp arg func
)))
2358 (let ((lterm (getfun-lt func
)))
2359 (setq funord
(e lterm
) fun-lc
(c lterm
)))
2361 (when (rczerop (or psarg
(setq psarg
(get-lexp arg
(rcone) () ))))
2362 (if (e> (rczero) funord
)
2363 (if (rczerop (setq psarg
(get-lexp arg
(rcone) 't
)))
2365 (go begin-expansion
))
2366 (return (cond ((setq temp
(assoc funame tay-pole-expand
:test
#'eq
))
2367 (funcall (cdr temp
) arg psarg func
))
2368 ((rczerop funord
) fun-lc
)
2370 (when (pscoefp psarg
) (setq psarg
(taylor2 arg
)))
2371 (when (pscoefp psarg
)
2373 (cond ((null (mfree (rdis psarg
) tvars
))
2374 (symbolic-expand arg psarg func
))
2375 ((setq temp
(assoc funame tay-pole-expand
:test
#'eq
))
2376 (funcall (cdr temp
) arg psarg func
))
2378 (if (atom func
) `((,func
) ,(rcdisrep psarg
))
2379 `((mqapply) ,func
,(rcdisrep psarg
)))))))))
2380 (when (e> (rczero) (setq argord
(ps-le psarg
)))
2381 (cond ((not (member funame
'(%atan %asin %asinh %atanh
) :test
#'eq
))
2382 (if (e> (rczero) (ps-le* (setq psarg
(get-lexp arg
(rcone) 't
))))
2384 (go begin-expansion
)))
2386 (if (and (eq funame
'%atan
)
2387 (eq (coef-sign arg
) '$neg
))
2388 (return (psplus (atrigh arg func
) (taylor2 (m- '$%pi
))))
2389 (return (atrigh arg func
))))))
2390 (setq temp
(t-o-var (gvar psarg
)))
2391 (when (e> (e* funord argord
) temp
) (return (rczero)))
2392 ;; the following form need not be executed if psarg is really exact.
2393 ;; The constant problem does not allow one to determine this now,
2394 ;; so we always have to execute this currently.
2395 ;; This really should be
2396 ;; (unless (infp (trunc-lvl psarg)) ... )
2397 ;; Likewise, the infp checks shouldn't be there; have to assume
2398 ;; nothing is exact until constant problem is fixed.
2399 (setq arg-trunc
(if (and (not (infp (trunc-lvl psarg
)))
2400 (e= funord
(rcone)))
2402 (e- temp
(e* (e1- funord
) argord
)))
2403 psarg
(let-pw (get-datum (get-key-var (gvar psarg
)))
2405 (if (or (infp (trunc-lvl psarg
))
2406 (e> arg-trunc
(trunc-lvl psarg
)))
2409 ;; We must recalculate argord since pstrunc may have "picked"
2410 ;; a coeff out of a constant monomial; e.g. this occurs in
2411 ;; taylor(sin(x+y),x,0,0,y,0,1) where psarg is (Y+...)*X^0+...
2412 ;; which truncates to Y+... of order 1.
2413 argord
(ps-le* psarg
))
2414 (if (rczerop argord
)
2415 (cond ((member funame
'(%atan %asin %asinh %atanh
) :test
#'eq
)
2416 (return (atrigh arg func
)))
2417 ((setq temp
(assoc funame const-exp-funs
:test
#'eq
))
2418 (return (funcall (cdr temp
) arg psarg func
)))
2419 ((rczerop (ps-le* (setq psarg
(get-lexp arg
(rcone) 't
))))
2420 (return () )) ; Don't know an addition formula
2421 (t (go begin-expansion
)))
2423 (if (mono-term?
(terms psarg
))
2424 (get-series func
(current-trunc
2425 (get-datum (get-key-var (gvar psarg
))))
2426 (gvar-o psarg
) (ps-le psarg
) (ps-lc psarg
))
2428 (setq temp
(get-series func
2429 (e// temp argord
) (gvar-o psarg
)
2431 (cond ((not (psp temp
)) temp
)
2432 (t (pscsubst1 psarg temp
)))))))))
2434 (defun symbolic-expand (ign psarg func
) ; should be much stronger
2435 (declare (ignore ign
))
2436 (prep1 (simplifya (if (atom func
)
2437 `((,func
) ,(rcdisrep psarg
))
2438 `((mqapply) ,func
,(rcdisrep psarg
)))
2441 (defun expand-sing-trig?
(arg func
)
2442 (cond ((member func
*pscirc
:test
#'eq
) (tay-exponentialize arg func
))
2443 ((member func
*psacirc
:test
#'eq
) (atrigh arg func
))
2444 (t (essen-sing-err))))
2446 (defun trig-const (a arg func
)
2447 (let ((const (ps-lc* arg
)) (temp (cdr (assoc func trigdisp
:test
#'eq
))))
2448 (cond ((and (pscoefp const
)
2449 (member func
'(%tan %cot
) :test
#'eq
)
2450 (multiple-%pi a
(srdis const
) func
)))
2451 (temp (funcall temp
(setq const
(psdisrep const
))
2453 (t (tsexpt `((,(get func
'recip
)) ,(srdis arg
)) -
1)))))
2455 (defun multiple-%pi
(a const func
)
2457 (and (equal ($hipow const
'$%pi
) 1)
2458 ($ratnump
(setq coef
($ratcoef const
'$%pi
1)))
2459 (cond ((numberp coef
) (expand (m- a const
) func
))
2460 ((equal (caddr coef
) 2)
2461 (psminus (expand (m- a const
)
2462 (cond ((eq func
'%tan
) '%cot
)
2463 ((eq func
'%cot
) '%tan
)
2464 (t (merror "MULTIPLE-%PI: internal error in Taylor expansion."))))))))))
2466 (setq *pscirc
'(%cot %tan %csc %sin %sec %cos %coth
2467 %tanh %csch %sinh %sech %cosh
)
2469 *psacirc
'(%acot %atan %acsc %asin %asec %acos %acoth
2470 %atanh %acsch %asinh %asech %acosh
))
2472 (setq const-exp-funs
2473 `((%gamma . gam-const
) ($psi . plygam-const
)
2474 .
,(mapcar #'(lambda (q) (cons q
'trig-const
)) *pscirc
))
2476 trigdisp
'((%sin . psina
+b
) (%cos . pscosa
+b
) (%tan . pstana
+b
)
2477 (%sinh . psinha
+b
) (%cosh . pscosha
+b
) (%tanh . pstanha
+b
))
2479 tay-pole-expand
'((%gamma . plygam-pole
) ($psi . plygam-pole
))
2481 tay-const-expand
; !these should be handled by symbolic-expand
2482 ; be sure to change this with tay-exponentialize!
2483 (append (mapcar #'(lambda (q) (cons q
'tay-exponentialize
)) *pscirc
)
2484 (mapcar #'(lambda (q) (cons q
'tay-exponentialize
)) *psacirc
)))
2486 (mapc #'(lambda (q) (putprop q
'atrig-trans
'tay-trans
))
2487 '(%acos %acot %asec %acsc %acosh %acoth %asech %acsch
))
2489 (defprop mfactorial factorial-trans tay-trans
)
2491 (defun factorial-trans (arg func
)
2492 (declare (ignore func
))
2493 (taylor2 `((%gamma
) ,(m1+ arg
))))
2495 (defprop %signum signum-trans tay-trans
)
2497 ;; signum(x) => x/abs(x)
2498 ;; this gives an error at x=0 and a derivative of 0 elsewhere
2499 (defun signum-trans (arg func
)
2500 (declare (ignore func
))
2501 (taylor2 `((mtimes) ,arg
((mexpt) ((mabs) ,arg
) -
1))))
2503 (defprop %gamma_incomplete gamma-upper-trans tay-trans
)
2504 (defprop $gamma_incomplete gamma-upper-trans tay-trans
)
2505 (defprop %gamma_incomplete_lower gamma-lower-trans tay-trans
)
2506 (defprop $gamma_incomplete_lower gamma-lower-trans tay-trans
)
2508 ;; for gamma_incomplete(s,z)
2509 ;; translate into gamma_incomplete_lower if s>0 and z=0
2511 ;; June 2022: To workaround the bug
2512 ;; integrate(x*exp(-x^2)*sin(x),x,minf,inf)
2513 ;; limit: variable must be a symbol or subscripted symbol; found: sin(x)
2514 ;; I (Barton Willis) surrounded the call to $limit with errcatch with
2515 ;; $errormsg set to nil. This change allows Maxima to find the correct
2516 ;; value of this definite integral. But almost surely there is a bug
2517 ;; somewhere else that calls gamma-upper-trans with faulty arguments.
2518 ;; The real bug should be fixed, but inserting errcatch here is
2520 (defun gamma-upper-trans (arg func
)
2524 (eq ($sign s
) '$pos
)
2525 (let (($errormsg nil
))
2526 (zerop1 (car (errcatch ($limit z
(caar tlist
) (exp-pt (car tlist
))))))))
2527 (taylor2 `((mplus) ((%gamma
) ,s
)
2528 ((mtimes) -
1 ((%gamma_incomplete_lower
) ,s
,z
))))
2529 (taylor2 (diff-expand `((,func
) .
,arg
)
2532 ;; for gamma_incomplete_lower(s,z)
2533 ;; if z=0, use A&S 6.5.29
2537 ;;; gamma_incomplete_lower(s,z) = z^s * > ------------
2541 (defun gamma-lower-trans (arg func
)
2544 (if (zerop1 ($limit z
(caar tlist
) (exp-pt (car tlist
))))
2549 ((mexpt) ((mtimes) -
1 ,z
) k
)
2550 ((mexpt) ((mtimes) ((mfactorial) k
)
2556 (taylor2 (diff-expand `((,func
) .
,arg
)
2560 ;;; Not done properly yet
2562 ;;; (defprop $BETA BETA-TRANS TAY-TRANS)
2564 (defun psina+b
(a b
)
2565 (psplus (pstimes (expand a
'%sin
) (expand b
'%cos
))
2566 (pstimes (expand a
'%cos
) (expand b
'%sin
))))
2568 (defun pscosa+b
(a b
)
2569 (psdiff (pstimes (expand a
'%cos
) (expand b
'%cos
))
2570 (pstimes (expand a
'%sin
) (expand b
'%sin
))))
2572 (defun pstana+b
(a b
)
2573 (setq a
(expand a
'%tan
) b
(expand b
'%tan
))
2574 (pstimes (psplus a b
)
2575 (psexpt (psdiff (rcone) (pstimes a b
))
2578 (defun psinha+b
(a b
)
2579 (psplus (pstimes (expand a
'%sinh
) (expand b
'%cosh
))
2580 (pstimes (expand a
'%cosh
) (expand b
'%sinh
))))
2582 (defun pscosha+b
(a b
)
2583 (psplus (pstimes (expand a
'%cosh
) (expand b
'%cosh
))
2584 (pstimes (expand a
'%sinh
) (expand b
'%sinh
))))
2586 (defun pstanha+b
(a b
)
2587 (setq a
(expand a
'%tanh
) b
(expand b
'%tanh
))
2588 (pstimes (psplus a b
)
2589 (psexpt (psplus (rcone) (pstimes a b
))
2592 (defun atrig-trans (arg func
)
2594 (cond ((eq func
'%acos
)
2595 `((mplus) ,half%pi
((mtimes) -
1 ((%asin
) ,arg
))))
2598 `((mtimes) -
1 $%i
((mplus) ,half%pi
((mtimes) -
1 ((%asin
) ,arg
)))))
2601 `((,(cdr (assoc func
'((%acsc . %asin
) (%asec . %acos
)
2602 (%acot . %atan
) (%acsch . %asinh
)
2603 (%asech . %acosh
) (%acoth . %atanh
)) :test
#'eq
)))
2606 (defun atrigh (arg func
)
2607 (let ((full-log t
) ($logarc t
) (log-1 '((mtimes) $%i $%pi
))
2608 (log%i
'((mtimes) ((rat) 1 2) $%i $%pi
)))
2609 (taylor2 (simplify `((,func
) ,arg
)))))
2611 (defun tay-exponentialize (arg fun
) ; !this should be in symbolic-expand!
2612 (let (($exponentialize t
) ($logarc t
))
2613 (setq arg
(meval `((,fun
) ,arg
))))
2617 (do ((l (cdr l
) (cdr l
))
2618 (ans (taylor2 (car l
))
2619 (psplus ans
(taylor2 (car l
)))))
2622 (defun ts-formula (form var pw
)
2623 (let ((datum (get-datum (get-key-var (car var
)))))
2625 (taylor2 (subst (get-inverse (car var
)) 'sp2var form
)))))
2627 (defmacro next-series
(l) `(cdadr ,l
))
2629 (defun tstimes-get-pw (l pw
)
2630 (do ((l l
(cdr l
)) (vect))
2632 (setq pw
(mapcar #'(lambda (pwq ple
) (e+ pwq ple
))
2633 pw
(setq vect
(ord-vector (cdar l
)))))
2634 (rplacd (car l
) (cons (cdar l
) vect
))))
2636 (defun tstimes-l-mult (a)
2637 (do ((l (cdr a
) (cdr l
)) ($maxtayorder t
)
2638 (a (car a
) (pstimes a
(car l
))))
2644 (or (zfree e
(car l
)) (return () ))))
2646 ;;; The lists posl, negl and zerl have the following format:
2648 ;;; ( (<expression> <expansion> <order vector>) . . . )
2651 (*bind
* ((funl) (expl) (negl) (zerl) (posl)
2652 (pw) (negfl) (temp) (fixl (rcone)))
2653 (dolist (fun l
) ;; find the exponentials
2655 (push (if (free (caddr fun
) (car tvars
)) fun
2656 `((mexpt) $%e
,(m* (caddr fun
)
2657 `((%log
) ,(cadr fun
)))))
2661 (setq expl
(tsexp-comb expl
)) ;; simplify exps
2662 (setq expl
(tsbase-comb expl
))) ;; and again
2663 (setq l
(nconc expl funl
)) ;; now try expanding
2664 (let ((trunc-constants?
() ))
2665 (setq expl
(cons 0 (mapcar #'(lambda (exp)
2666 (cons exp
(taylor2 exp
)))
2668 ;; EXPL is now of the form (0 ( <form> . <taylor2(form)> ) ...)
2669 ;; L points behind the cons considered for destructive updating.
2670 (do ((l expl
) (tem))
2672 (cond ((rczerop (cdadr l
))
2673 ;; Consider taylor((a+1/x)*1/x,x,0,-2). Each factor will be on
2674 ;; zerl. Each factor will also appear to have le = 0 since its
2675 ;; series is 0, which would fool the get-pw routines below if
2676 ;; they tried to handle this case. The easiest fix for now
2677 ;; appears to be to always call get-lexp here, killing this:
2678 (cond ;((null $maxtayorder)
2679 ; (setq zerl (cons (cadr l) zerl))
2680 ; (rplacd l (cddr l)))
2681 ((rczerop (setq tem
(get-lexp (caadr l
) (rcone) ())))
2682 (return (setq zerl
0)))
2683 ('t
(setq posl
(cons (cons (caadr l
) tem
) posl
))
2684 (rplacd l
(cddr l
)))))
2685 ((pscoefp (cdadr l
))
2686 (cond ((mzfree (caadr l
) tvars
) ;must be zfree to permit ratfuns
2687 (setq fixl
(pstimes (cdadr l
) fixl
))
2688 (rplacd l
(cddr l
)))
2689 ((setq zerl
(cons (cadr l
) zerl
))
2690 (rplacd l
(cddr l
)))))
2691 ((rczerop (ps-le (cdadr l
)))
2692 (setq zerl
(cons (cadr l
) zerl
))
2693 (rplacd l
(cddr l
)))
2694 ((e> (ps-le (cdadr l
)) (rczero))
2695 (setq posl
(cons (cadr l
) posl
))
2696 (rplacd l
(cddr l
)))
2697 ('t
(setq l
(cdr l
)))))
2698 (when (equal zerl
0) (return (rczero)))
2699 (setq negl
(cdr expl
) temp
(ord-vector fixl
))
2700 (mapcar #'(lambda (x) (and (e> (rczero) x
) (setq negfl t
))) temp
)
2701 (tstimes-get-pw zerl temp
)
2702 (setq pw
(tstimes-get-pw posl
(tstimes-get-pw negl temp
)))
2705 (mapcar #'(lambda (x)
2706 (prog2 (mapcar #'(lambda (datum lel pwl
)
2708 (e+ (current-trunc datum
)
2712 (mapcar #'(lambda (datum) (pop-pw datum
))
2714 (nconc posl zerl negl
)))
2715 (setq posl
(nconc (mapcar 'cadr posl
) (mapcar 'cadr zerl
)
2716 (mapcar 'cadr negl
))))
2717 (setq posl
(tstimes-l-mult posl
))
2718 (let ((ans (cond ((null fixl
) posl
)
2720 ('t
(pstimes fixl posl
)))))
2721 (if $maxtayorder ans
(pstrunc ans
)))))
2723 ;;; This routine transforms a list of exponentials as follows:
2725 ;;; a^c*b^(n*c) ===> (a*b^n)^c, where n is free of var.
2727 ;;; This transformation is only applicable when c is not free of var.
2729 (defun tsexp-comb (l) ;; ***** clobbers l *****
2730 (setq l
(cons '* l
))
2731 (do ((a l
) (e)) ;; updated by a rplacd or cdr.
2732 ((null (cddr a
)) (cdr l
)) ;; get rid of the *
2733 (rplaca (cddadr a
) (setq e
($ratsimp
(caddr (cadr a
)))))
2734 ;; Must delete e^0 lest we divide by the 0 below. RWG's byzero bug
2735 ;; of 3/1/78 used to cause this.
2736 (if (equal e
0) (rplacd a
(cddr a
))
2737 (if (mfree (caddr (cadr a
)) tvars
) (pop a
)
2738 (do ((b (cddr a
) (cdr b
)) (n))
2739 ((null b
) (setq a
(cdr a
)))
2740 (when (mfree (setq n
($ratsimp
(m// (caddar b
)
2743 (rplaca b
(list '(mexpt simp
)
2745 (m^
(cadar b
) n
)) ;; b^n
2747 (rplacd a
(cddr a
)) ;; delete a^c
2750 ;;; This routine transforms a list of exponentials as follows:
2752 ;;; a^b*a^c ===> a^(b+c),
2754 ;;; this is only necessary when b and c depend on "var."
2756 (defun tsbase-comb (l) ;;; *******clobbers l********
2757 (setq l
(cons '* l
))
2758 (do ((a l
)) ;;; updated by a rplacd or cdr
2759 ((null (cddr a
)) (cdr l
))
2760 (do ((b (cddr a
) (cdr b
)))
2761 ((null b
) (pop a
)) ;;; did not return early so pop.
2762 (when (alike1 (cadar b
) (cadadr a
))
2763 (rplaca b
(m^
(cadar b
) (m+ (caddar b
) (caddr (cadr a
)))))
2768 (cond ((and (atom b
) (mnump e
)
2770 (not (eq (exp-pt (get-datum b
)) '$minf
)))
2771 ;; one could remove this clause and let this case be handled by tsexpt1
2772 (var-expand b e
() ))
2773 ((mfree e tvars
) (tsexpt1 b e
))
2774 ((eq '$%e b
) (tsexpt-red (list e
)))
2775 (t (tsexpt-red (list (list '(%log
) b
) e
)))))
2777 (defun tsexpt-red (l)
2778 (*bind
* ((free) (nfree) (full-log) ($logarc t
) (expt) (ps)
2779 (log-1 '((mtimes) $%i $%pi
))
2780 (log%i
'((mtimes) ((rat) 1 2) $%i $%pi
)))
2781 a
(do ((l l
(cdr l
)))
2783 (cond ((mtimesp (car l
)) (setq l
(append l
(cdar l
))))
2784 ((mfree (car l
) tvars
) (push (car l
) free
))
2785 (t (push (car l
) nfree
))))
2786 (cond ((or (cdr nfree
) (atom (car nfree
))) )
2787 ((eq (caaar nfree
) '%log
)
2788 (return (tsexpt1 (cadar nfree
) (m*l free
))))
2789 ((member (caaar nfree
) *psacirc
:test
#'eq
)
2790 (setq l
(ncons (simplifya ;; simplify after removing simp flag
2791 (cons (ncons (caaar nfree
)) (cdar nfree
))
2795 ;; Must have truncs > 0 so that logs in the expt aren't trunc'd.
2796 ;; E.g, consider taylor(x^(x-1),x,0,-1).
2797 (tlist-mapc d
(push-pw d
(emax (current-trunc d
) (rcone))))
2798 (setq ps
(taylor2 (setq expt
(m*l
(append nfree free
)))))
2799 (tlist-mapc d
(pop-pw d
))
2800 ;; Here we must account for the truncation gain or lossage that
2801 ;; is encountered in exp(c*log(x)+y) -> x^c*exp(y).
2802 (let ((c0 (if (pscoefp ps
) ps
(psterm (terms ps
) (rczero))))
2804 (unless (rczerop c0
)
2805 (setq ord-e^c0
(ord-vector (setq e^c0
(psexpt-fn c0
))))
2806 ;; Must emax with 0 so that new singular kernels won't be trunc'd
2807 ;; e.g exp(1/x+...) to degree -2 should be exp(-1/x)+...
2808 ;; Also try taylor(screwa,x,0,-2).
2809 (mapc #'(lambda (d o
) (push-pw d
(emax (e- (current-trunc d
) o
)
2812 (setq ps
(psdiff (taylor2 expt
) c0
)))
2813 (setq ps
(psexpt-fn ps
))
2815 (tlist-mapc d
(pop-pw d
))
2816 (setq ps
(pstimes e^c0 ps
)))
2819 ;; Taylor's b^e, where e is independent of tvars.
2821 (defun tsexpt1 (b e
)
2823 (setq e
(let ((modulus () )) ; Don't mod exponents! See WGM's bug
2824 (prog2 (mapcar ; of 3/6/83 for an example.
2827 (emax (current-trunc datum
) (rczero))))
2830 (mapcar #'(lambda (datum) (pop-pw datum
)) tlist
)))
2833 pw
(if (psp tb
) (current-trunc (get-datum
2834 (get-key-var (gvar tb
))))
2835 ;; Constant problem kludge.
2836 (if (rczerop tb
) (current-trunc (car tlist
)) (rczero))))
2837 (if (floatp (car s
))
2838 (setq s
(maxima-rationalize (quot (car s
) (cdr s
)))))
2839 ;; We must ensure that the lc is non-zero since it will be inverted in
2841 (setq tb
(strip-zeroes tb
't
))
2843 (when (or ;; When 1 > s we need more terms since -le*(s-1) > 0.
2845 (and (e> (rczero) pw
) (e> s
(rcone))))
2846 (setq tb
(get-lexp b
() 't
)))
2847 (setq le
(ps-le* tb
)))
2848 ((psp tb
) (setq le
(ps-le tb
)))
2849 (t (return (rcexpt tb e
))))
2850 (and (e> (e* s le
) pw
) (null $maxtayorder
) (return (rczero)))
2851 (setq s
(e- pw
(e* le
(e1- s
))))
2852 ;(setq le (increment-truncs tb))
2857 (pstrunc1 tb
(list (cons (gvar tb
) s
))))
2858 ;; because of constants not retaining info, have to
2859 ;; just keep the constant here
2860 (cond ((not (psp tb
)) tb
)
2861 (t (let-pw (get-datum (get-key-var (gvar tb
))) s
(strip-zeroes (taylor2 b
) 't
)))))
2864 ;;; the method of calculating truncation levels below is incorrect.
2865 ;;; (i.e. increment-truncs & decrement-truncs, also used above)
2866 ;;; Examples which exhibit this incorrectness are:
2867 ;;; taylor(log(sin(y)+x),x,0,2,y,0,1) is missing a y/6*x and -1/6*x^2
2868 ;;; taylor(log(sin(z)+sin(y)+x),x,0,f1,y,0,3,z,0,5) misses a z^5*y^3 term.
2870 ;;; TSLOG must find the lowest degree term in the expansion of the
2871 ;;; log arg, then expand with the orders of all var's in this low term
2872 ;;; incremented by their order in this low term. Note that this is
2873 ;;; only necessary for var's with ord > 0, since otherwise we have
2874 ;;; already expanded to a higher ord than required. Also we must
2875 ;;; not do this for var's with trunc < 0, since this may incorrectly
2876 ;;; truncate terms which should end up as logs.
2878 (defun increment-truncs (ps)
2879 (do ((ps ps
(ps-lc ps
)) (trunc (t-o-var (gvar ps
))) (data () ))
2881 (when (e> (ps-le ps
) (rczero))
2882 (push (assoc (get-key-var (gvar ps
)) tlist
:test
#'eq
) data
)
2883 (push-pw (car data
) (e+ (e* (e+ trunc
(rctwo)) (ps-le ps
))
2884 (current-trunc (car data
))))
2885 (setq trunc
(e+ trunc
(current-trunc (car data
))))
2888 (defun decrement-truncs (data)
2889 (mapc #'(lambda (data) (pop-pw data
)) data
))
2892 (let ((psarg (taylor2 arg
)) datum
)
2893 (when (rczerop psarg
) (setq psarg
(get-lexp arg
() 't
)))
2894 ;; We must ensure that the lc is non-zero since it will be inverted in pslog
2895 (setq psarg
(strip-zeroes psarg
't
))
2896 (do ((ps psarg
(ps-lc ps
)) (shift (rcone) (e* shift
(rctwo))))
2899 (when (rczerop (setq psarg
(taylor2 arg
)))
2900 (setq psarg
(get-lexp arg
() 't
)))
2901 (mapc #'(lambda (data) (pop-pw data
)) datum
))
2903 (push (get-datum (get-key-var (gvar ps
))) datum
)
2904 (if (and (e> (ps-le ps
) (rczero))
2905 (e> (current-trunc (car datum
)) (rczero)))
2906 (push-pw (car datum
) (e+ (e* shift
(ps-le ps
))
2907 (current-trunc (car datum
))))
2910 ;; When e-start is non-null we start expanding at order e-start, ... , 2^m,
2911 ;; then 2^m*pow, instead of the normal sequence pow, ... , 2^m*pow
2912 ;; (where m = $taylordepth, pow = ord of var). This is done because it is
2913 ;; usually much more efficient for large, non-trivial expansions when we only
2914 ;; want the lowest order term.
2916 (defun get-lexp (exp e-start zerocheck?
)
2922 (tlist-mapc d
(push-pw d
(or e-start
(emax (orig-trunc d
) (rcone)))))
2923 (do ((psexp) (i (1+ $taylordepth
) (1- i
)))
2925 (tlist-mapc d
(pop-pw d
))
2929 (unless silent-taylor-flag
(zero-warn exp
))
2931 (declare (fixnum i
))
2932 (cond ((and (rczerop (setq psexp
(if zerocheck?
2933 (strip-zeroes (taylor2 exp
) 't
)
2935 (not (member exp zerolist
:test
#'eq
))) )
2936 ;; Info not needed yet.
2937 ;; ((and lexp-non0 (rczerop (le (terms psexp)))
2938 ;; (mono-term? (terms psexp))))
2939 (t (tlist-mapc d
(pop-pw d
))
2941 (cond ((and (= i
1) e-start
)
2942 (setq e-start
() i
2)
2943 (tlist-mapc d
(push-pw d
(prog1 (e* (orig-trunc d
) (current-trunc d
))
2945 (t (tlist-mapc d
(push-pw d
(prog1 (e* (rctwo) (current-trunc d
))
2949 (or (equal x
1) (equal x
1.0)))
2951 (defun [max-trunc
] ()
2952 (do ((l tlist
(cdr l
)) (emax (rczero)))
2953 ((null l
) (1+ (truncate (car emax
) (cdr emax
))))
2954 (when (e> (current-trunc (car l
)) emax
)
2955 (setq emax
(orig-trunc (car l
))))))
2957 (defun tsprsum (f l type
)
2958 (if (mfree f tvars
) (newsym f
)
2959 (let ((li (ncons (car l
))) (hi (caddr l
)) (lv (ncons (cadr l
))) a aa
2960 ($maxtayorder
() ));; needed to determine when terms are 0
2961 (if (and (numberp (car lv
)) (numberp hi
) (> (car lv
) hi
))
2962 (if (eq type
'%sum
) (taylor2 0) (taylor2 1))
2964 (if (eq type
'%sum
) (setq type
'()))
2965 (do ((m (* ([max-trunc
]) (ash 1 $taylordepth
)))
2967 (ans (taylor2 (maxima-substitute (car lv
) (car li
) f
))))
2968 ((equal hi
(car lv
)) ans
)
2969 (rplaca lv
(m1+ (car lv
)))
2970 ;; A cheap heuristic to catch infinite recursion when
2971 ;; possible, should be improved in the future
2972 (if (> k m
) (exp-pt-err)
2973 (setq a
;(mlet li lv (taylor2 (setq aa (meval f))))
2974 (taylor2 (maxima-substitute (car lv
) (car li
) f
))))
2976 (if (and (1p (car a
)) (1p (cdr a
)) (not (1p aa
)))
2978 (setq ans
(pstimes a ans
)))
2979 (if (and (rczerop a
) (not (signp e aa
)))
2981 (setq ans
(psplus ans a
))))))))))
2983 (defun tsdiff (e l check
)
2984 (*bind
* ((n) (v) (u))
2985 (do ((l l
(cddr l
)))
2987 (if (and (atom (car l
)) (numberp (cadr l
))
2988 (assoc (car l
) tlist
:test
#'eq
))
2989 (setq n
(cons (cadr l
) n
) v
(cons (car l
) v
))
2990 (setq u
(cons (car l
) (cons (cadr l
) u
)))))
2991 ;; If N is null, it means that the derivative is wrt to some variable
2992 ;; other than one of the taylor variables.
2993 ;; If so, give up; let the caller handle it.
2995 (if u
(setq e
(meval (cons '($diff
) (cons e l
)))))
2996 (setq l
(mapcar #'(lambda (x) (get-datum x
)) v
))
2997 (mapcar #'(lambda (datum pw
)
2998 (push-pw datum
(e+ (current-trunc datum
) (prep1 pw
))))
3000 (setq e
(taylor2 e
))
3001 (mapc #'(lambda (datum) (pop-pw datum
)) l
)
3002 (do ((vl v
(cdr vl
))
3007 (mapc #'(lambda (a b
)
3008 (putprop a
(prep1 (sdiff b
(car v
)))
3011 (setq e
(psdp e
))))))
3014 (defun no-sing-err (x) ;; try to catch all singularities
3016 (let ((ans (catch 'errorsw
(eval x
))))
3017 (if (eq ans t
) (unfam-sing-err) ans
))))
3019 ;; evaluate deriv at location var=pt
3020 (defun eval-deriv (deriv var pt
)
3022 (let ((ans (no-sing-err `(meval '(($at
) ,deriv
((mequal) ,var
,pt
))))))
3025 (defun check-inf-sing (pt-list) ; don't know behavior of random fun's @ inf
3026 (and (or (member '$inf pt-list
:test
#'eq
) (member '$minf pt-list
:test
#'eq
))
3029 (defun diff-expand (exp l
) ;l is tlist
3030 (check-inf-sing (mapcar (function caddr
) l
))
3033 (setq exp
(diff-expand exp
(cdr l
)))
3034 (do ((deriv (sdiff exp
(caar l
)) (sdiff deriv var
))
3036 (coef 1 (* coef
(1+ cnt
)))
3038 (pt (exp-pt (car l
)))
3039 (lim (rcdisrep (current-trunc (car l
))))
3040 (ans (list (no-sing-err `(meval '(($at
) ,exp
((mequal) ,(caar l
) ,(exp-pt (car l
)))))))
3041 (cons `((mtimes) ((rat simp
) 1 ,coef
)
3042 ,(eval-deriv deriv var pt
)
3043 ((mexpt) ,(sub* var pt
) ,cnt
))
3045 ((or (great cnt lim
) (equal deriv
0)) (cons '(mplus) ans
))))))
3047 ;;; subtitle disreping routines
3050 (if (= (cdr e
) 1) (car e
) (list '(rat) (car e
) (cdr e
))))
3052 (defun striptimes (a)
3053 (if (mtimesp a
) (cdr a
) (ncons a
)))
3056 (let (($psexpand
() )) ; Called only internally, no need to expand.
3058 (cons (list 'mrat
'simp varlist genvar tlist
'trunc
)
3062 (let ((varlist (mrat-varlist r
)) (genvar (mrat-genvar r
)))
3063 (mapc #'(lambda (exp genv
) (putprop genv exp
'disrep
))
3065 (setup-multivar-disrep r
)
3066 ;; This used to return 0 if psdisrep returned () but this is wrong
3067 ;; since taylor(false,x,0,0) would lose. If psdisrep really wants to
3068 ;; return () for 0 then we will probably find out soon.
3069 (if (eq $psexpand
'$multi
) (psdisexpand (cdr r
))
3070 (psdisrep (cdr r
)))))
3074 (psdisrep+ (psdisrep2 (terms p
) (getdisrep (gvar-o p
)) (trunc-lvl p
))
3075 (if (or $psexpand
(trunc-lvl p
)) '(mplus trunc
)
3079 (defun psdisrep^
(n var
)
3080 ;; If var = () then it is an internal var generated in a multivariate
3081 ;; expansion so it shouldn't be displayed. If var = 1 then it probably
3082 ;; resulted from the substitution in srdisrep, so it depends on an
3083 ;; internal var and likewise shouldn't be displayed.
3084 (cond ((or (rczerop n
) (null var
) (equal var
1)) 1)
3085 ((equal n
(rcone)) var
)
3086 ((and (mexptp var
) (equal (caddr var
) -
1))
3087 ;; This used to be conditioned on ps-bmt-disrep which was
3088 ;; always true. But I (rtoy) can't find a case where it
3089 ;; would have made a difference. The testsuite doesn't have
3091 (psdisrep^
(e- n
) (cadr var
)))
3092 ('t
`((mexpt ratsimp
) ,var
,(edisrep n
)))))
3094 ;;; There used to be a hack below that would print a series consisting
3095 ;;; of merely one term as exact polynomial (i.e. no trailing "..."'s).
3096 ;;; This is, of course, wrong but the problem with the fix is that
3097 ;;; now exact things like taylor(y*x,x,0,f1,y,0,1) will display like
3098 ;;; (y+...) x+... because of the problem with $MAXTAYORDER being internally
3099 ;;; bound to ()---which causes exact things to look inexact, such as
3100 ;;; x and y above. See the comment above taylor* for the $MAXTAYORDER problem.
3102 (defun psdisrep+ (p plush
&aux lowest-degree-term
)
3103 (if;; An exact sum of one arg is just that arg.
3104 (and (null (cdr p
)) (eq (cadr plush
) 'exact
))
3107 ;; Since the DISPLAY package prints trunc'd sum's arguments
3108 ;; from right to left we must put the terms of any constant term
3109 ;; in decreasing order. Note that only a constant (wrt to the
3110 ;; mainvar) term can be a term which is a sum.
3111 (when (mplusp (setq lowest-degree-term
(car (last p
))))
3112 (rplacd lowest-degree-term
(nreverse (cdr lowest-degree-term
))))
3115 (defun psdisrep* (a b
)
3116 (cond ((equal a
1) b
)
3118 (t (cons '(mtimes ratsimp
)
3119 (nconc (striptimes a
) (striptimes b
))))))
3121 (defun psdisrep2 (p var trunc
)
3122 (if (or $ratexpand $psexpand
) (psdisrep2expand p var
)
3123 (do ((a () (cons (psdisrep* (psdisrep (lc p
)) (psdisrep^
(le p
) var
))
3126 ((or (null p
) (e> (le p
) trunc
)) a
))))
3128 (defun psdisrep2expand (p var
)
3130 (l () (nconc (psdisrep*expand
(psdisrep (lc p
)) (psdisrep^
(le p
) var
))
3134 (defun psdisrep*expand
(a b
)
3135 (cond ((equal a
1) (list b
))
3136 ((equal b
1) (list a
))
3138 (list (cons '(mtimes ratimes
) (nconc (striptimes a
) (striptimes b
)))))
3139 ('t
(mapcar #'(lambda (z) (psdisrep* z b
))
3143 (defun psdisexpand (p)
3144 (let ((ans (ncons ())))
3145 (declare (special ans
)) ;used in pans-add
3146 (psdisexcnt p
() (rczero))
3149 (mapcar #'(lambda (x) (cond ((not (cddr x
)) (cadr x
))
3150 (t (cons '(mplus trunc
) (cdr x
)))))
3152 (cond ((not (cdr ans
)) (car ans
))
3153 (t (cons '(mplus trunc
) ans
)))))
3155 (defun psdisexcnt (p l n
)
3157 (do ((var (getdisrep (gvar-o p
))) (ll (terms p
) (n-term ll
)))
3159 (if (rczerop (le ll
)) (psdisexcnt (lc ll
) l n
)
3161 (cons (psdisrep^
(le ll
) var
) l
)
3163 (psans-add (cond ((not l
) (rcdisrep p
))
3164 (t (psdisrep* (rcdisrep p
)
3165 (cond ((not (cdr l
)) (car l
))
3166 (t (cons '(mtimes trunc
) l
))))))
3169 (defun psans-add (exp n
)
3170 (declare (special ans
)) ;bound in psdisexpand
3171 (do ((l ans
(cdr l
)))
3172 ((cond ((null (cdr l
)) (rplacd l
(ncons (list n exp
))))
3173 ((e= (caadr l
) n
) (rplacd (cadr l
) (cons exp
(cdadr l
))))
3174 ((e> (caadr l
) n
) (rplacd l
(cons (list n exp
) (cdr l
))))))))
3176 (defun srconvert (r)
3177 (cond ((not (atom (caadr (cdddar r
))))
3178 (cons (car r
) (psdisextend (cdr r
))))
3180 (*bind
* ((trunclist (cadr (cdddar r
)))
3185 (gens (cadddr (car r
))))
3186 (setq gps
(mapcar #'cons gens vs
))
3187 (do ((tl (cdr trunclist
) (cddr tl
)))
3188 ((null tl
) (cons (list 'mrat
'simp vs gens tlist
'trunc
) (srconvert1 (cdr r
))))
3189 (setq temp
(cdr (assoc (car tl
) gps
:test
#'eq
)))
3190 (cond ((null (member (car tl
) (cdr trunclist
) :test
#'eq
)))
3191 ((mplusp temp
) (merror "SRCONVERT: internal error."))
3194 (cons (list* temp
(tay-order (cadr tl
)) 0 nil
3195 (cons (car tl
) (symbol-value (car tl
))))
3198 (defun srconvert1 (p)
3199 (cond ((not (member (car p
) genvar
:test
#'eq
)) p
)
3201 (do ((l (cdr p
) (cddr l
))
3202 (a nil
(cons (term (prep1 (car l
)) (srconvert1 (cadr l
))) a
)))
3204 (make-ps (cons (car p
) (symbol-value (car p
)))
3205 (tay-order (zl-get trunclist
(car p
))) a
))))))
3207 ;;; subtitle error handling
3209 (defun tay-error (msg exp
)
3210 (if silent-taylor-flag
(throw 'taylor-catch
())
3212 (merror "taylor: ~A~%~M" msg exp
)
3213 (merror "taylor: ~A" msg
))))
3215 (defun exp-pt-err ()
3216 (tay-err (intl:gettext
"unable to expand at a point specified in:")))
3218 (defun essen-sing-err ()
3219 (tay-err (intl:gettext
"encountered an essential singularity in:")))
3221 (defun unfam-sing-err ()
3222 (tay-err (intl:gettext
"encountered an unfamiliar singularity in:")))
3224 (defun infin-ord-err ()
3225 (tay-err (intl:gettext
"expansion to infinite order?")))
3227 (defun tay-depth-err ()
3228 (tay-err (intl:gettext
"'taylordepth' exceeded while expanding:")))
3230 ;;; Subtitle TAYLORINFO
3232 (defun taylor-trunc (q)
3233 (setq q
(current-trunc q
))
3234 (cond ((null q
) '$inf
)
3235 ((equal (cdr q
) 1) (car q
))
3236 (t `((rat) ,(car q
) ,(cdr q
)))))
3238 (defun taylor-info (q)
3239 (let ((acc-var nil
) (acc-pt nil
) (acc-ord nil
) (qk) (acc))
3240 (cond ((null q
) nil
)
3243 (cond ((and (fourth qk
) (consp (fourth qk
)) (eq (caar (fourth qk
)) 'multivar
)) nil
)
3244 ((and (fourth qk
) (consp (fourth qk
)) (eq (caar (fourth qk
)) 'multi
))
3245 (while (and (fourth qk
) (consp (fourth qk
)) (eq (caar (fourth qk
)) 'multi
))
3247 (push (taylor-trunc qk
) acc-ord
)
3248 (push (exp-pt qk
) acc-pt
)
3249 (push (datum-var qk
) acc-var
)
3251 (push '(mlist) acc-ord
)
3252 (push '(mlist) acc-pt
)
3253 (push '(mlist) acc-var
)
3254 (setq q
(taylor-info q
))
3255 (if (null q
) (list acc-var acc-pt acc-ord
) (append q
(list acc-var acc-pt acc-ord
))))
3258 (setq acc
(if (and (fourth qk
) (consp (fourth qk
)) (eq '$asympt
(caar (fourth qk
))))
3259 (list '$asympt
) nil
))
3260 (push (taylor-trunc qk
) acc
)
3261 (push (exp-pt qk
) acc
)
3262 (push (datum-var qk
) acc
)
3264 (setq q
(taylor-info q
))
3265 (if (null q
) (list acc
) (append q
(list acc
)))))))))
3267 (defmfun $taylorinfo
(x)
3268 (if (and (consp x
) (member 'trunc
(first x
) :test
#'eq
))
3269 (cons '(mlist) (taylor-info (mrat-tlist x
)))
3274 ;;; Lisp let-pw Indent:2
