src/rat3c.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 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 (defmfun 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 (defmfun 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 (defmfun 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 (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)) (exch 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 (defmfun 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 (defmfun 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 (apply #'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 (defmfun 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 (defun rtzerl2 (n)
 426   (cond ((zerop n) 0)
 427         (t (do ((n n (ash n -2)))
 428                ((not (zerop (haipart n -2))) n)))))
 429
 430 (defmfun $jacobi (p q)
 431   (cond ((null (and (integerp p) (integerp q)))
 432          (list '($jacobi) p q))
 433         ((zerop q) (merror (intl:gettext "jacobi: zero denominator.")))
 434         ((minusp q) ($jacobi p (- q)))
 435         ((and (evenp (setq q (rtzerl2 q)))
 436               (setq q (ash q -1))
 437               (evenp p)) 0)
 438         ((equal q 1) 1)
 439         ((minusp (setq p (rem p q)))
 440          (jacobi (rtzerl2 (+ p q)) q))
 441         (t (jacobi (rtzerl2 p) q))))
 442
 443 (defun jacobi (p q)
 444   (do ((r1 p (rtzerl2 (rem r2 r1)))
 445        (r2 q r1)
 446        (bit2 (haipart q -2))
 447        (odd 0 (boole boole-xor odd (boole boole-and bit2 (setq bit2 (haipart r1 -2))))))
 448       ((zerop r1) 0)
 449     (cond ((evenp r1)
 450            (setq r1 (ash r1 -1))
 451            (setq odd (boole boole-xor odd (ash (expt (haipart r2 -4) 2) -2)))))
 452     (and (equal r1 1) (return (expt -1 (boole  boole-and 1 (ash odd -1)))))))
 453
 454 ;; it is convenient to have the *bigprimes* be actually less than
 455 ;; half the size of the most positive fixnum, so that arithmetic is easier
 456
 457 (defvar *bigprimes* (loop with p = (ash most-positive-fixnum -1) repeat 20 do
 458                          (setq p (next-prime (1- p) -1))
 459                        collect p))
 460
 461 (defmvar *alpha (car *bigprimes*))
 462
 463 (defun newprime (p)
 464   (do ((pl *bigprimes* (cdr pl)))
 465       ((null pl)
 466        (setq p (next-prime (1- p) -1))
 467        (setq *bigprimes* (nconc *bigprimes* (list p)))
 468        p)
 469     (when (< (car pl) p)
 470       (return (car pl)))))
 471
 472 (defun leadcoefficient (p)
 473   (if (pcoefp p) p (leadcoefficient (caddr p))))
 474
 475 (defun maxcoefficient (p)
 476   (if (pcoefp p) (abs p) (maxcoef1 (cdr p))))
 477
 478 (defun maxcoef1 (p)
 479   (if (null p) 0 (max (maxcoefficient (cadr p)) (maxcoef1 (cddr p)))))
 480
 481 (defun maxnorm (poly)
 482   (if (null poly) 0 (max (norm (cadr poly)) (maxnorm (cddr poly)))))
 483
 484 (defun norm (poly)
 485   (cond ((null poly) 0)
 486         ((pcoefp poly) (abs poly))
 487         (t (+ (norm (caddr poly)) (norm1 (cdddr poly)) )) ))
 488
 489 (defun norm1 (poly)
 490   (if (null poly) 0 (+ (norm (cadr poly)) (norm1 (cddr poly)) )) )
 491
 492 (defmfun pdegree (p var)
 493   (cond ((pcoefp p) 0)
 494         ((eq var (p-var p)) (p-le p))
 495         ((pointergp var (p-var p)) 0)
 496         (t (do ((l (p-red p) (pt-red l))
 497                 (e (pdegree (p-lc p) var) (max e (pdegree (pt-lc l) var))))
 498                ((null l) e)))))
 499
 500 (defun poly-in-var (p v)
 501   (cond ((or (pcoefp p) (pointergp v (p-var p))) (list 0 p))
 502         ((eq (p-var p) v) (p-terms p))
 503         ((loop with ans
 504                 for (exp coef) on (p-terms p) by #'cddr
 505                 do (setq ans (ptptplus ans
 506                                        (everysubst2 (poly-in-var coef v)
 507                                                     (list (p-var p) exp 1))))
 508                 finally (return ans)))))
 509
 510 (defun univar (x)
 511   (or (null x) (and (pcoefp (pt-lc x)) (univar (pt-red x)))))
 512
 513 ;;**THE CHINESE REMAINDER ALGORITHM IS A SPECIAL CASE OF LAGRANGE INTERPOLATION
 514
 515 (defun lagrange3 (u uk p qk)
 516   (set-modulus p)
 517   (setq uk (pdifference uk (pmod u)))
 518   (cond ((pzerop uk) (setq modulus nil) u)
 519         (t (setq uk (pctimes (crecip (cmod qk)) uk))
 520            (setq modulus nil)
 521            (pplus u (pctimes qk uk)))))
 522
 523
 524 (defun lagrange33 (u uk qk xk)
 525   (declare (special xv))
 526   (setq uk (pdifference uk (pcsubst u xk xv)))
 527   (cond ((pzerop uk) u)
 528         (t (pplus u (ptimes
 529                      (pctimes (crecip (pcsubst qk xk xv)) uk)
 530                      qk)))))
 531
 532
 533 ;;;*************************************************************
 534
 535 ;;      THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 3.
 536 ;;      IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS