1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module residu
)
15 (load-macsyma-macros rzmac
)
17 (declare-top (special $noprincipal
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
))
110 ;; Like POLELIST, but dependency on VAR is explicit. Use this instead
112 (defun polelist-var (var1 d region region1
)
113 (prog (roots $breakup r rr ss r1 s pole wflag cf
)
115 (setq *leadcoef
* (polyinx d var1
'leadcoef
))
116 (setq roots
(solvecase-var d var1
))
117 (if (eq roots
'failure
) (return ()))
118 ;; Loop over all the roots. SOLVECASE returns the roots in the form
125 (cond ((and *semirat
*
126 (> (+ (length s
) (length r
))
127 (+ (length ss
) (length rr
))))
128 ;; Return CF, repeated roots (*semirat*), simple
129 ;; roots (*semirat*), roots in region 1.
130 (return (list cf rr ss r1
)))
132 ;; Return CF, repeated roots, simple roots, roots in region 1.
133 (return (list cf r s r1
)))))
135 ;; Set POLE to the actual root and set D to the
136 ;; multiplicity of the root.
137 (setq pole
(caddar roots
))
138 (setq d
(cadr roots
))
140 ;; Is it possible for LEADCOEF to be NIL ever?
142 ;; Push (pole (x - pole)^d) onto the list CF.
145 (m^
(m+ var1
(m* -
1 pole
))
148 ;; Don't change the order of clauses here. We want to call REGION and then REGION1.
149 (cond ((funcall region pole
)
150 ;; The pole is in REGION
152 ;; A simple pole, so just push the pole onto the list S.
155 ;; A multiple pole, so push (pole d) onto the list R.
156 (push (list pole d
) r
))))
157 ((funcall region1 pole
)
158 ;; The pole is in REGION1
159 (cond ((not $noprincipal
)
160 ;; Put the pole onto the R1 list. (Don't know what
164 ;; Return NIL if we get here.
167 ;; (What does *SEMIRAT* mean?) Anyway if we're here, the
168 ;; pole is not in REGION or REGION1, so push the pole onto
169 ;; SS or RR depending if the pole is repeated or not.
172 (t (push (list pole d
) rr
)))))
173 ;; Pop this root and multiplicity and move on.
174 (setq roots
(cddr roots
))
178 (cond ((not (among var e
)) nil
)
179 (t (let (*failures
*roots
)
181 (cond (*failures
'failure
)
185 (defun solvecase-var (e var1
)
186 (cond ((not (among var1 e
)) nil
)
187 (t (let (*failures
*roots
)
189 (cond (*failures
'failure
)
193 ;; Compute the sum of the residues of n/d.
194 (defun res (n d region region1
)
195 (let ((pl (polelist d region region1
))
196 dp a b c factors
*leadcoef
*)
200 (setq factors
(car pl
))
202 ;; PL now contains the list of the roots in region, roots in
203 ;; region1, and everything else.
206 (setq dp
(sdiff d var
))))
208 ;; Compute the sum of the residues of n/d for the
209 ;; multiple roots in REGION.
210 (setq a
(m+l
(residue n
(cond (*leadcoef
* factors
)
215 ;; Compute the sum of the residues of n/d for the simple
216 ;; roots in REGION1. Since the roots are simple, we can
217 ;; use RES1 to compute the residues instead of the more
218 ;; complicated $RESIDUE. (This works around a bug
219 ;; described in bug 1073338.)
221 (setq b
(m+l
(mapcar #'(lambda (pole)
222 ($residue
(m// n d
) var pole
))
224 (setq b
(m+l
(res1 n dp
(cadr pl
)))))
227 ;; Compute the sum of the residues of n/d for the roots
228 ;; not in REGION nor REGION1.
229 (setq c
(m+l
(mapcar #'(lambda (pole)
230 ($residue
(m// n d
) var pole
))
233 ;; Return the sum of the residues in the two regions and the
234 ;; sum of the residues outside the two regions.
235 (list (m+ a b
) c
)))))
237 (defun res-var (var1 n d region region1
)
238 (let ((pl (polelist-var var1 d region region1
))
239 dp a b c factors
*leadcoef
*)
243 (setq factors
(car pl
))
245 ;; PL now contains the list of the roots in region, roots in
246 ;; region1, and everything else.
249 (setq dp
(sdiff d var1
))))
251 ;; Compute the sum of the residues of n/d for the
252 ;; multiple roots in REGION.
253 (setq a
(m+l
(let ((var var1
))
254 ;; This binding is very important for
255 ;; defint. It prevents the plog
256 ;; simplifier from simplifying plog until
257 ;; the pole has been substituted in when
258 ;; computing the residue.
259 (declare (special var
))
260 (residue-var var1 n
(cond (*leadcoef
* factors
)
265 ;; Compute the sum of the residues of n/d for the simple
266 ;; roots in REGION1. Since the roots are simple, we can
267 ;; use RES1 to compute the residues instead of the more
268 ;; complicated $RESIDUE. (This works around a bug
269 ;; described in bug 1073338.)
271 (setq b
(m+l
(mapcar #'(lambda (pole)
272 ($residue
(m// n d
) var1 pole
))
274 (setq b
(m+l
(res1-var var1 n dp
(cadr pl
)))))
277 ;; Compute the sum of the residues of n/d for the roots
278 ;; not in REGION nor REGION1.
279 (setq c
(m+l
(mapcar #'(lambda (pole)
280 ($residue
(m// n d
) var1 pole
))
283 ;; Return the sum of the residues in the two regions and the
284 ;; sum of the residues outside the two regions.
285 (list (m+ a b
) c
)))))
287 (defun residue (zn factors pl
)
289 (mapcar #'(lambda (j)
290 (destructuring-let (((factor1 factor2
) (remfactor factors
(car j
) zn
)))
291 (resm0 factor1 factor2
(car j
) (cadr j
))))
293 (t (mapcar #'(lambda (j)
294 (resm1 (div* zn factors
) (car j
)))
297 (defun residue-var (var1 zn factors pl
)
299 (mapcar #'(lambda (j)
300 (destructuring-let (((factor1 factor2
)
301 (remfactor factors
(car j
) zn
)))
302 (resm0-var var1 factor1 factor2
(car j
) (cadr j
))))
304 (t (mapcar #'(lambda (j)
305 (resm1-var var1
(div* zn factors
) (car j
)))
308 ;; Compute the residues of zn/d for the simple poles in the list PL1.
309 ;; The poles must be simple because ZD must be the derivative of
310 ;; denominator. For simple poles the residue can be computed as
311 ;; limit(n(z)/d(z)*(z-a),z,a). Since the pole is simple we have the
312 ;; indeterminate form (z-a)/d(z). Use L'Hospital's rule to get
313 ;; 1/d'(z). Hence, the residue is easily computed as n(a)/d'(a).
314 (defun res1 (zn zd pl1
)
315 (setq zd
(div* zn zd
))
316 (mapcar #'(lambda (j)
317 ;; In case the pole is messy, call $RECTFORM. This
318 ;; works around some issues with gcd bugs in certain
319 ;; cases. (See bug 1073338.)
320 ($rectform
($expand
(subin ($rectform j
) zd
))))
323 (defun res1-var (var1 zn zd pl1
)
324 (setq zd
(div* zn zd
))
325 (mapcar #'(lambda (j)
326 ;; In case the pole is messy, call $RECTFORM. This
327 ;; works around some issues with gcd bugs in certain
328 ;; cases. (See bug 1073338.)
329 ($rectform
($expand
(subin-var ($rectform j
) zd var1
))))
332 (defun resprog0 (f g n n2
)
334 (setq a
(resprog f g
))
335 (setq b
(cadr a
) c
(ptimes (cddr a
) n2
) a
(caar a
))
336 (setq a
(ptimes n a
) b
(ptimes n b
))
337 (setq r
(pdivide a g
))
338 (setq a
(cadr r
) r
(car r
))
339 (setq b
(cons (pplus (ptimes (car r
) f
) (ptimes (cdr r
) b
))
341 (return (cons (cons (car a
) (ptimes (cdr a
) c
))
342 (cons (car b
) (ptimes (cdr b
) c
))))))
344 (defun resprog0-var (var1 f g n n2
)
346 (setq a
(resprog-var var1 f g
))
347 (setq b
(cadr a
) c
(ptimes (cddr a
) n2
) a
(caar a
))
348 (setq a
(ptimes n a
) b
(ptimes n b
))
349 (setq r
(pdivide a g
))
350 (setq a
(cadr r
) r
(car r
))
351 (setq b
(cons (pplus (ptimes (car r
) f
) (ptimes (cdr r
) b
))
353 (return (cons (cons (car a
) (ptimes (cdr a
) c
))
354 (cons (car b
) (ptimes (cdr b
) c
))))))
357 (defun resm0 (e n pole m
)
359 (setq e
($diff e var
(1- m
)))
360 (setq e
($rectform
($expand
(subin pole e
))))
361 (div* e
(simplify `((mfactorial) ,(1- m
)))))
363 (defun resm0-var (var1 e n pole m
)
365 (setq e
($diff e var1
(1- m
)))
366 (setq e
($rectform
($expand
(subin-var pole e var1
))))
367 (div* e
(simplify `((mfactorial) ,(1- m
)))))
369 (defun remfactor (l p n
)
372 (return (list (m*l
(cons *leadcoef
* g
)) n
)))
373 ((equal p
(car l
)) (setq f
(cadr l
)))
374 (t (setq g
(cons (cadr l
) g
))))
378 (defun resprog (p1b p2b
)
379 (prog (temp coef1r coef2r fac coef1s coef2s zeropolb f1 f2
)
380 (setq coef2r
(setq coef1s
0))
381 (setq coef2s
(setq coef1r
1))
382 b1
(cond ((not (< (pdegree p1b var
) (pdegree p2b var
))) (go b2
)))
392 b2
(cond ((zerop (pdegree p2b var
))
393 (return (cons (cons coef2r p2b
) (cons coef2s p2b
)))))
394 (setq zeropolb
(psimp var
395 (list (- (pdegree p1b var
) (pdegree p2b var
))
397 (setq fac
(pgcd (caddr p1b
) (caddr p2b
)))
398 (setq f1
(pquotient (caddr p1b
) fac
))
399 (setq f2
(pquotient (caddr p2b
) fac
))
400 (setq p1b
(pdifference (ptimes f2
(psimp (car p1b
) (cdddr p1b
)))
405 (setq coef1r
(pdifference (ptimes f2 coef1r
)
406 (ptimes (ptimes f1 coef2r
) zeropolb
)))
407 (setq coef1s
(pdifference (ptimes f2 coef1s
)
408 (ptimes (ptimes f1 coef2s
) zeropolb
)))
411 (defun resprog-var (var1 p1b p2b
)
412 (prog (temp coef1r coef2r fac coef1s coef2s zeropolb f1 f2
)
413 (setq coef2r
(setq coef1s
0))
414 (setq coef2s
(setq coef1r
1))
415 b1
(cond ((not (< (pdegree p1b var1
) (pdegree p2b var1
))) (go b2
)))
425 b2
(cond ((zerop (pdegree p2b var1
))
426 (return (cons (cons coef2r p2b
) (cons coef2s p2b
)))))
427 (setq zeropolb
(psimp var1
428 (list (- (pdegree p1b var1
) (pdegree p2b var1
))
430 (setq fac
(pgcd (caddr p1b
) (caddr p2b
)))
431 (setq f1
(pquotient (caddr p1b
) fac
))
432 (setq f2
(pquotient (caddr p2b
) fac
))
433 (setq p1b
(pdifference (ptimes f2
(psimp (car p1b
) (cdddr p1b
)))
438 (setq coef1r
(pdifference (ptimes f2 coef1r
)
439 (ptimes (ptimes f1 coef2r
) zeropolb
)))
440 (setq coef1s
(pdifference (ptimes f2 coef1s
)
441 (ptimes (ptimes f1 coef2s
) zeropolb
)))
444 ;;;Looks for polynomials. puts polys^(pos-num) in sn* polys^(neg-num) in sd*.
448 (setq sn
* (cons e sn
*)))
450 (integerp (caddr e
)))
451 (cond ((polyinx (cadr e
) var nil
)
452 (cond ((minusp (caddr e
))
453 (setq sd
* (cons (cond ((equal (caddr e
) -
1) (cadr e
))
457 (t (setq sn
* (cons e sn
*)))))))
459 (setq sn
* (cons e sn
*)))))
461 ;; Like SNUMDEN, but dependency on VAR is explicit. Use this instead
463 (defun snumden-var (e var1 sn sd
)
466 (setq sn
(cons e sn
))
469 (integerp (caddr e
)))
470 (cond ((polyinx (cadr e
) var1 nil
)
471 (cond ((minusp (caddr e
))
472 (setq sd
(cons (cond ((equal (caddr e
) -
1) (cadr e
))
478 (setq sn
(cons e sn
))
479 (values t sn sd
))))))
480 ((polyinx e var1 nil
)
481 (setq sn
(cons e sn
))
484 (setq sn
* nil sd
* nil
)
486 (defmfun $residue
(e var p
)
489 (list '(%residue
) e var p
))
492 (if (and (mtimesp e
) (andmapcar #'snumden
(cdr e
)))
493 (setq nn
* (m*l sn
*) dn
* (m*l sd
*))
495 (resm1 (div* nn
* dn
*) p
))))
497 (defun resm1 (e pole
)
498 ;; Call taylor with silent-taylor-flag t and catch an error.
499 (if (setq e
(catch 'taylor-catch
500 (let ((silent-taylor-flag t
))
501 ;; Things like residue(s/(s^2-a^2),s,a) fails if use -1.
502 ($taylor e var pole
1))))
503 (coeff (ratdisrep e
) (m^
(m+ (m* -
1 pole
) var
) -
1) 1)
504 (merror (intl:gettext
"residue: taylor failed."))))
506 (defun resm1 (e pole
)
507 ;; Call taylor with silent-taylor-flag t and catch an error.
508 (if (setq e
(catch 'taylor-catch
509 (let ((silent-taylor-flag t
))
510 ;; Things like residue(s/(s^2-a^2),s,a) fails if use -1.
511 ($taylor e var pole
1))))
512 (coeff (ratdisrep e
) (m^
(m+ (m* -
1 pole
) var
) -
1) 1)
513 (merror (intl:gettext
"residue: taylor failed."))))
515 (defun resm1-var (var1 e pole
)
516 ;; Call taylor with silent-taylor-flag t and catch an error.
517 (if (setq e
(catch 'taylor-catch
518 (let ((silent-taylor-flag t
))
519 ;; Things like residue(s/(s^2-a^2),s,a) fails if use -1.
520 ($taylor e var1 pole
1))))
521 (coeff (ratdisrep e
) (m^
(m+ (m* -
1 pole
) var1
) -
1) 1)
522 (merror (intl:gettext
"residue: taylor failed."))))