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 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module residu
)
15 (load-macsyma-macros rzmac
)
17 (declare-top (special $breakup $noprincipal varlist
18 leadcoef var
*roots
*failures wflag nn
*
19 sn
* sd
* genvar dn
* zn
))
22 ;; Compute the poles (roots) of the polynomial D and return them.
23 ;; Divide these roots into several parts: Those in REGION, REGION1,
24 ;; and everywhere else. These are returned in a list. (In a more
25 ;; modern style, we'd probably return them in 4 different values.)
27 ;; The regions are determined by functions REGION and REGION1, which
28 ;; should return non-NIL if the root is in the given region.
30 ;; The description below applies if *SEMIRAT* is NIL. If *SEMIRAT* is
31 ;; non-NIL, somewhat different results are returned. I (rtoy) am not
32 ;; exactly sure what *SEMIRAT* is intended to mean.
34 ;; The first part of the list of the form ((r1 (x - r1)^d1) (r2 (x -
35 ;; r2)^d2) ...) where r1, r2 are the roots, d1, d2 are the
36 ;; multiplicity of each root, and x is the variable.
38 ;; The second part is a list of the repeated roots in REGION. Each
39 ;; element of the list is of the form (r d) where r is the root and d
40 ;; is the multiplicity.
42 ;; The third part is a list of the simple roots in REGION.
44 ;; Finally, the fourth part is NIL, unless *semirat* is T.
45 (defun polelist (d region region1
)
46 (prog (roots $breakup r rr ss r1 s pole wflag cf
)
48 (setq leadcoef
(polyinx d var
'leadcoef
))
49 (setq roots
(solvecase d
))
50 (if (eq roots
'failure
) (return ()))
51 ;; Loop over all the roots. SOLVECASE returns the roots in the form
59 (> (+ (length s
) (length r
))
60 (+ (length ss
) (length rr
))))
61 ;; Return CF, repeated roots (*semirat*), simple
62 ;; roots (*semirat*), roots in region 1.
63 (return (list cf rr ss r1
)))
65 ;; Return CF, repeated roots, simple roots, roots in region 1.
66 (return (list cf r s r1
)))))
68 ;; Set POLE to the actual root and set D to the
69 ;; multiplicity of the root.
70 (setq pole
(caddar roots
))
73 ;; Is it possible for LEADCOEF to be NIL ever?
75 ;; Push (pole (x - pole)^d) onto the list CF.
78 (m^
(m+ var
(m* -
1 pole
))
81 ;; Don't change the order of clauses here. We want to call REGION and then REGION1.
82 (cond ((funcall region pole
)
83 ;; The pole is in REGION
85 ;; A simple pole, so just push the pole onto the list S.
88 ;; A multiple pole, so push (pole d) onto the list R.
89 (push (list pole d
) r
))))
90 ((funcall region1 pole
)
91 ;; The pole is in REGION1
92 (cond ((not $noprincipal
)
93 ;; Put the pole onto the R1 list. (Don't know what
97 ;; Return NIL if we get here.
100 ;; (What does *SEMIRAT* mean?) Anyway if we're here, the
101 ;; pole is not in REGION or REGION1, so push the pole onto
102 ;; SS or RR depending if the pole is repeated or not.
105 (t (push (list pole d
) rr
)))))
106 ;; Pop this root and multiplicity and move on.
107 (setq roots
(cddr roots
))
111 (cond ((not (among var e
)) nil
)
112 (t (let (*failures
*roots
)
114 (cond (*failures
'failure
)
118 ;; Compute the sum of the residues of n/d.
119 (defun res (n d region region1
)
120 (let ((pl (polelist d region region1
))
121 dp a b c factors leadcoef
)
125 (setq factors
(car pl
))
127 ;; PL now contains the list of the roots in region, roots in
128 ;; region1, and everything else.
131 (setq dp
(sdiff d var
))))
133 ;; Compute the sum of the residues of n/d for the
134 ;; multiple roots in REGION.
135 (setq a
(m+l
(residue n
(cond (leadcoef factors
)
140 ;; Compute the sum of the residues of n/d for the simple
141 ;; roots in REGION1. Since the roots are simple, we can
142 ;; use RES1 to compute the residues instead of the more
143 ;; complicated $RESIDUE. (This works around a bug
144 ;; described in bug 1073338.)
146 (setq b
(m+l
(mapcar #'(lambda (pole)
147 ($residue
(m// n d
) var pole
))
149 (setq b
(m+l
(res1 n dp
(cadr pl
)))))
152 ;; Compute the sum of the residues of n/d for the roots
153 ;; not in REGION nor REGION1.
154 (setq c
(m+l
(mapcar #'(lambda (pole)
155 ($residue
(m// n d
) var pole
))
158 ;; Return the sum of the residues in the two regions and the
159 ;; sum of the residues outside the two regions.
160 (list (m+ a b
) c
)))))
162 (defun residue (zn factors pl
)
164 (mapcar #'(lambda (j)
165 (destructuring-let (((factor1 factor2
) (remfactor factors
(car j
) zn
)))
166 (resm0 factor1 factor2
(car j
) (cadr j
))))
168 (t (mapcar #'(lambda (j)
169 (resm1 (div* zn factors
) (car j
)))
172 ;; Compute the residues of zn/d for the simple poles in the list PL1.
173 ;; The poles must be simple because ZD must be the derivative of
174 ;; denominator. For simple poles the residue can be computed as
175 ;; limit(n(z)/d(z)*(z-a),z,a). Since the pole is simple we have the
176 ;; indeterminate form (z-a)/d(z). Use L'Hospital's rule to get
177 ;; 1/d'(z). Hence, the residue is easily computed as n(a)/d'(a).
178 (defun res1 (zn zd pl1
)
179 (setq zd
(div* zn zd
))
180 (mapcar #'(lambda (j)
181 ;; In case the pole is messy, call $RECTFORM. This
182 ;; works around some issues with gcd bugs in certain
183 ;; cases. (See bug 1073338.)
184 ($rectform
($expand
(subin ($rectform j
) zd
))))
187 (defun resprog0 (f g n n2
)
189 (setq a
(resprog f g
))
190 (setq b
(cadr a
) c
(ptimes (cddr a
) n2
) a
(caar a
))
191 (setq a
(ptimes n a
) b
(ptimes n b
))
192 (setq r
(pdivide a g
))
193 (setq a
(cadr r
) r
(car r
))
194 (setq b
(cons (pplus (ptimes (car r
) f
) (ptimes (cdr r
) b
))
196 (return (cons (cons (car a
) (ptimes (cdr a
) c
))
197 (cons (car b
) (ptimes (cdr b
) c
))))))
200 (defun resm0 (e n pole m
)
202 (setq e
($diff e var
(1- m
)))
203 (setq e
($rectform
($expand
(subin pole e
))))
204 (div* e
(simplify `((mfactorial) ,(1- m
)))))
206 (defun remfactor (l p n
)
209 (return (list (m*l
(cons leadcoef g
)) n
)))
210 ((equal p
(car l
)) (setq f
(cadr l
)))
211 (t (setq g
(cons (cadr l
) g
))))
215 (defun resprog (p1b p2b
)
216 (prog (temp coef1r coef2r fac coef1s coef2s zeropolb f1 f2
)
217 (setq coef2r
(setq coef1s
0))
218 (setq coef2s
(setq coef1r
1))
219 b1
(cond ((not (< (pdegree p1b var
) (pdegree p2b var
))) (go b2
)))
229 b2
(cond ((zerop (pdegree p2b var
))
230 (return (cons (cons coef2r p2b
) (cons coef2s p2b
)))))
231 (setq zeropolb
(psimp var
232 (list (- (pdegree p1b var
) (pdegree p2b var
))
234 (setq fac
(pgcd (caddr p1b
) (caddr p2b
)))
235 (setq f1
(pquotient (caddr p1b
) fac
))
236 (setq f2
(pquotient (caddr p2b
) fac
))
237 (setq p1b
(pdifference (ptimes f2
(psimp (car p1b
) (cdddr p1b
)))
242 (setq coef1r
(pdifference (ptimes f2 coef1r
)
243 (ptimes (ptimes f1 coef2r
) zeropolb
)))
244 (setq coef1s
(pdifference (ptimes f2 coef1s
)
245 (ptimes (ptimes f1 coef2s
) zeropolb
)))
248 ;;;Looks for polynomials. puts polys^(pos-num) in sn* polys^(neg-num) in sd*.
252 (setq sn
* (cons e sn
*)))
254 (integerp (caddr e
)))
255 (cond ((polyinx (cadr e
) var nil
)
256 (cond ((minusp (caddr e
))
257 (setq sd
* (cons (cond ((equal (caddr e
) -
1) (cadr e
))
261 (t (setq sn
* (cons e sn
*)))))))
263 (setq sn
* (cons e sn
*)))))
265 (setq sn
* nil sd
* nil
)
267 (defmfun $residue
(e var p
)
270 (list '(%residue
) e var p
))
273 (if (and (mtimesp e
) (andmapcar #'snumden
(cdr e
)))
274 (setq nn
* (m*l sn
*) dn
* (m*l sd
*))
276 (resm1 (div* nn
* dn
*) p
))))
278 (defun resm1 (e pole
)
279 ;; Call taylor with silent-taylor-flag t and catch an error.
280 (if (setq e
(catch 'taylor-catch
281 (let ((silent-taylor-flag t
))
282 (declare (special silent-taylor-flag
))
283 ;; Things like residue(s/(s^2-a^2),s,a) fails if use -1.
284 ($taylor e var pole
1))))
285 (coeff (ratdisrep e
) (m^
(m+ (m* -
1 pole
) var
) -
1) 1)
286 (merror (intl:gettext
"residue: taylor failed."))))