1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 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;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module residu) 14 15(load-macsyma-macros rzmac) 16 17(declare-top (special $breakup $noprincipal varlist 18 leadcoef var *roots *failures wflag nn* 19 sn* sd* genvar dn* zn)) 20 21 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.) 26;; 27;; The regions are determined by functions REGION and REGION1, which 28;; should return non-NIL if the root is in the given region. 29;; 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. 33;; 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. 37;; 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. 41;; 42;; The third part is a list of the simple roots in REGION. 43;; 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 ;; ...) 55 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))) 64 (t 65 ;; Return CF, repeated roots, simple roots, roots in region 1. 66 (return (list cf r s r1))))) 67 (t 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? 74 ;; 75 ;; Push (pole (x - pole)^d) onto the list CF. 76 (setq cf (cons pole 77 (cons 78 (m^ (m+ var (m* -1 pole)) 79 d) 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)) 87 (t 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)) 96 (t 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))) 109 110(defun solvecase (e) 111 (cond ((not (among var e)) nil) 112 (t (let (*failures *roots) 113 (solve e var 1) 114 (cond (*failures 'failure) 115 ((null *roots) ()) 116 (t *roots)))))) 117 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) 122 (cond 123 ((null pl) nil) 124 (t 125 (setq factors (car pl)) 126 (setq pl (cdr pl)) 127 ;; PL now contains the list of the roots in region, roots in 128 ;; region1, and everything else. 129 (cond ((or (cadr pl) 130 (caddr pl)) 131 (setq dp (sdiff d var)))) 132 (cond ((car pl) 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) 136 (t d)) 137 (car pl))))) 138 (t (setq a 0))) 139 (cond ((cadr pl) 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.) 145 #+nil 146 (setq b (m+l (mapcar #'(lambda (pole) 147 ($residue (m// n d) var pole)) 148 (cadr pl)))) 149 (setq b (m+l (res1 n dp (cadr pl))))) 150 (t (setq b 0))) 151 (cond ((caddr 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)) 156 (caddr pl))))) 157 (t (setq c ()))) 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))))) 161 162(defun residue (zn factors pl) 163 (cond (leadcoef 164 (mapcar #'(lambda (j) 165 (destructuring-let (((factor1 factor2) (remfactor factors (car j) zn))) 166 (resm0 factor1 factor2 (car j) (cadr j)))) 167 pl)) 168 (t (mapcar #'(lambda (j) 169 (resm1 (div* zn factors) (car j))) 170 pl)))) 171 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)))) 185 pl1)) 186 187(defun resprog0 (f g n n2) 188 (prog (a b c r) 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)) 195 (cdr r))) 196 (return (cons (cons (car a) (ptimes (cdr a) c)) 197 (cons (car b) (ptimes (cdr b) c)))))) 198 199 200(defun resm0 (e n pole m) 201 (setq e (div* n e)) 202 (setq e ($diff e var (1- m))) 203 (setq e ($rectform ($expand (subin pole e)))) 204 (div* e (simplify `((mfactorial) ,(1- m))))) 205 206(defun remfactor (l p n) 207 (prog (f g) 208 loop (cond ((null l) 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)))) 212 (setq l (cddr l)) 213 (go loop))) 214 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))) 220 (setq temp p2b) 221 (setq p2b p1b) 222 (setq p1b temp) 223 (setq temp coef2r) 224 (setq coef2r coef1r) 225 (setq coef1r temp) 226 (setq temp coef2s) 227 (setq coef2s coef1s) 228 (setq coef1s temp) 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)) 233 1))) 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))) 238 (ptimes f1 239 (ptimes zeropolb 240 (psimp (car p2b) 241 (cdddr p2b)))))) 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))) 246 (go b1))) 247 248;;;Looks for polynomials. puts polys^(pos-num) in sn* polys^(neg-num) in sd*. 249(defun snumden (e) 250 (cond ((or (atom e) 251 (mnump e)) 252 (setq sn* (cons e sn*))) 253 ((and (mexptp e) 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)) 258 (t (m^ (cadr e) 259 (- (caddr e))))) 260 sd*))) 261 (t (setq sn* (cons e sn*))))))) 262 ((polyinx e var nil) 263 (setq sn* (cons e sn*))))) 264 265(setq sn* nil sd* nil) 266 267(defmfun $residue (e var p) 268 (cond (($unknown e) 269 ($nounify '$residue) 270 (list '(%residue) e var p)) 271 (t 272 (let (sn* sd*) 273 (if (and (mtimesp e) (andmapcar #'snumden (cdr e))) 274 (setq nn* (m*l sn*) dn* (m*l sd*)) 275 (numden e))) 276 (resm1 (div* nn* dn*) p)))) 277 278(defun resm1 (e pole) 279 ;; Call $ratsimp to simplify pole. Not sure if this is really 280 ;; necessary or desired, but it fixes bug 1504505. It would be 281 ;; better to fix $taylor, but that seems much harder. 282 (let ((pole ($ratsimp ($rectform pole)))) 283 ;; Call taylor with silent-taylor-flag t and catch an error. 284 (if (setq e (catch 'taylor-catch 285 (let ((silent-taylor-flag t)) 286 (declare (special silent-taylor-flag)) 287 ;; Things like residue(s/(s^2-a^2),s,a) fails if use -1. 288 ($taylor e var pole 1)))) 289 (coeff (ratdisrep e) (m^ (m+ (m* -1 pole) var) -1) 1) 290 (merror (intl:gettext "residue: taylor failed."))))) 291