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