src/factor.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 factor)
  14
  15 ;;; This is the FACTOR package.
  16
  17 ;;;      THIS IS THE NEW FACTORING PACKAGE. THE FUNCTION
  18 ;;;      FACTOR72 TAKES A PRIMITIVE SQUARE-FREE POLY AS INPUT THE OUTPUT IS A
  19 ;;;      LIST OF FACTORS THE FUNCTION FACTOR1972 IS ABOVE FACTOR72 AND IT
  20 ;;;      TAKES CARE OF REPEATED FACTORS OVER THE GAUSSIAN INTEGERS BEFORE
  21 ;;;      CALLING FACTOR72 THE FUNCTION Z1 TAKES TWO FACTORS IN ONE VARIABLE
  22 ;;;      AND ONE POLY IN SEVERAL VARIABLES AS INPUT Z1 TAKES THESE FACTORS IN
  23 ;;;      ONE VARIABLES AND BUILDS OUT OF THEM TWO FACTORS OF THE GIVEN POLY IN
  24 ;;;      SEVERAL VARIABLES
  25
  26 (load-macsyma-macros ratmac)
  27
  28 (declare-top (special *stop* trl* *xn sharpcont subvar1 anotype invc fctc
  29                       subval1 var mcflag alcinv *ab* monic* intbs*
  30                       *prime *g* modulu* plim listelm many* *inl3
  31                       *sharpa *sharpb limk split* alc ind p l dosimp *odr*
  32                       *i* mcflag elm ne res fact1 fact2 subvar
  33                       subval ovarlist valist dlp nn* df1 df2 dn* fcs* uu*))
  34
  35 (defvar *afixn*)
  36 (defvar *fctcfixn*)
  37 (defvar *invcfixn*)
  38
  39 (defvar *negflag*)
  40
  41 ;; Internal specials
  42
  43 (defmvar gauss nil)
  44 (defmvar *min* nil)
  45 (defmvar *mx* nil)
  46 (defmvar minpoly* nil)
  47 (defmvar mplc* nil)
  48 (defmvar mm* 1)
  49 (defmvar alpha nil)
  50 ;(defmvar smallprimes '(3 5 7 11. 13. 17. 19. 23. 29. 31. 37. 41. 43. 47. 53. 59. 61.))
  51
  52 ;; External specials
  53
  54 (defmvar $nalgfac t "If t use bmt's algebraic factoring algorithm")
  55
  56 (defun primcyclo (n)
  57   (let ((*g* (gensym "PRIMCYCLO-"))
  58         (nl (loop for (c e) on (cfactorw n) by #'cddr
  59                nconc (make-list e :initial-element c))))
  60     (setf (symbol-value *g*) 0)
  61     (let ((res (cyclotomic (list n nl))))
  62       (cond ((consp res) (p-terms res))
  63             ((eql 0 res) nil)
  64             (t (list 0 res))))))
  65
  66 (defun factxn+-1 (p)
  67   (let ((*g* (car p))
  68         ($factorflag t))
  69     (cond ((equal 1 (cadr p)) (list p))
  70           ((equal (cddr p) '(1 0 1))
  71            (factxn+1 (cadr p)))
  72           ((equal (cddr p) '(1 0 -1))
  73            (factxn-1 (cadr p))))))
  74
  75 (defun cfactorw (n) (let (($factorflag t)) (cfactor n)))
  76
  77 (defun factxn-1 (n)
  78   (cond ((evenp n)
  79          (append (factxn-1 (ash n -1)) (factxn+1 (ash n -1))))
  80         (t (mapcar #'cyclotomic (divisors (cfactor n))))))
  81
  82 (defun factxn+1 (n)
  83   (cond (gauss
  84          (let* ((gauss nil) (facl (factxn+1 n)))
  85            (cond ((oddp n) facl)
  86                  (t (let (($gcd '$subres)
  87                           (pfac (list *g* (ash n -1) 1 0 alpha)))
  88                       (mapcan #'(lambda (q) (subseq (pgcdcofacts q pfac) 0 2)) facl))))))
  89         (t (let ((m 1) (nl (reverse (cfactor n))))
  90              (when (equal 2 (cadr nl))
  91                (setq m (expt 2 (car nl)))
  92                (setq nl (cddr nl)))
  93              (setq m (list *g* m -1))
  94              (if (null nl) (ncons (list *g* n 1 0 1))
  95                  (mapcar #'(lambda (p) (pabs (pcsubst p m (car p))))
  96                          (mapcar #'cyclotomic (divisors (reverse nl)))))))))
  97
  98
  99 (defun cyclp (n ind)
 100   (loop for i downfrom (1- n) to 0
 101          nconc (list (* ind i) 1)))
 102
 103 (defun csf (l)
 104   (cond ((null l) nil) (t (list* (car l) 1 (csf (cdr l))))))
 105
 106 (defun condense (l)
 107   (cond ((null (cdr l)) l)
 108         ((eq (car l) (cadr l)) (condense (cdr l)))
 109         (t (cons (car l) (condense (cdr l))))))
 110
 111 (defun cyclotomic (nl)
 112   (prog (n dp dpl num den p)
 113      (cond ((equal 1 (car nl)) (return (list *g* 1 1 0 -1)))
 114            ((null (cdr (setq p (condense (cadr nl)))))
 115             (return (cons *g*
 116                           (cyclp (car p)
 117                                  (expt (car p) (1- (length (cadr nl)))))))))
 118      (setq num 1 den 1 n (car nl) dpl (divisors (csf p)))
 119      loop (cond ((null dpl) (return (pquotient num den))))
 120      (setq dp (car dpl))
 121      (setq dpl (cdr dpl))
 122      (setq p (list *g* (quotient n (car dp)) 1 0 -1))
 123      (cond ((or (evenp (length (cadr dp))) (equal (car dp) 1))
 124             (setq num (ptimes p num)))
 125            (t (setq den (ptimes p den))))
 126      (go loop)))
 127
 128 (defun divisors (l)
 129   (if (equal l '(1 1)) (setq l nil))
 130   (do ((ans (list '(1 ()) ))
 131        (l l (cddr l)))
 132       ((null l) ans)
 133     (do ((u ans)
 134          (factor (car l))
 135          (mult (cadr l) (1- mult)))
 136         ((zerop mult))
 137       (setq u (mapcar #'(lambda (q) (list (* factor (car q))
 138                                           (cons factor (cadr q))))
 139                       u))
 140       (setq ans (nconc ans u)))))
 141
 142
 143 (defun estcheck2 (d lc c)
 144   (prog (p)
 145    loop (cond ((null d) (return nil)))
 146    (setq p (car d) d (cdr d))
 147    (cond ((or (and (not (equal (rem c p) 0))
 148                    (not (equal (rem lc (* p p)) 0)))
 149               (and (not (equal (rem lc p) 0))
 150                    (not (equal (rem c (* p p)) 0))))
 151           (return t)))
 152    (go loop)))
 153
 154 (defun estcheck (p)
 155   (prog (lc c d)
 156      (cond ((or (atom p) (null (cddr p)) (equal (ptterm p 0) 0))
 157             (return nil)))
 158      (setq lc (cadr p))
 159      (setq p (nreverse (cdr (oddelm (cdr p)))))
 160      (setq c (car p))
 161      (setq d (cgcdlist p))
 162      (cond ((equal 1 d) (return nil)))
 163      (setq d (oddelm (cfactorw d)))
 164      (return (estcheck2 d lc c))))
 165
 166
 167 (defun cgcdlist (l)
 168   (cond ((null l) nil)
 169         ((null (cdr l)) (abs (car l)))
 170         ((or (member 1 l) (member -1 l)) 1)
 171         ((null (cddr l)) (gcd (car l) (cadr l)))
 172         (t (cgcdlist (cons (gcd (car l) (cadr l)) (cddr l))))))
 173
 174 (defun dropterms (p)
 175   (prog (ans c)
 176      (cond ((atom p) (return p))
 177            ((not (eq (car p) var)) (return (kterms p dlp))))
 178      (setq ans (cons (car p) ans) p (cdr p))
 179      loop (cond ((null p) (return (cond ((cdr ans) (nreverse ans)) (t 0)))))
 180      (setq c (kterms (cadr p) dlp))
 181      (cond ((not (equal c 0)) (setq ans (cons c (cons (car p) ans)))))
 182      (setq p (cddr p))
 183      (go loop)))
 184
 185
 186 (defun restorelc (l lc)
 187   (prog (h r ans var c d deg)
 188      (cond ((equal 1 lc)
 189             (cond ((and (not many*) algfac* (not (equal intbs* 1)))
 190                    (return (mapcar #'intbasehk l)))
 191                   (t (return (reverse l))))))
 192      (setq r (lcprodl l) h 1)
 193      loop (cond ((null l) (return ans)))
 194      (setq d (car l) l (cdr l) var (car d) deg (cadr d) c (caddr d))
 195      (setq d (ptimes (ptimes h (car r)) (psimp var (cdddr d))))
 196      (cond (many* (setq d (dropterms d))))
 197      (setq d (pplus (list var deg lc)d))
 198      (cond ((and (not many*) algfac* (not (equal intbs* 1)))
 199             (setq d (intbasehk d))))
 200      (let ((modulus))
 201        (setq ans (cons (cadr (oldcontent d)) ans)))
 202      (setq h (ptimes h c) r (cdr r))
 203      (go loop)))
 204
 205 (defun iredup (p)
 206   (let ((mm* 1) (algfac*))
 207     (cond ((sqfrp p(car p))
 208            (setq p (catch 'splt (cpber1 p)))
 209            (and (null (car p)) (null (cdadr p)))))))
 210
 211 (defun zerolp (a)
 212   (every #'zerop1 a))
 213
 214 ;;; TESTDIVIDE
 215 ;;;
 216 ;;; Check if y divides x, assuming that x and y are polynomials in the
 217 ;;; same variable.
 218 (defun testdivide (x y)
 219   (if algfac*
 220       (algtestd x y)
 221       (eztestdivide x y)))
 222
 223 ;;; TESTDIVIDE*
 224 ;;;
 225 ;;; Check if y divides x, for general polynomials x and y (not
 226 ;;; necessarily in the same variable).
 227 (defun testdivide* (x y)
 228   (if algfac*
 229       (algtestd x y)
 230       (ignore-rat-err (pquotient x y))))
 231
 232 (defun algtestd (x y)
 233   (and (div-deg-chk (nreverse (pdegreevector x)) (nreverse (pdegreevector y))
 234                     (reverse genvar))
 235        (cond ((setq x (ignore-rat-err (rquotient x y)))
 236               (setq adn* (* adn* (cdr x)))
 237               (car x)) )))
 238
 239 (defun div-deg-chk (xl yl gl)
 240   (cond ((or (null gl) (algv (car gl))) t)
 241         ((> (car yl) (car xl)) nil)
 242         (t (div-deg-chk (cdr xl) (cdr yl) (cdr gl)))))
 243
 244 ;; FUU is used by systems programmers such as BMT and PAULW while debugging.
 245 (defun fuu nil
 246   (setq tellratlist nil varlist nil genvar nil genpairs nil))
 247
 248 (defun linout (u)
 249   (prog (m linfac x y)
 250      (setq y (list (setq x (car u)) 1 1) m modulus)
 251      loop (decf m)
 252      (cond ((< m 0) (return (list u linfac)))
 253            ((equal (cadr u) 1) (return (list 1 (cons u linfac))))
 254            ((zerop (pcsubsty (cmod m) x u))
 255             (setq linfac
 256                   (cons (append y
 257                                 (cond ((zerop m) nil)
 258                                       (t (list 0 (cmod (- m))))))
 259                         linfac))
 260             (setq u (car (pmodquo u (car linfac))))))
 261      (go loop)))
 262
 263 (defun onevarp (p)
 264   (if algfac* (every #'pacoefp (cdr p))
 265       (every #'numberp (cdr p))))
 266
 267 (defun putodr (l)
 268   (do ((l l (cdr l))
 269        (i 1 (1+ i))
 270        (ans))
 271       ((null l) ans)
 272     (push (cons (car l) i) ans)))
 273
 274 (defun kterms (p k)
 275   (cond ((pacoefp p) p)
 276         ((= k 0) (consta p))
 277         (t (prog (v ans w)
 278               (setq v (car p))
 279               (setq p (cdr p))
 280               loop (cond ((null p) (return 0))
 281                          ((> (car p) k) (setq p (cddr p)) (go loop)))
 282               ag   (cond ((null p)
 283                           (return (psimp v ans))))
 284               (setq w (kterms (cadr p) (- k (car p))))
 285               (cond ((not (pzerop w))
 286                      (setq ans (nconc ans (list (car p) w)))))
 287               (setq p (cddr p))
 288               (go ag)))))
 289
 290 (defun consta (p)
 291   (cond ((or (pcoefp p) (alg p)) p)
 292         (t (consta (ptterm (cdr p) 0)))))
 293
 294 (defun constacl (p)                     ;NO LONGER USED?
 295   (cond ((atom p)
 296          (cond ((equal p 1) (throw 'cnt 1))
 297                (t (list p))))
 298         ((every #'numberp (cdr p))
 299          (setq p (oddelm p))
 300          (cond ((member 1 p) (throw 'cnt 1))
 301                (t (cdr p))))
 302         (t (apply #'append (mapcar #'constacl (cdr (oddelm p)))))))
 303
 304 (defun z1 (poly fact1 fact2)
 305   (prog (res hsteps steps kterm a b c d *ab* m df1 df2 dlr step *sharpa *sharpb)
 306      (let ((modulus))
 307        (set-modulus *prime)
 308        (setq *sharpb (fact20 fact1 fact2 limk)))
 309      (setq *sharpa (car *sharpb))
 310      (setq *sharpb (cadr *sharpb))
 311      (setq *ab* (list (list 0 *sharpa *sharpb)))
 312      (setq steps dlp
 313            hsteps (ash steps -1))
 314      (setq res (pdifference (ptimes (pmod fact1) (pmod fact2)) (pmod poly)))
 315      (setq poly nil)
 316      (setq step 0)
 317      (setq df1 fact1)
 318      (setq df2 fact2)
 319      loop (cond ((equal res 0) (go out)))
 320      (incf step)
 321      (cond ((> step steps) (go out)))
 322      (cond ((eq (car res) var) (setq c (cdr res)))
 323            (t (setq c (list 0 res))))
 324      (setq a 0 b 0)
 325      nextm (cond ((null c) (z2 a b step hsteps) (go loop)))
 326      (setq m (car c) dlr (cadr c))
 327      (setq c (cddr c))
 328      (setq kterm (kterms dlr step) dlr nil)
 329      (cond ((equal 0 kterm) (go nextm)))
 330      (setq d (obtainabm m))
 331      (setq b (pplus b (ptimes (car d) kterm))
 332            a (pplus a (ptimes (cadr d) kterm))
 333            kterm nil)
 334      (go nextm)
 335      out  (return (list df1 df2))))
 336
 337 (defun z2 (a b step hsteps)
 338   (unless (and (equal a 0) (equal b 0))
 339     (setq step
 340           (pdifference
 341            (pdifference (cond ((not (< step hsteps))
 342                                (dropterms (ptimes a b)))
 343                               (t (ptimes a b)))
 344                         (cond ((not (< step hsteps))
 345                                (dropterms (ptimes df1 b)))
 346                               (t (ptimes df1 b))))
 347            (cond ((not (< step hsteps))
 348                   (dropterms (ptimes df2 a)))
 349                  (t (ptimes df2 a)))))
 350     (setq res (pplus res step))
 351     (setq df1 (pdifference df1 a))
 352     (setq df2 (pdifference df2 b))))
 353
 354 (defun obtainabm (m)
 355   (prog (ans)
 356      (cond ((setq ans (cdr (assoc m *ab* :test #'equal))) (return ans)))
 357      (setq ans (obtainab (list var m 1)))
 358      (setq *ab* (cons (cons m ans) *ab*))
 359      (return ans)))
 360
 361 (defun fact20 (f1 g1 limk)
 362   (prog (f g a pk b reml qlp h k b1)
 363      (setq k 0)
 364      (setq reml (ppprog (pmod f1) (pmod g1)))
 365      (setq a (car reml))
 366      (setq b (cadr reml))
 367      sharp   (cond ((> k limk) (return (list a b))))
 368      (setq pk modulus)
 369      (set-modulus (* modulus modulus))
 370      (setq f(pmod f1) g (pmod g1))
 371      (setq h (pquo (pmod (pdifference (pplus (ptimes a f) (ptimes b g))
 372                                       1))
 373                    pk))
 374      (setq qlp (pmodquo (ptimes a h) g))
 375      (setq b1 (pplus (ptimes b h) (ptimes (car qlp) f)))
 376      (setq a (pdifference a (pmod (pctimes pk (cdr qlp)))))
 377      (setq b (pdifference b (pmod (pctimes pk b1))))
 378      (incf k)
 379      (go sharp)))
 380
 381 (defun baselist (n)
 382   (setq *i* n)
 383   (completevector nil 0 n elm))
 384
 385 (defun inlist3 (l)
 386   (cond ((null l) (setq *inl3 nil))
 387         ((zerop (car l)) (cons 1 (cdr l)))
 388         ((equal (car l) 1) (cons -1 (cdr l)))
 389         (t (cons 0 (inlist3 (cdr l))))))
 390
 391 (defun newrep (p)
 392   (let ((modulus))
 393     (if subvar (pcsubsty (mapcar #'(lambda (a b) (list a 1 1 0 b))
 394                                  subvar subval)
 395                          subvar
 396                          p)
 397         p)))
 398
 399 (defun oldrep (p)
 400   (let ((modulus))
 401     (if subvar (pcsubsty (mapcar #'(lambda (a b) (list a 1 1 0 (- b)))
 402                                  subvar subval)
 403                          subvar
 404                          p)
 405         p)))
 406
 407 (defun completevector (l n m v)
 408   (do ((i m (1- i)))
 409       ((= i n) l)
 410     (push v l)))
 411
 412 ;; special variable *odr* contains order of nesting of variables in c
 413 (defun degvector (l n c)
 414   (prog (lf ans j)
 415    bk (cond ((numberp c)
 416              (return  (list (completevector l n nn* 0)))))
 417    (setq j (cdr (assoc (car c) *odr* :test #'equal)))
 418              ;;; IN CASE (CAR C) IS ALGEBRAIC
 419    (cond ((null j) (setq c 0)(go bk)))
 420    (setq c (cdr c))
 421    (setq lf (completevector l n j 0))
 422    loop (cond ((null c) (return ans)))
 423    (setq ans
 424          (nconc (degvector (cons (car c) lf) (1+ j) (cadr c)) ans))
 425    (cond (*mx* (setq ans (ncons (maxlist ans))))
 426          (*min* (setq ans (ncons (minlist ans)))))
 427    (setq c (cddr c))
 428    (go loop)))
 429
 430 (defun union1 (a b)
 431   (do ((a a (cdr a))
 432        (ans b))
 433       ((null a) ans)
 434     (or (member (car a) ans :test #'equal)
 435         (setq ans (cons (car a) ans)))))
 436
 437 (defun obtainab (u)
 438   (prog (c ql)
 439      (setq c (pmod u))
 440      (setq ql (pmodquo (ptimes *sharpa c) fact2))
 441      (return (list (cdr ql) (pmod (pplus (ptimes (car ql) fact1)
 442                                          (ptimes *sharpb c)))))))
 443
 444 (defun pcdifconc (v j)
 445   (do ((l v (cddr l)))
 446       ((null (cdr l))
 447        (or (= j 0)
 448            (rplacd l (list 0 j)))
 449        v)
 450     (cond ((= (cadr l) 0)
 451            (cond ((= j 0)
 452                   (rplacd l nil))
 453                  ((rplaca (cddr l) j)))
 454            (return v)))))
 455
 456 (defun orde (a l)
 457   (cond ((null l) (list a))
 458         (t (cond ((< a (car l)) (cons a l))
 459                  (t (cons (car l) (orde a (cdr l))))))))
 460
 461 (defun pquo (x y)
 462   (let (modulus)
 463     (pquotient x y)))
 464
 465 ;; maintains order of list x
 466 (defun intersect (x y)
 467   (if x (if (member (car x) y :test #'equal)
 468             (cons (car x) (intersect (cdr x) y))
 469             (intersect (cdr x) y))))
 470
 471 ;; Like APL IOTA function.
 472 (defun index* (k)
 473   (declare (fixnum k))
 474   (if (< k 2) (list 1) (cons k (index* (f1- k)))))
 475
 476 ; sets plim to a power of p1 likely to be a large enough modulus for polynomial u
 477 ; returns limk, where plim = p1^(2^(limk+1))
 478 (defun klim (u p1)
 479   (prog (bcoef)
 480      (setq bcoef (* 5 (maxcoefficient u)))
 481      (when algfac* (setq bcoef (* bcoef intbs*)))
 482      (when (< bcoef 10000.) (setq bcoef 20000.))
 483      (setq limk 0)
 484      test (setq p1 (* p1 p1))
 485      (when (> p1 bcoef)
 486        (setq plim p1)
 487        (return limk))
 488      (incf limk)
 489      (go test)))
 490
 491 (declare-top (special b b2))
 492
 493 (defun cpberl (u)
 494   (prog (ql d)
 495      (setq ql (catch 'splt (cpber1 u)) u (caddr ql))
 496      (setq d (car ql) ql (cadr ql))
 497      (cond ((null ql)(return d))
 498            ((null (cdr ql)) (return (cons u d))))
 499      (return (append d
 500                      (cond ((or alpha (> modulus 70.))
 501                             (cpbgzass ql (pmod u) (length ql)))
 502                            (t (cpbg ql (pmod u) (length ql))))))))
 503
 504 ;; Returns a list of monomials in G of degree less than N.
 505 (defun powrs (g n &aux (ans (ncons 1)))
 506   (declare (fixnum n))
 507   (do ((i 1 (1+ i))) ((= i n) ans)
 508     (declare (fixnum i))
 509     (push (make-poly g i 1) ans)))
 510
 511 ;; Finds polynomials A and B such that A*F+B*G=1 when MODULUS
 512 ;; is non-NIL.  Same algorithm as INVMOD.
 513 (defun ppprog (f g)
 514   (prog (a1 a2 b1 b2 r1 r2 ql ans ap bp g1 f1 s)
 515      (cond ((> (cadr g) (cadr f)) (setq g1 g) (setq f1 f))
 516            (t (setq g1 f) (setq f1 g) (setq s t)))
 517      (setq ql (pmodquo g1 f1))
 518      (setq a1 1)
 519      (setq b1 0)
 520      (setq a2 (pminus (car ql)))
 521      (setq b2 1)
 522      (setq r1 f1)
 523      (setq r2 (cdr ql))
 524      test (cond ((or (numberp r2) (and alpha (alg r2))) (go end)))
 525      (setq ql (pmodquo r1 r2))
 526      (setq ap (pdifference a1 (ptimes (car ql) a2)))
 527      (setq bp (pdifference b1 (ptimes (car ql) b2)))
 528      (setq r1 r2)
 529      (setq r2 (cdr ql))
 530      (setq a1 a2)
 531      (setq b1 b2)
 532      (setq a2 ap)
 533      (setq b2 bp)
 534      (go test)
 535      end  (cond ((pzerop r2)
 536                  (cond ((equal 1 (setq ans (caddr r1)))
 537                         (setq ans (list b1 a1)))
 538                        (t (setq ans (list (car (pmodquo b1 ans))
 539                                           (car (pmodquo a1 ans))))))
 540                  (go out)))
 541      (setq ans (list (car (pmodquo b2 r2)) (car (pmodquo a2 r2))))
 542      out  (cond ((not s) (return (reverse ans))) (t (return ans)))))
 543
 544
 545 (defun zff (v f g)
 546   (if many*
 547       (z1 v f g)
 548       (fact2z v f g limk)))
 549
 550 (defun zfact (u fl limk many*)
 551   (prog (fcs* prodl)
 552      (when many*
 553        (set-modulus plim)
 554        (setq dlp (reduce #'max (mapcar #'multideg (cdr (oddelm u))))))
 555      (when (= (length fl) 1) (return (list u)))
 556      (setq prodl (fsplit fl 'v))
 557      (zfactsplit prodl u)
 558      (return fcs*)))
 559
 560 (defun zfactsplit (fl v)
 561   (prog (d)
 562      (cond ((null (cdr fl)) (return (setq fcs* (cons v fcs*))))
 563            ((null (cddr fl))
 564             (setq fl (cadr fl))
 565             (return (setq fcs* (nconc (zff v (car fl) (cadr fl)) fcs*))))
 566            (t (setq fl (cdr fl))
 567               (setq d (zff v (caar fl) (caadr fl)))
 568               (setq v nil)
 569               (zfactsplit (car fl) (car d))
 570               (return (zfactsplit (cadr fl) (cadr d)))))))
 571
 572 (defun split2 (l)
 573   (prog (s)
 574      (setq s (nthcdr (1- (ash (length l) -1)) l))
 575      (setq dn* (copy-list (cdr s)))
 576      (rplacd s nil)
 577      (setq nn* l)))
 578
 579 (defun fsplit (l ind)
 580   (prog (nn* dn*)
 581      (cond ((null (cdr l)) (return l))
 582            ((null (cddr l))
 583             (return (list (apply #'ptimes l) l))))
 584      (split2 l)
 585      (setq nn* (fsplit nn* nil))
 586      (setq dn* (fsplit dn* nil))
 587      (return (list (cond (ind ind) (t (ptimes (car nn*) (car dn*))))
 588                    nn*
 589                    dn*))))
 590
 591 ;; this page contains routines changed for non-monic hack
 592
 593 (defun pexptmod (p n q)
 594   (prog (u x)
 595      (when (pcoefp p) (return (cexpt p n)))
 596      (setq q (cdr q) x (car p))
 597      (cond   ((oddp n)
 598               (setq p (setq u (pgcd1 (cdr p) q)))
 599               (go b))
 600              (t (setq u '(0 1))))
 601      (setq p (cdr p))
 602      a    (setq p (pgcd1 p q))
 603      b    (setq n (ash n -1))
 604      (when (zerop n) (return (cons x u)))
 605      (setq p (ptimes1 p p))
 606      (when (oddp n) (setq u (pgcd1 (ptimes1 u p) q)))
 607      (go a)))
 608
 609 (defun sqfrp (u var)
 610   (cond ((and (equal 0 (ptterm (cdr u) 0)) (equal 0 (ptterm (cdr u) 1)))
 611          nil)
 612         ((onevarp u)
 613          (setq u (pgcd u (pderivative u var)))
 614          (or (numberp u) (alg u)))
 615         (t (quick-sqfr-check u var))))
 616
 617 (declare-top (special p))
 618
 619 (defun fixvl0 (l1 l2 ov)
 620   (prog (a b c)
 621    loop (cond ((null ov) (setq subvar a subval b valist c) (return nil))
 622               ((member (car ov) l1 :test #'eq)
 623                (setq a (cons (car ov) a)
 624                      b (cons (assso (car ov) l1 l2) b)
 625                      c (cons (car b) c)))
 626               (t (setq c (cons 0 c))))
 627    (setq ov (cdr ov))
 628    (go loop)))
 629
 630 (defun assso (a l1 l2)
 631   (prog ()
 632    loop (cond ((null l1) (return nil))
 633               ((eq (car l1) a) (return (car l2))))
 634    (setq l1 (cdr l1) l2 (cdr l2))
 635    (go loop)))
 636
 637 (defun zerohk (l)
 638   (prog (ans i)
 639      (when (null l) (return nil))
 640      ag   (setq ans (car l)
 641                 i (count 0 ans))
 642      loop (setq l (cdr l))
 643      (cond ((null l) (return ans))
 644            ((> (count 0 (car l)) i) (go ag)))
 645      (go loop)))
 646
 647 (defun multfact (poly)
 648   (prog (*inl3 *i* *min* *mx* nn* *odr* lc elm listelm plim origenvar ne var valist val1
 649          ovarlist p subvar subvar1 subval1 subval dlp)
 650      ;;      (declare (special p))
 651      (setq var (car poly) elm (listovars poly)
 652            origenvar genvar
 653            genvar (intersect genvar (if algfac*
 654                                         (delete (car alpha) elm :test #'equal)
 655                                         elm))
 656            ovarlist (butlast genvar)    ;this depends on the order of the above intersection!
 657            nn* (1+ (length ovarlist)))
 658      (setq listelm 0)
 659      (setq lc (caddr poly))
 660      (setq elm 1 *i* 1 ne 1)
 661      (setq subval (reverse poly))
 662      (setq *odr*(putodr (reverse ovarlist)))
 663      (setq val1 (zerohk (nconc (degvector nil 1 lc)
 664                                (cond ((or (> (cadr subval) 0)
 665                                           (> (cadddr subval) 1))
 666                                       (degvector nil 1 (car subval)))))))
 667      (setq subval nil)
 668      (setq p poly)
 669      (when (null val1)
 670        (setq subvar1 ovarlist)
 671        (setq subval1 (polysubst (make-list (length subvar1) :initial-element 0) subvar1))
 672        (go tag))
 673      (fixvl val1 ovarlist)
 674      (fixvl1 val1 ovarlist)
 675      (when subval1 (setq subval1 (polysubst subval1 subvar1)))
 676      (setq subval (polysubst (completevector nil 0 (length subval) 1) subvar))
 677      tag  (fixvl subval1 subvar1)
 678      (setq subval1 nil subvar1 nil)
 679      (fixvl0 subvar subval (reverse ovarlist))
 680      (when algfac* (push (car alpha) genvar))
 681      (setq poly (cpber3 poly p))
 682      (setq genvar origenvar)
 683      (return poly)))
 684
 685 (defun polysubst (a b)
 686   (prog (lc *inl3 d n modulus)
 687      (when modulu* (setq modulus modulu*))
 688      (setq *inl3 t lc (caddr p) n (length a))
 689      loop (setq d (pcsubsty a b lc))
 690      (when (equal 0 d) (go inl))
 691      (let ((modulus nil))
 692        (setq d (pcsubsty a b p)))
 693      (when (sqfrp (pmod d) (car d))
 694        (setq p d) (return a))
 695      inl  (setq a (increaselist a n))
 696      (go loop)))
 697
 698 (declare-top (unspecial p))
 699
 700 (defun fixvl1 (l r)
 701   (prog nil
 702    loop (cond ((null l)
 703                (setq subval1 (nreverse subval1) subvar1 (nreverse subvar1))
 704                (return nil))
 705               ((zerop (car l))
 706                (setq subval1 (cons (car l) subval1))
 707                (setq subvar1 (cons (car r) subvar1))))
 708    (setq l (cdr l))
 709    (setq r (cdr r))
 710    (go loop)))
 711
 712 (defun fixvl (l r)
 713   (prog nil
 714    loop (cond ((null l)
 715                (setq subval (nreverse subval) subvar (nreverse subvar))
 716                (return nil))
 717               ((not (zerop (car l)))
 718                (setq subval (cons (car l) subval))
 719                (setq subvar (cons (car r) subvar))))
 720    (setq l (cdr l))
 721    (setq r (cdr r))
 722    (go loop)))
 723
 724 (defun maxcoef (p)
 725   (maxcoefficient p))
 726
 727 (defun incrlimk (p)
 728   (prog (v min-plim)
 729      (cond (modulu* (setq plim modulu* *prime modulu* limk -1) (return nil))
 730            ((null limk)(setq plim *alpha *prime *alpha limk -1)(return nil)))
 731      (setq v (butlast (pdegreevector p)))
 732      (setq v
 733            (apply
 734             '*
 735             (mapcar
 736              #'(lambda (a b)
 737                  (cond ((equal b 0) 1)
 738                        (t (max (* (simpbinocoef (list '(%binocoef)
 739                                                       a
 740                                                       (ash a -1))
 741                                                 1
 742                                                 t)
 743                                   (expt b (ash a -1)))
 744                                (expt b a)))))
 745              v
 746              valist)))
 747      (setq min-plim (* (max (maxcoef p) plim) v))
 748      loop (cond ((< min-plim plim) (return nil)))
 749      (incf limk)
 750      (setq plim (* plim plim))
 751      (go loop)))
 752
 753
 754
 755 (defun increaselist (l n)
 756   (cond (*inl3 (setq l (inlist3 l))))
 757   (cond (*inl3 l)
 758         (t (cond ((equal elm 2)
 759                   (cond (modulu*
 760                          (merror (intl:gettext "factor: not enough choices for substitution.")))
 761                         (t (rand n 13.))))
 762                  ((equal ne n)
 763                   (incf elm)
 764                   (setq ne 1)
 765                   (completevector (baselist ne) ne n listelm))
 766                  (t (cond ((equal *i* n)
 767                            (incf ne)
 768                            (completevector (baselist ne) ne n listelm))
 769                           (t (incf *i*)
 770                              (butlast (cons listelm l)))))))))
 771
 772
 773 ;; Returns a list of N random numbers.  If MODULUS is set, then the
 774 ;; numbers will be modulo MODULUS.   Otherwise, between 0 and 1000.
 775 (defun rand (n modulus)
 776   (declare (fixnum n))
 777   (do ((i n (1- i)) (l))
 778       ((= i 0) (cond (modulus (mapcar #'cmod l))
 779                      (t l)))
 780     (declare (fixnum i))
 781     (push (random 1000.) l)))
 782
 783 (defun trufac (v lp olfact many* modulus)
 784   (prog (ans olc lc af qnt factor lfunct)
 785      (setq lc 1 olc 1)
 786      (set-modulus modulus)
 787      (setq lfunct (setq olfact (cons nil olfact)))
 788      test (cond
 789             ((equal v 1) (setq ans factor) (go out))
 790             ((null lp)
 791              (setq
 792               ans
 793               (cond ((< (length olfact) 4) (cons v factor))
 794                     (t (nconc factor
 795                               (nprod lc
 796                                      v
 797                                      (cons (let ((modulus plim))
 798                                              (ptimes olc (cadr olfact)))
 799                                            (cddr olfact)))))))
 800              (go out))
 801             ((and (null (cdr lp)) (or (null (cdr olfact)) (null (cddr olfact))))
 802              (setq ans (cons v factor))
 803              (go out)))
 804      (setq af (car lp))
 805      (cond ((setq qnt (let ((modulus modulu*))
 806                         (testdivide v af)))
 807             (setq factor (cons af factor))
 808             (setq lc (ptimes lc (caddr af)))
 809             (setq v qnt)
 810             (let ((modulus plim))
 811               (setq olc (ptimes (caddr (cadr lfunct)) olc)))
 812             (rplacd lfunct (cddr lfunct)))
 813            (t (setq lfunct (cdr lfunct))))
 814      (setq lp (cdr lp))
 815      (go test)
 816      out  (return ans)))
 817
 818 (defun multideg (p)
 819   (prog (m d)
 820      (cond ((numberp p) (return 0)) ((onevarp p) (return (cadr p))))
 821      (setq p (cdr p) m (car p))
 822      loop (cond ((null p) (return m)))
 823      (setq d (+ (car p) (multideg (cadr p))) p (cddr p) m (max d m))
 824      (go loop)))
 825
 826 (defun oddelm (list)
 827   "return a list of the first. third etc. elements of list"
 828   (loop for el in list by #'cddr
 829         collect el))
 830
 831 (defun cpber3 (v u)
 832   (prog (factz alcinv lc plim monic* sharpcont limk var vfact)
 833      (setq var (car u))
 834      (cond ((and algfac* (not (atom (caddr u))))
 835             (setq alc (caddr u))
 836             (setq u (ptimes u (car(setq alcinv(rainv alc))) ))
 837             (setq v (ptimes v (car alcinv)))
 838             (setq adn* (* adn* (cdr alcinv)))))
 839      (setq u (oldcontent u))
 840      (setq sharpcont (car u) u (cadr u))
 841      (setq lc (caddr v))
 842      (cond ((equal lc 1) (setq monic* t)))
 843      (setq factz (fact5 u))
 844      ;; this is the barry trick
 845      (cond (*stop* (setq *stop* plim) (return (cons (car subval)  factz))))
 846      (setq u nil)
 847      (cond ((null (cdr factz)) (return (list v)))
 848            ((and algfac* (not (equal adn* 1)))
 849             (setq v (pctimes adn* v) lc (pctimes adn* lc))))
 850      (incrlimk v)
 851      (setq modulus plim)
 852      (setq u v v (newrep v))
 853      (cond ((numberp (car factz))
 854             (setq sharpcont (ptimes sharpcont (car factz)) factz (cdr factz))))
 855      (cond ((not (equal sharpcont 1))
 856             (setq factz (cons (ptimes sharpcont (car factz)) (cdr factz)))))
 857      (setq vfact (zfact v factz limk t))
 858
 859      (setq factz (cond (monic* (reverse vfact))
 860                        (t (restorelc vfact (newrep lc)))))
 861      (cond ((and algfac* (not (equal adn* 1)))
 862             (setq v (pctimes (crecip adn*) v))(setq adn* 1)))
 863      (setq vfact (trufac v factz (nreverse vfact) t modulu*))
 864      (setq factz nil)
 865      (cond ((null (cdr vfact)) (return (list u)))
 866            (t (return (mapcar #'oldrep vfact))))))
 867
 868
 869
 870 (defun nprod (lc u lfunct)
 871   (prog (stage v d2 af0 r lcindex factor llc ltuple lprod lindex qnt af
 872          funct tuple ltemp lpr f l li lf modulus)
 873      (setq lpr (copy-tree (setq ltemp (cons nil nil))))
 874      (setq lprod (cons nil lfunct))
 875      (setq d2 (ash (cadr u) -1))
 876      (remov0 lprod d2)
 877      (setq lfunct (cdr lprod))
 878      (setq lindex (index* (setq r (length lfunct))))
 879      (cond ((not monic*)
 880             (setq llc (mapcar #'caddr lfunct))
 881             (setq lcindex (copy-list lindex))
 882             (remov3 llc lcindex)
 883             (setq v (ptimes lc (ptimes (caddr u) u))))
 884            (t (setq v u)))
 885      (setq ltuple (cons nil (mapcar  #'list lindex)))
 886      (setq stage 1)
 887      (setq lindex (cons nil lindex))
 888      (setq lfunct (copy-list lprod))
 889      tloop (incf stage)
 890      nextuple (cond ((or (> stage d2) (> stage (1- r))
 891                          (null ltuple) (null (cdr ltuple)))
 892                      (return (cons u factor))))
 893      (setq li (cdr lindex))
 894      (setq lf (cdr lfunct))
 895      (setq tuple (cadr ltuple))
 896      (setq funct (cadr lprod))
 897      (rplacd ltuple (cddr ltuple))
 898      (rplacd lprod (cddr lprod))
 899      iloop(setq l (car li))
 900      (setq f (car lf))
 901      (setq li (cdr li))
 902      (setq lf (cdr lf))
 903      (cond ((and (not (member l tuple :test #'equal))
 904                  (not (> (+ (cadr f) (cadr funct)) d2))
 905                  (not (member (setq l (orde l tuple)) ltemp :test #'equal)))
 906             (set-modulus plim)
 907             (setq af0 (setq af (ptimes(pmod f) (pmod funct))))
 908             (cond (llc (setq af (ptimes (pmod (lchk llc lcindex l)) af))))
 909             (cond (many* (setq af (dropterms af)))
 910                   ((and algfac* (not (equal intbs* 1)))(setq af (intbasehk af))))
 911             (set-modulus nil)
 912             (cond ((setq qnt (testdivide v af))
 913                    (cond (llc (setq af (oldcontent af))
 914                               (setq v (ptimes (car af) qnt)af (cadr af))
 915                               (setq u (cond (algfac*(car (ignore-rat-err
 916                                                            (rquotient u af))))
 917                                             (t (pquotient u af)))))
 918                          (t (setq u qnt v qnt)))
 919                    (setq factor (cons af factor))
 920                    (cond ((equal u 1) (return factor)))
 921                    (setq d2 (ash (cadr u) -1))
 922                    (cond ((< d2 stage) (return (cons u factor))))
 923                    (remov1 l ltuple lprod d2)
 924                    (remov1 l ltemp lpr d2)
 925                    (remov2 l lindex lfunct d2)
 926                    (setq r (- r stage))
 927                    (setq li nil))       ; exit iloop
 928                   (t (setq ltemp (nconc ltemp (list l)))
 929                      (setq lpr (nconc lpr (list af0)))))))
 930      (cond (li (go iloop)) ((cdr ltuple) (go nextuple)))
 931      (setq ltuple ltemp lprod lpr ltemp nil lpr nil)
 932      (go tloop)))
 933
 934 (defun remov2 (a b c d2)
 935   (prog nil
 936    tag1 (cond ((null (cdr b)) (return nil))
 937               ((or (member (cadr b) a :test #'equal) (> (cadadr c) d2))
 938                (rplacd b (cddr b))
 939                (rplacd c (cddr c))
 940                (go tag1)))
 941    (setq b (cdr b))
 942    (setq c (cdr c))
 943    (go tag1)))
 944
 945 (defun remov1 (a lt1 lp1 d2)
 946   (prog nil
 947    tag1 (cond ((null (cdr lt1)) (return nil))
 948               ((and (not (> (cadadr lp1) d2))
 949                     (null (intersection a (cadr lt1) :test #'equal)))
 950                (setq lt1 (cdr lt1))
 951                (setq lp1 (cdr lp1))
 952                (go tag1)))
 953    (rplacd lt1 (cddr lt1))
 954    (rplacd lp1 (cddr lp1))
 955    (go tag1)))
 956
 957 (defun remov0 (lf d2)
 958   (prog (d)(setq d lf)
 959         tag  (cond ((null (cdr lf)) (return nil))
 960                    ((> (cadadr lf) d2)(setq d2 (caddr (cadr lf))) (rplacd lf (cddr lf))
 961                     (cond ((equal d2 1) nil)(t (rplacd d (cons (ptimes d2 (cadr d)) (cddr d)))))
 962                     (return nil)))
 963         (setq lf (cdr lf))
 964         (go tag)))
 965
 966 (defun remov3 (a b)
 967   (prog nil
 968    loop      (cond ((null (cdr a)) (return nil))
 969                    ((equal (cadr a) 1)
 970                     (rplacd a (cddr a))
 971                     (rplacd b (cddr b))(go loop)))
 972    (setq a (cdr a) b (cdr b))(go loop)))
 973
 974 (defun lchk (a b c)
 975   (prog (ans)
 976      (setq ans 1)
 977      loop (cond ((null a) (return ans))
 978                 ((not (member (car b) c :test #'equal)) (setq ans (ptimes ans (car a)))))
 979      (setq a (cdr a) b (cdr b))
 980      (go loop)))
 981
 982
 983
 984 (defun lcprodl (l)
 985   (prog (ans d)
 986      (setq d 1 l (reverse l) ans '(1))
 987      loop (cond ((null (cdr l)) (return ans)))
 988      (setq d (ptimes d (caddar l)))
 989      (setq l (cdr l))
 990      (setq ans (cons d ans))
 991      (go loop)))
 992
 993
 994 (defun fact5 (poly)
 995   (prog (ql trl* linfac uu* lc deg factp factz modulus monic*  split* var
 996          anotype fctc invc *afixn* *fctcfixn* *invcfixn*)
 997      (setq var (car poly))
 998      (cond ((null (cdddr poly)) (return (list poly))))
 999      (cond((and algfac* (not (atom (caddr poly))))
1000            (setq alc (caddr poly))
1001            (setq poly (rattimes (cons poly 1) (setq alcinv(rainv alc)) t))
1002            (setq adn*(* adn* (cdr poly)))
1003            (setq poly (car poly))))
1004      (cond((and algfac* minpoly* (or $nalgfac (equal (cdr minpoly*) '(4 1 0 1))))
1005            (setq ql 'splitcase) (go tag0)))
1006      (setq uu* poly)
1007      (cond ((equal (setq lc (caddr uu*)) 1) (setq monic* t)))
1008      (setq deg (cadr poly))
1009      (cond ((not algfac*)
1010             (setq *fctcfixn* (make-array deg :initial-element 0)
1011                   *invcfixn* (make-array deg :initial-element 0)
1012                   *afixn* (make-array (list deg deg) :initial-element 0)))
1013            (t (setq fctc (make-array deg)
1014                     invc (make-array deg)
1015                     anotype (make-array (list deg deg))
1016                     *fctcfixn* (make-array mm* :initial-element 0)
1017                     *invcfixn* (make-array mm* :initial-element 0)
1018                     *afixn* (make-array (list mm* mm*) :initial-element 0))))
1019      (cond (modulu* (return (fact5mod poly))))
1020      (cond ((not (atom (setq ql (choozp uu*))))
1021             (setq linfac (car ql) uu* (caddr ql) ql (cadr ql))))
1022      (setq *prime modulus)
1023      tag0            (cond ((eq ql 'splitcase)
1024                             (setq poly(nalgfac poly (cons (car alpha) (cdr minpoly*))))
1025                             (setq plim *alpha *prime plim limk -1)
1026                             (return poly))
1027                            ((null (cdr (append linfac ql)))
1028                             (setq poly (list poly))
1029                             (go out))
1030                            ((equal uu* 1) (setq factp nil) (go on)))
1031      (cond  (algfac* (setq factp (cpbgzass ql  uu* (length ql))))
1032             ((not (equal uu* 1))
1033              (setq factp (cpbg ql  uu* (length ql)))))
1034      (setq uu* nil)
1035      on   (setq factp (nconc factp linfac)
1036                 linfac nil
1037                 factp (cons (pctimes (pmod lc) (car factp)) (cdr factp)))
1038      (setq limk (klim poly modulus))
1039      (setq factz (zfact poly factp limk nil)factp nil)
1040      (setq poly (trufac poly
1041                         (let ((modulus plim))
1042                           (restorelc factz lc))
1043                         (nreverse factz)
1044                         nil
1045                         nil))
1046      (setq modulus nil)
1047      out  (return poly)))
1048
1049 (defun fact5mod (u)
1050   (prog (lc poly)
1051      (setq poly (copy-list u))
1052      (set-modulus modulu*)
1053      (setq poly (pmod poly))
1054      (setq lc (caddr poly))
1055      (pmonicize (cdr poly))
1056      (setq poly (cpberl poly))
1057      (cond ((null (cdr poly))
1058             (return (list u)))
1059            (t (return (if (= lc 1)
1060                           poly
1061                           (cons lc poly)))))))
1062
1063 (defun cpbg (qlist v m)
1064   (declare (fixnum m))
1065   (prog (y vj factors u w (j 0)
1066          (p1  (ash modulus -1))
1067          (p2 1)
1068          fnj fnq oldfac)
1069      (declare (fixnum j p1 p2))
1070      (when (= m 1) (return (list v)))
1071      (setq p1 (ash modulus -1))
1072      (setq p2 1)
1073      (setq qlist (cdr (nreverse qlist)))
1074      (setq oldfac (list nil v))
1075      (setq v nil)
1076      tag3 (setq vj (nconc (car qlist) (list 0 0)))
1077      (setq qlist (cdr qlist))
1078      (setq j (- p1))
1079      (setq oldfac (nconc oldfac fnq))
1080      (setq fnq nil)
1081      incrj(setq factors (nconc oldfac fnj))
1082      (setq fnj nil)
1083      (pcdifconc vj j)
1084      tag2 (setq u (cadr factors))
1085      (setq w (pgcdu vj u))
1086      (cond ((or (numberp w) (= (cadr w) (cadr u))) (go agg)))
1087      (setq y (car (pmodquo u w)))
1088      (setq fnq (cons (copy-list w) fnq))
1089      (setq fnj (cons y fnj))
1090      (incf p2)
1091      (rplacd factors (cddr factors))
1092      (cond ((equal p2 m) (go out)) (t (go tag1)))
1093      agg  (setq factors (cdr factors))
1094      tag1 (cond ((cdr factors) (go tag2))
1095                 ((< j p1) (incf j) (go incrj))
1096                 (qlist (go tag3)))
1097      out  (return (nconc fnq fnj (cdr oldfac)))))
1098
1099
1100
1101
1102 (defun fact2z (u f g limk)
1103   (prog (a a1 w pk mpk b c r p ql qlp h (k 0) b1)
1104      (declare (fixnum k))
1105      (setq p modulus)
1106      (setq r (ppprog f g))
1107      (setq a (car r))
1108      (setq b (cadr r))
1109      (let ((modulus nil))
1110        (setq r (pdifference (ptimes f g) u)))
1111      sharp (cond ((or (equal r 0) (> k limk)) (go on)))
1112      (setq pk modulus mpk (- pk))
1113      (setq modulus (* modulus modulus))
1114      (setq w (pmod r))
1115      (cond ((equal w 0) (go tag1)))
1116      (setq c (npquo w pk))(setq w nil)
1117      (setq ql (pmodquo (ptimes a c) g))
1118      (setq a1 (npctimes mpk
1119                         (pplus (ptimes (car ql) f)
1120                                (ptimes b c))))
1121      (setq b1 (npctimes mpk (cdr ql)))
1122      (let ((modulus plim))
1123        (setq r (pplus (pplus r (ptimes a1 b1))
1124                       (pplus (ptimes a1 g) (ptimes b1 f))))
1125        (setq f (pplus f a1))
1126        (setq g (pplus g b1)))
1127      (setq a1 nil b1 nil)
1128      tag1         (cond ((or (equal r 0)(> (incf k) limk)) (go on)))
1129      (setq h (npquo (pplus (pplus (ptimes a f)
1130                                   (ptimes b g))
1131                            -1)
1132                     pk))
1133      (setq qlp (pmodquo (ptimes a h) g))
1134      (setq b1 (pplus (ptimes b h) (ptimes (car qlp) f)))
1135      (setq a (pplus a (npctimes mpk (cdr qlp))))
1136      (setq b (pplus b (npctimes mpk b1)))
1137      (setq h nil b1 nil qlp nil)
1138      (go sharp)
1139      on   (set-modulus p)
1140      (return (list f g))))
1141
1142
1143
1144 (defun npctimes (c p)
1145   (setq p (npctimes1 c p))
1146   (if (and (not (atom p)) (null (cdr p)))
1147       0
1148       p))
1149
1150 (defun npquo (p c)
1151   (prog (u modulus)
1152      (cond ((equal c 1) (return p))
1153            ((pcoefp p) (return (quotient p c))))
1154      (setq u p)
1155      loop (cond ((null (cdr u)) (return p)))
1156      (setq u (cddr u))
1157      (rplaca u (cond ((pcoefp (car u))
1158                       (quotient (car u) c))
1159                      (t (npquo (copy-list (car u)) c))))
1160      (go loop)))
1161
1162 (defun npctimes1 (c p)
1163   (prog (u a)
1164      (cond((equal c 1)(return p))
1165           ((pcoefp p)(return (ctimes c p))))
1166      (setq u p)
1167      loop (cond ((null (cdr u))(return p)))
1168      (setq a (cond ((pcoefp (caddr u)) (ctimes c (caddr u)))
1169                    (t (npctimes c (copy-list (caddr u))))))
1170      (cond ((equal a 0) (rplacd u (cdddr u)))
1171            (t (setq u (cddr u))
1172               (rplaca u a)))
1173      (go loop)))
1174
1175 (defun x**q1 (term u m p)
1176   (declare (fixnum m))
1177   (prog ((i 1))
1178      (declare (fixnum i))
1179      (setq trl* (list term))
1180      loop (when (= i m) (return (pexptmod term p u)))
1181      (setq term (pexptmod term p u))
1182      (setq trl* (cons term trl*))
1183      (incf i)
1184      (go loop)))
1185
1186 (defun cptomf (p u n)
1187   (declare (fixnum n p))
1188   (prog (l s *xn (j 0) (i 0) ind (n-1 (1- n)) )
1189      (declare (fixnum i j))
1190      loop (incf j)
1191      (cond ((= j n) (return nil))
1192            (ind (go sa))
1193            ((> (* p j) n-1)
1194             (setq *xn (mapcar #'- (p2cpol (cddr u) n-1))
1195                   s (copy-tree *xn)
1196                   ind t)
1197             (setq i (- (* p j) n))
1198             (go sa1)))
1199      (setq s (p2cpol (list var (* p j) 1) n-1))
1200      (go st)
1201      sa   (setq i p)
1202      sa1  (cond ((= i 0) (go st)))
1203      (cptimesx s)
1204      (decf i)
1205      (go sa1)
1206      st (cond ((and (= j 1)
1207                     (equal '(1 0) (last s 2))
1208                     (= 1 (apply #'+ s)))
1209                (return (setq split* t))))
1210      (setq l s)
1211      (setq i n-1)
1212      sharp2  (when (null l) (go on))
1213      (setf (aref *afixn* j i) (car l))
1214      (setq l (cdr l))
1215      (decf i)
1216      (go sharp2)
1217      on   (decf (aref *afixn* j j))
1218      (go loop)))
1219
1220 (defun p2cpol (p n)
1221   (declare (fixnum n))
1222   (prog (l)
1223      (setq p (cdr p))
1224      loop (cond ((= n -1) (return (nreverse l)))
1225                 ((or (null p) (> n (car p))) (setq l (cons 0 l)))
1226                 ((= n (car p))
1227                  (setq l (cons (cadr p) l))
1228                  (setq p (cddr p))))
1229      (decf n)
1230      (go loop)))
1231
1232 (defun cptimesx (p)
1233   (prog (xn q lc)
1234      (setq xn *xn q p lc (car p))
1235      loop (cond ((cdr q)
1236                  (rplaca q (cplus (cadr q) (ctimes lc (car xn))))
1237                  (setq q (cdr q) xn (cdr xn)))
1238                 (t (rplaca q (ctimes lc (car xn))) (return p)))
1239      (go loop)))
1240
1241
1242 (defun cmnullf (n)
1243   (declare (fixnum n))
1244   (prog (nullsp mone (k 1) (j 0) s (n-1 (1- n)) nullv vj m aks)
1245      (declare (fixnum k j n-1))
1246      (setq mone (cmod -1))
1247      (do ((i 0 (1+ i))) ((> i n-1))
1248        (setf (aref *fctcfixn* i) -1)
1249        (setf (aref *invcfixn* i) -1))
1250      (setq  nullsp (list 1))
1251      n2   (when (> k n-1) (return nullsp))
1252      (setq j 0)
1253      n3a  (cond ((> j n-1) (go null))
1254                 ((or (= (aref *afixn* k j) 0) (> (aref *fctcfixn* j) -1))
1255                  (incf j)
1256                  (go n3a)))
1257      (setf (aref *invcfixn* k) j)
1258      (setf (aref *fctcfixn* j) k)
1259      (setq m (aref *afixn* k j))
1260      (setq m (crecip (ctimes mone m)))
1261      (do ((s k (1+ s))) ((> s n-1))
1262        (setf (aref *afixn* s j) (ctimes m (aref *afixn* s j))))
1263      ;; go through columns
1264      (setq s 0)
1265      sharp2  (when (> s n-1) (go nextk))
1266      ;; go through rows in each column
1267      (cond ((= s j) nil)
1268            (t (setq aks (aref *afixn* k s))
1269               (do ((i k (1+ i)))
1270                   ((> i n-1))
1271                 (setf (aref *afixn* i s)
1272                        (cplus (aref *afixn* i s)
1273                               (ctimes (aref *afixn* i j) aks))))))
1274      (incf s)
1275      (go sharp2)
1276      null (setq nullv nil)
1277      (do ((s 0 (1+ s))) ((> s n-1))
1278        (cond ((= s k) (setq nullv (cons s (cons 1 nullv))))
1279              ((> (aref *invcfixn* s) -1)
1280               (setq vj (aref *afixn* k (aref *invcfixn* s)))
1281               (cond ((= vj 0) nil)
1282                     (t (setq nullv (cons s (cons vj nullv))))))))
1283      (cond ((equal (car nullv) 0) (setq nullv (cadr nullv)))
1284            ((setq nullv (cons var nullv))))
1285      (setq nullsp (cons nullv nullsp))
1286      nextk (incf k)
1287      (go n2)))
1288
1289
1290 (defun choozp (v)
1291   (prog (lchar1 u tr n (ncont 1) bmod b1 b mincont (lmin 0) (nf 0)
1292          (deg  (cadr v)) (algcont 0))
1293      (declare (special ncont lmin nf deg algcont))
1294      (setq nf (integer-length deg))
1295      (setq lchar1 (if gauss '(3 7 11. 19. 23. 29. 31. 37.) (cdr *small-primes*)))
1296      test (setq modulus (car lchar1))
1297      (setq u (pmod v))
1298      (cond ((or (zerop (rem sharpcont modulus))
1299                 (and (not monic*)
1300                      (or (pcoefp u)
1301                          (> deg (cadr u)))))
1302             (go nextp)))
1303      (cond ((or (null (sqfrp u var))
1304                 (and algfac*
1305                      (not gauss)
1306                      (not (iredup (pmod minpoly*)))))
1307             (incf algcont)
1308             (go nextp)))
1309      (pmonicize (cdr u))
1310      (setq b1 (catch 'splt (cpber1 u)))(setq algcont 0)
1311      (incf ncont)
1312      (setq n (+ (length (car b1)) (length (cadr b1))))
1313      (cond ((or (zerop lmin) (< n lmin))
1314             (setq lmin n  mincont 1 bmod modulus b b1)
1315             (cond (algfac* (setq tr trl*))))
1316            ((= n lmin) (incf mincont)))
1317      (cond ((or (> ncont nf) (not(> n nf)) (= mincont 3)) (go out)))
1318      nextp (setq lchar1 (cdr lchar1))
1319      (cond ((null lchar1)
1320             (cond ((not (zerop lmin)) (go out))
1321                   (t (merror (intl:gettext "factor: ran out of primes.")))))
1322            ((and algfac* minpoly* (> algcont 6))
1323             (cond ((ziredup minpoly*)(setq trl* tr)(setq modulus nil)
1324                    (return 'splitcase))
1325                   (t (merror (intl:gettext "factor: the minimal polynomial must be irreducible over the integers."))))))
1326      (go test)
1327      out  (setq modulus bmod trl* tr)
1328      (return b)))
1329
1330
1331 (defun cpbq1 (a n)
1332   (declare (fixnum  n ))
1333   (prog ()
1334      (setq split* nil)
1335      (cond ((not (integerp modulus)) (*array 'a t n n)))
1336      (cond ((or algfac* (not (integerp modulus)))
1337             (cptom modulus mm* a n))
1338            (t (cptomf modulus a n)))
1339      (cond (split*
1340             (return (powrs (car a) (cadr a)))))
1341      (return (cond ((or algfac* (not (integerp modulus)))
1342                     (cmnull n))
1343                    (t (cmnullf n))))))
1344
1345
1346 (defun cpber1 (u)
1347   (prog (linfac)
1348      (setq var (car u))
1349      (setq linfac
1350            (linout u)
1351            u
1352            (car linfac)
1353            linfac
1354            (cadr linfac))
1355      (cond ((equal u 1) (return (list linfac nil u))))
1356      (return (list linfac (cpbq1 u (cadr u)) u))))
1357
1358
1359 (defun factor1972 (p)
1360   (let ((modulu* modulus) many* *stop* modulus mcflag *negflag*)
1361     (if (or (atom p) (numberp p) (and algfac* (alg p)))
1362         (list p)
1363         (factor72 p))))
1364
1365 (defun factor72 (p)
1366   (let ((sharpcont 1) plim)
1367     (setq p (cond ((onevarp p)
1368                    (mapcar #'posize (fact5 p)))
1369                   (t
1370                    (setq many* t)
1371                    (multfact p))))
1372     (if *negflag*
1373         (cons (pminus (car p)) (cdr p))
1374         p)))
1375
1376 (defun posize (p)
1377   (cond ((pminusp p)
1378          (setq *negflag* (not *negflag*))
1379          (pminus p))
1380         (t p)))