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