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 result) 14 15(declare-top (special varlist genvar $ratfac $keepfloat modulus *alpha xv)) 16 17(load-macsyma-macros ratmac) 18 19(defmfun $poly_discriminant (poly var) 20 (let* ((varlist (list var)) 21 ($ratfac nil) 22 (genvar ()) 23 (rform (rform poly)) 24 (rvar (car (last genvar))) 25 (n (pdegree (setq poly (car rform)) rvar))) 26 27 (cond ((= n 1) 1) 28 ((or (= n 0) (not (atom (cdr rform)))) 29 (merror (intl:gettext "poly_discriminant: first argument must be a polynomial in ~:M; found: ~M") var poly)) 30 (t (pdis (presign 31 (ash (* n (1- n)) -1) 32 (pquotient (resultant poly (pderivative poly rvar)) 33 (p-lc poly)))))))) 34 35(defmfun $resultant (a b mainvar) 36 (let ((varlist (list mainvar)) (genvar nil) 37 ($ratfac t) ($keepfloat nil) 38 formflag res (ans 1)) 39 (when ($ratp a) (setf a ($ratdisrep a) formflag t)) 40 (when ($ratp b) (setf b ($ratdisrep b) formflag t)) 41 (newvar a) 42 (newvar b) 43 (setq a (lmake2 (cadr (ratrep* a)) nil)) 44 (setq b (lmake2 (cadr (ratrep* b)) nil)) 45 (setq mainvar (caadr (ratrep* mainvar))) 46 47 (dolist (a-term a) 48 (dolist (b-term b) 49 (setq res (result1 (car a-term) (car b-term) mainvar)) 50 (setq ans (ptimes ans 51 (pexpt 52 (if (zerop (third res)) 53 (first res) 54 (ptimeschk (first res) 55 (pexpt (makprod (second res) nil) 56 (third res)))) 57 (* (cdr a-term) (cdr b-term))))))) 58 59 (if formflag (pdis* ans) (pdis ans)))) 60 61(defun result1 (p1 p2 var) 62 (cond ((or (pcoefp p1) (pointergp var (car p1))) 63 (list 1 p1 (pdegree p2 var))) 64 ((or (pcoefp p2) (pointergp var (car p2))) 65 (list 1 p2 (pdegree p1 var))) 66 ((null (cdddr p1)) 67 (cond ((null (cdddr p2)) (list 0 0 1)) 68 (t (list (pexpt (caddr p1) (cadr p2)) 69 (pcsubsty 0 var p2) 70 (cadr p1))))) 71 ((null (cdddr p2)) 72 (list (pexpt (caddr p2) (cadr p1)) 73 (pcsubsty 0 var p1) 74 (cadr p2))) 75 ((> (setq var (gcd (pgcdexpon p1) (pgcdexpon p2))) 1) 76 (list 1 (resultant (pexpon*// p1 var nil) 77 (pexpon*// p2 var nil)) var)) 78 (t (list 1 (resultant p1 p2) 1)))) 79 80(defmvar $resultant '$subres "Designates which resultant algorithm") 81 82(defvar *resultlist '($subres $mod $red)) 83 84(defun resultant (p1 p2) ;assumes same main var 85 (if (> (p-le p2) (p-le p1)) 86 (presign (* (p-le p1) (p-le p2)) (resultant p2 p1)) 87 (case $resultant 88 ($subres (subresult p1 p2)) 89 #+broken ($mod (modresult p1 p2)) 90 ($red (redresult p1 p2)) 91 (t (merror (intl:gettext "resultant: no such algorithm: ~M") $resultant))))) 92 93(defun presign (n p) 94 (if (oddp n) (pminus p) p)) 95 96;; Computes resultant using subresultant PRS. See Brown, "The Subresultant PRS 97;; Algorithm" (TOMS Sept. 1978). This is Algorithm 1, as found on page 241. 98;; 99;; This routine will not work if the coefficients contain hidden copies of the 100;; main polynomial variable. For example, if P is x*sqrt(x^2-1) + 2 then, when 101;; encoded as a polynomial in x, it appears as x*c + 2 for some (opaque) 102;; c. While doing the PRS calculation, the SUBRESULT code will square c, causing 103;; extra x's to appear and making the degree of polynomials derived from P 104;; behave erratically. The code might error or, worse, return a wrong result. 105;; 106;; Write G[1] = P, G[2] = Q and write g[i] for the leading coefficient of 107;; G[i]. On the k'th iteration of the loop, we are computing G[k+2], which is 108;; given by the formula 109;; 110;; G[k+2] = (-1)^(delta[k]+1) prem(G[k], G[k+1]) / (g[k] h[k]^delta[k]) 111;; 112;; except we set the denominator to 1 when k=1. Here, h[2] = g[2]^delta[1] and, 113;; for i >= 3, satisfies the recurrence: 114;; 115;; h[i] = g[i]^delta[i-1] h[i-1]^(1 - delta[i-1]) 116;; 117;; Here, delta[i] = deg(G[i]) - deg(G[i+1]), which is non-negative. 118;; 119;; Dictionary between program variables and values computed by the algorithm: 120;; 121;; - g is set to g[i-2] 122;; - h is set to g[i-2]^d / "h^1-d" 123;; 124;; Since d and h^1-d haven't yet been set on this iteration, they get their 125;; previous values, which turn out to give: 126;; 127;; - h is set to g[i-2]^delta[i-3] h[i-3]^[1-delta[i-3]] which, substituting 128;; above, is h[i-2]. 129;; 130;; Continuing: 131;; 132;; - degq is set to deg(G[i-1]) 133;; - d is set to delta[i-2] 134;; - h^1-d is (confusingly!) set to h[i-2]^(delta[i-2] - 1). 135 136(defun subresult (p q) 137 (loop for g = 1 then (p-lc p) 138 for h = 1 then (pquotient (pexpt g d) h^1-d) 139 for degq = (pdegree q (p-var p)) 140 for d = (- (p-le p) degq) 141 for h^1-d = (if (equal h 1) 1 (pexpt h (1- d))) 142 if (zerop degq) return (if (pzerop q) q (pquotient (pexpt q d) h^1-d)) 143 do (psetq p q 144 q (presign (1+ d) (pquotient (prem p q) 145 (ptimes g (ptimes h h^1-d))))))) 146 147;; PACKAGE FOR CALCULATING MULTIVARIATE POLYNOMIAL RESULTANTS 148;; USING MODIFIED REDUCED P.R.S. 149 150(defun redresult (u v) 151 (prog (a r sigma c) 152 (setq a 1) 153 (setq sigma 0) 154 (setq c 1) 155 a (if (pzerop (setq r (prem u v))) (return (pzero))) 156 (setq c (ptimeschk c (pexpt (p-lc v) 157 (* (- (p-le u) (p-le v)) 158 (- (p-le v) (pdegree r (p-var u)) 159 1))))) 160 (setq sigma (+ sigma (* (p-le u) (p-le v)))) 161 (if (zerop (pdegree r (p-var u))) 162 (return 163 (presign sigma 164 (pquotient (pexpt (pquotientchk r a) (p-le v)) c)))) 165 (psetq u v 166 v (pquotientchk r a) 167 a (pexpt (p-lc v) (+ (p-le u) 1 (- (p-le v))))) 168 (go a))) 169 170 171;; PACKAGE FOR CALCULATING MULTIVARIATE POLYNOMIAL RESULTANTS 172;; USING MODULAR AND EVALUATION HOMOMORPHISMS. 173;; modresultant fails on the following example 174;;RESULTANT(((-4)*Z)^4+(Y+8*Z)^4+(X-5*Z)^4-1, 175;; ((-4)*Z)^4-(X-5*Z)^3*((-4)*Z)^3+(Y+8*Z)^3*((-4)*Z)^2 176;; +(-2)*(Y+8*Z)^4+((-4)*Z)^4+1,Z) 177 178#+broken 179(progn 180 (defun modresult (a b) 181 (modresult1 a b (sort (union* (listovars a) (listovars b)) 182 (function pointergp)))) 183 184 (defun modresult1 (x y varl) 185 (cond ((null modulus) (pres x y (car varl) (cdr varl))) 186 (t (cpres x y (car varl) (cdr varl))) )) 187 188 (defun pres (a b xr1 varl) 189 (prog (m n f a* b* c* p q c modulus) 190 (setq m (cadr a)) 191 (setq n (cadr b)) 192 (setq f (coefbound m n (maxnorm (cdr a)) (maxnorm (cdr b)) )) 193 (setq q 1) 194 (setq c 0) 195 (setq p *alpha) 196 (go step3) 197 step2 (setq p (newprime p)) 198 step3 (set-modulus p) 199 (setq a* (pmod a)) 200 (setq b* (pmod b)) 201 (cond ((or (reject a* m xr1) (reject b* n xr1)) (go step2))) 202 (setq c* (cpres a* b* xr1 varl)) 203 (set-modulus nil) 204 (setq c (lagrange3 c c* p q)) 205 (setq q (* p q)) 206 (cond ((> q f) (return c)) 207 (t (go step2)) ) )) 208 209 (defun reject (a m xv) 210 (not (eqn (pdegree a xv) m))) 211 212 (defun coefbound (m n d e) 213 (* 2 (expt (1+ m) (ash n -1)) 214 (expt (1+ n) (ash m -1)) 215 (cond ((oddp n) (1+ ($isqrt (1+ m)))) 216 (t 1)) 217 (cond ((oddp m) (1+ ($isqrt (1+ n)))) 218 (t 1)) 219 ;; (FACTORIAL (PLUS M N)) USED TO REPLACE PREV. 4 LINES. KNU II P. 375 220 (expt d n) 221 (expt e m) )) 222 223 (defun main2 (a var exp tot) 224 (cond ((null a) (cons exp tot)) 225 (t (main2 (cddr a) var 226 (max (setq var (pdegree (cadr a) var)) exp) 227 (max (+ (car a) var) tot))) )) 228 229 (defun cpres (a b xr1 varl) ;XR1 IS MAIN VAR WHICH 230 (cond ((null varl) (cpres1 (cdr a) (cdr b))) ;RESULTANT ELIMINATES 231 (t (prog ( m2 ( m1 (cadr a)) 232 ( n1 (cadr b)) n2 (k 0) c d a* b* c* bp xv) ;XV IS INTERPOLATED VAR 233 (declare (fixnum m1 n1 k)) 234 235 step2 236 (setq xv (car varl)) 237 (setq varl (cdr varl)) 238 (setq m2 (main2 (cdr a) xv 0 0)) ;<XV DEG . TOTAL DEG> 239 (setq n2 (main2 (cdr b) xv 0 0)) 240 (cond ((zerop (+ (car m2) (car n2))) 241 (cond ((null varl) (return (cpres1 (cdr a) (cdr b)))) 242 (t (go step2)) ) )) 243 (setq k (1+ (min (+ (* m1 (car n2)) (* n1 (car m2))) 244 (+ (* m1 (cdr n2)) (* n1 (cdr m2)) 245 (- (* m1 n1))) ))) 246 (setq c 0) 247 (setq d 1) 248 (setq m2 (car m2) n2 (car n2)) 249 (setq bp (- 1)) 250 step3 251 (cond ((equal (setq bp (1+ bp)) modulus) 252 (merror "CPRES: resultant primes too small.")) 253 ((zerop m2) (setq a* a)) 254 (t (setq a* (pcsubst a bp xv)) 255 (cond ((reject a* m1 xr1)(go step3)) )) ) 256 (cond ((zerop n2) (setq b* b)) 257 (t (setq b* (pcsubst b bp xv)) 258 (cond ((reject b* n1 xr1) (go step3))) )) 259 (setq c* (cpres a* b* xr1 varl)) 260 (setq c (lagrange33 c c* d bp)) 261 (setq d (ptimeschk d (list xv 1 1 0 (cminus bp)))) 262 (cond ((> (cadr d) k) (return c)) 263 (t (go step3)))))))) 264 265 266;; *** NOTE THAT MATRIX PRODUCED IS ALWAYS SYMETRIC 267;; *** ABOUT THE MINOR DIAGONAL. 268 269(defmfun $bezout (p q var) 270 (let ((varlist (list var)) genvar) 271 (newvar p) 272 (newvar q) 273 (setq p (cadr (ratrep* p)) 274 q (cadr (ratrep* q))) 275 (setq p (cond ((> (cadr q) (cadr p)) (bezout q p)) 276 (t (bezout p q)))) 277 (cons '($matrix) 278 (mapcar #'(lambda (l) (cons '(mlist) (mapcar 'pdis l))) 279 p)))) 280 281(defun vmake (poly n *l) 282 (do ((i (1- n) (1- i))) ((minusp i)) 283 (cond ((or (null poly) (< (car poly) i)) 284 (setq *l (cons 0 *l))) 285 (t (setq *l (cons (cadr poly) *l)) 286 (setq poly (cddr poly))))) 287 (nreverse *l)) 288 289(defun bezout (p q) 290 (let* ((n (1+ (p-le p))) 291 (n2 (- n (p-le q))) 292 (a (vmake (p-terms p) n nil)) 293 (b (vmake (p-terms q) n nil)) 294 (ar (reverse (nthcdr n2 a))) 295 (br (reverse (nthcdr n2 b))) 296 (l (make-list n :initial-element 0))) 297 (rplacd (nthcdr (1- (p-le p)) a) nil) 298 (rplacd (nthcdr (1- (p-le p)) b) nil) 299 (nconc 300 (mapcar 301 #'(lambda (ar br) 302 (setq l (mapcar #'(lambda (a b l) 303 (ppluschk l (pdifference 304 (ptimes br a) (ptimes ar b)))) 305 a b (cons 0 l)))) 306 ar br) 307 (and (pzerop (car b)) 308 (do ((b (vmake (cdr q) (cadr p) nil) (rot* b)) 309 (m nil (cons b m))) 310 ((not (pzerop (car b))) (cons b m))))) )) 311 312(defun rot* (b) 313 (setq b (copy-list b)) 314 (prog2 315 (nconc b b) 316 (cdr b) 317 (rplacd b nil))) 318 319 320(defun ppluschk (p q) 321 (cond ((pzerop p) q) 322 (t (pplus p q)))) 323