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 1981 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module rat3c) 14 15;; THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 3. 16;; IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS 17 18(declare-top (special $keepfloat $algebraic $ratfac genvar)) 19 20;; List of GCD algorithms. Default one is first. 21(defmvar *gcdl* '($spmod $subres $ez $red $mod $algebraic)) 22 23(defmvar $gcd (car *gcdl*)) ;Sparse Modular 24 25(defun cgcd (a b) 26 (cond (modulus 1) 27 ((and $keepfloat (or (floatp a) (floatp b))) 1) 28 (t (gcd a b)))) 29 30(defun pquotientchk (a b) 31 (if (equal b 1) a (pquotient a b))) 32 33;; divides polynomial x by polynomial y 34;; avoids error "quotient by polynomial of higher degree" 35;; (returns nil in this case) 36(defun pquotientchk-safe (x y) 37 (ignore-rat-err (pquotientchk x y))) 38 39(defun ptimeschk (a b) 40 (cond ((equal a 1) b) 41 ((equal b 1) a) 42 (t (ptimes a b)))) 43 44(defun pfloatp (x) 45 (catch 'float (if (pcoefp x) (floatp x) (pfloatp1 x)))) 46 47(defun pfloatp1 (x) 48 (mapc #'(lambda (q) (cond ((pcoefp q) (when (floatp q) (throw 'float t))) 49 ((pfloatp1 q)))) 50 (cdr x)) 51 nil) 52 53(defun pgcd (x y) 54 (setq x (car (pgcda x y nil))) 55 (cond ((pminusp x) (pminus x)) 56 (modulus (monize x)) 57 (t x))) 58 59(defun plcm (x y) 60 (setq x (pgcdcofacts x y)) 61 (ptimes (car x) (ptimes (cadr x) (caddr x)))) 62 63(defun plcmcofacts (x y) 64 (setq x (pgcdcofacts x y)) 65 (list (ptimes (car x) (ptimes (cadr x) (caddr x))) 66 (caddr x) (cadr x))) 67 68; returns list (gcd xx yy alg) 69; where x * y = gcd^2 * xx * yy / alg^2 70; and alg is non-nil only when $algebraic is true 71(defun pgcdcofacts (x y) 72 (let ((a (pgcda x y t))) 73 (cond ((cdr a) a) 74 ((equal (setq a (car a)) 1) (list 1 x y)) 75 ((and $algebraic (not (pcoefp a))) 76 (cons a (prog2 (setq x (rquotient x a) 77 y (rquotient y a) 78 a (pgcdcofacts (cdr x) (cdr y))) 79 (list (ptimes (car x) (caddr a)) 80 (ptimes (car y) (cadr a)) 81 (ptimes (cadr a) (cdr y)))))) 82 ((eq a x) (list x 1 (pquotient y x))) 83 ((eq a y) (list a (pquotient x y) 1)) 84 (t (list a (pquotient x a) (pquotient y a)))))) 85 86(defun pgcda (x y cofac? &aux a c) 87 (cond ((not $gcd) (list 1 x y)) 88 ((and $keepfloat (or (pfloatp x) (pfloatp y))) 89 (cond ((or (pcoefp x) (pcoefp y) 90 (pcoefp (setq a (car (ptermcont x)))) 91 (pcoefp (setq a (pgcd a (car (ptermcont y)))))) 92 (list 1 x y)) 93 (t (list a)))) 94 ((pcoefp x) 95 (cond ((pcoefp y) 96 (cons (setq a (cgcd x y)) 97 (and cofac? 98 (list (cquotient x a) ;(CQUOTIENT 0 0) = 0 99 (cquotient y a))))) 100 ((zerop x) (list y x 1)) 101 (t (list (pcontent1 (cdr y) x))))) 102 ((pcoefp y) (cond ((zerop y) (list x 1 y)) 103 (t (list (pcontent1 (cdr x) y))))) 104 ((equal x y) (list x 1 1)) 105 ($ratfac (multiple-value-list (fpgcdco x y))) 106 ((not (eq (p-var x) (p-var y))) 107 (list (if (pointergp (p-var x) (p-var y)) 108 (oldcontent1 (p-terms x) y) 109 (oldcontent1 (p-terms y) x)))) 110 ((progn (desetq (a x) (ptermcont x)) 111 (desetq (c y) (ptermcont y)) 112 (not (and (equal a 1) (equal c 1)))) 113 (mapcar #'ptimes (monomgcdco a c cofac?) (pgcda x y cofac?))) 114 ((and (not $algebraic) (not modulus) 115 (desetq (a . c) (lin-var-find (nreverse (pdegreevector x)) 116 (nreverse (pdegreevector y)) 117 (reverse genvar)))) 118 (cond ((= a 1) (linhack x y (car c) (cadr c) cofac?)) 119 (t (setq a (linhack y x a (cadr c) cofac?)) 120 (if (cdr a) (rplacd a (nreverse (cdr a)))) 121 a))) 122 ((eq $gcd '$spmod) (list (zgcd x y))) 123 ((eq $gcd '$subres) (list (oldgcd x y))) 124 ((eq $gcd '$algebraic) 125 (if (or (palgp x) (palgp y)) 126 (let (($gcd '$subres)) (list (oldgcd x y))) 127 (let (($gcd '$spmod)) (list (zgcd x y))))) 128 ((eq $gcd '$ez) (ezgcd2 x y)) 129 ((eq $gcd '$red) (list (oldgcd x y))) 130 ((eq $gcd '$mod) (newgcd x y modulus)) 131 ((not (member $gcd *gcdl* :test #'eq)) 132 (merror (intl:gettext "gcd: 'gcd' variable must be one of ~M; found: ~M") *gcdl* $gcd)) 133 (t (list 1 x y)))) 134 135(defun monomgcdco (p q cofac?) 136 (let ((gcd (monomgcd p q))) 137 (cons gcd (if cofac? (list (pquotient p gcd) (pquotient q gcd)) ())))) 138 139(defun monomgcd (p q) 140 (cond ((or (pcoefp p) (pcoefp q)) 1) 141 ((eq (p-var p) (p-var q)) 142 (make-poly (p-var p) (min (p-le p) (p-le q)) 143 (monomgcd (p-lc p) (p-lc q)))) 144 ((pointergp (car p) (car q)) (monomgcd (p-lc p) q)) 145 (t (monomgcd p (p-lc q))))) 146 147(defun linhack (pol1 pol2 nonlindeg var cofac?) 148 (prog (coeff11 coeff12 gcdab rpol1 rpol2 gcdcd gcdcoef) 149 (desetq (coeff11 . coeff12) (bothprodcoef (make-poly var) pol1)) 150 (setq gcdab (if (pzerop coeff12) coeff11 151 (pgcd coeff11 coeff12))) 152 (cond ((equal gcdab 1) 153 (cond ((setq coeff11 (testdivide pol2 pol1)) 154 (return (list pol1 1 coeff11))) 155 (t (return (list 1 pol1 pol2)))))) 156 (setq rpol1 (pquotient pol1 gcdab)) 157 (desetq (gcdcd rpol2) (linhackcontent var pol2 nonlindeg)) 158 (cond ((equal gcdcd 1) 159 (cond ((setq coeff12 (testdivide rpol2 rpol1)) 160 (return (list rpol1 gcdab coeff12))) 161 (t (return (list 1 pol1 pol2)))))) 162 (cond (cofac? (desetq (gcdcoef coeff11 coeff12) 163 (pgcdcofacts gcdab gcdcd)) 164 (cond ((setq gcdcd (testdivide rpol2 rpol1)) 165 (return (list (ptimes gcdcoef rpol1) 166 coeff11 167 (ptimes coeff12 gcdcd)))) 168 (t (return (list gcdcoef 169 (ptimes coeff11 rpol1) 170 (ptimes coeff12 rpol2)))))) 171 (t (setq gcdcoef (pgcd gcdcd gcdab)) 172 (cond ((testdivide rpol2 rpol1) 173 (return (list (ptimes gcdcoef rpol1)))) 174 (t (return (list gcdcoef)))))))) 175 176(defun lin-var-find (a b c) 177 (do ((varl c (cdr varl)) 178 (degl1 a (cdr degl1)) 179 (degl2 b (cdr degl2))) 180 ((or (null degl1) (null degl2)) nil) 181 (if (equal (min (car degl1) (car degl2)) 1) 182 (return (list (car degl1) (car degl2) (car varl)))))) 183 184(defun linhackcontent (var pol nonlindeg &aux (npol pol) coef gcd) 185 (do ((i nonlindeg (1- i))) 186 ((= i 0) (list (setq gcd (pgcd gcd npol)) (pquotient pol gcd))) 187 (desetq (coef . npol) (bothprodcoef (make-poly var i 1) npol)) 188 (unless (pzerop coef) 189 (setq gcd (if (null gcd) coef (pgcd coef gcd))) 190 (if (equal gcd 1) (return (list 1 pol)))))) 191 192;;*** THIS IS THE REDUCED POLYNOMIAL REMAINDER SEQUENCE GCD (COLLINS') 193 194(defun oldgcd (x y &aux u v s egcd) ;only called from pgcda 195 (desetq (x u) (oldcontent x)) 196 (desetq (y v) (oldcontent y)) 197 (setq egcd (gcd (pgcdexpon u) (pgcdexpon v))) 198 (if (> egcd 1) 199 (setq u (pexpon*// u egcd nil) 200 v (pexpon*// v egcd nil))) 201 (if (> (p-le v) (p-le u)) (rotatef u v)) 202 (setq s (case $gcd 203 ($red (redgcd u v)) 204 ($subres (subresgcd u v)) 205 (t (merror "OLDGCD: found gcd = ~M; how did that happen?" $gcd)))) 206 ;; Check for gcd that simplifies to 0. SourceForge bugs 831445 and 1313987 207 (unless (ignore-rat-err (rainv s)) 208 (setq s 1)) 209 (unless (equal s 1) 210 (setq s (pexpon*// (primpart 211 (if $algebraic s 212 (pquotient s (pquotient (p-lc s) 213 (pgcd (p-lc u) (p-lc v)))))) 214 egcd t))) 215 (setq s (ptimeschk s (pgcd x y))) 216 (and $algebraic (not (pcoefp (setq u (leadalgcoef s)))) 217 (not (equal u s)) (setq s (algnormal s))) 218 (cond (modulus (monize s)) 219 ((pminusp s) (pminus s)) 220 (t s))) 221 222(defun pgcdexpon (p) 223 (if (pcoefp p) 0 224 (do ((d (cadr p) (gcd d (car l))) 225 (l (cdddr p) (cddr l))) 226 ((or (null l) (= d 1)) d)))) 227 228(defun pexpon*// (p n *?) 229 (if (or (pcoefp p) (= n 1)) p 230 (do ((ans (list (car p)) 231 (cons (cadr l) 232 (cons (if *? (* (car l) n) 233 (truncate (car l) n)) 234 ans))) 235 (l (cdr p) (cddr l))) 236 ((null l) (nreverse ans))))) 237 238;;polynomial gcd using reduced prs 239 240(defun redgcd (p q &aux (d 0)) 241 (loop until (zerop (pdegree q (p-var p))) 242 do (psetq p q 243 q (pquotientchk-safe (prem p q) (pexpt (p-lc p) d)) 244 d (+ (p-le p) 1 (- (p-le q)))) 245 (if (< d 1) (return 1)) 246 finally (return (if (pzerop q) p 1)))) 247 248;;computes gcd's using subresultant prs 249;;ACM Transactions On Mathematical Software Sept. 1978 250 251(defun subresgcd (p q) 252 (loop for g = 1 then (p-lc p) 253 for h = 1 then (pquotientchk-safe (pexpt g d) h^1-d) 254 for d = (- (p-le p) (p-le q)) 255 for h^1-d = 1 then (if (< d 1) 256 (return 1) 257 (pexpt h (1- d))) 258 do (psetq p q 259 q (pquotientchk-safe (prem p q) (ptimes g (ptimes h h^1-d)))) 260 if (zerop (pdegree q (p-var p))) return (if (pzerop q) p 1))) 261 262;;*** THIS COMPUTES PSEUDO REMAINDERS 263 264(defun psquorem1 (u v quop) 265 (prog (k (m 0) lcu lcv quo lc) 266 (declare (special lcu lcv)) 267 (setq lcv (pt-lc v)) 268 (setq k (- (pt-le u) (pt-le v))) 269 (cond ((minusp k) (return (list 1 '(0 0) u)))) 270 (if quop (setq lc (pexpt (pt-lc v) (1+ k)))) 271 a (setq lcu (pminus (pt-lc u))) 272 (if quop (setq quo (cons (ptimes (pt-lc u) (pexpt (pt-lc v) k)) 273 (cons k quo)))) 274 (cond ((null (setq u (pgcd2 (pt-red u) (pt-red v) k))) 275 (return (list lc (nreverse quo) '(0 0)))) 276 ((minusp (setq m (- (pt-le u) (pt-le v)))) 277 (setq u (cond ((zerop k) u) 278 (t (pctimes1 (pexpt lcv k) u)))) 279 (return (list lc (nreverse quo) u))) 280 ((> (1- k) m) 281 (setq u (pctimes1 (pexpt lcv (- (1- k) m)) u)))) 282 (setq k m) 283 (go a))) 284 285(defun prem (p q) 286 (cond ((pcoefp p) 287 (if (pcoefp q) 288 (if (or modulus (floatp p) (floatp q)) 289 0 290 (rem p q)) 291 p)) 292 ((pcoefp q) (pzero)) 293 (t (psimp (p-var p) (pgcd1 (p-terms p) (p-terms q)))))) 294 295(defun pgcd1 (u v) (caddr (psquorem1 u v nil))) 296 297(defun pgcd2 (u v k &aux (i 0)) 298 (declare (special lcu lcv) (fixnum k i)) 299 (cond ((null u) (pcetimes1 v k lcu)) 300 ((null v) (pctimes1 lcv u)) 301 ((zerop (setq i (+ (pt-le u) (- k) (- (car v))))) 302 (pcoefadd (pt-le u) (pplus (ptimes lcv (pt-lc u)) 303 (ptimes lcu (pt-lc v))) 304 (pgcd2 (pt-red u) (pt-red v) k))) 305 ((minusp i) 306 (list* (+ (pt-le v) k) (ptimes lcu (pt-lc v)) (pgcd2 u (pt-red v) k))) 307 (t (list* (pt-le u) (ptimes lcv (pt-lc u)) (pgcd2 (pt-red u) v k))))) 308 309;;;*** OLDCONTENT REMOVES ALL BUT MAIN VARIABLE AND PUTS THAT IN CONTENT 310;;;*** OLDCONTENT OF 3*A*X IS 3*A (WITH MAINVAR=X) 311 312(defun rcontent (p) ;RETURNS RAT-FORMS 313 (let ((q (oldcontenta p))) 314 (list (cons q 1) (cond ($algebraic (rquotient p q)) 315 (t (cons (pquotient p q) 1)))))) 316 317(defun oldcontenta (x) 318 (cond ((pcoefp x) x) 319 (t (setq x (contsort (cdr x))) 320 (oldcontent2 (cdr x) (car x))))) 321 322(defun oldcontent (x) 323 (cond ((pcoefp x) (list x 1)) 324 ((null (p-red x)) 325 (list (p-lc x) (make-poly (p-var x) (p-le x) 1))) 326 (t (let ((u (contsort (cdr x))) v) 327 (setq u (oldcontent2 (cdr u) (car u)) 328 v (cond ($algebraic (car (rquotient x u))) 329 (t (pcquotient x u)))) 330 (cond ((pminusp v) (list (pminus u) (pminus v))) 331 (t (list u v))))))) 332 333(defun oldcontent1 (x gcd) 334 (cond ((equal gcd 1) 1) 335 ((null x) gcd) 336 (t (oldcontent2 (contsort x) gcd)))) 337 338(defun oldcontent2 (x gcd) 339 (do ((x x (cdr x)) 340 (gcd gcd (pgcd (car x) gcd))) 341 ((or (null x) (equal gcd 1)) gcd))) 342 343(defun contsort (x) 344 (setq x (coefl x)) 345 (cond ((member 1 x) '(1)) 346 ((null (cdr x)) x) 347 (t (sort x #'contodr)))) 348 349(defun coefl (x) 350 (do ((x x (cddr x)) 351 (ans nil (cons (cadr x) ans))) 352 ((null x) ans))) 353 354(defun contodr (a b) 355 (cond ((pcoefp a) t) 356 ((pcoefp b) nil) 357 ((eq (car a) (car b)) (not (> (cadr a) (cadr b)))) 358 (t (pointergp (car b)(car a))))) 359 360;;;*** PCONTENT COMPUTES INTEGER CONTENT 361;;;*** PCONTENT OF 3*A*X IS 3 IF MODULUS = NIL 1 OTHERWISE 362 363(defun pcontent (x) 364 (cond ((pcoefp x) (list x 1)) 365 (t (let ((u (pcontentz x))) 366 (if (equal u 1) (list 1 x) 367 (list u (pcquotient x u))))))) 368 369(defun pcontent1 (x gcd) 370 (do ((x x (cddr x)) 371 (gcd gcd (cgcd gcd (pcontentz (cadr x))))) 372 ((or (null x) (equal gcd 1)) gcd))) 373 374(defun pcontentz (p) 375 (cond ((pcoefp p) p) 376 (t (pcontent1 (p-red p) (pcontentz (p-lc p)))))) 377 378(defun ucontent (p) ;CONTENT OF UNIV. POLY 379 (cond ((pcoefp p) (abs p)) 380 (t (setq p (mapcar #'abs (coefl (cdr p)))) 381 (let ((m (reduce #'min p))) 382 (oldcontent2 (delete m p :test #'equal) m))))) 383 384;;*** PGCDU CORRESPONDS TO BROWN'S ALGORITHM U 385 386;;;PGCDU IS NOT NOW IN RAT;UFACT > 387 388(defun pgcdu (p q) 389 (do () ((pzerop q) (monize p)) 390 (psetq p q q (pmodrem p q)))) 391 392(defun pmodrem (x y) 393 (cond ((null modulus) 394 (merror "PMODREM: null modulus; how did that happen?")) 395 ((pacoefp y) (if (pzerop y) x 0)) 396 ((pacoefp x) x) 397 ((eq (p-var x) (p-var y)) 398 (psimp (car x) (pgcdu1 (p-terms x) (p-terms y) nil))) 399 (t (merror "PMODREM: I can't handle this; x = ~M, y = ~M" x y)))) 400 401(defun pmodquo (u v &aux quo) 402 (declare (special quo)) 403 (cond ((null modulus) 404 (merror "PMODQUO: null modulus; how did that happen?")) 405 ((pcoefp v) (cons (ptimes (crecip v) u) 0)) 406 ((alg v) (cons (ptimes (painvmod v) u) 0)) 407 ((pacoefp u) (cons 0 u)) 408 ((not (eq (p-var u) (p-var v))) 409 (merror "PMODQUO: arguments have different variables; how did that happen?")) 410 (t (xcons (psimp (car u) (pgcdu1 (cdr u) (cdr v) t)) 411 (psimp (car u) quo))))) 412 413 414(defun pgcdu1 (u v pquo*) 415 (let ((invv (painvmod (pt-lc v))) (k 0) q*) 416 (declare (special k quo q*) (fixnum k)) 417 (loop until (minusp (setq k (- (pt-le u) (pt-le v)))) 418 do (setq q* (ptimes invv (pt-lc u))) 419 if pquo* do (setq quo (nconc quo (list k q*))) 420 when (ptzerop (setq u (ptpt-subtract-powered-product 421 (pt-red u) (pt-red v) q* k))) 422 return (ptzero) 423 finally (return u)))) 424 425;; it is convenient to have the *bigprimes* be actually less than 426;; half the size of the most positive fixnum, so that arithmetic is easier 427 428(defvar *bigprimes* (loop with p = (ash most-positive-fixnum -1) repeat 20 do 429 (setq p (next-prime (1- p) -1)) 430 collect p)) 431 432(defmvar *alpha (car *bigprimes*)) 433 434(defun newprime (p) 435 (do ((pl *bigprimes* (cdr pl))) 436 ((null pl) 437 (setq p (next-prime (1- p) -1)) 438 (setq *bigprimes* (nconc *bigprimes* (list p))) 439 p) 440 (when (< (car pl) p) 441 (return (car pl))))) 442 443(defun leadcoefficient (p) 444 (if (pcoefp p) p (leadcoefficient (caddr p)))) 445 446(defun maxcoefficient (p) 447 (if (pcoefp p) (abs p) (maxcoef1 (cdr p)))) 448 449(defun maxcoef1 (p) 450 (if (null p) 0 (max (maxcoefficient (cadr p)) (maxcoef1 (cddr p))))) 451 452(defun maxnorm (poly) 453 (if (null poly) 0 (max (norm (cadr poly)) (maxnorm (cddr poly))))) 454 455(defun norm (poly) 456 (cond ((null poly) 0) 457 ((pcoefp poly) (abs poly)) 458 (t (+ (norm (caddr poly)) (norm1 (cdddr poly)) )) )) 459 460(defun norm1 (poly) 461 (if (null poly) 0 (+ (norm (cadr poly)) (norm1 (cddr poly)) )) ) 462 463(defun pdegree (p var) 464 (cond ((pcoefp p) 0) 465 ((eq var (p-var p)) (p-le p)) 466 ((pointergp var (p-var p)) 0) 467 (t (do ((l (p-red p) (pt-red l)) 468 (e (pdegree (p-lc p) var) (max e (pdegree (pt-lc l) var)))) 469 ((null l) e))))) 470 471(defun poly-in-var (p v) 472 (cond ((or (pcoefp p) (pointergp v (p-var p))) (list 0 p)) 473 ((eq (p-var p) v) (p-terms p)) 474 ((loop with ans 475 for (exp coef) on (p-terms p) by #'cddr 476 do (setq ans (ptptplus ans 477 (everysubst2 (poly-in-var coef v) 478 (list (p-var p) exp 1)))) 479 finally (return ans))))) 480 481(defun univar (x) 482 (or (null x) (and (pcoefp (pt-lc x)) (univar (pt-red x))))) 483 484;;**THE CHINESE REMAINDER ALGORITHM IS A SPECIAL CASE OF LAGRANGE INTERPOLATION 485 486(defun lagrange3 (u uk p qk) 487 (set-modulus p) 488 (setq uk (pdifference uk (pmod u))) 489 (cond ((pzerop uk) (setq modulus nil) u) 490 (t (setq uk (pctimes (crecip (cmod qk)) uk)) 491 (setq modulus nil) 492 (pplus u (pctimes qk uk))))) 493 494 495(defun lagrange33 (u uk qk xk) 496 (declare (special xv)) 497 (setq uk (pdifference uk (pcsubst u xk xv))) 498 (cond ((pzerop uk) u) 499 (t (pplus u (ptimes 500 (pctimes (crecip (pcsubst qk xk xv)) uk) 501 qk))))) 502 503 504;;;************************************************************* 505 506;; THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 3. 507;; IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS 508