1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; **************************************************************
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 varlist
;List of all the variables occurring in a power
119 ;series, the power series variables at the end
120 genvar
;The list of gensyms corresponding to varlist
122 *a
* ;Temporary special
123 silent-taylor-flag
;If true indicates that errors will be
124 ;returned via a throw to TAY-ERR
125 tlist
;An association list which contains the
126 ;relevant information for the expansion which
127 ;is passed in at toplevel invocation.
128 $float
;Indicates whether to convert rational numbers
129 ;to floating point numbers.
130 $keepfloat
;When true retains floatin point numbers
133 log-1
;What log(-1) should be log(-1) or pi*i.
134 log%i
;Similarly for log(i)
135 exact-poly
;Inicates whether polynomials are to be
136 ;considered exact or not. True within SRF,
137 ;false within TAYLOR.
139 half%pi
;Has pi/2 to save space.
141 tay-const-expand
;For rediculousness like csch(log(x))
142 $exponentialize
;which we do by exponentiation.
145 last-exp
;last-expression through taylor2
148 genpairs
;List of dotted pairs
150 ivars
;Pairlist if gensym and disreped version
151 key-vars
;Pairlist of gensym and key var (for searching
159 *within-srf?
* ;flag for in srf
161 least_term?
; If non-null then the addition routines
162 ; are adding or subtracting coeff's of the
163 ; least term of a sum so they should do
164 ; zero checking on it if it is desired.
165 taylor_simplifier
; This is set by taylor1 to the object
166 ; which will be funcalled whenever
167 ; coefficient simplification is desired.
168 zerolist
; A list of constant expressions which have
169 ; been verified to be zero by a call to
170 ; $TAYLOR_SIMPLIFIER in taylor2. It is used to
171 ; suppress the message that TAYLOR is assumming
172 ; an expression to be zero.
173 ; 0p-funord lexp-non0 ; referenced only in commented-out code, so comment out here too
175 ) ;Don't want to see closed compilation notes.
177 (defmvar $psexpand
()
178 "When TRUE extended rational function expressions will be displayed fully
179 expanded. (RATEXPAND will also cause this.) If FALSE, multivariate
180 expressions will be displayed just as in the rational function package.
181 If PSEXPAND:MULTI, then terms with the same total degree in the variables
182 are grouped together.")
184 (defmvar $maxtayorder t
185 "When true TAYLOR retains as many terms as are certain to be correct
186 during power series arithmetic. Otherwise, truncation is controlled
187 by the arguments specified to TAYLOR.")
189 (defmvar $taylor_truncate_polynomials t
190 "When FALSE polynomials input to TAYLOR are considered to have infinite
191 precison; otherwise (the default) they are truncated based upon the input
194 (defmvar $taylor_logexpand t
195 "Unless FALSE log's of products will be expanded fully in TAYLOR (the default)
196 to avoid identically-zero constant terms which involve log's. When FALSE,
197 only expansions necessary to produce a formal series will be executed.")
199 (defmvar $taylor_simplifier
'simplify
200 "A function of one argument which TAYLOR uses to simplify coefficients
203 (defvar taylor_simplifier nil
)
205 ;;; Subtitle General Macsyma Free Predicates
208 (cond ((equal e x
) () )
211 (null (member x
(cdr ($listofvars e
)) :test
#'equal
)))
212 ('t
(do ((l (cdr e
) (cdr l
))) ((null l
) 't
)
213 (or (zfree (car l
) x
) (return () ))))))
215 (defun mfree (exp varl
)
216 (declare (special dummy-variable-operators
))
217 (cond ((atom exp
) (not (member exp varl
:test
#'eq
)))
218 ((eq (caar exp
) 'mrat
)
219 (do ((l (mrat-varlist exp
) (cdr l
)))
221 (unless (mfree (car l
) varl
) (return () ))))
222 ((or (member (caar exp
) dummy-variable-operators
:test
#'eq
)
223 (member 'array
(cdar exp
) :test
#'eq
))
224 (do ((vars varl
(cdr vars
)))
226 (unless (freeof (car vars
) exp
) (return () ))))
227 ('t
(and (mfree (caar exp
) varl
) (mfreel (cdr exp
) varl
)))))
229 (defun mfreel (l varl
)
230 (or (null l
) (and (mfree (car l
) varl
) (mfreel (cdr l
) varl
))))
232 ;;; Subtitle Coefficient Arithmetic
235 (cond ((equal x
(rcone)) (rcone))
236 ((rczerop y
) (rcone))
237 ((and (equal (cdr y
) 1) (fixnump (car y
)))
239 ((and $radexpand
(numberp (car y
)) (numberp (cdr y
)))
241 (setq y
(maxima-rationalize (quot (car y
) (cdr y
)))))
242 (ratexpt (rcquo (rcexpt1 (car x
) (cdr y
))
243 (rcexpt1 (cdr x
) (cdr y
)))
246 (prep1 (m^
(rcdisrep x
) (rcdisrep y
)))))))
249 (cond ((equal p
1) (rcone))
250 ((pcoefp p
) (prep1 (m^
(pdis p
) (*red
1 n
))))
251 ;; psfr does a square-free decom on p yielding (p1 e1 p2 e2 ... pn en)
252 ;; where p = p1^e1 p2^e2 ... pn^en, the pi being square-free
253 (t (do ((l (psqfr p
) (cddr l
))
256 (if (not (equal (rem (cadr l
) n
) 0))
257 (setq ans
(rctimes ans
(prep1 (m^
(pdis (car l
))
258 (*red
(cadr l
) n
)))))
259 ;; If pi<0, n=2m and n|ei then ei=2e and
260 ;; (pi^ei)^(1/(2m)) = (-pi)^(e/m)
262 (when (and (evenp n
) (eq ($sign
(pdis (car l
))) '$neg
))
263 (rplaca l
(pminus (car l
))))
264 (setq ans
(rctimes ans
(ratexpt (cons (car l
) 1)
265 (truncate (cadr l
) n
))))))))))
267 (defun rccoefp (e) ;a sure check, but expensive
270 (member (caar e
) genvar
:test
#'eq
))
272 (member (cadr e
) genvar
:test
#'eq
))))
274 ;;; Subtitle Exponent arithmetic
277 (and (not (infp x
)) (signp e
(car x
))))
280 (cond ((or (infp x
) (infp y
)) (inf))
281 ((and (equal (cdr x
) 1) (equal (cdr y
) 1))
282 (cons (+ (car x
) (car y
)) 1))
283 (t (ereduce (+ (* (car x
) (cdr y
)) (* (cdr x
) (car y
)))
284 (* (cdr x
) (cdr y
))))))
287 (cond ((infp x
) (inf))
288 ((and (equal (cdr x
) 1) (equal (cdr y
) 1))
289 (cons (- (car x
) (car y
)) 1))
290 (t (ereduce (- (* (car x
) (cdr y
)) (* (cdr x
) (car y
)))
291 (* (cdr x
) (cdr y
))))))
296 ((equal (cdr x
) (cdr y
)) (cons (min (car x
) (car y
)) (cdr x
)))
297 ((< (* (car x
) (cdr y
)) (* (cdr x
) (car y
))) x
)
301 (cond ((or (infp x
) (infp y
)) (inf))
302 ((equal (cdr x
) (cdr y
)) (cons (max (car x
) (car y
)) (cdr x
)))
303 ((> (* (car x
) (cdr y
)) (* (cdr x
) (car y
))) x
)
307 (cond ((or (infp x
) (infp y
)) (inf))
308 ((and (equal (cdr x
) 1) (equal (cdr y
) 1))
309 (cons (* (car x
) (car y
)) 1))
310 (t (ereduce (* (car x
) (car y
)) (* (cdr x
) (cdr y
))))))
314 (cons (- (cdr e
)) (- (car e
)))
315 (cons (cdr e
) (car e
))))
318 (cond ((infp x
) (inf))
320 (t (ereduce (* (car x
) (cdr y
))
321 (* (cdr x
) (car y
))))))
324 (cond ((infp x
) (inf))
325 ((= (cdr x
) 1) (cons (1+ (car x
)) 1))
326 (t (cons (+ (cdr x
) (car x
)) (cdr x
)))))
329 (cond ((infp x
) (inf))
330 ((equal (cdr x
) 1) (cons (1- (car x
)) 1))
331 (t (cons (- (car x
) (cdr x
)) (cdr x
)))))
336 ((equal (cdr x
) (cdr y
)) (> (car x
) (car y
)))
337 (t (> (* (car x
) (cdr y
)) (* (car y
) (cdr x
))))))
341 ((or (null e1
) (null e2
)) ())
342 (t (and (equal (car e1
) (car e2
))
343 (equal (cdr e1
) (cdr e2
))))))
346 (if (signp l d
) (setq d
(- d
) n
(- n
)))
347 (if (zerop n
) (rczero)
348 (let ((gcd (gcd n d
)))
349 (cons (/ n gcd
) (/ d gcd
)))))
352 (let ((xn (abs (car x
))) (xd (cdr x
))
353 (yn (abs (car y
))) (yd (cdr y
)))
354 (cons (gcd xn yn
) (* xd
(/ yd
(gcd xd yd
))))))
356 ;;; Subtitle polynomial arithmetic
358 (declare-top (special vars
))
360 (defun ord-vector (p)
361 (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum
))) tlist
)))
362 (declare (special vars
))
363 (cond ((not (cdr vars
)) (ncons (ps-le* p
)))
364 (t (ord-vect1 p
) (mapcar #'(lambda (x) (or (cdr x
) (rczero))) vars
)))))
367 (declare (special vars
))
369 (let ((data (assoc (gvar p
) vars
:test
#'eq
))
371 (rplacd data
(cond ((not (cdr data
)) le
)
372 (t (emin (cdr data
) le
))))
373 (mapl #'(lambda (l) (ord-vect1 (lc l
))) (terms p
)))))
375 (defun trunc-vector (p min?
)
376 (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum
))) tlist
)))
377 (declare (special vars
))
378 (if (null (cdr vars
)) (ncons (if (psp p
) (trunc-lvl p
) () ))
381 (mapcar 'cdr vars
)))))
383 (defun trunc-vect1 (p min?
)
384 (declare (special vars
))
386 (let ((data (assoc (gvar p
) vars
:test
#'eq
))
387 (trunc (trunc-lvl p
)))
389 (rplacd data
(if (null (cdr data
)) trunc
390 (if min?
(emin (cdr data
) trunc
)
391 (emax (cdr data
) trunc
))))))
392 (dolist (term (terms p
))
393 (trunc-vect1 (c term
) min?
))))
395 (declare-top (unspecial vars
))
399 (cond ((pscoefp y
) (rcplus x y
))
402 ((pscoefp y
) (if (rczerop y
) x
(pscplus y x
)))
403 ((eqgvar (gvar-o x
) (gvar-o y
)) (psplus1 x y
))
404 ((pointerp (gvar-o x
) (gvar-o y
)) (pscplus y x
))
408 (if (not (and least_term? taylor_simplifier
)) (rcplus x y
)
409 (prep1 (funcall taylor_simplifier
(m+ (rcdisrep x
) (rcdisrep y
))))))
412 (cond ((pscoefp x
) (cond ((pscoefp y
) (rcdiff x y
))
413 ((rczerop x
) (pstimes (rcmone) y
))
414 (t (pscdiff x y
() ))))
415 ((pscoefp y
) (if (rczerop y
) x
(pscdiff y x t
)))
416 ((eqgvar (gvar-o x
) (gvar-o y
)) (psdiff1 x y
))
417 ((pointerp (gvar-o x
) (gvar-o y
)) (pscdiff y x t
))
418 (t (pscdiff x y
() ))))
421 (if (not (and least_term? taylor_simplifier
)) (rcdiff x y
)
422 (prep1 (funcall taylor_simplifier
(m- (rcdisrep x
) (rcdisrep y
))))))
425 (let ((ans (cons () () )))
426 (psplus2 (gvar-o x
) (emin (trunc-lvl x
) (trunc-lvl y
))
427 (cons 0 (terms x
)) (cons 0 (terms y
)) ans ans
)))
430 (if (e> (rczero) (trunc-lvl p
)) p
431 (pscheck (gvar-o p
) (poly-data p
) (pscplus1 c
(terms p
)))))
433 (defun pscdiff (c p fl
)
434 (if (e> (rczero) (trunc-lvl p
))
435 (if fl p
(psminus p
))
436 (pscheck (gvar-o p
) (poly-data p
)
437 (cond ((not fl
) (pscplus1 c
(psminus-terms (terms p
))))
438 (t (pscplus1 (psminus c
) (terms p
)))))))
440 (defun strip-zeroes (terms ps?
)
441 (cond ((or (null terms
) (null taylor_simplifier
)) terms
)
443 (do ((terms terms
(n-term terms
)))
445 (change-lc terms
(strip-zeroes (lc terms
) 't
))
446 (unless (rczerop (lc terms
)) (return terms
))))
448 (if (null taylor_simplifier
) terms
449 (let ((exp (rcdisrep terms
)))
450 ;; If a pscoeff is not free of tvars then the ps is a
451 ;; multivar series and we can't handle a recursive
452 ;; call to taylor (as opposed to a call to prep1, as below)
453 ;; because this would be circuler (e.g. try
454 ;; taylor(x/ (x^2+1),[x],%i,-1) ). Besides, in this case
455 ;; the pscoeff contains a tvar hence should not be 0.
456 (if (not (mfree exp tvars
)) terms
457 (prep1 (funcall taylor_simplifier exp
))))))
458 (t (pscheck (gvar-o terms
) (poly-data terms
)
459 (strip-zeroes (terms terms
) () )))))
461 (defun pscplus1 (c l
)
462 (cond ((null l
) (list (term (rczero) c
)))
463 ((rczerop (le l
)) (setq c
(psplus c
(lc l
)))
464 (if (rczerop c
) (strip-zeroes (n-term l
) () )
465 (cons (term (rczero) c
) (n-term l
))))
466 ((e> (le l
) (rczero)) (cons (term (rczero) c
) l
))
467 (t (cons (lt l
) (let ((least_term?
)) (pscplus1 c
(n-term l
)))))))
469 ;;; Both here and in psdiff2 xx and yy point one before where one
470 ;;; might think they should point so that extensions will be retained.
472 (defun psplus2 (varh trunc xx yy ans a
)
474 a
(cond ((mono-term? xx
)
475 (if (mono-term? yy
) (go end
) (go null
)))
476 ((mono-term? yy
) (setq yy xx
) (go null
)))
477 (cond ((equal (le (n-term xx
)) (le (n-term yy
)))
478 (setq xx
(n-term xx
) yy
(n-term yy
))
479 (setq c
(let ((least_term?
(null (n-term ans
))))
480 (psplus (lc xx
) (lc yy
))))
481 (if (rczerop c
) (go a
) (add-term a
(le xx
) c
)))
482 ((e> (le (n-term xx
)) (le (n-term yy
)))
483 (setq yy
(n-term yy
))
484 (add-term a
(lt yy
)))
485 (t (setq xx
(n-term xx
))
486 (add-term a
(lt xx
))))
489 null
(if (or (mono-term? yy
) (e> (le (n-term yy
)) trunc
))
492 (setq yy
(n-term yy
))
493 (add-term-&-pop a
(lt yy
))
495 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
498 (let ((ans (cons () () )))
499 (psdiff2 (gvar-o x
) (emin (trunc-lvl x
) (trunc-lvl y
))
500 (cons 0 (terms x
)) (cons 0 (terms y
)) ans ans
)))
502 (defun psdiff2 (varh trunc xx yy ans a
)
504 a
(cond ((mono-term? xx
)
509 (cons 0 (mapcar #'(lambda (q)
510 (term (e q
) (psminus (c q
))))
514 (setq yy xx
) (go null
)))
515 (cond ((equal (le (n-term xx
)) (le (n-term yy
)))
516 (setq xx
(n-term xx
) yy
(n-term yy
))
517 (setq c
(let ((least_term?
(null (n-term ans
))))
518 (psdiff (lc xx
) (lc yy
))))
519 (if (rczerop c
) (go a
)
520 (add-term a
(le xx
) c
)))
521 ((e> (le (n-term xx
)) (le (n-term yy
)))
522 (setq yy
(n-term yy
))
523 (add-term a
(le yy
) (psminus (lc yy
))))
524 (t (setq xx
(n-term xx
))
525 (add-term a
(lt xx
))))
528 null
(if (or (mono-term? yy
) (e> (le (n-term yy
)) trunc
))
531 (setq yy
(n-term yy
))
532 (add-term-&-pop a
(le yy
) (lc yy
))
534 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
537 (if (psp x
) (make-ps x
(psminus-terms (terms x
)))
540 (defun psminus-terms (terms)
541 (let ((ans (cons () () )))
542 (do ((q terms
(n-term q
))
545 (add-term a
(le q
) (psminus (lc q
))))))
547 (defun pscheck (a b terms
)
548 (cond ((null terms
) (rczero))
549 ((and (mono-term? terms
) (rczerop (le terms
)))
551 (t (make-ps a b terms
))))
553 (defun pstrim-terms (terms e
)
555 (cond ((null terms
) (return () ))
556 ((null (e> e
(le terms
))) (return terms
))
557 (t (setq terms
(n-term terms
))))))
559 (defun psterm (terms e
)
560 (psterm1 (pstrim-terms terms e
) e
))
563 (cond ((null l
) (rczero))
564 ((e= (le l
) e
) (lc l
))
567 (defun pscoeff1 (a b c
) ;a is an mrat!!!
568 (let ((tlist (mrat-tlist a
)))
569 (cons (nconc (list 'mrat
'simp
(mrat-varlist a
) (mrat-genvar a
))
570 (do ((l (mrat-tlist a
) (cdr l
))
571 (ans () (cons (car l
) ans
)))
573 (when (alike1 (caar l
) b
)
575 (and (or ans
(cdr l
))
576 (list (nreconc ans
(cdr l
)) 'trunc
))))))
577 (pscoef (mrat-ps a
) (int-gvar (get-datum b
)) (prep1 c
)))))
579 (defun pscoef (a b c
)
580 (cond ((pscoefp a
) (if (rczerop c
) a
(rczero)))
581 ((eq b
(gvar a
)) (psterm (terms a
) c
))
582 (t (do ((gvar-o (gvar-o a
))
583 (poly-data (poly-data a
))
585 (terms (terms a
) (n-term terms
))
588 (unless (rczerop (setq temp
(pscoef (lc terms
) b c
)))
589 (setq ans
(psplus ans
590 (make-ps gvar-o poly-data
591 (ncons (term (le terms
)
594 (defun psdisextend (p)
595 (cond ((not (psp p
)) p
)
596 (t (make-ps p
(mapcar #'(lambda (q) (cons (car q
) (psdisextend (cdr q
))))
600 (if (psp p
) (psfloat1 p
(trunc-lvl p
) (terms p
) (ncons 0))
601 (rctimes (rcfone) p
)))
603 (defun psfloat1 (p trunc l ans
)
605 (a (last ans
) (n-term a
)))
606 ((or (null l
) (e> (le l
) trunc
))
607 (pscheck (gvar-o p
) (poly-data p
) (cdr ans
)))
608 (add-term a
(le l
) (psfloat (lc l
)))
609 (setq l
(n-term l
))))
612 (pstrunc1 p
(mapcar #'(lambda (q) (cons (int-gvar q
) (current-trunc q
)))
615 (defun pstrunc1 (p trlist
)
619 (let ((trnc (cdr (assoc (gvar p
) trlist
:test
#'eq
))) (trunc-ps) (a nil
))
620 (do ((l (terms p
) (n-term l
)))
621 ((null l
) (pscheck (gvar-o p
) (ncons (trunc-lvl p
)) (nreverse a
)))
622 (when (e> (le l
) trnc
)
623 (return (pscheck (gvar-o p
) (ncons trnc
) (nreverse a
))))
624 (unless (rczerop (setq trunc-ps
(pstrunc1 (lc l
) trlist
)))
625 (push (term (le l
) trunc-ps
) a
)))))))
628 (cond ((or (rczerop x
) (rczerop y
)) (rczero))
629 ((pscoefp x
) (cond ((pscoefp y
) (rctimes x y
))
630 ((equal x
(rcone)) y
)
631 (t (psctimes* x y
))))
632 ((pscoefp y
) (if (equal y
(rcone)) x
(psctimes* y x
)))
633 ((eqgvar (gvar-o x
) (gvar-o y
)) (pstimes*1 x y
))
634 ((pointerp (gvar-o x
) (gvar-o y
)) (psctimes* y x
))
635 (t (psctimes* x y
))))
637 (defun psctimes* (c p
)
638 (make-ps p
(maplist #'(lambda (l)
639 (term (le l
) (pstimes c
(lc l
))))
642 (defun pstimes*1 (xa ya
)
643 (let ((ans (cons () () ))
644 (trunc (let ((lex (ps-le xa
)) (ley (ps-le ya
)))
645 (e+ (emin (e- (trunc-lvl xa
) lex
) (e- (trunc-lvl ya
) ley
))
648 (setq trunc
(emin trunc
(t-o-var (gvar xa
)))))
649 (pstimes*2 xa ya trunc ans
)))
651 (defun pstimes*2 (xa ya trunc ans
)
653 (setq x
(terms xa
) y
(setq yy
(terms ya
)) a ans
)
654 a
(cond ((or (null y
) (e> (setq e
(e+ (le x
) (le y
))) trunc
))
656 ((not (rczerop (setq c
(pstimes (lc x
) (lc y
)))))
657 (add-term-&-pop a e c
)))
660 b
(unless (setq x
(n-term x
))
661 (return (pscheck (gvar-o xa
) (list trunc
) (cdr ans
))))
663 c
(when (or (null y
) (e> (setq e
(e+ (le x
) (le y
))) trunc
))
665 (setq c
(pstimes (lc x
) (lc y
)))
666 d
(cond ((or (mono-term? a
) (e> (le (n-term a
)) e
))
667 (add-term-&-pop a e c
))
668 ((e> e
(le (n-term a
)))
671 (t (setq c
(psplus c
(lc (n-term a
))))
673 (rplacd a
(n-term (n-term a
)))
675 (change-lc (n-term a
) c
)
676 (setq a
(n-term a
))))))
680 (defun pscsubst (c v p
)
681 (cond ((pscoefp p
) p
)
682 ((eq v
(gvar p
)) (pscsubst1 c p
))
683 ((pointerp v
(gvar p
)) p
)
684 (t (make-ps p
(maplist
685 #'(lambda (q) (term (le q
)
686 (pscsubst c v
(lc q
))))
689 (defun pscsubst1 (v u
)
691 (ul (terms u
) (n-term ul
)))
693 (setq a
(psplus a
(pstimes (lc ul
) (psexpt v
(le ul
)))))))
695 (defun get-series (func trunc var e c
)
696 (let ((pw (e// trunc e
)))
697 (setq e
(if (and (equal e
(rcone)) (equal c
(rcone)))
698 (getexp-fun func var pw
)
699 (psmonsubst (getexp-fun func var pw
) trunc e c
)))
700 (if (and $float $keepfloat
) (psfloat e
) e
)))
702 (defun psmonsubst (p trunc e c
)
704 (psmonsubst1 p trunc e c
705 `(() .
,(terms p
)) (ncons () ) (rcone) (rczero))
709 (defun psmonsubst1 (p trunc e c l ans cc el
)
710 ;; We set $MAXTAYORDER to () here so that the calls to psexpt below
711 ;; won't do needless extra work, e.g. see rwg's complaint of 9/7/82.
712 (prog (a ee varh $maxtayorder
)
713 (setq a ans varh
(gvar-o p
))
714 a
(cond ((or (mono-term? l
)
715 (e> (setq ee
(e* e
(le (n-term l
)))) trunc
))
719 (psexpt c
(e- (le (setq l
(n-term l
)))
722 (add-term a ee
(pstimes cc
(lc l
)))))
723 (setq a
(n-term a
) el
(le l
))
725 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
727 (defun psexpon-gcd (terms)
728 (do ((gcd (le terms
) (egcd (le l
) gcd
))
729 (l (n-term terms
) (n-term l
)))
733 (if (psp p
) (psfind-s (psterm (terms p
) (rczero)))
737 (cond ((null (atom (cdr r
))) (rczero))
739 (t (do ((p (ptterm (cdar r
) 0) (ptterm (cdr p
) 0)))
740 ((atom p
) (cons p
(cdr r
)))))))
743 (cond ((rczerop n
) ;; p^0
744 (if (rczerop p
) ;; 0^0
745 (merror (intl:gettext
"taylor: 0^0 is undefined."))
746 (rcone))) ;; Otherwise can let p^0 = 1
747 ((or (equal n
(rcone)) (equal n
(rcfone))) p
) ;; p^1 cases
748 ((pscoefp p
) (rcexpt p n
))
749 ((mono-term?
(terms p
)) ;; A monomial to a power
750 (let ((s (psfind-s n
)) (n-s) (x) (l (terms p
)))
751 ;; s is the numeric part of the exponent
752 (if (floatp (car s
)) ;; Perhaps we souldn't
753 ;; rationalize if $keepfloat is true?
754 (setq s
(maxima-rationalize (quot (car s
) (cdr s
)))))
755 (setq n-s
(psdiff n s
) ;; the non-numeric part of exponent
756 x
(e* s
(le l
))) ;; the degree of the lowest term
757 (setq x
(if (and (null $maxtayorder
) ;; if not getting all terms
758 (e> x
(t-o-var (gvar p
))))
759 ;; and result is of high order
760 (rczero) ;; then zero is enough
761 (pscheck (gvar-o p
) ;; otherwise
762 (ncons (e+ (trunc-lvl p
) ;; new trunc-level
763 (e- x
(le l
)))) ;; kick exponent
764 (ncons (term x
(psexpt (lc l
) n
))))))
766 (if (or (rczerop n-s
) (rczerop x
)) ;; is that good enough?
767 x
;; yes! The rest is bletcherous.
768 (pstimes x
(psexpt (prep1 (m^
(get-inverse (gvar p
))
771 (t (prog (l lc le inc trunc s lt mr lim lcinv ans
)
772 (setq lc
(lc (setq l
(terms p
)))
773 le
(le l
) lt
(lt l
) trunc
(trunc-lvl p
)
774 inc
(psexpon-gcd l
) s
(psfind-s n
))
775 (when (floatp (car s
))
776 (setq s
(maxima-rationalize (quot (car s
) (cdr s
)))))
777 (setq ans
(psexpt (setq lt
(pscheck (gvar-o p
) (list trunc
)
779 lcinv
(psexpt lc
(rcmone))
780 mr
(e+ inc
(e* s le
))
781 lim
(if (and (infp trunc
) (not (e> s
(rczero))))
783 ;; See the comment in PSEXPT1 below which tells
784 ;; why we don't allow inf. trunc's here.
785 (e+ (if (and (infp trunc
) (not (rcintegerp s
)))
786 (if (infp (setq lim
(t-o-var (gvar p
))))
792 (if (or (pscoefp ans
) (null (eq (gvar p
) (gvar ans
))))
793 (list 0 (term (rczero) ans
))
794 (cons 0 (terms ans
))))
795 (and (null $maxtayorder
)
798 (e> (e* s
(le (last l
))) (t-o-var (gvar p
))))
799 (setq lim
(emin lim
(t-o-var (gvar p
)))))
800 ;;(and (infp lim) (n-term l) (e> (rczero) n)
802 (return (psexpt1 (gvar-o p
)
803 lim l n s inc
1 mr ans le lcinv
))))))
805 (defun psexpt1 (varh trunc l n s inc m mr ans r linv
)
806 ;; n is the power we are raising the series to
807 ;; inc is the exponent increment
808 ;; mr is the current exponent
809 ;; tr is the truncation level desired
810 ;; m is the term index
812 (prog (a (k 0) ak cm-k c ma0 sum kr tr
)
814 ;; truly unfortunate that we need so many variables in this hack
815 (setq a
(last ans
) tr trunc
)
816 ;; I don't see what's wrong with truncating exact series when
817 ;; raising them to fractional powers so we'll allow it for now.
818 ;; This is accomplished above in PSEXPT (see the comment). Thus,
819 ;; presumably, this check should never be needed anymore.
820 ;; Bugs catching this clause were sqrt(1-x)*taylor(f1,x,0,0)
821 ;; and sqrt(taylor(x+x^2,x,0,2)),taylor_truncate_polynomials=false.
824 (setq tr
(e* s
(le (last l
))))
825 (merror (intl:gettext
"taylor: expected an integer, instead found: ~:M") s
)))
826 (when (infp tr
) (setq tr
(t-o-var (car varh
))))
827 b
(and (e> mr tr
) (go end
))
828 (setq kr inc ak l ma0
(pstimes (cons 1 m
) linv
)
830 a
(if (or (> k m
) (null (setq cm-k
(psterm (cdr ans
) (e- mr kr
)))))
832 (setq ak
(or (pstrim-terms ak
(e+ kr r
)) (go add-term
))
833 c
(pstimes (psdiff (pstimes (cons k
1) n
)
835 (pstimes (if (e= (e+ kr r
) (le ak
))
839 (setq sum
(psplus sum c
)
840 k
(1+ k
) kr
(e+ kr inc
))
843 (and (null (rczerop sum
))
844 (add-term-&-pop a mr
(pstimes ma0 sum
)))
845 (setq m
(1+ m
) mr
(e+ mr inc
))
847 end
(return (pscheck varh
(list trunc
) (cdr ans
)))))
849 (defun psderivative (p v
)
850 (cond ((pscoefp p
) (rcderiv p v
))
852 (if (prog1 (rczerop (ps-le p
))
853 (setq p
(psderiv1 (gvar-o p
)
854 (trunc-lvl p
) (cons 0 (terms p
)) (list 0))))
855 (strip-zeroes p
't
) p
))
856 (t (psderiv2 (gvar-o p
)
857 (trunc-lvl p
) v
(cons 0 (terms p
)) (list 0)))))
859 (defun psderiv1 (varh trunc l ans
)
861 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
862 (pscheck varh
(list (e1- trunc
)) (cdr ans
)))
864 (when (not (rczerop (le l
)))
865 (add-term-&-pop a
(e1- (le l
)) (pstimes (le l
) (lc l
))))))
867 (defun psderiv2 (varh trunc v l ans
)
868 (do ((a (last ans
) (n-term a
)) (c))
869 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
870 (pscheck varh
(list trunc
) (cdr ans
)))
872 (or (rczerop (setq c
(psderivative (lc l
) v
)))
873 (add-term a
(le l
) c
))))
877 (cond ((pscoefp p
) (rcderivx p
))
878 ((or (rczerop (setq temp
(getdiff (gvar-o p
))))
879 (eq (car temp
) 'multi
))
880 (setq temp2
(psdp2 (gvar-o p
) (trunc-lvl p
)
881 (cons 0 (terms p
)) (list 0)))
882 (if (eq (car temp
) 'multi
)
884 (make-ps (gvar-o p
) (ncons (inf))
885 (list (term (cdr temp
) (rcone)))))
888 (trunc-lvl p
) (cons 0 (terms p
))
891 (defun psdp1 (varh trunc l ans dx
)
892 (do ((a (last ans
)) (c (rczero)))
893 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
894 (psplus c
(pscheck varh
(list (e1- trunc
)) (cdr ans
))))
896 (if (rczerop (le l
)) (setq c
(psdp (lc l
)))
898 a
(e1- (le l
)) (pstimes (le l
) (pstimes dx
(lc l
)))))))
900 (defun psdp2 (varh trunc l ans
)
901 (do ((a (last ans
)) (c))
902 ((or (mono-term? l
) (e> (le (n-term l
)) trunc
))
903 (pscheck varh
(list trunc
) (cdr ans
)))
905 (when (null (rczerop (setq c
(psdp (lc l
)))))
906 (add-term-&-pop a
(le l
) c
))))
910 ;;; (defun psintegrate (p v)
911 ;;; (cond ((rczerop p) (rczero))
913 ;;; (pstimes p (taylor2 (get-inverse (car v)))))
914 ;;; ((eqgvar v (gvar-o p))
915 ;;; (psinteg1 (gvar-o p)
916 ;;; (trunc-lvl p) (cons 0 (terms p)) (list 0)))
917 ;;; (t (psinteg2 (gvar-o p)
918 ;;; (trunc-lvl p) v (cons 0 (terms p)) (list 0)))))
920 ;;; (defun psinteg1 (varh trunc l ans)
922 ;;; (setq a (last ans))
923 ;;; a (if (or (null (n-term l)) (e> (le (n-term l)) trunc))
925 ;;; (add-term a (e1+ (le (setq l (n-term l))))
927 ;;; (if (e= (le l) (rcmone))
928 ;;; (prep1 (list '(%LOG)
932 ;;; (setq a (n-term a)))
934 ;;; end (return (pscheck varh (list (e1+ trunc)) (cdr ans)))))
936 ;;; (defun psinteg2 (varh trunc v l ans)
938 ;;; (setq a (last ans))
939 ;;; a (if (or (null (n-term l)) (e> (le (n-term l)) trunc))
941 ;;; (add-term a (le l)
942 ;;; (psintegrate (lc (setq l (n-term l))) v))
943 ;;; (setq a (n-term a)))
945 ;;; end (return (pscheck varh (list trunc) (cdr ans)))))
947 (defun psexpt-log-ord (p)
948 (cond ((null $maxtayorder
) (emin (trunc-lvl p
) (t-o-var (gvar p
))))
949 ((infp (trunc-lvl p
)) (t-o-var (gvar p
)))
958 (cond ((pscoefp p
) (psexpt-fn2 (rcdisrep p
)))
959 ((ps-lim-infp p
) (psexpt-fn-sing p
))
960 ((prog2 (setq ord
<0?
(e> (rczero) (ps-le p
)))
961 (null (n-term (terms p
))))
962 (setq ans
(get-series '%exp
(psexpt-log-ord p
) (gvar-o p
)
963 (if ord
<0?
(e- (ps-le p
)) (ps-le p
))
965 (if ord
<0?
(ps-invert-var ans
) ans
))
967 (when (e= (rczero) (e (setq ans
(ps-gt p
))))
968 (pstimes (psexpt-fn (pscheck (gvar-o p
) (list (trunc-lvl p
))
969 (delete ans
(terms p
) :test
#'eq
)))
970 (psexpt-fn2 (srdis (c ans
)))))
971 (when (e= (rczero) (ps-le p
))
972 (pstimes (psexpt-fn2 (srdis (lc (terms p
))))
973 (psexpt-fn (pscheck (gvar-o p
) (list (trunc-lvl p
))
974 (n-term (terms p
))))))) )
975 (t (prog (l inc trunc ea0 ans
)
978 ;(return (ps-invert-var (psexpt-fn (ps-invert-var p))))
979 (setq l
(invert-terms l
)))
980 (setq trunc
(trunc-lvl p
)
981 inc
(psexpon-gcd l
) ea0
(rcone))
982 (unless (e> (le l
) (rczero))
983 ;; MEANING OF FOLLOWING MESSAGE IS OBSCURE
984 (merror "PSEXPT-FN: unreachable point."))
986 (if (or (pscoefp ea0
) (null (eq (gvar p
) (gvar ea0
))))
987 (list 0 (term (rczero) ea0
))
988 (cons 0 (terms ea0
))))
990 (setq trunc
(emin trunc
(t-o-var (gvar p
)))))
991 (when (infp trunc
) (setq trunc
(t-o-var (gvar p
))))
992 (setq ans
(psexpt-fn1 (gvar-o p
) trunc l inc
1 inc ans
))
993 (return (if ord
<0?
(ps-invert-var ans
) ans
)))))))
995 (defun psexpt-fn-sing (p)
996 (let ((inf-var?
(member (gvar-lim (gvar p
)) '($inf $minf
) :test
#'eq
))
997 (c*logs
(c*logs
(lt-poly p
))) c strongest-term
)
998 ;; Must pull out out logs here: exp(ci*log(ui)+x) -> ui^ci*exp(x)
999 ;; since its much harder for adjoin-tvar to do this transformation
1000 ;; below after things have been disrepped.
1001 (setq c
(exp-c*logs c
*logs
) p
(psdiff p
(sum-c*logs c
*logs
)))
1002 (if (not (ps-lim-infp p
))
1003 ;; Here we just subtracted the only infinite term, e.g.
1004 ;; p = 1/2*log(x)+1/log(x)+...
1005 (pstimes c
(psexpt-fn p
))
1007 (setq strongest-term
(if inf-var?
(ps-gt p
) (ps-lt p
)))
1008 ;; If the strongest term has degree 0 in the mainvar then the singular
1009 ;; terms occur in some other weaker var. There may be terms in this
1010 ;; coef which arent singular (e.g. 1 in (1/x+1+...)+exp(-1/x)+...) so
1011 ;; we must recursively psexpt-fn this term to get only what we need.
1012 (if (rczerop (e strongest-term
))
1013 (setq c
(pstimes c
(psexpt-fn (c strongest-term
))))
1014 (dolist (exp (expand-and-disrep strongest-term p
))
1015 (setq c
(pstimes c
(adjoin-tvar (m^
'$%e exp
))))))
1016 (pstimes c
(psexpt-fn (pscheck (gvar-o p
) (list (trunc-lvl p
))
1018 (delete strongest-term
(terms p
) :test
#'eq
)
1019 (n-term (terms p
))))))))))
1021 (defun gvar-logp (gvar)
1022 (let ((var (get-inverse gvar
)))
1023 (and (consp var
) (eq (caar var
) 'mexpt
) (equal (caddr var
) -
1)
1024 (consp (setq var
(cadr var
))) (eq (caar var
) '%log
)
1029 (let ((log (gvar-logp (gvar p
))) c
)
1033 (setq c
(psconst (psterm (terms p
) (rcmone))))
1034 ;; We don't want e.g. exp(x^a*log(x)) -> x^x^a
1035 (if (not (mfree (rcdisrep c
) tvars
)) ()
1036 (cons (cons c
(cons log p
))
1037 (c*logs
(psterm (terms p
) (rczero))))))))))
1040 (if (pscoefp p
) p
(psconst (psterm (terms p
) (rczero)))))
1042 (defun exp-c*logs
(c*logs
)
1043 (if (null c
*logs
) (rcone)
1044 (pstimes (taylor2 `((mexpt) ,(cadr (cadr (car c
*logs
)))
1045 ,(rcdisrep (caar c
*logs
))))
1046 (exp-c*logs
(cdr c
*logs
)))))
1048 (defun sum-c*logs
(c*logs
)
1049 (if (null c
*logs
) (rczero)
1050 (let ((ps (cddr (car c
*logs
))))
1051 (psplus (make-ps ps
(ncons (term (ps-le ps
) (caar c
*logs
))))
1052 (sum-c*logs
(cdr c
*logs
))))))
1054 ;; Calculatest the limit of a series at the expansion point. Returns one of
1055 ;; {$zeroa, $zerob, $pos, $neg, $inf, $minf}.
1057 (defvar tvar-limits
()
1058 "A list of the form ((gvar . limit(gvar)) ...)")
1060 (defun ps-lim-infp (ps)
1062 ;; Assume taylor vars at 0+ for now. Should handle the cases when
1063 ;; the expansion point is INF, MINF,etc.
1064 (let* ((lim (gvar-lim (gvar ps
)))
1066 (if (member lim
'($inf $minf
) :test
#'eq
) (ps-gt ps
) (ps-lt ps
))))
1067 (if (ezerop (e strongest-term
))
1068 (ps-lim-infp (c strongest-term
))
1070 (setq lim
(lim-power lim
(e strongest-term
)))
1071 (and (lim-infp lim
) (not (eq lim
'$infinity
))))))))
1073 (defun lim-zerop (lim)
1074 (member lim
'($zeroa $zerob $zeroim
) :test
#'eq
))
1076 (defun lim-plusp (lim)
1077 (member lim
'($zeroa $pos $inf $finite
) :test
#'eq
))
1079 (defun lim-finitep (lim)
1080 (member lim
'($pos $neg $im $finite
) :test
#'eq
))
1082 (defun lim-infp (lim)
1083 (member lim
'($inf $minf $infinity
) :test
#'eq
))
1085 (defun lim-imagp (lim)
1086 (member lim
'($im $infinity
) :test
#'eq
))
1088 (defun lim-minus (lim)
1089 (cdr (assoc lim
'(($zeroa . $zerob
) ($zerob . $zeroa
) ($pos . $neg
) ($zero . $zero
)
1090 ($neg . $pos
) ($inf . $minf
) ($minf . $inf
)
1091 ($im . $im
) ($infinity . $infinity
) ($finite . $finite
)) :test
#'eq
)))
1092 (defun lim-abs (lim)
1093 (or (cdr (assoc lim
'(($zerob . $zeroa
) ($neg . $pos
) ($minf . $inf
)) :test
#'eq
))
1096 (defun lim-times (lim1 lim2
)
1098 (cond ((or (eq lim1
'$zero
) (eq lim2
'$zero
)) (setq lim
'$zero
))
1099 ((and (lim-infp lim1
) (lim-infp lim2
)) (setq lim
'$inf
))
1100 ((and (lim-zerop lim1
) (lim-zerop lim2
)) (setq lim
'$pos
))
1101 ((or (when (lim-finitep lim2
) (rotatef lim1 lim2
) 't
)
1103 (when (and (eq lim1
'$finite
) (lim-infp lim1
))
1104 (break "Undefined finite*inf in lim-times"))
1105 (setq lim
(lim-abs lim2
)))
1106 (t (break "Undefined limit product ~A * ~A in lim-times" lim1 lim2
)))
1107 (if (or (lim-imagp lim1
) (lim-imagp lim2
))
1108 (if (lim-infp lim
) '$infinity
'$im
)
1109 (if (and (lim-plusp lim1
) (lim-plusp lim2
)) lim
(lim-minus lim
)))))
1111 (defun lim-power (lim power
)
1112 (cond ((ezerop power
) '$pos
)
1113 ((e> (rczero) power
) (lim-recip (lim-power lim
(e- power
))))
1114 ((not (oddp (car power
)))
1115 (if (lim-plusp lim
) lim
(lim-minus lim
)))
1118 (defun lim-recip (lim)
1119 (or (cdr (assoc lim
'(($zeroa . $inf
) ($zerob . $minf
)
1120 ($inf . $zeroa
) ($minf . $zerob
)) :test
#'eq
))
1121 (if (eq lim
'$finite
) (break "inverting $finite?")
1124 (defun lim-exp (lim)
1126 (($zeroa $zerob $zero $pos $neg $minf
) '$zeroa
)
1127 (($inf $finite
) lim
)
1128 ($infinity
'$infinity
) ; actually only if Re lim = inf
1129 (t (break "Unhandled limit in lim-exp"))))
1131 (defun lim-log (lim)
1137 (t (break "Unhandled limit in lim-log"))))
1139 (defun expand-and-disrep (term p
)
1140 (let ((x^n
(list '(mexpt) (get-inverse (gvar p
)) (edisrep (e term
))))
1142 (if (pscoefp a
) (ncons (m* (srdis a
) x^n
))
1143 (mapcar #'(lambda (subterm)
1144 (m* (cons '(mtimes) (expand-and-disrep subterm a
)) x^n
))
1147 (defun adjoin-sing-datum (d)
1148 (let ((r (prep1 (datum-var d
))) (g (gensym)) (kernel (datum-var d
))
1149 (no (1+ (cdr (int-var (car (last tlist
)))))))
1150 (unless (and (equal (car r
) 1) (equal (cddr r
) '(1 1)))
1151 (break "bad singular datum"))
1152 (putprop g kernel
'disrep
)
1153 (rplacd (cdddr d
) (cons g no
))
1155 (push (cons (cadr r
) kernel
) key-vars
)
1156 (push (cons g kernel
) key-vars
)
1157 (push (car key-vars
) ivars
)
1158 ;(push (cons kernel (cons (pget g) 1)) genpairs)
1159 (push (cons g
(exp-pt d
)) tvar-limits
)))
1161 (defun adjoin-tvar (exp) (rat->ps
(prep1 exp
)))
1163 (defun rat->ps
(rat)
1164 (pstimes (poly->ps
(car rat
))
1165 (psexpt (poly->ps
(cdr rat
)) (rcmone))))
1167 (defun poly->ps
(poly)
1168 (if (or (pcoefp poly
) (mfree (pdis poly
) tvars
)) (prep1 poly
)
1169 (let ((g (p-var poly
)) datum
(pow (rcone)))
1170 (if (setq datum
(key-var-pow g
)) (desetq (g . pow
) datum
)
1171 (desetq (g . pow
) (adjoin-pvar g
)))
1172 (if (and (not (atom g
)) (psp g
))
1175 (setq datum
(gvar-data g
))
1176 (do ((po-terms (p-terms poly
) (p-red po-terms
))
1178 (push (term (e* pow
(prep1 (pt-le po-terms
)))
1179 (poly->ps
(pt-lc po-terms
)))
1182 ;; This must be exact so that when we invert in rat-ps above
1183 ;; we don't lose any terms. E.g. try
1184 ;; taylor(log(1+exp(-1/x)),x,0,5). When taylor2'ing exp(-1/x),
1185 ;; if you used current trunc here this would return exp(1/x)...5
1186 ;; which would then be trunc'd to degree 3 by psexpt.
1187 (make-ps (int-var datum
)
1188 (ncons (current-trunc datum
))
1189 (if (eq g
(data-gvar datum
)) ps-terms
1190 (invert-terms ps-terms
))))))))))
1192 (defun key-var-pow (g)
1193 (let ((var (get-key-var g
)) datum
)
1195 (setq datum
(get-datum var
))
1196 (if (eq g
(setq g
(data-gvar datum
))) (cons g
(rcone))
1197 (cons g
(rcmone))))))
1199 (defun adjoin-pvar (g)
1200 (let ((kernel (get g
'disrep
)) g
* lim datum ans
1201 (no (1+ (cdr (int-var (car (last tlist
)))))) (pow (rcone)) expt
)
1202 (when (assol kernel tlist
) (break "bad1"))
1203 (if (and (eq (caar kernel
) 'mexpt
) (eq (cadr kernel
) '$%e
)
1204 (not (atom (setq expt
(caddr kernel
))))
1205 (eq (caar expt
) 'mtimes
) (not (mfree expt
(ncons '$%i
))))
1206 (destructuring-let (((rpart . ipart
) (trisplit expt
)))
1207 (cons (pstimes (prep1 (m^
'$%e rpart
))
1208 (psplus (adjoin-tvar `((%cos
) ,ipart
))
1209 (pstimes (prep1 '$%i
)
1210 (adjoin-tvar `((%sin
) ,ipart
)))))
1213 (when (eq (caar kernel
) 'mexpt
)
1214 (when (and (not (atom (setq expt
(caddr kernel
))))
1215 (eq (caar expt
) 'mtimes
)
1216 ($ratnump
(cadr expt
)))
1217 (setq pow
(cadr expt
) kernel
(m^ kernel
(m// pow
))
1218 g
(prep1 kernel
) pow
(prep1 pow
))
1219 (unless (and (equal (cdr g
) 1) (equal (cdar g
) '(1 1)))
1220 (break "Illegal kernel in `adjoin-pvar'"))
1221 (setq g
(caar g
) kernel
(get g
'disrep
))))
1222 (if (setq ans
(key-var-pow g
))
1223 (cons (car ans
) (e* pow
(cdr ans
)))
1225 (when (lim-infp (or lim
(setq lim
(tvar-lim kernel
))))
1226 (setq g
* g g
(gensym) kernel
(m// kernel
)
1227 lim
(lim-recip lim
) pow
(e* (rcmone) pow
))
1228 (putprop g kernel
'disrep
)
1229 ;(push g genvar) (push kernel varlist)
1230 (push (cons g
* kernel
) key-vars
))
1231 (when (assol kernel tlist
) (break "bad2"))
1232 (setq datum
(list* kernel
1233 ;(mapcar #'(lambda (e) (emax e (rczero)))
1234 ; (trunc-stack (car tlist)))
1235 (copy-list (trunc-stack (car tlist
)))
1237 ;(setq tlist (nconc tlist (ncons datum)))
1238 (adjoin-datum datum
)
1239 (push (cons g kernel
) key-vars
)
1240 (push (car key-vars
) ivars
)
1241 ;(push (cons kernel (cons (pget g) 1)) genpairs)
1242 (push (cons g lim
) tvar-limits
)
1245 (defun adjoin-datum (datum)
1246 (do ((tlist* tlist
(cdr tlist
*))
1247 (tlist** () tlist
*))
1248 ((null tlist
*) (setq tlist
(nconc tlist
(ncons datum
))))
1249 (when (stronger-var?
(datum-var (car tlist
*)) (datum-var datum
))
1250 (return (if (null tlist
**)
1253 (renumber-tlist tlist
))
1255 (rplacd tlist
** (cons datum tlist
*))
1256 (renumber-tlist (cdr tlist
**))))))))
1258 ;; Maybe this should just permute the numbering in case it isn't sequential?
1260 (defun renumber-tlist (tlist)
1261 (rplacd (data-gvar-o (car tlist
)) (cdr (data-gvar-o (cadr tlist
))))
1262 (do ((tlist* (cdr tlist
) (cdr tlist
*)))
1264 (rplacd (data-gvar-o (car tlist
*))
1265 (1+ (cdr (data-gvar-o (car tlist
*)))))))
1267 (defun tvar?
(var) (or (atom var
) (member 'array
(cdar var
) :test
#'eq
)))
1269 ;; Needs to be extended to handle singular tvars in > 1 var's.
1271 (defun stronger-var?
(v1 v2
)
1272 (let ((e1 (rcone)) (e2 (rcone)) reverse? ans
)
1273 (when (alike1 v1 v2
)
1274 (tay-err (intl:gettext
"taylor: stronger-var? called on equal vars.")))
1275 (when (and (mexptp v1
) ($ratnump
(caddr v1
)))
1276 (setq e1
(prep1 (caddr v1
)) v1
(cadr v1
)))
1277 (when (and (mexptp v2
) ($ratnump
(caddr v2
)))
1278 (setq e2
(prep1 (caddr v2
)) v2
(cadr v2
)))
1282 (intl:gettext
"taylor: stronger-var? called on equal vars."))
1285 (when (eq (tvar-lim v2
) '$finite
)
1288 (setq reverse?
(not reverse?
)))
1289 (if (eq (tvar-lim v1
) '$finite
)
1290 (if (eq (tvar-lim v2
) '$finite
)
1291 (great v1 v2
) reverse?
)
1296 (setq reverse?
(not reverse?
)))
1299 (stronger-vars?
(order-vars-by-strength (cdr v1
))
1300 (order-vars-by-strength (if (mtimesp v2
) (cdr v2
)
1301 (ncons (m^ v2
(edisrep e2
))))))
1306 (setq reverse?
(not reverse?
)))
1309 (let ((n1 (cdr (data-gvar-o (get-datum v1 t
))))
1310 (n2 (cdr (data-gvar-o (get-datum v2 t
)))))
1312 ((mfree v2
(ncons v1
))
1314 (intl:gettext
"taylor: Unhandled multivar datum comparison.")))
1315 ((eq (caar v2
) '%log
) 't
)
1316 ((and (eq (caar v2
) 'mexpt
) (eq (cadr v2
) '$%e
))
1317 (stronger-var?
`((%log
) ,v1
) (caddr v2
)))
1318 (t (tay-err (intl:gettext
"taylor: Unhandled var in stronger-var?."))))
1320 (when (eq (caar v2
) '%log
)
1323 (setq reverse?
(not reverse?
)))
1324 (if (eq (caar v1
) '%log
)
1325 (cond ((eq (caar v2
) '%log
)
1326 (stronger-var?
(cadr v1
) (cadr v2
)))
1327 ((and (eq (caar v2
) 'mexpt
) (eq (cadr v2
) '$%e
))
1328 (stronger-var?
`((%log
) ,v1
) (caddr v2
)))
1329 (t (tay-err (intl:gettext
"taylor: Unhandled var in stronger-var?"))))
1330 (if (and (eq (caar v1
) 'mexpt
) (eq (cadr v1
) '$%e
)
1331 (eq (caar v2
) 'mexpt
) (eq (cadr v2
) '$%e
))
1332 (stronger-var?
(caddr v1
) (caddr v2
))
1333 (tay-err (intl:gettext
"taylor: Unhandled var in stronger-var?")))))))))
1334 (if reverse?
(not ans
) ans
)))))))
1336 (defun neg-monom?
(exp)
1337 (and (mtimesp exp
) (equal (cadr exp
) -
1) (null (cdddr exp
))
1340 (defun order-vars-by-strength (vars)
1341 (do ((vars* vars
(cdr vars
*)) (ordvars () ))
1342 ((null vars
*) ordvars
)
1343 (unless (mfree (car vars
*) tvars
) ; ignore constants
1344 (do ((ordvars* ordvars
(cdr ordvars
*)))
1346 (if (null ordvars
) (setq ordvars
(ncons (car vars
*)))
1347 (rplacd (last ordvars
) (ncons (car vars
*)))))
1348 (when (stronger-var?
(car vars
*) (car ordvars
*))
1349 (rplacd ordvars
* (cons (car ordvars
*) (cdr ordvars
*)))
1350 (rplaca ordvars
* (car vars
*))
1353 (defun stronger-vars?
(vars1 vars2
)
1354 (do ((vars1* vars1
(cdr vars1
*))
1355 (vars2* vars2
(cdr vars2
*)))
1357 (cond ((null vars1
*)
1359 ;; two equal vars generated
1361 (let ((lim (tvar-lim (car vars2
*))))
1363 (cond ((lim-infp lim
) ())
1364 ((lim-zerop lim
) 't
)
1365 (t (break "var with non-zero finite lim?")))))))
1367 (let ((lim (tvar-lim (car vars1
*))))
1369 (cond ((lim-infp lim
) 't
)
1370 ((lim-zerop lim
) ())
1371 (t (break "var with non-zero finite lim?"))))))
1372 ((alike1 (car vars1
*) (car vars2
*)) )
1373 ((return (stronger-var?
(car vars1
*) (car vars2
*)))))))
1375 (defun stronger-datum?
(d1 d2
)
1376 (setq d1
(datum-var d1
) d2
(datum-var d2
))
1377 (do ((end-flag) (answer))
1378 (end-flag (member answer
'($yes $y
) :test
#'eq
))
1379 (setq answer
(retrieve `((mtext) |Is |
,d1 | stronger than |
,d2 |?|
)
1381 (if (member answer
'($yes $y $no $n
) :test
#'eq
) (setq end-flag
't
)
1382 (mtell "~%Acceptable answers are: yes, y, no, n~%"))))
1384 (defun datum-lim (datum)
1385 (if (not (tvar?
(datum-var datum
)))
1387 (let ((pt (exp-pt datum
)))
1388 (if (member pt
'($inf $minf
) :test
#'eq
) pt
'$zeroa
))))
1390 (defun tvar-lim (kernel)
1391 (if (mfree kernel tvars
) (coef-sign kernel
)
1392 (let ((datum (get-datum kernel t
)) lim
)
1393 (or (and datum
(datum-lim datum
))
1394 (and (setq datum
(get-datum (m// kernel
) t
))
1395 (setq lim
(datum-lim datum
))
1399 (cond ((eq (caar kernel
) 'mexpt
)
1400 (cond ((and (setq datum
(get-datum (cadr kernel
) t
))
1401 ($ratnump
(caddr kernel
)))
1402 (lim-power (datum-lim datum
)
1403 (prep1 (caddr kernel
))))
1404 (($ratnump
(caddr kernel
))
1405 (lim-power (tvar-lim (cadr kernel
))
1406 (prep1 (caddr kernel
))))
1407 ((eq (cadr kernel
) '$%e
)
1408 (lim-exp (tvar-lim (caddr kernel
))))
1409 (t (tay-error "Unhandled case in tvar-lim" kernel
))))
1410 ((eq (caar kernel
) 'mtimes
)
1411 (do ((ans (tvar-lim (cadr kernel
))
1412 (lim-times ans
(tvar-lim (car facs
))))
1413 (facs (cddr kernel
) (cdr facs
)))
1415 ((eq (caar kernel
) '%log
)
1416 (lim-log (datum-lim (get-datum (cadr kernel
) t
))))
1417 ((member (caar kernel
) '(%sin %cos
) :test
#'eq
)
1418 (unless (lim-infp (tvar-lim (cadr kernel
)))
1419 (tay-error "Invalid trig kernel in tvar-lim" kernel
))
1421 (t (tay-error "Unhandled kernel in tvar-lim" kernel
))))
1424 (defun coef-sign (coef)
1425 (if (not ($freeof
'$%i
($rectform coef
)))
1429 (defun gvar-lim (gvar)
1430 (or (cdr (assoc gvar tvar-limits
:test
#'eq
))
1431 (if (member (gvar->var gvar
) tvars
:test
#'eq
) '$zeroa
; user tvars assumed 0+ now
1432 (break "Invalid gvar"))))
1434 (defun psexpt-fn1 (varh trunc l inc m mr ans
)
1435 (declare (fixnum m
))
1436 (prog (a (k 0) ak cm-k c sum kr lim
)
1437 (declare (fixnum k
))
1438 ;; truly unfortunate that we need so many variables in this hack
1440 b
(and (e> mr trunc
) (go end
))
1441 (setq kr inc ak l k
1 sum
(rczero) lim m
)
1442 a
(cond ((or (> k lim
)
1443 (null (setq cm-k
(psterm (cdr ans
) (e- mr kr
)))))
1445 (setq ak
(or (pstrim-terms ak kr
)
1447 c
(pstimes (ereduce k m
)
1448 (pstimes (psterm1 ak kr
) cm-k
))
1450 (setq k
(1+ k
) kr
(e+ kr inc
))
1453 (unless (rczerop sum
) (add-term-&-pop a mr sum
))
1454 (setq m
(1+ m
) mr
(e+ mr inc
))
1457 (return (pscheck varh
(list trunc
) (cdr ans
)))))
1459 ;;; PSEXPT-FN2 and RED-MONO-LOG are needed to reduce exponentials of logs.
1461 (defun psexpt-fn2 (p)
1462 (cond ((atom p
) (if (get-datum p
)
1463 (psexpt-fn (taylor2 p
))
1464 (prep1 `((mexpt) $%e
,p
))))
1465 ((eq (caar p
) '%log
)
1466 (if (get-datum (cadr p
)) (taylor2 (cadr p
)) (prep1 (cadr p
))))
1467 ((or (eq (caar p
) 'mplus
) (eq (caar p
) 'mtimes
))
1468 (let ((e ($ratexpand p
)) temp
)
1469 (cond ((not (and (consp e
) (member (caar e
) '(mplus mtimes
) :test
#'eq
)))
1472 (if (eq (caar e
) 'mplus
)
1473 (do ((sumnds (cdr e
) (cdr sumnds
)) (log-facs) (l))
1475 (cond ((not log-facs
) (tsexpt '$%e p
))
1476 (t (tstimes (cons (m^t
'$%e
(m+l l
)) log-facs
)))))
1477 (if (setq temp
(red-mono-log (car sumnds
)))
1478 (push temp log-facs
)
1479 (push (car sumnds
) l
)))
1481 (setq temp
(red-mono-log e
))
1484 (prep1 (power '$%e p
)))))))))
1485 (t (prep1 (power '$%e p
)))))
1487 (defun red-mono-log (e)
1489 ((eq (caar e
) '%log
) (cadr e
))
1491 (do ((facs (cdr e
) (cdr facs
)) (log-term))
1494 (m^t
(cadr log-term
) (m*l
(remove log-term
(cdr e
) :test
#'eq
)))))
1495 (if (and (null (atom (car facs
))) (eq (caaar facs
) '%log
))
1496 (if log-term
(return ()) (setq log-term
(car facs
)))
1497 (unless (mfree (car facs
) tvars
) (return nil
)))))
1501 (if (pscoefp p
) (pslog2 (rcdisrep p
))
1502 (let ((terms (terms p
)))
1503 (cond ((mono-term? terms
) ; log(c x^n) = log(c) + n log(x)
1504 ;; do this always for now
1505 (if 't
;$TAYLOR_LOGEXPAND
1506 ;(psplus (pslog (lc terms))
1507 ; (pstimes (le terms) (pslog-of-gvar (gvar p))))
1509 ;(prep1 `((%LOG) ,(term-disrep (lt terms) p)))
1511 ;; expand log(1+ax^n) directly by series substitution
1512 ((not (or (n-term (setq terms
(terms (psplus p
(rcmone)))))
1513 ;(e> (rczero) (le terms))
1515 (setq p
(get-series '%log
(psexpt-log-ord p
) (gvar-o p
)
1516 (if (e> (rczero) (le terms
)) (e- (le terms
))
1519 (if (e> (rczero) (le terms
)) (ps-invert-var p
) p
))
1520 (t (prog (l inc trunc lt ans lterm $maxtayorder gvar-lim gt
)
1521 ;; log(lt+y) = log(lt) + log(1 + y/lt) = lterm + p
1522 (setq trunc
(trunc-lvl p
))
1523 (if (not (member (setq gvar-lim
(gvar-lim (gvar p
)))
1524 '($zeroa $zerob $inf $minf
) :test
#'eq
))
1525 (tay-error "bad gvar lim" gvar-lim
)
1526 (if (member gvar-lim
'($inf $minf
) :test
#'eq
)
1527 (setq lt
(ps-gt p
) gt lt
)
1528 (setq lt
(ps-lt p
) gt
() )))
1530 (setq lt
(pscheck (gvar-o p
)
1533 p
(pstimes p
(let (($maxtayorder
't
))
1534 (psexpt lt
(rcmone)))))
1535 (when (and (member gvar-lim
'($inf $minf
) :test
#'eq
)
1536 (e> (le terms
) (rczero)))
1537 (return (psplus lterm
(pslog p
))))
1539 (unless (equal p
(rcone))
1540 (merror "PSLOG: internal error."))
1542 (setq l
(terms p
) inc
(psexpon-gcd l
))
1543 (if gt
(setq l
(delete (last l
) l
:test
#'equal
))
1544 (setq l
(n-term l
)))
1545 (setq ans
(ncons 0))
1546 (unless $maxtayorder
1547 (setq trunc
(emin trunc
(t-o-var (gvar p
)))))
1548 ;; When we've divided by the greatest term, all terms
1549 ;; have non-positive exponents and we must perform the
1550 ;; transformation x -> 1/x befor calling pslog1 and then
1551 ;; perform the inverse afterwards.
1552 (when gt
(setq l
(invert-terms l
)))
1553 (when (e> (rczero) inc
) (setq inc
(e- inc
)))
1554 (setq ans
(psplus lterm
1555 (pslog1 (gvar-o p
) trunc l inc
1 inc ans
)))
1557 (if (and gt
(psp ans
) (eq (gvar ans
) (gvar p
)))
1561 (defun invert-terms (terms)
1562 (nreverse (mapc #'(lambda (x) (rplaca x
(e- (e x
)))) terms
)))
1564 (defun ps-invert-var (ps)
1565 (when (psp ps
) (rplacd (cddr ps
) (invert-terms (terms ps
))))
1569 (if (pscoefp ps
) (term (rczero) ps
)
1570 (lt (last (terms ps
)))))
1572 (defun pslog1 (varh trunc l inc m mr ans
)
1573 (declare (fixnum m
))
1574 (prog (a (k 0) ak cm-k c sum kr m-kr
)
1575 (declare (fixnum k
))
1576 ;; truly unfortunate that we need so many variables in this hack
1579 b
(and (e> mr trunc
) (go end
))
1580 (setq kr inc ak l k
1 sum
(rczero))
1581 a
(cond ((or (= k m
)
1582 (null (setq cm-k
(psterm (cdr ans
)
1583 (setq m-kr
(e- mr kr
))))))
1585 (setq ak
(or (pstrim-terms ak kr
)
1587 c
(pstimes m-kr
(pstimes (psterm1 ak kr
) cm-k
))
1589 k
(1+ k
) kr
(e+ kr inc
))
1592 (cond ((setq c
(pstrim-terms ak mr
))
1593 (setq c
(psterm1 c mr
)))
1594 ((setq c
(rczero))))
1595 (setq sum
(psdiff c
(pstimes sum
(e// mr
))))
1596 (unless (rczerop sum
) (add-term-&-pop a mr sum
))
1597 (setq m
(1+ m
) mr
(e+ mr inc
))
1600 (return (pscheck varh
(list trunc
) (cdr ans
)))))
1602 ;; Computes log(monom), where monom = c x^n. Is extra careful trying to keep
1603 ;; singular logs going to INF and not generating log(-1)'s unless it is
1604 ;; necessary to transform a log at MINF to INF.
1606 (defun pslog-monom (monom)
1607 (let* ((gvar (gvar monom
))
1608 (datum (gvar-data gvar
)) var pt logvar c
)
1609 (if (switch 'multivar datum
)
1610 (pslog (ps-lc monom
))
1612 (setq var
(datum-var datum
))
1614 (if (not (member (setq pt
(exp-pt datum
)) '($inf $minf
) :test
#'eq
))
1615 (setq logvar
(adjoin-tvar `((%log
) ,(m- var pt
))))
1617 ;; At x = inf: log(c (1/x)^n) -> log(c) - n log(x)
1618 ;; At x = minf: log(c (-1/x)^n) -> log(c (-1)^n) - n log(x)
1619 (setq logvar
(psminus (adjoin-tvar `((%log
) ,var
))))
1620 (when (eq pt
'$minf
)
1621 (setq c
(rcexpt (rcmone) (ps-le monom
))))))
1622 (if (eq (caar var
) 'mexpt
)
1623 (if (equal (caddr var
) -
1);; var must be 1/log(y)
1624 ;; Try to keep inf's real. Here we want
1625 ;; log(c (1/log(x))^n) -> log(c (-1)^n) - n log(-log(x))
1626 (if (equal (tvar-lim (cadr var
)) '$minf
)
1627 (setq c
(rcexpt (rcmone) (ps-le monom
))
1629 (psminus (adjoin-tvar
1630 `((%log
) ,(m- (cadr var
))))))
1631 (setq logvar
(psminus
1632 (adjoin-tvar `((%log
) ,(cadr var
))))))
1633 (if (equal (cadr var
) '$%e
)
1634 (setq logvar
(taylor2 (caddr var
)))
1635 (break "Unhandled gvar in `pslog-of-gvar'")))))
1636 (psplus (pslog (if c
(pstimes c
(ps-lc monom
)) (ps-lc monom
)))
1637 (pstimes (ps-le monom
) logvar
))))))
1639 ;; Computes log(p), where p is an rcdisrep'd pscoef.
1641 (defun pslog2 (p) (let ($logarc
) (pslog3 p
)))
1645 (prep1 (cond ((equal p
1) 0)
1646 ((equal p -
1) log-1
)
1650 (merror (intl:gettext
"taylor: log(0) encountered while processing ~:M") last-exp
))
1653 (prep1 (cond ((not $taylor_logexpand
) `((%log
) ,p
))
1654 (t (m- `((%log
) ,(cadr p
)) `((%log
) ,(caddr p
)))))))
1655 ((and full-log
(not (free p
'$%i
)))
1656 (let ((full-log () )) (pslog3 ($polarform p
))))
1657 ((eq (caar p
) 'mexpt
)
1658 ;; Must handle things like x^a, %e^(a*x), etc. which are pscoef's.
1659 (pstimes (taylor2 (caddr p
)) (pslog (taylor2 (cadr p
)))))
1660 ((and (eq (caar p
) 'mtimes
) $taylor_logexpand
)
1661 (do ((l (cddr p
) (cdr l
))
1662 (ans (pslog3 (cadr p
)) (psplus ans
(pslog3 (car l
)))))
1664 (t (prep1 `((%log
) ,p
)))))
1666 ;;; Subtitle Extending Routines
1668 (defun getfun-lt (fun)
1669 (let ((exp-datum (get (oper-name fun
) 'exp-form
)))
1671 ;; Info not needed yet.
1672 ;; (or (atom (car exp-datum))
1673 ;; (setq 0p-funord (copy-tree (cdar exp-datum))))
1674 (exp-datum-lt fun exp-datum
))
1675 ((setq exp-datum
(get (oper-name fun
) 'sp2
))
1676 (setq exp-datum
(get-lexp (subst (dummy-var) 'sp2var exp-datum
)
1678 ;; Info not needed yet; need to bind lexp-non0 to T when
1679 ;; this is used though so n-term will be there.
1680 ;; (and (rczerop (le exp-datum))
1681 ;; (setq 0p-funord (le (n-term exp-datum))))
1682 (if (psp exp-datum
) (ps-lt exp-datum
)
1683 (term (rczero) exp-datum
)))
1684 (t (merror "GETFUN-LT: unknown function ~A" fun
)))))
1686 (declare-top (special var
))
1688 (defun getexp-fun (fun var pw
)
1689 (declare (special var
))
1690 (let ((exp-datum (copy-tree (get (oper-name fun
) 'exp-form
))))
1691 (cond ((infp pw
) (infin-ord-err))
1693 (if (null (setq exp-datum
1694 (get-ps-form (if (atom fun
) fun
(caar fun
)))))
1695 (merror (intl:gettext
"taylor: power series unavailable for function ~A") fun
)
1698 (do ((subvals (cdr fun
) (cdr subvals
))
1699 (subs (safe-get (caar fun
) 'sp2subs
) (cdr subs
)))
1700 ((or (null subvals
) (null subs
))
1701 (when (or subvals subs
)
1702 (merror (intl:gettext
"taylor: incorrect number of subscripts to the deftaylor'd function ~A") (caar fun
))))
1703 (setq exp-datum
(maxima-substitute (car subvals
) (car subs
)
1705 (ts-formula exp-datum var pw
))))
1706 ((e> (exp-datum-le fun exp-datum
) pw
) (pszero var pw
))
1708 (apply (exp-fun exp-datum
)
1709 (if (atom fun
) (cons pw
(cdr exp-datum
))
1710 (cons pw
(cons (cdr fun
) (cdr exp-datum
))))))
1711 (cond ((null exp-datum
) (pszero var pw
))
1712 ((psp exp-datum
) exp-datum
)
1713 (t (make-ps var
(ncons pw
) exp-datum
)))))))
1715 (declare-top (unspecial var
))
1717 (defun expexp-funs (pw l sign chng inc
)
1719 (setq e
(e l
) lt-l
(setq l
(ncons l
)))
1720 a
(cond ((e> (setq e
(e+ e inc
)) pw
) (return l
))
1725 (cond ((e= inc
(rcone)) e
)
1727 (cons 1 (cdr (lc lt-l
)))))
1728 (setq sign
(e* sign chng
))))
1731 ;; returns series expansion of %expintegral_si
1732 ;; from term x^i to term x^max
1733 ;; sign is sign of term x^i (1 or -1)
1734 ;; ifac is i! (to avoid repeated computation)
1735 (defun expsi_series (i sign ifac max
)
1737 nil
; we have all powers up to max
1738 (cons (cons (cons i
1) ; this is i'th term
1739 (cons sign
(* i ifac
))); i'th coefficient
1740 (expsi_series (+ 2 i
) ; advance to next even/odd power
1741 (* -
1 sign
) ; flip sign
1742 (* (+ 2 i
) (+ 1 i
) ifac
) ; update factorial
1743 max
)))) ; pass through the max
1745 ;; returns series expansion of %expintegral_si up to term with x^pw
1746 (defun exp_%expintegral_si
(pw l
)
1750 (/ (float (car pw
)) (float (cdr pw
)))))
1752 (defun explog-funs (pw l sign chng inc
)
1754 (setq e
(e l
) lt-l
(setq l
(ncons l
)))
1755 a
(cond ((e> (setq e
(e+ e inc
)) pw
) (return l
))
1756 (t (add-term lt-l e
(e// sign e
))
1757 (setq lt-l
(n-term lt-l
)
1758 sign
(e* sign chng
))))
1761 (defun exptan-funs (pw l chng
)
1762 (prog (e lt-l sign fact pow
)
1763 (setq e
(e l
) lt-l
(setq l
(ncons l
))
1764 sign
(rcone) fact
'(1 .
2) pow
'(4 .
1))
1765 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1766 (t (setq fact
(e// fact
(e* e
(e1+ e
)))
1767 pow
(e* '(4 .
1) pow
)
1768 sign
(e* chng sign
))
1769 (add-term lt-l e
(e* (e* sign fact
)
1771 ($bern
(rcdisrep (e1+ e
))))
1772 (e* pow
(e1- pow
)))))
1773 (setq lt-l
(n-term lt-l
))))
1776 (defun expcot-funs (pw l sign chng plus
)
1777 (prog (e lt-l fact pow
)
1778 (setq e
(e l
) lt-l
(setq l
(ncons l
))
1779 fact
(rcone) pow
(rcone))
1780 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1781 (t (setq fact
(e// fact
(e* e
(e1+ e
)))
1782 pow
(e* '(4 .
1) pow
)
1783 sign
(e* chng sign
))
1784 (add-term lt-l e
(e* (e* sign fact
)
1786 ($bern
(rcdisrep (e1+ e
))))
1788 (setq lt-l
(n-term lt-l
))))
1791 (defun expsec-funs (pw l chng
)
1792 (prog (e lt-l sign fact
)
1793 (setq e
(e l
) lt-l
(setq l
(ncons l
))
1794 sign
(rcone) fact
(rcone))
1795 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1796 (t (setq fact
(e// fact
(e* e
(e1- e
)))
1797 sign
(e* chng sign
))
1798 (add-term lt-l e
(e* (e* sign fact
)
1799 (prep1 ($euler
(rcdisrep e
)))))
1800 (setq lt-l
(n-term lt-l
))))
1803 (defun expasin-funs (pw l chng
)
1804 (prog (e lt-l sign n d
)
1805 (setq e
(e l
) lt-l
(setq l
(ncons l
)) sign
1 n
1 d
1)
1806 a
(cond ((e> (setq e
(e+ (rctwo) e
)) pw
) (return l
))
1807 (t (setq n
(* n
(car (e- e
(rctwo))))
1808 d
(* d
(car (e1- e
)))
1810 (add-term lt-l e
; need to reduce here ? - check this.
1811 (let ((x (*red
(* n sign
)
1814 (cons (cadr x
) (caddr x
)))))
1815 (setq lt-l
(n-term lt-l
))))
1818 ;;; This is the table of expansion data for known functions.
1819 ;;; The format of the EXP-FORM property is as follows:
1820 ;;; (<name of the expanding routine for the function or
1821 ;;; (name . le of n-term) if expansion is of order 0>
1822 ;;; <first term in the expansion or the name of a routine which
1823 ;;; computes the order when it may depend on parameters (e.g subsripts)>
1824 ;;; <data for the expanding routine>)
1827 (loop for
(fun exp
) on
1828 '(%exp
((expexp-funs 1 .
1) ((0 .
1) 1 .
1) (1 .
1) (1 .
1) (1 .
1))
1829 %sin
(expexp-funs ((1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (2 .
1))
1830 %cos
((expexp-funs 2 .
1) ((0 .
1) 1 .
1) (-1 .
1) (-1 .
1) (2 .
1))
1831 %sinh
(expexp-funs ((1 .
1) 1 .
1) (1 .
1) (1 .
1) (2 .
1))
1832 %cosh
((expexp-funs 2 .
1) ((0 .
1) 1 .
1) (1 .
1) (1 .
1) (2 .
1))
1833 %log
(explog-funs ((1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (1 .
1))
1834 %atan
(explog-funs ((1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (2 .
1))
1835 %atanh
(explog-funs ((1 .
1) 1 .
1) (1 .
1) (1 .
1) (2 .
1))
1836 %cot
(expcot-funs ((-1 .
1) 1 .
1) (1 .
1) (-1 .
1) (0 .
1))
1837 %csc
(expcot-funs ((-1 .
1) 1 .
1) (-1 .
1) (-1 .
1) (-2 .
1))
1838 %csch
(expcot-funs ((-1 .
1) 1 .
1) (-1 .
1) (1 .
1) (-2 .
1))
1839 %coth
(expcot-funs ((-1 .
1) 1 .
1) (1 .
1) (1 .
1) (0 .
1))
1840 %tan
(exptan-funs ((1 .
1) 1 .
1) (-1 .
1))
1841 %tanh
(exptan-funs ((1 .
1) 1 .
1) (1 .
1))
1842 %sec
((expsec-funs 2 .
1) ((0 .
1) 1 .
1) (-1 .
1))
1843 %sech
((expsec-funs 2 .
1) ((0 .
1) 1 .
1) (1 .
1))
1844 %asin
(expasin-funs ((1 .
1) 1 .
1) 1)
1845 %asinh
(expasin-funs ((1 .
1) 1 .
1) -
1)
1846 %gamma
(expgam-fun ((-1 .
1) 1 .
1))
1847 $li
(exp$li-fun li-ord
)
1848 %expintegral_si
(exp_%expintegral_si
((1 .
1) 1 .
1))
1849 $psi
(expplygam-funs plygam-ord
))
1851 do
(putprop fun exp
'exp-form
))
1854 (defun known-ps (fun)
1855 (getl fun
'(exp-form sp2
)))
1857 ;;; Autoload Properties
1859 ;;; Taylor series expansion routines
1861 ;;; SRF is only called externally; by RATF and SIMPEXPT.
1864 (let ((exact-poly t
) (tlist) (*within-srf?
* 't
))
1865 (setq x
(taylor1 x
()) tlist
(mrat-tlist x
))
1866 ;; Set trunc levels in the local tlist to correspond to the maximum
1867 ;; level occurring in any series.
1868 (do ((data tlist
(cdr data
))
1869 (truncs (trunc-vector (mrat-ps x
) () )))
1871 (when (and (car truncs
) (e> (car truncs
) (current-trunc (car data
))))
1872 (setf (current-trunc (car data
)) (car truncs
))))
1875 ;;; [var, pt, order, asymp]
1877 (defmfun $taylor
(e &rest args
)
1878 (when (not ($taylorp e
))
1879 ;; Not a taylor expression. Remove the special representation.
1880 (setq e
(specrepcheck e
)))
1883 (defun taylor* (arg l
)
1884 ;; We must bind $MAXTAYORDER to () below because of the problem of constants
1885 ;; not retaining their truncation level. This means that when we add a
1886 ;; series which has more terms than the user-specified truncation to a
1887 ;; constant we must truncate the series with more terms down to the user
1888 ;; specified level because, in the worst case, the constant could be a
1889 ;; series no better than to the user-specified level. Hence $MAXTAYORDER
1890 ;; is essentially useless until the constant problem is fixed. If one
1891 ;; decides to not bind $MAXTAYORDER below then the sum routines must
1892 ;; be updated to truncate series with more terms than the user-specified
1893 ;; level down to that level---taylor(sin(x)^2-cos(x)^2-1,x,0,1) would
1894 ;; give x^2+... in this case if the sum routines weren't updated.
1895 ;; Also, batch(mquery,160,aljabr) for another truncation bug which crops
1896 ;; up when $maxtayorder isn't bound here. Similarly, loadfile(taybad,rl,
1897 ;; aljabr) and see tomh's bug note of 4/15/81.
1898 (let ((tlist () ) ($maxtayorder
() ) (*within-srf?
* () )
1899 (exact-poly (if l
(not $taylor_truncate_polynomials
) 'user-specified
)))
1900 (declare (special *within-srf?
*))
1903 (taylor1 arg
(ncons tlist
))))
1905 (defun tay-order (n)
1906 (let (($float
) (modulus))
1907 (cond ((eq n
'$inf
) (ncons (inf)))
1908 ((null n
) (wna-err '$taylor
))
1910 (merror (intl:gettext
"taylor: expansion order must be a number; found: ~:M") n
))
1911 (t (ncons (prep1 n
))))))
1913 (defun re-erat (head exp
)
1914 (taylor1 exp
(list (cadddr (cdr head
)))))
1916 (defun parse-tay-args (l)
1919 (merror (intl:gettext
"taylor: variable of expansion cannot be a number: ~M") (car l
)))
1920 ((or (symbolp (car l
)) (not (eq (caaar l
) 'mlist
)))
1921 (parse-tay-args1 (list (car l
) ($ratdisrep
(cadr l
)) (caddr l
)))
1922 (parse-tay-args (cdddr l
)))
1923 ((do ((l (cddar l
) (cdr l
)))
1925 (and (or (mnump (car l
)) (eq (car l
) '$inf
))
1927 (parse-tay-args1 (cdar l
))
1928 (parse-tay-args (cdr l
)))
1929 (t (parse-tay-args2 (list (car l
) (cadr l
) (caddr l
)))
1930 (parse-tay-args (cdddr l
)))))
1932 (defun parse-tay-args1 (l)
1933 (if ($listp
(car l
)) (parse-tay-args2 l
)
1935 (pt ($ratdisrep
(cadr l
)))
1936 (ord (tay-order (caddr l
)))
1937 (switches (make-switch-list (cdddr l
))))
1938 (push (list v ord pt switches
) tlist
))))
1940 (defun parse-tay-args2 (l)
1941 (let ((label (gensym))
1943 (pts (make-long-list (if ($listp
(cadr l
))
1944 (copy-list (cdadr l
))
1945 (ncons (ratdisrep (cadr l
))))))
1947 (switches (make-switch-list (cdddr l
)))
1951 (setq lcm ord max ord ord
(make-long-list (ncons ord
)))
1953 (l (cdr ord
) (cdr l
)))
1954 ((null a
) (setq ord
(cdr ord
)))
1955 (cond ((not l
) (merror "PARSE-TAY-ARGS2: ran out of truncation levels."))
1956 (t (setq lcm
(lcm lcm
(car l
)) max
(max max
(car l
)))))))
1957 (push (list label
(tay-order max
) 0
1958 (ncons (list 'multivar lcm vs
)))
1960 (do ((vl vs
(cdr vl
))
1961 (ordl ord
(cdr ordl
))
1962 (ptl pts
(cdr ptl
)))
1964 (cond ((not ptl
) (merror "PARSE-TAY-ARGS2: ran out of matching points of expansion."))
1967 (list (car vl
) (tay-order (car ordl
)) (car ptl
)
1968 (cons (list 'multi label
(timesk lcm
(expta (car ordl
) -
1))) switches
))
1971 (defun make-switch-list (l)
1972 (mapcar #'(lambda (q) (cons q t
)) l
))
1974 (defun make-long-list (q)
1977 ;;; This checks to ensure that there isn't more than one set of multi-
1978 ;;; dependent variables with different orders of expansion, e.g.
1979 ;;; taylor(u+v+x+y,[u,v],0,[2,3],[x,y],0,[5,7]) is one.
1981 (defun ratwtsetup (l)
1982 (do ((l l
(cdr l
)) (a) (sw))
1984 (when (setq a
(switch 'multivar
(car l
)))
1985 (do ((ll (cadr a
) (cdr ll
)))
1987 (cond ((equal (cadr (switch 'multi
(get-datum (car ll
)))) 1) )
1988 (sw (merror (intl:gettext
"taylor: multiple dependent variables must all have the same order of expansion.")))
1989 ('t
(setq sw
't
) (return 't
)))))))
1991 (defmvar $taylor_order_coefficients t
1992 "When `true', coefficients of taylor series will be ordered canonically.")
1994 (defun taylor1 (e tlist
)
1995 (declare (special *within-srf?
* ))
1996 (setq tlist
(tlist-merge (nconc (find-tlists e
) tlist
)))
1997 (prog ($zerobern $simp $algebraic genpairs varlist tvars sing-tvars
1998 log-1 log%i ivars key-vars ans full-log last-exp
1999 mainvar-datum zerolist taylor_simplifier least_term? tvar-limits
2001 (setq tlist
(mapcan #'(lambda (d)
2002 (if (tvar?
(datum-var d
))
2008 (setq $zerobern t $simp t $algebraic t last-exp e least_term?
't
2009 log-1
'((%log simp
) -
1) log%i
'((%log simp
) $%i
)
2010 tvars
(mapcar 'car tlist
) varlist
(copy-list tvars
))
2011 (when $taylor_simplifier
2012 ; This symbolp/fboundp check is presumably for efficiency (so it
2013 ; can be directly funcalled).
2014 (setq taylor_simplifier
2015 (if (and (symbolp $taylor_simplifier
)
2016 (fboundp $taylor_simplifier
))
2018 'taylor_simplifier_caller
)))
2019 ;; Ensure that the expansion points don't depend on the expansion vars.
2020 ;; This could cause an infinite loop, e.g. taylor(x,x,x,1).
2021 (do ((tl tlist
(cdr tl
)))
2023 (unless (mfree (exp-pt (car tl
)) tvars
)
2024 (merror (intl:gettext
"taylor: attempt to expand ~M at a point depending on ~M.") e
(caar tl
))))
2025 ;; This drastic initialization ensures that ALGEBRAIC, TELLRAT, DISREP,
2026 ;; etc. prop's are removed from our gensyms. RATSETUP does not appear
2027 ;; to do this correctly, e.g. see ASB's bug of 1/10/83 (MQUERY 17).
2028 (mapc #'(lambda (g) (setf (symbol-plist g
) nil
)) genvar
)
2029 (ratsetup varlist genvar
)
2030 (when (and $taylor_order_coefficients
(not *within-srf?
*)) (newvar e
))
2031 (orderpointer varlist
)
2032 (maplist #'(lambda (q g
)
2033 (push (cons (car g
) (car q
)) key-vars
)
2034 (let ((data (get-datum (car q
))))
2035 (rplaca q
(transform-tvar (car q
) data
))
2036 (push (cons (car g
) (car q
)) ivars
)
2037 ;(setf (data-gvar-o data)
2038 ; (cons (car g) (valget (car g))))
2039 (rplacd (cdddr data
)
2040 (cons (car g
) (valget (car g
))))))
2041 (do ((v varlist
(cdr v
)))
2042 ((eq (car v
) (car tvars
)) v
))
2043 (do ((v varlist
(cdr v
))
2045 ((eq (car v
) (car tvars
)) g
)))
2046 (setq genpairs
(mapcar #'(lambda (y z
)
2047 (putprop z y
'disrep
)
2048 (cons y
(cons (pget z
) 1)))
2051 (setup-multivar-disrep () )
2052 (setq mainvar-datum
(car (last tlist
)))
2053 (mapc #'(lambda (d) (adjoin-sing-datum d
)) sing-tvars
)
2054 (setq ans
(catch 'tay-err
(taylor3 e
)))
2056 (if (atom (car ans
)) (tay-error (car ans
) (cadr ans
)) ans
))))
2058 (defun transform-tvar (var data
)
2059 (if (not (tvar? var
)) var
2060 (cond ((and (signp e
(exp-pt data
)) (null (switch '$asymp data
)))
2062 ((member (exp-pt data
) '($inf infinity
) :test
#'eq
)
2064 ((eq (exp-pt data
) '$minf
)
2066 ((let ((temp (m- var
(exp-pt data
))))
2067 (if (switch '$asymp data
) (m^ temp -
1) temp
))))))
2069 (defun taylor_simplifier_caller (exp)
2070 (mcall $taylor_simplifier exp
))
2072 (defun taylor_simplify_recurse (ps)
2073 (if (pscoefp ps
) (taylor2 (funcall taylor_simplifier
(rcdisrep ps
)))
2074 (let ((datum (ps-data ps
)) (var () ))
2075 ;; We must treat multivars like 1, since they'll reappear again
2076 ;; when we call taylor2 on their disrep'd coeff's.
2077 (if (switch 'multivar datum
)
2080 (setq var
(getdisrep (gvar-o ps
)))
2081 ;; Don't push pw's < 0, else constants will be truncated
2082 (push-pw datum
(emax (trunc-lvl ps
) (rczero)))))
2083 (do ((terms (terms ps
) (n-term terms
))
2084 (ans (rczero) (psplus (if (null datum
)
2085 (taylor_simplify_recurse (lc terms
))
2086 (pstimes (taylor_simplify_recurse
2089 ;; (taylor2 (funcall taylor_simplifier
2090 ;; (m^ var (edisrep (le terms)))))
2091 ;; causes terms to be lost when inverting. E.g.
2092 ;; taylor(log(1+exp(-1/x)),x,0,5) calls psexpt(<exp(1/x)^3>...3,-1)
2093 ;; which must return a series good to 3+3(-1-1)=-3 which, when added
2094 ;; to other terms will truncate them to degree -3 also.
2095 (if (ezerop (le terms
)) (rcone)
2098 (term (le terms
) (rcone)))))))
2102 (when datum
(pop-pw datum
))
2105 (defun push-pw (datum pw
)
2106 (push pw
(trunc-stack datum
))
2107 ;; When changing the truncation on an internal multivar we must also
2108 ;; propagate the change to all var's which depend upon it. See WGD's
2109 ;; bug report of 9/15/82 which exhibits the necessity of this.
2110 (when (setq datum
(switch 'multivar datum
))
2111 (do ((vars (cadr datum
) (cdr vars
)))
2113 (push pw
(trunc-stack (get-datum (car vars
)))))))
2115 (defun pop-pw (datum)
2116 (pop (trunc-stack datum
))
2117 ;; See the comment above in push-pw; here we must undo the propagation.
2118 (when (setq datum
(switch 'multivar datum
))
2119 (do ((vars (cadr datum
) (cdr vars
)))
2121 (pop (trunc-stack (get-datum (car vars
)))))))
2123 (defun setup-multivar-disrep (mrat?
)
2124 (let ((varlist varlist
) (genvar genvar
) (multivars () ))
2126 (setq varlist
(mrat-varlist mrat?
) genvar
(mrat-genvar mrat?
)))
2127 (mapc #'(lambda (datum)
2128 (when (switch 'multivar datum
)
2129 (push (car datum
) multivars
)
2130 ;; All genvar's corresponding to internally generated
2131 ;; multivars must "disappear" when disrep'd. If this
2132 ;; were not done then disrep'ing gx*gt would give x*t
2133 ;; which, upon, re-tayloring would give (gx*gt)*gt,
2134 ;; where t is the internal multivar for x, and gt, gx
2135 ;; are their genvars. An example where this would occur is
2136 ;; taylor(sin(x+y),[x],0,f1,[y],0,1).
2137 (putprop (int-gvar datum
) 1 'disrep
)))
2138 (if mrat?
(mrat-tlist mrat?
) tlist
))
2139 ;; Here we must substitute 1 for any genvars which depend on multivars.
2140 ;; For example, taylor(x^n,[x],0,0) generates a multivar^n.
2142 (do ((expl varlist
(cdr expl
))
2143 (gvarl genvar
(cdr gvarl
)))
2145 (unless (mfree (car expl
) multivars
)
2146 (putprop (car gvarl
) 1 'disrep
))))))
2148 ;; An example showing why this flag is need is given by
2149 ;; taylor(-exp(exp(-1/x)+2/x),x,0,-1). Without it, tstimes and
2150 ;; taylor_simplify_recurse end up trunc'ing the -1.
2152 (defvar trunc-constants?
't
)
2155 (cond ((mbagp e
) (cons (car e
) (mapcar #'taylor3
(cdr e
))))
2156 ((and (null tlist
) (not (eq exact-poly
'user-specified
)))
2158 (list 'mrat
'simp varlist genvar
)))
2159 (t (xcons (if (null taylor_simplifier
)
2162 (setq e
(taylor2 e
))
2163 (let ((exact-poly () ) (trunc-constants?
() ))
2164 (taylor_simplify_recurse e
))))
2165 (list 'mrat
'simp varlist genvar tlist
'trunc
)))))
2167 (defun find-tlists (e) (let (*a
*) (findtl1 e
) *a
*))
2170 (cond ((or (atom e
) (mnump e
)) )
2171 ((eq (caar e
) 'mrat
)
2172 (when (member 'trunc
(car e
) :test
#'eq
)
2173 (push (mapcar #'copy-tree
(mrat-tlist e
)) *a
*)))
2174 (t (mapc #'findtl1
(cdr e
)))))
2176 (defun tlist-merge (tlists)
2177 (do ((tlists tlists
(cdr tlists
)) (tlist () ))
2178 ((null tlists
) (nreverse tlist
))
2179 (do ((a_tlist (car tlists
) (cdr a_tlist
)) (temp nil
))
2181 (if (null (setq temp
(get-datum (datum-var (car a_tlist
)) t
)))
2182 (if (prog2 (setq temp
(car a_tlist
))
2183 (or (tvar?
(datum-var temp
))
2184 (member (caar (datum-var temp
)) '(mexpt %log
) :test
#'eq
)))
2185 (push (list (datum-var temp
) (trunc-stack temp
)
2186 (exp-pt temp
) (switches temp
))
2188 (merror (intl:gettext
"taylor: ~M cannot be a variable.") (datum-var temp
)))
2191 ;; We must take the max truncation level when $maxtayorder
2192 ;; is T, cf. JPG's bug of 9/15/82.
2193 (when (e> (current-trunc (car a_tlist
)) (current-trunc temp
))
2194 (setf (current-trunc temp
) (current-trunc (car a_tlist
))))
2195 (unless (e> (current-trunc (car a_tlist
)) (current-trunc temp
))
2196 (setf (current-trunc temp
) (current-trunc (car a_tlist
)))))
2197 (unless (alike1 (exp-pt temp
) (exp-pt (car a_tlist
)))
2198 (merror (intl:gettext
"taylor: cannot combine expressions expanded at different points.")))
2199 (setf (switches temp
)
2200 (union* (switches temp
) (switches (car a_tlist
)))))))))
2202 (defun compattlist (list)
2203 (do ((l list
(cdr l
)))
2205 (or (alike1 (exp-pt (get-datum (datum-var (car l
)))) (exp-pt (car l
)))
2209 (let ((last-exp e
)) ;; lexp-non0 should be bound here when needed
2210 (cond ((assolike e tlist
) (var-expand e
1 () ))
2211 ((or (mnump e
) (atom e
) (mfree e tvars
))
2212 (if (or (e> (rczero) (current-trunc mainvar-datum
))
2214 (pszero (data-gvar-o mainvar-datum
)
2215 (current-trunc mainvar-datum
))
2216 (if (and taylor_simplifier
(not (atom e
)))
2217 (let ((e-simp (prep1 (funcall taylor_simplifier e
))))
2218 (when (and (rczerop e-simp
) (not (member e-simp zerolist
:test
#'eq
)))
2222 ((null (atom (caar e
))) (merror "TAYLOR2: internal error."))
2224 (if (and (compatvarlist varlist
(mrat-varlist e
)
2225 genvar
(mrat-genvar e
))
2226 (compattlist (mrat-tlist e
)))
2228 (let ((exact-poly () )) (re-taylor e
))))
2229 ((eq (caar e
) 'mplus
) (tsplus (cdr e
)))
2230 ((eq (caar e
) 'mtimes
) (tstimes (cdr e
)))
2231 ((eq (caar e
) 'mexpt
) (tsexpt (cadr e
) (caddr e
)))
2232 ((eq (caar e
) '%log
) (tslog (cadr e
)))
2233 ((and (or (known-ps (caar e
)) (get (caar e
) 'tay-trans
))
2234 (not (member 'array
(cdar e
) :test
#'eq
))
2235 (try-expansion (if (cddr e
) (cdr e
) (cadr e
))
2238 (cond ((get (subfunname e
) 'spec-trans
)
2239 (funcall (get (subfunname e
) 'spec-trans
) e
))
2240 ((known-ps (subfunname e
))
2241 (try-expansion (caddr e
) (cadr e
))))) )
2242 ((and (member (caar e
) '(%sum %product
) :test
#'eq
)
2243 (mfreel (cddr e
) tvars
))
2244 (tsprsum (cadr e
) (cddr e
) (caar e
)))
2245 ;; If TSDIFF returns NIL, it means that the derivative is wrt to some variable
2246 ;; other than one of the taylor variables.
2247 ;; If so, keep going and handle E as if it were a general expression.
2248 ((and (eq (caar e
) '%derivative
) (tsdiff (cadr e
) (cddr e
) e
)))
2249 ((or (eq (caar e
) '%at
)
2250 (do ((l (mapcar 'car tlist
) (cdr l
)))
2252 (or (free e
(car l
)) (return ()))))
2254 (t (let ((exact-poly () )) ; Taylor series aren't exact
2255 (let ((ee (diff-expand e tlist
)))
2256 (cond ((equal (sratsimp e
) (sratsimp ee
)) (return-from taylor2
(prep1 ee
)))
2257 (t (taylor2 ee
)))))))))
2259 (defun compatvarlist (a b c d
)
2261 ((or (null b
) (null c
) (null d
)) () )
2262 ((alike1 (car a
) (car b
))
2263 (if (not (eq (car c
) (car d
))) ()
2264 (compatvarlist (cdr a
) (cdr b
) (cdr c
) (cdr d
))))
2265 (t (compatvarlist a
(cdr b
) c
(cdr d
)))))
2268 (defun re-taylor (mrat)
2269 (let ((old-tlist (mrat-tlist mrat
)) (old-varlist (mrat-varlist mrat
))
2270 (old-genvar (mrat-genvar mrat
)) old-ivars
)
2271 (declare (special old-tlist old-ivars
))
2272 ;; Put back the old disrpes so rcdisrep's will work correctly.
2273 (mapc #'(lambda (g v
) (putprop g v
'disrep
)) old-genvar old-varlist
)
2274 (setup-multivar-disrep mrat
)
2275 (setq old-ivars
(mapcar #'(lambda (g v
) (cons g v
))
2276 old-genvar old-varlist
))
2277 (prog1 (re-taylor-recurse (mrat-ps mrat
))
2278 ;; Restore the correct disreps.
2279 (mapc #'(lambda (g v
) (putprop g v
'disrep
)) genvar varlist
)
2280 (setup-multivar-disrep () ))))
2282 (defun re-taylor-recurse (ps)
2283 (declare (special old-tlist old-ivars
))
2284 (if (not (psp ps
)) (taylor2 (rcdisrep ps
))
2285 (let (var (datum () ))
2286 (setq var
(cdr (assoc (gvar ps
) old-ivars
:test
#'eq
)))
2287 ;; We must treat multivars like 1, since they'll reappear again
2288 ;; when we call taylor2 or var-expand below.
2289 (if (switch 'multivar
(assoc var old-tlist
:test
#'equal
))
2291 (when (setq datum
(var-data var
))
2292 (push-pw datum
(trunc-lvl ps
))))
2294 (do ((terms (terms ps
) (n-term terms
))
2296 (psplus (if (null var
) (re-taylor-recurse (lc terms
))
2297 (pstimes (re-taylor-recurse (lc terms
))
2299 (var-expand (car datum
)
2300 (edisrep (le terms
))
2303 (m^t var
(edisrep (le terms
)))))))
2306 (when datum
(pop-pw datum
))))))
2308 (defun var-expand (var exp dont-truncate?
)
2309 (let (($keepfloat
) ($float
) (modulus))
2310 (setq exp
(prep1 exp
))) ;; exp must be a rational integer
2311 (let ((temp (get-datum var
't
)))
2312 (cond ((null temp
) (merror "VAR-EXPAND: invalid call."))
2313 ((member (exp-pt temp
) '($inf $minf $infinity
) :test
#'eq
)
2314 (cond ((switch '$asymp temp
)
2315 (merror (intl:gettext
"taylor: cannot create an asymptotic expansion at infinity.")))
2316 ((e> (setq exp
(rcminus exp
)) (current-trunc temp
))
2318 (t (make-ps (int-var temp
)
2319 (ncons (if exact-poly
(inf) (current-trunc temp
)))
2321 (if (eq (exp-pt temp
) '$minf
)
2324 ;; multivar expansion does not work at infinity, so
2325 ;; expansion at infinity is handled by above clause even if doing multivar.
2326 ((switch 'multi temp
) ;; multivar expansion
2328 ;; The reason we call var-expand below instead of taylor2
2329 ;; is that we must be sure the call is not truncated to
2330 ;; 0 which would cause an error in psexpt if exp < 0.
2331 ;; For example, this occurred in TAYLOR(X^2/Y,[X,Y],0,2).
2333 ;; Must ensure that we get back a series truncated
2334 ;; to at least what is specified by tlist. This means
2335 ;; we'll have to push-pw unless exp>0 since psexpt'n
2336 ;; kills (exp-1) terms. The bug that discovered this
2337 ;; is taylor(li[2](x+1/2)/x,[x],0,0) missing 2*log(2).
2338 (if (not (e> exp
(rczero)))
2339 (let-pw (get-datum (car (switch 'multi temp
)))
2340 (e+ (current-trunc temp
) (e- (e1- exp
)))
2341 (var-expand (car (switch 'multi temp
)) 1 't
))
2342 (var-expand (car (switch 'multi temp
)) 1 't
))
2343 (cons (list (int-gvar temp
) 1 1) 1))
2344 (taylor2 (exp-pt temp
)))
2346 ((signp e
(exp-pt temp
))
2347 (let ((exp>trunc?
() ))
2348 (if (and (e> exp
(current-trunc temp
)) (setq exp
>trunc?
't
)
2349 (not dont-truncate?
))
2351 (make-ps (int-var temp
)
2352 (ncons (if exact-poly
(inf)
2353 (if exp
>trunc? exp
(current-trunc temp
))))
2354 (ncons (term (if (switch '$asymp temp
) (rcminus exp
)
2358 (make-ps (int-var temp
)
2359 (ncons (if exact-poly
(inf) (current-trunc temp
)))
2360 (ncons (term (if (switch '$asymp temp
)
2364 (taylor2 (exp-pt temp
)))
2367 (defun expand (arg func
)
2368 (or (try-expansion arg func
) (exp-pt-err)))
2370 (defun try-expansion (arg func
)
2371 (prog (funame funord fun-lc argord psarg arg-trunc temp exact-poly
)
2372 ;; We bind exact-poly to () since we don't want psexpt retaining
2373 ;; higher order terms when subst'ing into series (which aren't exact).
2374 ;; Try diff-expanding unknown subsripted functions.
2375 (unless (or (atom func
) (known-ps (caar func
)))
2376 (taylor2 (diff-expand `((mqapply) ,func
,arg
) tlist
)))
2377 (when (setq temp
(get (setq funame
(oper-name func
)) 'tay-trans
))
2378 (return (funcall temp arg func
)))
2379 (let ((lterm (getfun-lt func
)))
2380 (setq funord
(e lterm
) fun-lc
(c lterm
)))
2382 (when (rczerop (or psarg
(setq psarg
(get-lexp arg
(rcone) () ))))
2383 (if (e> (rczero) funord
)
2384 (if (rczerop (setq psarg
(get-lexp arg
(rcone) 't
)))
2386 (go begin-expansion
))
2387 (return (cond ((setq temp
(assoc funame tay-pole-expand
:test
#'eq
))
2388 (funcall (cdr temp
) arg psarg func
))
2389 ((rczerop funord
) fun-lc
)
2391 (when (pscoefp psarg
) (setq psarg
(taylor2 arg
)))
2392 (when (pscoefp psarg
)
2394 (cond ((null (mfree (rdis psarg
) tvars
))
2395 (symbolic-expand arg psarg func
))
2396 ((setq temp
(assoc funame tay-pole-expand
:test
#'eq
))
2397 (funcall (cdr temp
) arg psarg func
))
2399 (if (atom func
) `((,func
) ,(rcdisrep psarg
))
2400 `((mqapply) ,func
,(rcdisrep psarg
)))))))))
2401 (when (e> (rczero) (setq argord
(ps-le psarg
)))
2402 (cond ((not (member funame
'(%atan %asin %asinh %atanh
) :test
#'eq
))
2403 (if (e> (rczero) (ps-le* (setq psarg
(get-lexp arg
(rcone) 't
))))
2405 (go begin-expansion
)))
2407 (if (and (eq funame
'%atan
)
2408 (eq (coef-sign arg
) '$neg
))
2409 (return (psplus (atrigh arg func
) (taylor2 (m- '$%pi
))))
2410 (return (atrigh arg func
))))))
2411 (setq temp
(t-o-var (gvar psarg
)))
2412 (when (e> (e* funord argord
) temp
) (return (rczero)))
2413 ;; the following form need not be executed if psarg is really exact.
2414 ;; The constant problem does not allow one to determine this now,
2415 ;; so we always have to execute this currently.
2416 ;; This really should be
2417 ;; (unless (infp (trunc-lvl psarg)) ... )
2418 ;; Likewise, the infp checks shouldn't be there; have to assume
2419 ;; nothing is exact until constant problem is fixed.
2420 (setq arg-trunc
(if (and (not (infp (trunc-lvl psarg
)))
2421 (e= funord
(rcone)))
2423 (e- temp
(e* (e1- funord
) argord
)))
2424 psarg
(let-pw (get-datum (get-key-var (gvar psarg
)))
2426 (if (or (infp (trunc-lvl psarg
))
2427 (e> arg-trunc
(trunc-lvl psarg
)))
2430 ;; We must recalculate argord since pstrunc may have "picked"
2431 ;; a coeff out of a constant monomial; e.g. this occurs in
2432 ;; taylor(sin(x+y),x,0,0,y,0,1) where psarg is (Y+...)*X^0+...
2433 ;; which truncates to Y+... of order 1.
2434 argord
(ps-le* psarg
))
2435 (if (rczerop argord
)
2436 (cond ((member funame
'(%atan %asin %asinh %atanh
) :test
#'eq
)
2437 (return (atrigh arg func
)))
2438 ((setq temp
(assoc funame const-exp-funs
:test
#'eq
))
2439 (return (funcall (cdr temp
) arg psarg func
)))
2440 ((rczerop (ps-le* (setq psarg
(get-lexp arg
(rcone) 't
))))
2441 (return () )) ; Don't know an addition formula
2442 (t (go begin-expansion
)))
2444 (if (mono-term?
(terms psarg
))
2445 (get-series func
(current-trunc
2446 (get-datum (get-key-var (gvar psarg
))))
2447 (gvar-o psarg
) (ps-le psarg
) (ps-lc psarg
))
2449 (setq temp
(get-series func
2450 (e// temp argord
) (gvar-o psarg
)
2452 (cond ((not (psp temp
)) temp
)
2453 (t (pscsubst1 psarg temp
)))))))))
2455 (defun symbolic-expand (ign psarg func
) ; should be much stronger
2456 (declare (ignore ign
))
2457 (prep1 (simplifya (if (atom func
)
2458 `((,func
) ,(rcdisrep psarg
))
2459 `((mqapply) ,func
,(rcdisrep psarg
)))
2462 (defun expand-sing-trig?
(arg func
)
2463 (cond ((member func
*pscirc
:test
#'eq
) (tay-exponentialize arg func
))
2464 ((member func
*psacirc
:test
#'eq
) (atrigh arg func
))
2465 (t (essen-sing-err))))
2467 (defun trig-const (a arg func
)
2468 (let ((const (ps-lc* arg
)) (temp (cdr (assoc func trigdisp
:test
#'eq
))))
2469 (cond ((and (pscoefp const
)
2470 (member func
'(%tan %cot
) :test
#'eq
)
2471 (multiple-%pi a
(srdis const
) func
)))
2472 (temp (funcall temp
(setq const
(psdisrep const
))
2474 (t (tsexpt `((,(get func
'recip
)) ,(srdis arg
)) -
1)))))
2476 (defun multiple-%pi
(a const func
)
2478 (and (equal ($hipow const
'$%pi
) 1)
2479 ($ratnump
(setq coef
($ratcoef const
'$%pi
1)))
2480 (cond ((numberp coef
) (expand (m- a const
) func
))
2481 ((equal (caddr coef
) 2)
2482 (psminus (expand (m- a const
)
2483 (cond ((eq func
'%tan
) '%cot
)
2484 ((eq func
'%cot
) '%tan
)
2485 (t (merror "MULTIPLE-%PI: internal error in Taylor expansion."))))))))))
2487 (setq *pscirc
'(%cot %tan %csc %sin %sec %cos %coth
2488 %tanh %csch %sinh %sech %cosh
)
2490 *psacirc
'(%acot %atan %acsc %asin %asec %acos %acoth
2491 %atanh %acsch %asinh %asech %acosh
))
2493 (setq const-exp-funs
2494 `((%gamma . gam-const
) ($psi . plygam-const
)
2495 .
,(mapcar #'(lambda (q) (cons q
'trig-const
)) *pscirc
))
2497 trigdisp
'((%sin . psina
+b
) (%cos . pscosa
+b
) (%tan . pstana
+b
)
2498 (%sinh . psinha
+b
) (%cosh . pscosha
+b
) (%tanh . pstanha
+b
))
2500 tay-pole-expand
'((%gamma . plygam-pole
) ($psi . plygam-pole
))
2502 tay-const-expand
; !these should be handled by symbolic-expand
2503 ; be sure to change this with tay-exponentialize!
2504 (append (mapcar #'(lambda (q) (cons q
'tay-exponentialize
)) *pscirc
)
2505 (mapcar #'(lambda (q) (cons q
'tay-exponentialize
)) *psacirc
)))
2507 (mapc #'(lambda (q) (putprop q
'atrig-trans
'tay-trans
))
2508 '(%acos %acot %asec %acsc %acosh %acoth %asech %acsch
))
2510 (defprop mfactorial factorial-trans tay-trans
)
2512 (defun factorial-trans (arg func
)
2513 (declare (ignore func
))
2514 (taylor2 `((%gamma
) ,(m1+ arg
))))
2516 (defprop %signum signum-trans tay-trans
)
2518 ;; signum(x) => x/abs(x)
2519 ;; this gives an error at x=0 and a derivative of 0 elsewhere
2520 (defun signum-trans (arg func
)
2521 (declare (ignore func
))
2522 (taylor2 `((mtimes) ,arg
((mexpt) ((mabs) ,arg
) -
1))))
2524 (defprop %gamma_incomplete gamma-upper-trans tay-trans
)
2525 (defprop $gamma_incomplete gamma-upper-trans tay-trans
)
2526 (defprop %gamma_incomplete_lower gamma-lower-trans tay-trans
)
2527 (defprop $gamma_incomplete_lower gamma-lower-trans tay-trans
)
2529 ;; for gamma_incomplete(s,z)
2530 ;; translate into gamma_incomplete_lower if s>0 and z=0
2532 ;; June 2022: To workaround the bug
2533 ;; integrate(x*exp(-x^2)*sin(x),x,minf,inf)
2534 ;; limit: variable must be a symbol or subscripted symbol; found: sin(x)
2535 ;; I (Barton Willis) surrounded the call to $limit with errcatch with
2536 ;; $errormsg set to nil. This change allows Maxima to find the correct
2537 ;; value of this definite integral. But almost surely there is a bug
2538 ;; somewhere else that calls gamma-upper-trans with faulty arguments.
2539 ;; The real bug should be fixed, but inserting errcatch here is
2541 (defun gamma-upper-trans (arg func
)
2545 (eq ($sign s
) '$pos
)
2546 (let (($errormsg nil
))
2547 (zerop1 (car (errcatch ($limit z
(caar tlist
) (exp-pt (car tlist
))))))))
2548 (taylor2 `((mplus) ((%gamma
) ,s
)
2549 ((mtimes) -
1 ((%gamma_incomplete_lower
) ,s
,z
))))
2550 (taylor2 (diff-expand `((,func
) .
,arg
)
2553 ;; for gamma_incomplete_lower(s,z)
2554 ;; if z=0, use A&S 6.5.29
2558 ;;; gamma_incomplete_lower(s,z) = z^s * > ------------
2562 (defun gamma-lower-trans (arg func
)
2565 (if (zerop1 ($limit z
(caar tlist
) (exp-pt (car tlist
))))
2570 ((mexpt) ((mtimes) -
1 ,z
) k
)
2571 ((mexpt) ((mtimes) ((mfactorial) k
)
2577 (taylor2 (diff-expand `((,func
) .
,arg
)
2581 ;;; Not done properly yet
2583 ;;; (defprop $BETA BETA-TRANS TAY-TRANS)
2585 (defun psina+b
(a b
)
2586 (psplus (pstimes (expand a
'%sin
) (expand b
'%cos
))
2587 (pstimes (expand a
'%cos
) (expand b
'%sin
))))
2589 (defun pscosa+b
(a b
)
2590 (psdiff (pstimes (expand a
'%cos
) (expand b
'%cos
))
2591 (pstimes (expand a
'%sin
) (expand b
'%sin
))))
2593 (defun pstana+b
(a b
)
2594 (setq a
(expand a
'%tan
) b
(expand b
'%tan
))
2595 (pstimes (psplus a b
)
2596 (psexpt (psdiff (rcone) (pstimes a b
))
2599 (defun psinha+b
(a b
)
2600 (psplus (pstimes (expand a
'%sinh
) (expand b
'%cosh
))
2601 (pstimes (expand a
'%cosh
) (expand b
'%sinh
))))
2603 (defun pscosha+b
(a b
)
2604 (psplus (pstimes (expand a
'%cosh
) (expand b
'%cosh
))
2605 (pstimes (expand a
'%sinh
) (expand b
'%sinh
))))
2607 (defun pstanha+b
(a b
)
2608 (setq a
(expand a
'%tanh
) b
(expand b
'%tanh
))
2609 (pstimes (psplus a b
)
2610 (psexpt (psplus (rcone) (pstimes a b
))
2613 (defun atrig-trans (arg func
)
2615 (cond ((eq func
'%acos
)
2616 `((mplus) ,half%pi
((mtimes) -
1 ((%asin
) ,arg
))))
2619 `((mtimes) -
1 $%i
((mplus) ,half%pi
((mtimes) -
1 ((%asin
) ,arg
)))))
2622 `((,(cdr (assoc func
'((%acsc . %asin
) (%asec . %acos
)
2623 (%acot . %atan
) (%acsch . %asinh
)
2624 (%asech . %acosh
) (%acoth . %atanh
)) :test
#'eq
)))
2627 (defun atrigh (arg func
)
2628 (let ((full-log t
) ($logarc t
) (log-1 '((mtimes) $%i $%pi
))
2629 (log%i
'((mtimes) ((rat) 1 2) $%i $%pi
)))
2630 (taylor2 (simplify `((,func
) ,arg
)))))
2632 (defun tay-exponentialize (arg fun
) ; !this should be in symbolic-expand!
2633 (let (($exponentialize t
) ($logarc t
))
2634 (setq arg
(meval `((,fun
) ,arg
))))
2638 (do ((l (cdr l
) (cdr l
))
2639 (ans (taylor2 (car l
))
2640 (psplus ans
(taylor2 (car l
)))))
2643 (defun ts-formula (form var pw
)
2644 (let ((datum (get-datum (get-key-var (car var
)))))
2646 (taylor2 (subst (get-inverse (car var
)) 'sp2var form
)))))
2648 (defmacro next-series
(l) `(cdadr ,l
))
2650 (defun tstimes-get-pw (l pw
)
2651 (do ((l l
(cdr l
)) (vect))
2653 (setq pw
(mapcar #'(lambda (pwq ple
) (e+ pwq ple
))
2654 pw
(setq vect
(ord-vector (cdar l
)))))
2655 (rplacd (car l
) (cons (cdar l
) vect
))))
2657 (defun tstimes-l-mult (a)
2658 (do ((l (cdr a
) (cdr l
)) ($maxtayorder t
)
2659 (a (car a
) (pstimes a
(car l
))))
2665 (or (zfree e
(car l
)) (return () ))))
2667 ;;; The lists posl, negl and zerl have the following format:
2669 ;;; ( (<expression> <expansion> <order vector>) . . . )
2672 (*bind
* ((funl) (expl) (negl) (zerl) (posl)
2673 (pw) (negfl) (temp) (fixl (rcone)))
2674 (dolist (fun l
) ;; find the exponentials
2676 (push (if (free (caddr fun
) (car tvars
)) fun
2677 `((mexpt) $%e
,(m* (caddr fun
)
2678 `((%log
) ,(cadr fun
)))))
2682 (setq expl
(tsexp-comb expl
)) ;; simplify exps
2683 (setq expl
(tsbase-comb expl
))) ;; and again
2684 (setq l
(nconc expl funl
)) ;; now try expanding
2685 (let ((trunc-constants?
() ))
2686 (setq expl
(cons 0 (mapcar #'(lambda (exp)
2687 (cons exp
(taylor2 exp
)))
2689 ;; EXPL is now of the form (0 ( <form> . <taylor2(form)> ) ...)
2690 ;; L points behind the cons considered for destructive updating.
2691 (do ((l expl
) (tem))
2693 (cond ((rczerop (cdadr l
))
2694 ;; Consider taylor((a+1/x)*1/x,x,0,-2). Each factor will be on
2695 ;; zerl. Each factor will also appear to have le = 0 since its
2696 ;; series is 0, which would fool the get-pw routines below if
2697 ;; they tried to handle this case. The easiest fix for now
2698 ;; appears to be to always call get-lexp here, killing this:
2699 (cond ;((null $maxtayorder)
2700 ; (setq zerl (cons (cadr l) zerl))
2701 ; (rplacd l (cddr l)))
2702 ((rczerop (setq tem
(get-lexp (caadr l
) (rcone) ())))
2703 (return (setq zerl
0)))
2704 ('t
(setq posl
(cons (cons (caadr l
) tem
) posl
))
2705 (rplacd l
(cddr l
)))))
2706 ((pscoefp (cdadr l
))
2707 (cond ((mzfree (caadr l
) tvars
) ;must be zfree to permit ratfuns
2708 (setq fixl
(pstimes (cdadr l
) fixl
))
2709 (rplacd l
(cddr l
)))
2710 ((setq zerl
(cons (cadr l
) zerl
))
2711 (rplacd l
(cddr l
)))))
2712 ((rczerop (ps-le (cdadr l
)))
2713 (setq zerl
(cons (cadr l
) zerl
))
2714 (rplacd l
(cddr l
)))
2715 ((e> (ps-le (cdadr l
)) (rczero))
2716 (setq posl
(cons (cadr l
) posl
))
2717 (rplacd l
(cddr l
)))
2718 ('t
(setq l
(cdr l
)))))
2719 (when (equal zerl
0) (return (rczero)))
2720 (setq negl
(cdr expl
) temp
(ord-vector fixl
))
2721 (mapcar #'(lambda (x) (and (e> (rczero) x
) (setq negfl t
))) temp
)
2722 (tstimes-get-pw zerl temp
)
2723 (setq pw
(tstimes-get-pw posl
(tstimes-get-pw negl temp
)))
2726 (mapcar #'(lambda (x)
2727 (prog2 (mapcar #'(lambda (datum lel pwl
)
2729 (e+ (current-trunc datum
)
2733 (mapcar #'(lambda (datum) (pop-pw datum
))
2735 (nconc posl zerl negl
)))
2736 (setq posl
(nconc (mapcar 'cadr posl
) (mapcar 'cadr zerl
)
2737 (mapcar 'cadr negl
))))
2738 (setq posl
(tstimes-l-mult posl
))
2739 (let ((ans (cond ((null fixl
) posl
)
2741 ('t
(pstimes fixl posl
)))))
2742 (if $maxtayorder ans
(pstrunc ans
)))))
2744 ;;; This routine transforms a list of exponentials as follows:
2746 ;;; a^c*b^(n*c) ===> (a*b^n)^c, where n is free of var.
2748 ;;; This transformation is only applicable when c is not free of var.
2750 (defun tsexp-comb (l) ;; ***** clobbers l *****
2751 (setq l
(cons '* l
))
2752 (do ((a l
) (e)) ;; updated by a rplacd or cdr.
2753 ((null (cddr a
)) (cdr l
)) ;; get rid of the *
2754 (rplaca (cddadr a
) (setq e
($ratsimp
(caddr (cadr a
)))))
2755 ;; Must delete e^0 lest we divide by the 0 below. RWG's byzero bug
2756 ;; of 3/1/78 used to cause this.
2757 (if (equal e
0) (rplacd a
(cddr a
))
2758 (if (mfree (caddr (cadr a
)) tvars
) (pop a
)
2759 (do ((b (cddr a
) (cdr b
)) (n))
2760 ((null b
) (setq a
(cdr a
)))
2761 (when (mfree (setq n
($ratsimp
(m// (caddar b
)
2764 (rplaca b
(list '(mexpt simp
)
2766 (m^
(cadar b
) n
)) ;; b^n
2768 (rplacd a
(cddr a
)) ;; delete a^c
2771 ;;; This routine transforms a list of exponentials as follows:
2773 ;;; a^b*a^c ===> a^(b+c),
2775 ;;; this is only necessary when b and c depend on "var."
2777 (defun tsbase-comb (l) ;;; *******clobbers l********
2778 (setq l
(cons '* l
))
2779 (do ((a l
)) ;;; updated by a rplacd or cdr
2780 ((null (cddr a
)) (cdr l
))
2781 (do ((b (cddr a
) (cdr b
)))
2782 ((null b
) (pop a
)) ;;; did not return early so pop.
2783 (when (alike1 (cadar b
) (cadadr a
))
2784 (rplaca b
(m^
(cadar b
) (m+ (caddar b
) (caddr (cadr a
)))))
2789 (cond ((and (atom b
) (mnump e
)
2791 (not (eq (exp-pt (get-datum b
)) '$minf
)))
2792 ;; one could remove this clause and let this case be handled by tsexpt1
2793 (var-expand b e
() ))
2794 ((mfree e tvars
) (tsexpt1 b e
))
2795 ((eq '$%e b
) (tsexpt-red (list e
)))
2796 (t (tsexpt-red (list (list '(%log
) b
) e
)))))
2798 (defun tsexpt-red (l)
2799 (*bind
* ((free) (nfree) (full-log) ($logarc t
) (expt) (ps) (e)
2800 (log-1 '((mtimes) $%i $%pi
))
2801 (log%i
'((mtimes) ((rat) 1 2) $%i $%pi
)))
2802 (declare (special e
))
2803 a
(do ((l l
(cdr l
)))
2805 (cond ((mtimesp (car l
)) (setq l
(append l
(cdar l
))))
2806 ((mfree (car l
) tvars
) (push (car l
) free
))
2807 (t (push (car l
) nfree
))))
2808 (cond ((or (cdr nfree
) (atom (car nfree
))) )
2809 ((eq (caaar nfree
) '%log
)
2810 (return (tsexpt1 (cadar nfree
) (m*l free
))))
2811 ((member (caaar nfree
) *psacirc
:test
#'eq
)
2812 (setq l
(ncons (simplifya ;; simplify after removing simp flag
2813 (cons (ncons (caaar nfree
)) (cdar nfree
))
2817 ;; Must have truncs > 0 so that logs in the expt aren't trunc'd.
2818 ;; E.g, consider taylor(x^(x-1),x,0,-1).
2819 (tlist-mapc d
(push-pw d
(emax (current-trunc d
) (rcone))))
2820 (setq ps
(taylor2 (setq expt
(m*l
(append nfree free
)))))
2821 (tlist-mapc d
(pop-pw d
))
2822 ;; Here we must account for the truncation gain or lossage that
2823 ;; is encountered in exp(c*log(x)+y) -> x^c*exp(y).
2824 (let ((c0 (if (pscoefp ps
) ps
(psterm (terms ps
) (rczero))))
2826 (unless (rczerop c0
)
2827 (setq ord-e^c0
(ord-vector (setq e^c0
(psexpt-fn c0
))))
2828 ;; Must emax with 0 so that new singular kernals won't be trunc'd
2829 ;; e.g exp(1/x+...) to degree -2 should be exp(-1/x)+...
2830 ;; Also try taylor(screwa,x,0,-2).
2831 (mapc #'(lambda (d o
) (push-pw d
(emax (e- (current-trunc d
) o
)
2834 (setq ps
(psdiff (taylor2 expt
) c0
)))
2835 (setq ps
(psexpt-fn ps
))
2837 (tlist-mapc d
(pop-pw d
))
2838 (setq ps
(pstimes e^c0 ps
)))
2841 ;; Taylor's b^e, where e is independent of tvars.
2843 (defun tsexpt1 (b e
)
2845 (setq e
(let ((modulus () )) ; Don't mod exponents! See WGM's bug
2846 (prog2 (mapcar ; of 3/6/83 for an example.
2849 (emax (current-trunc datum
) (rczero))))
2852 (mapcar #'(lambda (datum) (pop-pw datum
)) tlist
)))
2855 pw
(if (psp tb
) (current-trunc (get-datum
2856 (get-key-var (gvar tb
))))
2857 ;; Constant problem kludge.
2858 (if (rczerop tb
) (current-trunc (car tlist
)) (rczero))))
2859 (if (floatp (car s
))
2860 (setq s
(maxima-rationalize (quot (car s
) (cdr s
)))))
2861 ;; We must ensure that the lc is non-zero since it will be inverted in
2863 (setq tb
(strip-zeroes tb
't
))
2865 (when (or ;; When 1 > s we need more terms since -le*(s-1) > 0.
2867 (and (e> (rczero) pw
) (e> s
(rcone))))
2868 (setq tb
(get-lexp b
() 't
)))
2869 (setq le
(ps-le* tb
)))
2870 ((psp tb
) (setq le
(ps-le tb
)))
2871 (t (return (rcexpt tb e
))))
2872 (and (e> (e* s le
) pw
) (null $maxtayorder
) (return (rczero)))
2873 (setq s
(e- pw
(e* le
(e1- s
))))
2874 ;(setq le (increment-truncs tb))
2879 (pstrunc1 tb
(list (cons (gvar tb
) s
))))
2880 ;; because of constants not retaining info, have to
2881 ;; just keep the constant here
2882 (cond ((not (psp tb
)) tb
)
2883 (t (let-pw (get-datum (get-key-var (gvar tb
))) s
(strip-zeroes (taylor2 b
) 't
)))))
2886 ;;; the method of calculating truncation levels below is incorrect.
2887 ;;; (i.e. increment-truncs & decrement-truncs, also used above)
2888 ;;; Examples which exhibit this incorrectness are:
2889 ;;; taylor(log(sin(y)+x),x,0,2,y,0,1) is missing a y/6*x and -1/6*x^2
2890 ;;; taylor(log(sin(z)+sin(y)+x),x,0,f1,y,0,3,z,0,5) misses a z^5*y^3 term.
2892 ;;; TSLOG must find the lowest degree term in the expansion of the
2893 ;;; log arg, then expand with the orders of all var's in this low term
2894 ;;; incremented by their order in this low term. Note that this is
2895 ;;; only necessary for var's with ord > 0, since otherwise we have
2896 ;;; already expanded to a higher ord than required. Also we must
2897 ;;; not do this for var's with trunc < 0, since this may incorrectly
2898 ;;; truncate terms which should end up as logs.
2900 (defun increment-truncs (ps)
2901 (do ((ps ps
(ps-lc ps
)) (trunc (t-o-var (gvar ps
))) (data () ))
2903 (when (e> (ps-le ps
) (rczero))
2904 (push (assoc (get-key-var (gvar ps
)) tlist
:test
#'eq
) data
)
2905 (push-pw (car data
) (e+ (e* (e+ trunc
(rctwo)) (ps-le ps
))
2906 (current-trunc (car data
))))
2907 (setq trunc
(e+ trunc
(current-trunc (car data
))))
2910 (defun decrement-truncs (data)
2911 (mapc #'(lambda (data) (pop-pw data
)) data
))
2914 (let ((psarg (taylor2 arg
)) datum
)
2915 (when (rczerop psarg
) (setq psarg
(get-lexp arg
() 't
)))
2916 ;; We must ensure that the lc is non-zero since it will be inverted in pslog
2917 (setq psarg
(strip-zeroes psarg
't
))
2918 (do ((ps psarg
(ps-lc ps
)) (shift (rcone) (e* shift
(rctwo))))
2921 (when (rczerop (setq psarg
(taylor2 arg
)))
2922 (setq psarg
(get-lexp arg
() 't
)))
2923 (mapc #'(lambda (data) (pop-pw data
)) datum
))
2925 (push (get-datum (get-key-var (gvar ps
))) datum
)
2926 (if (and (e> (ps-le ps
) (rczero))
2927 (e> (current-trunc (car datum
)) (rczero)))
2928 (push-pw (car datum
) (e+ (e* shift
(ps-le ps
))
2929 (current-trunc (car datum
))))
2932 ;; When e-start is non-null we start expanding at order e-start, ... , 2^m,
2933 ;; then 2^m*pow, instead of the normal sequence pow, ... , 2^m*pow
2934 ;; (where m = $taylordepth, pow = ord of var). This is done because it is
2935 ;; usually much more efficient for large, non-trivial expansions when we only
2936 ;; want the lowest order term.
2938 (defun get-lexp (exp e-start zerocheck?
)
2944 (tlist-mapc d
(push-pw d
(or e-start
(emax (orig-trunc d
) (rcone)))))
2945 (do ((psexp) (i (1+ $taylordepth
) (1- i
)))
2947 (tlist-mapc d
(pop-pw d
))
2951 (unless silent-taylor-flag
(zero-warn exp
))
2953 (declare (fixnum i
))
2954 (cond ((and (rczerop (setq psexp
(if zerocheck?
2955 (strip-zeroes (taylor2 exp
) 't
)
2957 (not (member exp zerolist
:test
#'eq
))) )
2958 ;; Info not needed yet.
2959 ;; ((and lexp-non0 (rczerop (le (terms psexp)))
2960 ;; (mono-term? (terms psexp))))
2961 (t (tlist-mapc d
(pop-pw d
))
2963 (cond ((and (= i
1) e-start
)
2964 (setq e-start
() i
2)
2965 (tlist-mapc d
(push-pw d
(prog1 (e* (orig-trunc d
) (current-trunc d
))
2967 (t (tlist-mapc d
(push-pw d
(prog1 (e* (rctwo) (current-trunc d
))
2971 (or (equal x
1) (equal x
1.0)))
2973 (defun [max-trunc
] ()
2974 (do ((l tlist
(cdr l
)) (emax (rczero)))
2975 ((null l
) (1+ (truncate (car emax
) (cdr emax
))))
2976 (when (e> (current-trunc (car l
)) emax
)
2977 (setq emax
(orig-trunc (car l
))))))
2979 (defun tsprsum (f l type
)
2980 (if (mfree f tvars
) (newsym f
)
2981 (let ((li (ncons (car l
))) (hi (caddr l
)) (lv (ncons (cadr l
))) a aa
2982 ($maxtayorder
() ));; needed to determine when terms are 0
2983 (if (and (numberp (car lv
)) (numberp hi
) (> (car lv
) hi
))
2984 (if (eq type
'%sum
) (taylor2 0) (taylor2 1))
2986 (if (eq type
'%sum
) (setq type
'()))
2987 (do ((m (* ([max-trunc
]) (ash 1 $taylordepth
)))
2989 (ans (taylor2 (maxima-substitute (car lv
) (car li
) f
))))
2990 ((equal hi
(car lv
)) ans
)
2991 (rplaca lv
(m1+ (car lv
)))
2992 ;; A cheap heuristic to catch infinite recursion when
2993 ;; possible, should be improved in the future
2994 (if (> k m
) (exp-pt-err)
2995 (setq a
;(mlet li lv (taylor2 (setq aa (meval f))))
2996 (taylor2 (maxima-substitute (car lv
) (car li
) f
))))
2998 (if (and (1p (car a
)) (1p (cdr a
)) (not (1p aa
)))
3000 (setq ans
(pstimes a ans
)))
3001 (if (and (rczerop a
) (not (signp e aa
)))
3003 (setq ans
(psplus ans a
))))))))))
3005 (defun tsdiff (e l check
)
3006 (*bind
* ((n) (v) (u))
3007 (do ((l l
(cddr l
)))
3009 (if (and (atom (car l
)) (numberp (cadr l
))
3010 (assoc (car l
) tlist
:test
#'eq
))
3011 (setq n
(cons (cadr l
) n
) v
(cons (car l
) v
))
3012 (setq u
(cons (car l
) (cons (cadr l
) u
)))))
3013 ;; If N is null, it means that the derivative is wrt to some variable
3014 ;; other than one of the taylor variables.
3015 ;; If so, give up; let the caller handle it.
3017 (if u
(setq e
(meval (cons '($diff
) (cons e l
)))))
3018 (setq l
(mapcar #'(lambda (x) (get-datum x
)) v
))
3019 (mapcar #'(lambda (datum pw
)
3020 (push-pw datum
(e+ (current-trunc datum
) (prep1 pw
))))
3022 (setq e
(taylor2 e
))
3023 (mapc #'(lambda (datum) (pop-pw datum
)) l
)
3024 (do ((vl v
(cdr vl
))
3029 (mapc #'(lambda (a b
)
3030 (putprop a
(prep1 (sdiff b
(car v
)))
3033 (setq e
(psdp e
))))))
3036 (defun no-sing-err (x) ;; try to catch all singularities
3038 (declare (special errorsw
))
3039 (let ((ans (catch 'errorsw
(eval x
))))
3040 (if (eq ans t
) (unfam-sing-err) ans
))))
3042 ;; evaluate deriv at location var=pt
3043 (defun eval-deriv (deriv var pt
)
3045 (declare (special errorsw
))
3046 (let ((ans (no-sing-err `(meval '(($at
) ,deriv
((mequal) ,var
,pt
))))))
3049 (defun check-inf-sing (pt-list) ; don't know behavior of random fun's @ inf
3050 (and (or (member '$inf pt-list
:test
#'eq
) (member '$minf pt-list
:test
#'eq
))
3053 (defun diff-expand (exp l
) ;l is tlist
3054 (check-inf-sing (mapcar (function caddr
) l
))
3057 (setq exp
(diff-expand exp
(cdr l
)))
3058 (do ((deriv (sdiff exp
(caar l
)) (sdiff deriv var
))
3060 (coef 1 (* coef
(1+ cnt
)))
3062 (pt (exp-pt (car l
)))
3063 (lim (rcdisrep (current-trunc (car l
))))
3064 (ans (list (no-sing-err `(meval '(($at
) ,exp
((mequal) ,(caar l
) ,(exp-pt (car l
)))))))
3065 (cons `((mtimes) ((rat simp
) 1 ,coef
)
3066 ,(eval-deriv deriv var pt
)
3067 ((mexpt) ,(sub* var pt
) ,cnt
))
3069 ((or (great cnt lim
) (equal deriv
0)) (cons '(mplus) ans
))))))
3071 ;;; subtitle disreping routines
3074 (if (= (cdr e
) 1) (car e
) (list '(rat) (car e
) (cdr e
))))
3076 (defun striptimes (a)
3077 (if (mtimesp a
) (cdr a
) (ncons a
)))
3080 (let (($psexpand
() )) ; Called only internally, no need to expand.
3082 (cons (list 'mrat
'simp varlist genvar tlist
'trunc
)
3086 (let ((varlist (mrat-varlist r
)) (genvar (mrat-genvar r
)))
3087 (mapc #'(lambda (exp genv
) (putprop genv exp
'disrep
))
3089 (setup-multivar-disrep r
)
3090 ;; This used to return 0 if psdisrep returned () but this is wrong
3091 ;; since taylor(false,x,0,0) would lose. If psdisrep really wants to
3092 ;; return () for 0 then we will probably find out soon.
3093 (if (eq $psexpand
'$multi
) (psdisexpand (cdr r
))
3094 (psdisrep (cdr r
)))))
3098 (psdisrep+ (psdisrep2 (terms p
) (getdisrep (gvar-o p
)) (trunc-lvl p
))
3099 (if (or $psexpand
(trunc-lvl p
)) '(mplus trunc
)
3103 (defun psdisrep^
(n var
)
3104 ;; If var = () then it is an internal var generated in a multivariate
3105 ;; expansion so it shouldn't be displayed. If var = 1 then it probably
3106 ;; resulted from the substitution in srdisrep, so it depends on an
3107 ;; internal var and likewise shouldn't be displayed.
3108 (cond ((or (rczerop n
) (null var
) (equal var
1)) 1)
3109 ((equal n
(rcone)) var
)
3110 ((and ps-bmt-disrep
(mexptp var
) (equal (caddr var
) -
1))
3111 (psdisrep^
(e- n
) (cadr var
)))
3112 ('t
`((mexpt ratsimp
) ,var
,(edisrep n
)))))
3114 ;;; There used to be a hack below that would print a series consisting
3115 ;;; of merely one term as exact polynomial (i.e. no trailing "..."'s).
3116 ;;; This is, of course, wrong but the problem with the fix is that
3117 ;;; now exact things like taylor(y*x,x,0,f1,y,0,1) will display like
3118 ;;; (y+...) x+... because of the problem with $MAXTAYORDER being internally
3119 ;;; bound to ()---which causes exact things to look inexact, such as
3120 ;;; x and y above. See the comment above taylor* for the $MAXTAYORDER problem.
3122 (defun psdisrep+ (p plush
&aux lowest-degree-term
)
3123 (if;; An exact sum of one arg is just that arg.
3124 (and (null (cdr p
)) (eq (cadr plush
) 'exact
))
3127 ;; Since the DISPLAY package prints trunc'd sum's arguments
3128 ;; from right to left we must put the terms of any constant term
3129 ;; in decreasing order. Note that only a constant (wrt to the
3130 ;; mainvar) term can be a term which is a sum.
3131 (when (mplusp (setq lowest-degree-term
(car (last p
))))
3132 (rplacd lowest-degree-term
(nreverse (cdr lowest-degree-term
))))
3135 (defun psdisrep* (a b
)
3136 (cond ((equal a
1) b
)
3138 (t (cons '(mtimes ratsimp
)
3139 (nconc (striptimes a
) (striptimes b
))))))
3141 (defun psdisrep2 (p var trunc
)
3142 (if (or $ratexpand $psexpand
) (psdisrep2expand p var
)
3143 (do ((a () (cons (psdisrep* (psdisrep (lc p
)) (psdisrep^
(le p
) var
))
3146 ((or (null p
) (e> (le p
) trunc
)) a
))))
3148 (defun psdisrep2expand (p var
)
3150 (l () (nconc (psdisrep*expand
(psdisrep (lc p
)) (psdisrep^
(le p
) var
))
3154 (defun psdisrep*expand
(a b
)
3155 (cond ((equal a
1) (list b
))
3156 ((equal b
1) (list a
))
3158 (list (cons '(mtimes ratimes
) (nconc (striptimes a
) (striptimes b
)))))
3159 ('t
(mapcar #'(lambda (z) (psdisrep* z b
))
3163 (defun psdisexpand (p)
3164 (let ((ans (ncons ())))
3165 (declare (special ans
)) ;used in pans-add
3166 (psdisexcnt p
() (rczero))
3169 (mapcar #'(lambda (x) (cond ((not (cddr x
)) (cadr x
))
3170 (t (cons '(mplus trunc
) (cdr x
)))))
3172 (cond ((not (cdr ans
)) (car ans
))
3173 (t (cons '(mplus trunc
) ans
)))))
3175 (defun psdisexcnt (p l n
)
3177 (do ((var (getdisrep (gvar-o p
))) (ll (terms p
) (n-term ll
)))
3179 (if (rczerop (le ll
)) (psdisexcnt (lc ll
) l n
)
3181 (cons (psdisrep^
(le ll
) var
) l
)
3183 (psans-add (cond ((not l
) (rcdisrep p
))
3184 (t (psdisrep* (rcdisrep p
)
3185 (cond ((not (cdr l
)) (car l
))
3186 (t (cons '(mtimes trunc
) l
))))))
3189 (defun psans-add (exp n
)
3190 (declare (special ans
)) ;bound in psdisexpand
3191 (do ((l ans
(cdr l
)))
3192 ((cond ((null (cdr l
)) (rplacd l
(ncons (list n exp
))))
3193 ((e= (caadr l
) n
) (rplacd (cadr l
) (cons exp
(cdadr l
))))
3194 ((e> (caadr l
) n
) (rplacd l
(cons (list n exp
) (cdr l
))))))))
3196 (defun srconvert (r)
3197 (cond ((not (atom (caadr (cdddar r
))))
3198 (cons (car r
) (psdisextend (cdr r
))))
3200 (*bind
* ((trunclist (cadr (cdddar r
)))
3205 (gens (cadddr (car r
))))
3206 (setq gps
(mapcar #'cons gens vs
))
3207 (do ((tl (cdr trunclist
) (cddr tl
)))
3208 ((null tl
) (cons (list 'mrat
'simp vs gens tlist
'trunc
) (srconvert1 (cdr r
))))
3209 (setq temp
(cdr (assoc (car tl
) gps
:test
#'eq
)))
3210 (cond ((null (member (car tl
) (cdr trunclist
) :test
#'eq
)))
3211 ((mplusp temp
) (merror "SRCONVERT: internal error."))
3214 (cons (list* temp
(tay-order (cadr tl
)) 0 nil
3215 (cons (car tl
) (symbol-value (car tl
))))
3218 (defun srconvert1 (p)
3219 (cond ((not (member (car p
) genvar
:test
#'eq
)) p
)
3221 (do ((l (cdr p
) (cddr l
))
3222 (a nil
(cons (term (prep1 (car l
)) (srconvert1 (cadr l
))) a
)))
3224 (make-ps (cons (car p
) (symbol-value (car p
)))
3225 (tay-order (zl-get trunclist
(car p
))) a
))))))
3227 ;;; subtitle error handling
3229 (defun tay-error (msg exp
)
3230 (if silent-taylor-flag
(throw 'taylor-catch
())
3232 (merror "taylor: ~A~%~M" msg exp
)
3233 (merror "taylor: ~A" msg
))))
3235 (defun exp-pt-err ()
3236 (tay-err (intl:gettext
"unable to expand at a point specified in:")))
3238 (defun essen-sing-err ()
3239 (tay-err (intl:gettext
"encountered an essential singularity in:")))
3241 (defun unfam-sing-err ()
3242 (tay-err (intl:gettext
"encountered an unfamiliar singularity in:")))
3244 (defun infin-ord-err ()
3245 (tay-err (intl:gettext
"expansion to infinite order?")))
3247 (defun tay-depth-err ()
3248 (tay-err (intl:gettext
"'taylordepth' exceeded while expanding:")))
3250 ;;; Subtitle TAYLORINFO
3252 (defun taylor-trunc (q)
3253 (setq q
(current-trunc q
))
3254 (cond ((null q
) '$inf
)
3255 ((equal (cdr q
) 1) (car q
))
3256 (t `((rat) ,(car q
) ,(cdr q
)))))
3258 (defun taylor-info (q)
3259 (let ((acc-var nil
) (acc-pt nil
) (acc-ord nil
) (qk) (acc))
3260 (cond ((null q
) nil
)
3263 (cond ((and (fourth qk
) (consp (fourth qk
)) (eq (caar (fourth qk
)) 'multivar
)) nil
)
3264 ((and (fourth qk
) (consp (fourth qk
)) (eq (caar (fourth qk
)) 'multi
))
3265 (while (and (fourth qk
) (consp (fourth qk
)) (eq (caar (fourth qk
)) 'multi
))
3267 (push (taylor-trunc qk
) acc-ord
)
3268 (push (exp-pt qk
) acc-pt
)
3269 (push (datum-var qk
) acc-var
)
3271 (push '(mlist) acc-ord
)
3272 (push '(mlist) acc-pt
)
3273 (push '(mlist) acc-var
)
3274 (setq q
(taylor-info q
))
3275 (if (null q
) (list acc-var acc-pt acc-ord
) (append q
(list acc-var acc-pt acc-ord
))))
3278 (setq acc
(if (and (fourth qk
) (consp (fourth qk
)) (eq '$asympt
(caar (fourth qk
))))
3279 (list '$asympt
) nil
))
3280 (push (taylor-trunc qk
) acc
)
3281 (push (exp-pt qk
) acc
)
3282 (push (datum-var qk
) acc
)
3284 (setq q
(taylor-info q
))
3285 (if (null q
) (list acc
) (append q
(list acc
)))))))))
3287 (defmfun $taylorinfo
(x)
3288 (if (and (consp x
) (member 'trunc
(first x
) :test
#'eq
))
3289 (cons '(mlist) (taylor-info (mrat-tlist x
)))
3294 ;;; Lisp let-pw Indent:2
