src/result.lisp

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