src/rat3d.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 1981 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module rat3d)
  14
  15 ;;      THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 4.
  16 ;;      IT INCLUDES THE POLYNOMIAL FACTORING ROUTINES.
  17
  18 (declare-top (special *odr* *checkagain))
  19
  20 (defmvar *irreds nil)
  21
  22 (defmvar $berlefact t)
  23
  24 (defmvar $factor_max_degree_print_warning t
  25   "Print a warning message when a polynomial is not factored because its
  26    degree is larger than $factor_max_degree?"
  27   boolean)
  28
  29 (defun listovars (q)
  30   (cond ((pcoefp q) nil)
  31         (t (let ((ans nil))
  32              (declare (special ans))
  33              (listovars0 q)))))
  34
  35 (defun listovars0 (q)
  36   (declare (special ans))
  37   (cond ((pcoefp q) ans)
  38         ((member (car q) ans :test #'eq) (listovars1 (cdr q)))
  39         (t (push (car q) ans)
  40            (listovars1 (cdr q)))))
  41
  42 (defun listovars1 (ql)
  43   (declare (special ans))
  44   (cond ((null ql) ans)
  45         (t (listovars0 (cadr ql)) (listovars1 (cddr ql)))))
  46
  47 (defun dontfactor (y)
  48   (cond ((or (null $dontfactor) (equal $dontfactor '((mlist)))) nil)
  49         ((memalike (pdis (make-poly y)) $dontfactor) t)))
  50
  51 (defun removealg (l)
  52   (loop for var in l
  53          unless (algv var) collect var))
  54
  55 (defun degvecdisrep (degl)
  56   (do ((l degl (cdr l))
  57        (gv genvar (cdr gv))
  58        (ans 1))
  59       ((null l) ans)
  60     (and (> (car l) 0)
  61          (setq ans (list (car gv) (car l) ans)))))
  62
  63 (defun ptermcont (p)
  64   (let ((tcont (degvecdisrep (pmindegvec p)))
  65         ($algebraic))
  66     (list tcont (pquotient p tcont))))
  67
  68 (defun pmindegvec (p)
  69   (minlist (let ((*odr* (putodr (reverse genvar)))
  70                  (nn* (1+ (length genvar)))
  71                  (*min* t))
  72              (degvector nil 1 p))))
  73
  74 (defun pdegreevector (p)
  75   (maxlist (let ((*odr* (putodr (reverse genvar)))
  76                  (nn* (1+ (length genvar)))
  77                  (*mx* t))
  78              (degvector nil 1 p))))
  79
  80 (defun maxlist(l) (maxminl l t))
  81
  82 (defun minlist(l) (maxminl l nil))
  83
  84 (defun maxminl (l switch)
  85   (do ((l1 (copy-list (car l)))
  86        (ll (cdr l) (cdr ll)))
  87       ((null ll) l1)
  88     (do ((v1 l1 (cdr v1))
  89          (v2 (car ll) (cdr v2)))
  90         ((null v1))
  91       (cond (switch
  92              (cond ((> (car v2) (car v1))
  93                     (rplaca v1 (car v2)))))
  94             (t (cond ((< (car v2) (car v1))
  95                       (rplaca v1 (car v2)))))))))
  96
  97 (defun quick-sqfr-check (p var)
  98   (let ((gv (delete var (listovars p) :test #'equal))
  99         (modulus (or modulus *alpha))
 100         (l) (p0))
 101     (if $algebraic (setq gv (removealg gv)))
 102     (and gv
 103          (not (pzerop (pcsubsty (setq l (rand (length gv) modulus))
 104                                 gv (pmod (p-lc p)))))
 105          (not (pcoefp (setq p0 (pcsubsty l gv (pmod p)))))
 106          (pcoefp (pgcd p0 (pderivative p0 (car p0))))
 107          (list l gv p0))))
 108
 109 (defun monom->facl (p)
 110   (cond ((pcoefp p) (if (equal p 1) nil (list p 1)))
 111         (t (list* (pget (car p)) (cadr p) (monom->facl (caddr p))))))
 112
 113 (defun psqfr (p)
 114   (prog (r varl var mult factors)
 115      (cond ((pcoefp p) (return (cfactor p)))
 116            ((pminusp p) (return (cons -1 (cons 1 (psqfr (pminus p)))))))
 117      (desetq (factors p) (ptermcont p))
 118      (setq factors (monom->facl factors))
 119      (cond ((pcoefp p) (go end)))
 120      (setq varl (sort (listovars p) 'pointergp))
 121      setvar
 122      (setq var (car varl) varl (cdr varl) mult 0)
 123      (cond ((pointergp var (car p)) (go nextvar))
 124            ((dontfactor var)
 125             (setq factors (cons p (cons 1 factors))
 126                   p 1)
 127             (go end)))
 128      (cond ((quick-sqfr-check p var)    ;QUICK SQFR CHECK BY SUBST.
 129             (setq r (oldcontent p))
 130             (setq p (car r) factors (cons (cadr r)
 131                                           (cons 1 factors)))
 132             (go nextvar)))
 133      (setq r (pderivative p var))
 134      (cond ((pzerop r) (go nextvar)))
 135      (cond ((and modulus (not (pcoefp r))) (pmonicize (cdr r))))
 136      (setq p (pgcdcofacts p r))
 137      (and algfac* (cadddr p) (setq adn* (ptimes adn* (cadddr p))))
 138      (setq r (cadr p)                   ; PRODUCT OF P[I]
 139            p (car p))
 140      a (setq r (pgcdcofacts r p)
 141              p (caddr r)
 142              mult (1+ mult))
 143      (and algfac* (cadddr r) (setq adn* (ptimes adn* (cadddr r))))
 144      (cond ((not (pcoefp (cadr r)))
 145             (setq factors
 146                   (cons (cadr r)
 147                         (cons mult factors)))))
 148      (cond ((not (pcoefp (setq r (car r)))) (go a)))
 149      nextvar
 150      (cond ((pcoefp p) (go end))
 151            (varl (go setvar))
 152            (modulus (setq factors (append (fixmult (psqfr (pmodroot p))
 153                                                    modulus)
 154                                           factors))
 155                     (setq p 1)))
 156      end  (setq p (cond ((equal 1 p) nil)
 157                         (t (cfactor p))))
 158      (return (append p factors))))
 159
 160 (defun fixmult (l n)
 161   (do ((l l (cddr l)))
 162       ((null l))
 163     (rplaca (cdr l) (* n (cadr l))))
 164   l)
 165
 166 (defun pmodroot (p)
 167   (cond ((pcoefp p) p)
 168         ((alg p) (pexpt p (expt modulus (1- (car (alg p))))))
 169         (t (cons (car p) (pmodroot1 (cdr p))))))
 170
 171 (defun pmodroot1 (x)
 172   (cond ((null x) x)
 173         (t (cons (truncate (car x) modulus)
 174                  (cons (pmodroot (cadr x))
 175                        (pmodroot1 (cddr x)))))))
 176
 177 (defun savefactors (l)
 178   (when $savefactors
 179     (savefactor1 (car l))
 180     (savefactor1 (cdr l)))
 181   l)
 182
 183 (defun savefactor1 (p)
 184   (unless (or (pcoefp p)
 185               (ptzerop (p-red p))
 186               (member p *checkfactors* :test #'equal))
 187     (push p *checkfactors*)))
 188
 189 (defun heurtrial1 (poly facs)
 190   (prog (h j)
 191      (setq h (pdegreevector poly))
 192      (cond ((or (member 1 h :test #'equal) (member 2 h :test #'equal)) (return (list poly))))
 193      (cond ((null facs) (return (list poly))))
 194      (setq h (pgcd poly (car facs)))
 195      (return (cond ((pcoefp h) (heurtrial1 poly (cdr facs)))
 196                    ((pcoefp (setq j (pquotient poly h)))
 197                     (heurtrial1 poly (cdr facs)))
 198                    (t (heurtrial (list h j) (cdr facs)))))))
 199
 200 (defun heurtrial (x facs)
 201   (cond ((null x) nil)
 202         (t (nconc (heurtrial1 (car x) facs)
 203                   (heurtrial (cdr x) facs)))))
 204
 205
 206 (defun pfactorquad (p)
 207   (prog (a b c d $dontfactor l v)
 208      (cond((or (onevarp p)(equal modulus 2))(return (list p))))
 209      (setq l (pdegreevector p))
 210      (cond ((not (member 2 l :test #'equal)) (return (list p))))
 211      (setq l (nreverse l) v (reverse genvar)) ;FIND MOST MAIN VAR
 212      loop (cond ((equal (car l) 2) (setq v (car v)))
 213                 (t (setq l (cdr l)) (setq v (cdr v)) (go loop)))
 214      (desetq (a . c) (bothprodcoef (make-poly v 2 1) p))
 215      (desetq (b . c) (bothprodcoef (make-poly v 1 1) c))
 216      (setq d (pgcd (pgcd a b) c))
 217      (cond ((pcoefp d) nil)
 218            (t (setq *irreds (nconc *irreds (pfactor1 d)))
 219               (return (pfactorquad (pquotient p d)))))
 220      (setq d (pplus (pexpt b 2) (ptimes -4 (ptimes a c))))
 221      (return
 222        (cond ((setq c (pnthrootp d 2))
 223               (setq d (ratreduce (pplus b c) (ptimes 2 a)))
 224               (setq d (pabs (pplus (ptimes (make-poly v) (cdr d))
 225                                    (car d))))
 226               (setq *irreds (nconc *irreds (list d (pquotient p d))))
 227               nil)
 228              (modulus (list p))    ;NEED TO TAKE SQRT(INT. MOD P) LCF.
 229              (t (setq *irreds (nconc *irreds (list p)))nil)))))
 230
 231 (defmfun $isqrt (x) ($inrt x 2))
 232
 233 (defmfun $inrt (x n)
 234   (cond ((not (integerp (setq x (mratcheck x))))
 235          (cond ((equal n 2) (list '($isqrt) x)) (t (list '($inrt) x n))))
 236         ((zerop x) x)
 237         ((not (integerp (setq n (mratcheck n)))) (list '($inrt) x n))
 238         (t (car (iroot (abs x) n)))))
 239
 240 (defun iroot (a n)   ; computes a^(1/n)  see Fitch, SIGSAM Bull Nov 74
 241   (cond ((< (integer-length a) n) (list 1 (1- a)))
 242         (t                              ;assumes integer a>0 n>=2
 243          (do ((x (expt 2 (1+ (truncate (integer-length a) n)))
 244                  (- x (truncate (+ n1 bk) n)))
 245               (n1 (1- n)) (xn) (bk))
 246              (nil)
 247            (cond ((signp le (setq bk (- x (truncate a (setq xn (expt x n1))))))
 248                   (return (list x (- a (* x xn))))))))))
 249
 250 (defmfun $nthroot (p n)
 251   (if (and (integerp n) (> n 0))
 252       (let ((k (pnthrootp (cadr ($rat p)) n)))
 253         (if k (pdis k) (merror (intl:gettext "nthroot: ~M is not a ~M power") p (format nil "~:r" n))))
 254     (merror (intl::gettext "nthroot: ~M is not a positive integer") n)))
 255
 256 (defun pnthrootp (p n)
 257   (ignore-rat-err (pnthroot p n)))
 258
 259 (defun pnthroot (poly n)
 260   (cond ((equal n 1) poly)
 261         ((pcoefp poly) (cnthroot poly n))
 262         (t (let* ((var (p-var poly))
 263                   (ans (make-poly var (cquotient (p-le poly) n)
 264                                   (pnthroot (p-lc poly) n)))
 265                   (ae (p-terms (pquotient (pctimes n (leadterm poly)) ans))))
 266              (do ((p (psimp var (p-red poly))
 267                      (pdifference poly (pexpt ans n))))
 268                  ((pzerop p) ans)
 269                (cond ((or (pcoefp p) (not (eq (p-var p) var))
 270                           (> (car ae) (p-le p)))
 271                       (rat-error "pnthroot error (should have been caught)")))
 272                (setq ans (nconc ans (ptptquotient (cdr (leadterm p)) ae)))
 273                )))))
 274
 275 (defun cnthroot(c n)
 276   (cond ((minusp c)
 277          (cond ((oddp n) (- (cnthroot (- c) n)))
 278                (t (rat-error "cnthroot error (should have been caught"))))
 279         ((zerop c) c)
 280         ((zerop (cadr (setq c (iroot c n)))) (car c))
 281         (t (rat-error "cnthroot error2 (should have been caught"))))
 282
 283
 284 (defun pabs (x) (cond ((pminusp x) (pminus x)) (t x)))
 285
 286 (defun pfactorlin (p l)
 287   (do ((degl l (cdr degl))
 288        (v genvar (cdr v))
 289        (a)(b))
 290       ((null degl) nil)
 291     (cond ((and (= (car degl) 1)
 292                 (not (algv (car v))))
 293            (desetq (a . b) (bothprodcoef (make-poly (car v)) p))
 294            (setq a (pgcd a b))
 295            (return (cons (pquotientchk p a)
 296                          (cond ((equal a 1) nil)
 297                                (t (pfactor1 a)))))))))
 298
 299
 300 (defun ffactor (l fn &aux (*alpha* *alpha*))
 301   ;;  (declare (special varlist))               ;i suppose...
 302   (prog (q)
 303      (cond ((and (null $factorflag) (mnump l)) (return l))
 304            ((or (atom l) algfac* modulus) nil)
 305            ((and (not gauss)(member 'irreducible (cdar l) :test #'eq))(return l))
 306            ((and gauss (member 'irreducibleg (cdar l) :test #'eq)) (return l))
 307            ((and (not gauss)(member 'factored (cdar l) :test #'eq))(return l))
 308            ((and gauss (member 'gfactored (cdar l) :test #'eq)) (return l)))
 309      (newvar l)
 310      (if algfac* (setq varlist (cons *alpha* (remove *alpha* varlist :test #'equal))))
 311      (setq q (ratrep* l))
 312      (when algfac*
 313        (setq *alpha* (cadr (ratrep* *alpha*)))
 314        (setq minpoly* (subst (car (last genvar))
 315                              (car minpoly*)
 316                              minpoly*)))
 317      (mapc #'(lambda (y z) (putprop y z (quote disrep)))
 318            genvar
 319            varlist)
 320      (return (retfactor (cdr q) fn l))))
 321
 322 (defun factorout1 (l p)
 323   (do ((gv genvar (cdr gv))
 324        (dl l (cdr dl))
 325        (ans))
 326       ((null dl) (list ans p))
 327     (cond ((zerop (car dl)))
 328           (t (setq ans (cons (pget (car gv)) (cons (car dl) ans))
 329                    p (pquotient p (list (car gv) (car dl) 1)))))))
 330
 331 (defun factorout (p)
 332   (cond ((and (pcoefp (ptterm (cdr p) 0))
 333               (not (zerop (ptterm (cdr p) 0))))
 334          (list nil p))
 335         (t (factorout1 (pmindegvec p) p))))
 336
 337 (defun pfactor (p &aux ($algebraic algfac*))
 338   (cond ((pcoefp p) (cfactor p))
 339         ($ratfac (pfacprod p))
 340         (t (setq p (factorout p))
 341            (cond ((equal (cadr p) 1) (car p))
 342                  ((numberp (cadr p)) (append (cfactor (cadr p)) (car p)))
 343                  (t (let ((cont (cond (modulus (list (leadalgcoef (cadr p)) (monize (cadr p))))
 344                                       (algfac* (algcontent (cadr p)))
 345                                       (t (pcontent (cadr p))))))
 346                       (nconc
 347                        (cond ((equal (car cont) 1) nil)
 348                              (algfac*
 349                               (cond (modulus (list (car cont) 1))
 350                                     ((equal (car cont) '(1 . 1)) nil)
 351                                     ((equal (cdar cont) 1) (list (caar cont) 1))
 352                                     (t (list (caar cont) 1 (cdar cont) -1))))
 353                              (t (cfactor (car cont))))
 354                        (pfactor11 (psqfr (cadr cont)))
 355                        (car p))))))))
 356
 357 (defun pfactor11 (p)
 358   (cond ((null p) nil)
 359         ((numberp (car p))
 360          (cons (car p) (cons (cadr p) (pfactor11 (cddr p)))))
 361         (t (let* ((adn* 1)
 362                   (f (pfactor1 (car p))))
 363              (nconc (if (equal adn* 1) nil
 364                         (list adn* (- (cadr p))))
 365                     (do ((l f (cdr l))
 366                          (ans nil (cons (car l) (cons (cadr p) ans))))
 367                         ((null l) ans))
 368                     (pfactor11 (cddr p)))))))
 369
 370 (defun pfactor1 (p)                     ;ASSUMES P SQFR
 371   (prog (factors *irreds *checkagain)
 372      (cond ((dontfactor (car p)) (return (list p)))
 373            ((and (not (zerop $factor_max_degree)) (> (apply 'max (pdegreevector p)) $factor_max_degree))
 374                  (when $factor_max_degree_print_warning
 375                    (mtell (intl:gettext "Refusing to factor polynomial ~M because its degree exceeds factor_max_degree (~M)~%") (pdis p) $factor_max_degree))
 376                  (return (list p)))
 377            ((onevarp p)
 378             (cond ((setq factors (factxn+-1 p))
 379                    (if (and (not modulus)
 380                             (or gauss (not algfac*)))
 381                        (setq *irreds factors
 382                              factors nil))
 383                    (go out))
 384                   ((and (not algfac*) (not modulus)
 385                         (not (equal (cadr p) 2)) (estcheck (cdr p)))
 386                    (return (list p))))))
 387      (and (setq factors (pfactorlin p (pdegreevector p)))
 388           (return factors))
 389      (setq factors(if (or algfac* modulus) (list p) ;SQRT(NUM. CONT OF DISC)
 390                       (pfactorquad p)))
 391      (cond ((null factors)(go out)))
 392      (when *checkfactors*
 393        (setq factors (heurtrial factors *checkfactors*))
 394        (setq *checkagain (cdr factors)))
 395      out (return (nconc *irreds (mapcan (function pfactorany) factors)))))
 396
 397 (defmvar $homog_hack nil)  ; If T tries to eliminate homogeneous vars.
 398
 399 (declare-top (special *hvar *hmat))
 400
 401 (defun pfactorany (p)
 402   (cond (*checkagain (let (*checkfactors*) (pfactor1 p)))
 403         ((and $homog_hack (not algfac*) (not (onevarp p)))
 404          (let ($homog_hack *hvar *hmat)
 405            (mapcar #'hexpand (pfactor (hreduce p)))))
 406         ($berlefact (factor1972 p))
 407         (t (pkroneck p))))
 408
 409
 410 (defun retfactor (x fn l &aux (a (ratfact x fn)))
 411   (prog ()
 412    b    (cond ((null (cddr a))
 413                (setq a (retfactor1 (car a) (cadr a)))
 414                (return (cond ((and scanmapp (not (atom a)) (not (atom l))
 415                                    (eq (caar a) (caar l)))
 416                               (tagirr l))
 417                              (t a))))
 418               ((equal (car a) 1) (setq a (cddr a)) (go b))
 419               (t (setq a (map2c #'retfactor1 a))
 420                  (return (cond ((member 0 a) 0)
 421                                (t (setq a (let (($expop 0) ($expon 0)
 422                                                 $negdistrib)
 423                                             (muln (sortgreat a) t)))
 424                                   (cond ((not (mtimesp a)) a)
 425                                         (t (cons '(mtimes simp factored)
 426                                                  (cdr a)))))))))))
 427
 428 ;;; FOR LISTS OF ARBITRARY EXPRESSIONS
 429 (defun retfactor1 (p e)
 430   (power (tagirr (simplify (pdisrep p))) e))
 431
 432 (defun tagirr (x)
 433   (cond ((or (atom x) (member 'irreducible (cdar x) :test #'eq)) x)
 434         (t (cons (append (car x) '(irreducible)) (cdr x)))))
 435
 436 (defun revsign (x)
 437   (cond ((null x) nil)
 438         (t (cons (car x)
 439                  (cons (- (cadr x)) (revsign (cddr x)))))))
 440
 441 ;;      THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 4