Rename tab to datatab in numeric/interpol.mac
[maxima.git] / src / residu.lisp
blobcfe98836eb60cabc0885f6f7aa2c8a3b11257eff
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
4 ;;; ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 (in-package :maxima)
13 (macsyma-module residu)
15 (load-macsyma-macros rzmac)
17 (declare-top (special $noprincipal
18 var *roots *failures
19 sn* sd* genvar 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)
47 (setq wflag t)
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
52 ;; ((x = r1) mult1
53 ;; (x = r2) mult2
54 ;; ...)
56 loop1
57 (cond ((null roots)
58 (cond ((and *semirat*
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))
71 (setq d (cadr roots))
72 (cond (*leadcoef*
73 ;; Is it possible for LEADCOEF to be NIL ever?
75 ;; Push (pole (x - pole)^d) onto the list CF.
76 (setq cf (cons pole
77 (cons
78 (m^ (m+ var (m* -1 pole))
80 cf)))))))
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
84 (cond ((equal d 1)
85 ;; A simple pole, so just push the pole onto the list S.
86 (push pole 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
94 ;; $NOPRINCIPAL is.)
95 (push pole r1))
97 ;; Return NIL if we get here.
98 (return nil))))
99 (*semirat*
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.
103 (cond ((equal d 1)
104 (push pole ss))
105 (t (push (list pole d) rr)))))
106 ;; Pop this root and multiplicity and move on.
107 (setq roots (cddr roots))
108 (go loop1)))
110 ;; Like POLELIST, but dependency on VAR is explicit. Use this instead
111 ;; when possible.
112 (defun polelist-var (var1 d region region1)
113 (prog (roots $breakup r rr ss r1 s pole wflag cf)
114 (setq wflag t)
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
119 ;; ((x = r1) mult1
120 ;; (x = r2) mult2
121 ;; ...)
123 loop1
124 (cond ((null roots)
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))
139 (cond (*leadcoef*
140 ;; Is it possible for LEADCOEF to be NIL ever?
142 ;; Push (pole (x - pole)^d) onto the list CF.
143 (setq cf (cons pole
144 (cons
145 (m^ (m+ var1 (m* -1 pole))
147 cf)))))))
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
151 (cond ((equal d 1)
152 ;; A simple pole, so just push the pole onto the list S.
153 (push pole 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
161 ;; $NOPRINCIPAL is.)
162 (push pole r1))
164 ;; Return NIL if we get here.
165 (return nil))))
166 (*semirat*
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.
170 (cond ((equal d 1)
171 (push pole ss))
172 (t (push (list pole d) rr)))))
173 ;; Pop this root and multiplicity and move on.
174 (setq roots (cddr roots))
175 (go loop1)))
177 (defun solvecase (e)
178 (cond ((not (among var e)) nil)
179 (t (let (*failures *roots)
180 (solve e var 1)
181 (cond (*failures 'failure)
182 ((null *roots) ())
183 (t *roots))))))
185 (defun solvecase-var (e var1)
186 (cond ((not (among var1 e)) nil)
187 (t (let (*failures *roots)
188 (solve e var1 1)
189 (cond (*failures 'failure)
190 ((null *roots) ())
191 (t *roots))))))
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*)
197 (cond
198 ((null pl) nil)
200 (setq factors (car pl))
201 (setq pl (cdr pl))
202 ;; PL now contains the list of the roots in region, roots in
203 ;; region1, and everything else.
204 (cond ((or (cadr pl)
205 (caddr pl))
206 (setq dp (sdiff d var))))
207 (cond ((car pl)
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)
211 (t d))
212 (car pl)))))
213 (t (setq a 0)))
214 (cond ((cadr pl)
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.)
220 #+nil
221 (setq b (m+l (mapcar #'(lambda (pole)
222 ($residue (m// n d) var pole))
223 (cadr pl))))
224 (setq b (m+l (res1 n dp (cadr pl)))))
225 (t (setq b 0)))
226 (cond ((caddr 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))
231 (caddr pl)))))
232 (t (setq c ())))
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*)
240 (cond
241 ((null pl) nil)
243 (setq factors (car pl))
244 (setq pl (cdr pl))
245 ;; PL now contains the list of the roots in region, roots in
246 ;; region1, and everything else.
247 (cond ((or (cadr pl)
248 (caddr pl))
249 (setq dp (sdiff d var1))))
250 (cond ((car pl)
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)
261 (t d))
262 (car pl))))))
263 (t (setq a 0)))
264 (cond ((cadr pl)
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.)
270 #+nil
271 (setq b (m+l (mapcar #'(lambda (pole)
272 ($residue (m// n d) var1 pole))
273 (cadr pl))))
274 (setq b (m+l (res1-var var1 n dp (cadr pl)))))
275 (t (setq b 0)))
276 (cond ((caddr 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))
281 (caddr pl)))))
282 (t (setq c ())))
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)
288 (cond (*leadcoef*
289 (mapcar #'(lambda (j)
290 (destructuring-let (((factor1 factor2) (remfactor factors (car j) zn)))
291 (resm0 factor1 factor2 (car j) (cadr j))))
292 pl))
293 (t (mapcar #'(lambda (j)
294 (resm1 (div* zn factors) (car j)))
295 pl))))
297 (defun residue-var (var1 zn factors pl)
298 (cond (*leadcoef*
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))))
303 pl))
304 (t (mapcar #'(lambda (j)
305 (resm1-var var1 (div* zn factors) (car j)))
306 pl))))
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))))
321 pl1))
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))))
330 pl1))
332 (defun resprog0 (f g n n2)
333 (prog (a b c r)
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))
340 (cdr r)))
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)
345 (prog (a b c r)
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))
352 (cdr r)))
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)
358 (setq e (div* n e))
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)
364 (setq e (div* n e))
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)
370 (prog (f g)
371 loop (cond ((null l)
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))))
375 (setq l (cddr l))
376 (go loop)))
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)))
383 (setq temp p2b)
384 (setq p2b p1b)
385 (setq p1b temp)
386 (setq temp coef2r)
387 (setq coef2r coef1r)
388 (setq coef1r temp)
389 (setq temp coef2s)
390 (setq coef2s coef1s)
391 (setq coef1s temp)
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))
396 1)))
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)))
401 (ptimes f1
402 (ptimes zeropolb
403 (psimp (car p2b)
404 (cdddr p2b))))))
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)))
409 (go b1)))
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)))
416 (setq temp p2b)
417 (setq p2b p1b)
418 (setq p1b temp)
419 (setq temp coef2r)
420 (setq coef2r coef1r)
421 (setq coef1r temp)
422 (setq temp coef2s)
423 (setq coef2s coef1s)
424 (setq coef1s temp)
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))
429 1)))
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)))
434 (ptimes f1
435 (ptimes zeropolb
436 (psimp (car p2b)
437 (cdddr p2b))))))
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)))
442 (go b1)))
444 ;;;Looks for polynomials. puts polys^(pos-num) in sn* polys^(neg-num) in sd*.
445 (defun snumden (e)
446 (cond ((or (atom e)
447 (mnump e))
448 (setq sn* (cons e sn*)))
449 ((and (mexptp e)
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))
454 (t (m^ (cadr e)
455 (- (caddr e)))))
456 sd*)))
457 (t (setq sn* (cons e sn*)))))))
458 ((polyinx e var nil)
459 (setq sn* (cons e sn*)))))
461 ;; Like SNUMDEN, but dependency on VAR is explicit. Use this instead
462 ;; when possible.
463 (defun snumden-var (e var1 sn sd)
464 (cond ((or (atom e)
465 (mnump e))
466 (setq sn (cons e sn))
467 (values t sn sd))
468 ((and (mexptp e)
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))
473 (t (m^ (cadr e)
474 (- (caddr e)))))
475 sd))
476 (values t sn sd))
478 (setq sn (cons e sn))
479 (values t sn sd))))))
480 ((polyinx e var1 nil)
481 (setq sn (cons e sn))
482 (values t sn sd))))
484 (setq sn* nil sd* nil)
486 (defmfun $residue (e var p)
487 (cond (($unknown e)
488 ($nounify '$residue)
489 (list '(%residue) e var p))
491 (let (sn* sd*)
492 (if (and (mtimesp e) (andmapcar #'snumden (cdr e)))
493 (setq nn* (m*l sn*) dn* (m*l sd*))
494 (numden e)))
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."))))