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 461