src/result.lisp

   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))))