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 (defun testdivide (x y)
 215   (if algfac*
 216       (algtestd x y)
 217       (eztestdivide x y)))
 218
 219 (defun algtestd (x y)
 220   (and (div-deg-chk (nreverse (pdegreevector x)) (nreverse (pdegreevector y))
 221                     (reverse genvar))
 222        (cond ((setq x (ignore-rat-err (rquotient x y)))
 223               (setq adn* (* adn* (cdr x)))
 224               (car x)) )))
 225
 226 (defun div-deg-chk (xl yl gl)
 227   (cond ((or (null gl) (algv (car gl))) t)
 228         ((> (car yl) (car xl)) nil)
 229         (t (div-deg-chk (cdr xl) (cdr yl) (cdr gl)))))
 230
 231 ;; FUU is used by systems programmers such as BMT and PAULW while debugging.
 232 (defun fuu nil
 233   (setq tellratlist nil varlist nil genvar nil genpairs nil))
 234
 235 (defun linout (u)
 236   (prog (m linfac x y)
 237      (setq y (list (setq x (car u)) 1 1) m modulus)
 238      loop (decf m)
 239      (cond ((< m 0) (return (list u linfac)))
 240            ((equal (cadr u) 1) (return (list 1 (cons u linfac))))
 241            ((zerop (pcsubsty (cmod m) x u))
 242             (setq linfac
 243                   (cons (append y
 244                                 (cond ((zerop m) nil)
 245                                       (t (list 0 (cmod (- m))))))
 246                         linfac))
 247             (setq u (car (pmodquo u (car linfac))))))
 248      (go loop)))
 249
 250 (defun onevarp (p)
 251   (if algfac* (every #'pacoefp (cdr p))
 252       (every #'numberp (cdr p))))
 253
 254 (defun putodr (l)
 255   (do ((l l (cdr l))
 256        (i 1 (1+ i))
 257        (ans))
 258       ((null l) ans)
 259     (push (cons (car l) i) ans)))
 260
 261 (defun kterms (p k)
 262   (cond ((pacoefp p) p)
 263         ((= k 0) (consta p))
 264         (t (prog (v ans w)
 265               (setq v (car p))
 266               (setq p (cdr p))
 267               loop (cond ((null p) (return 0))
 268                          ((> (car p) k) (setq p (cddr p)) (go loop)))
 269               ag   (cond ((null p)
 270                           (return (psimp v ans))))
 271               (setq w (kterms (cadr p) (- k (car p))))
 272               (cond ((not (pzerop w))
 273                      (setq ans (nconc ans (list (car p) w)))))
 274               (setq p (cddr p))
 275               (go ag)))))
 276
 277 (defun consta (p)
 278   (cond ((or (pcoefp p) (alg p)) p)
 279         (t (consta (ptterm (cdr p) 0)))))
 280
 281 (defun constacl (p)                     ;NO LONGER USED?
 282   (cond ((atom p)
 283          (cond ((equal p 1) (throw 'cnt 1))
 284                (t (list p))))
 285         ((every #'numberp (cdr p))
 286          (setq p (oddelm p))
 287          (cond ((member 1 p) (throw 'cnt 1))
 288                (t (cdr p))))
 289         (t (apply #'append (mapcar #'constacl (cdr (oddelm p)))))))
 290
 291 (defun z1 (poly fact1 fact2)
 292   (prog (res hsteps steps kterm a b c d *ab* m df1 df2 dlr step *sharpa *sharpb)
 293      (let ((modulus))
 294        (set-modulus *prime)
 295        (setq *sharpb (fact20 fact1 fact2 limk)))
 296      (setq *sharpa (car *sharpb))
 297      (setq *sharpb (cadr *sharpb))
 298      (setq *ab* (list (list 0 *sharpa *sharpb)))
 299      (setq steps dlp
 300            hsteps (ash steps -1))
 301      (setq res (pdifference (ptimes (pmod fact1) (pmod fact2)) (pmod poly)))
 302      (setq poly nil)
 303      (setq step 0)
 304      (setq df1 fact1)
 305      (setq df2 fact2)
 306      loop (cond ((equal res 0) (go out)))
 307      (incf step)
 308      (cond ((> step steps) (go out)))
 309      (cond ((eq (car res) var) (setq c (cdr res)))
 310            (t (setq c (list 0 res))))
 311      (setq a 0 b 0)
 312      nextm (cond ((null c) (z2 a b step hsteps) (go loop)))
 313      (setq m (car c) dlr (cadr c))
 314      (setq c (cddr c))
 315      (setq kterm (kterms dlr step) dlr nil)
 316      (cond ((equal 0 kterm) (go nextm)))
 317      (setq d (obtainabm m))
 318      (setq b (pplus b (ptimes (car d) kterm))
 319            a (pplus a (ptimes (cadr d) kterm))
 320            kterm nil)
 321      (go nextm)
 322      out  (return (list df1 df2))))
 323
 324 (defun z2 (a b step hsteps)
 325   (unless (and (equal a 0) (equal b 0))
 326     (setq step
 327           (pdifference
 328            (pdifference (cond ((not (< step hsteps))
 329                                (dropterms (ptimes a b)))
 330                               (t (ptimes a b)))
 331                         (cond ((not (< step hsteps))
 332                                (dropterms (ptimes df1 b)))
 333                               (t (ptimes df1 b))))
 334            (cond ((not (< step hsteps))
 335                   (dropterms (ptimes df2 a)))
 336                  (t (ptimes df2 a)))))
 337     (setq res (pplus res step))
 338     (setq df1 (pdifference df1 a))
 339     (setq df2 (pdifference df2 b))))
 340
 341 (defun obtainabm (m)
 342   (prog (ans)
 343      (cond ((setq ans (cdr (assoc m *ab* :test #'equal))) (return ans)))
 344      (setq ans (obtainab (list var m 1)))
 345      (setq *ab* (cons (cons m ans) *ab*))
 346      (return ans)))
 347
 348 (defun fact20 (f1 g1 limk)
 349   (prog (f g a pk b reml qlp h k b1)
 350      (setq k 0)
 351      (setq reml (ppprog (pmod f1) (pmod g1)))
 352      (setq a (car reml))
 353      (setq b (cadr reml))
 354      sharp   (cond ((> k limk) (return (list a b))))
 355      (setq pk modulus)
 356      (set-modulus (* modulus modulus))
 357      (setq f(pmod f1) g (pmod g1))
 358      (setq h (pquo (pmod (pdifference (pplus (ptimes a f) (ptimes b g))
 359                                       1))
 360                    pk))
 361      (setq qlp (pmodquo (ptimes a h) g))
 362      (setq b1 (pplus (ptimes b h) (ptimes (car qlp) f)))
 363      (setq a (pdifference a (pmod (pctimes pk (cdr qlp)))))
 364      (setq b (pdifference b (pmod (pctimes pk b1))))
 365      (incf k)
 366      (go sharp)))
 367
 368 (defun baselist (n)
 369   (setq *i* n)
 370   (completevector nil 0 n elm))
 371
 372 (defun inlist3 (l)
 373   (cond ((null l) (setq *inl3 nil))
 374         ((zerop (car l)) (cons 1 (cdr l)))
 375         ((equal (car l) 1) (cons -1 (cdr l)))
 376         (t (cons 0 (inlist3 (cdr l))))))
 377
 378 (defun newrep (p)
 379   (let ((modulus))
 380     (if subvar (pcsubsty (mapcar #'(lambda (a b) (list a 1 1 0 b))
 381                                  subvar subval)
 382                          subvar
 383                          p)
 384         p)))
 385
 386 (defun oldrep (p)
 387   (let ((modulus))
 388     (if subvar (pcsubsty (mapcar #'(lambda (a b) (list a 1 1 0 (- b)))
 389                                  subvar subval)
 390                          subvar
 391                          p)
 392         p)))
 393
 394 (defun completevector (l n m v)
 395   (do ((i m (1- i)))
 396       ((= i n) l)
 397     (push v l)))
 398
 399 ;; special variable *odr* contains order of nesting of variables in c
 400 (defun degvector (l n c)
 401   (prog (lf ans j)
 402    bk (cond ((numberp c)
 403              (return  (list (completevector l n nn* 0)))))
 404    (setq j (cdr (assoc (car c) *odr* :test #'equal)))
 405              ;;; IN CASE (CAR C) IS ALGEBRAIC
 406    (cond ((null j) (setq c 0)(go bk)))
 407    (setq c (cdr c))
 408    (setq lf (completevector l n j 0))
 409    loop (cond ((null c) (return ans)))
 410    (setq ans
 411          (nconc (degvector (cons (car c) lf) (1+ j) (cadr c)) ans))
 412    (cond (*mx* (setq ans (ncons (maxlist ans))))
 413          (*min* (setq ans (ncons (minlist ans)))))
 414    (setq c (cddr c))
 415    (go loop)))
 416
 417 (defun union1 (a b)
 418   (do ((a a (cdr a))
 419        (ans b))
 420       ((null a) ans)
 421     (or (member (car a) ans :test #'equal)
 422         (setq ans (cons (car a) ans)))))
 423
 424 (defun obtainab (u)
 425   (prog (c ql)
 426      (setq c (pmod u))
 427      (setq ql (pmodquo (ptimes *sharpa c) fact2))
 428      (return (list (cdr ql) (pmod (pplus (ptimes (car ql) fact1)
 429                                          (ptimes *sharpb c)))))))
 430
 431 (defun pcdifconc (v j)
 432   (do ((l v (cddr l)))
 433       ((null (cdr l))
 434        (or (= j 0)
 435            (rplacd l (list 0 j)))
 436        v)
 437     (cond ((= (cadr l) 0)
 438            (cond ((= j 0)
 439                   (rplacd l nil))
 440                  ((rplaca (cddr l) j)))
 441            (return v)))))
 442
 443 (defun orde (a l)
 444   (cond ((null l) (list a))
 445         (t (cond ((< a (car l)) (cons a l))
 446                  (t (cons (car l) (orde a (cdr l))))))))
 447
 448 (defun pquo (x y)
 449   (let (modulus)
 450     (pquotient x y)))
 451
 452 ;; maintains order of list x
 453 (defun intersect (x y)
 454   (if x (if (member (car x) y :test #'equal)
 455             (cons (car x) (intersect (cdr x) y))
 456             (intersect (cdr x) y))))
 457
 458 ;; Like APL IOTA function.
 459 (defun index* (k)
 460   (declare (fixnum k))
 461   (if (< k 2) (list 1) (cons k (index* (f1- k)))))
 462
 463 ; sets plim to a power of p1 likely to be a large enough modulus for polynomial u
 464 ; returns limk, where plim = p1^(2^(limk+1))
 465 (defun klim (u p1)
 466   (prog (bcoef)
 467      (setq bcoef (* 5 (maxcoefficient u)))
 468      (when algfac* (setq bcoef (* bcoef intbs*)))
 469      (when (< bcoef 10000.) (setq bcoef 20000.))
 470      (setq limk 0)
 471      test (setq p1 (* p1 p1))
 472      (when (> p1 bcoef)
 473        (setq plim p1)
 474        (return limk))
 475      (incf limk)
 476      (go test)))
 477
 478 (declare-top (special b b2))
 479
 480 (defun cpberl (u)
 481   (prog (ql d)
 482      (setq ql (catch 'splt (cpber1 u)) u (caddr ql))
 483      (setq d (car ql) ql (cadr ql))
 484      (cond ((null ql)(return d))
 485            ((null (cdr ql)) (return (cons u d))))
 486      (return (append d
 487                      (cond ((or alpha (> modulus 70.))
 488                             (cpbgzass ql (pmod u) (length ql)))
 489                            (t (cpbg ql (pmod u) (length ql))))))))
 490
 491 ;; Returns a list of monomials in G of degree less than N.
 492 (defun powrs (g n &aux (ans (ncons 1)))
 493   (declare (fixnum n))
 494   (do ((i 1 (1+ i))) ((= i n) ans)
 495     (declare (fixnum i))
 496     (push (make-poly g i 1) ans)))
 497
 498 ;; Finds polynomials A and B such that A*F+B*G=1 when MODULUS
 499 ;; is non-NIL.  Same algorithm as INVMOD.
 500 (defun ppprog (f g)
 501   (prog (a1 a2 b1 b2 r1 r2 ql ans ap bp g1 f1 s)
 502      (cond ((> (cadr g) (cadr f)) (setq g1 g) (setq f1 f))
 503            (t (setq g1 f) (setq f1 g) (setq s t)))
 504      (setq ql (pmodquo g1 f1))
 505      (setq a1 1)
 506      (setq b1 0)
 507      (setq a2 (pminus (car ql)))
 508      (setq b2 1)
 509      (setq r1 f1)
 510      (setq r2 (cdr ql))
 511      test (cond ((or (numberp r2) (and alpha (alg r2))) (go end)))
 512      (setq ql (pmodquo r1 r2))
 513      (setq ap (pdifference a1 (ptimes (car ql) a2)))
 514      (setq bp (pdifference b1 (ptimes (car ql) b2)))
 515      (setq r1 r2)
 516      (setq r2 (cdr ql))
 517      (setq a1 a2)
 518      (setq b1 b2)
 519      (setq a2 ap)
 520      (setq b2 bp)
 521      (go test)
 522      end  (cond ((pzerop r2)
 523                  (cond ((equal 1 (setq ans (caddr r1)))
 524                         (setq ans (list b1 a1)))
 525                        (t (setq ans (list (car (pmodquo b1 ans))
 526                                           (car (pmodquo a1 ans))))))
 527                  (go out)))
 528      (setq ans (list (car (pmodquo b2 r2)) (car (pmodquo a2 r2))))
 529      out  (cond ((not s) (return (reverse ans))) (t (return ans)))))
 530
 531
 532 (defun zff (v f g)
 533   (if many*
 534       (z1 v f g)
 535       (fact2z v f g limk)))
 536
 537 (defun zfact (u fl limk many*)
 538   (prog (fcs* prodl)
 539      (when many*
 540        (set-modulus plim)
 541        (setq dlp (reduce #'max (mapcar #'multideg (cdr (oddelm u))))))
 542      (when (= (length fl) 1) (return (list u)))
 543      (setq prodl (fsplit fl 'v))
 544      (zfactsplit prodl u)
 545      (return fcs*)))
 546
 547 (defun zfactsplit (fl v)
 548   (prog (d)
 549      (cond ((null (cdr fl)) (return (setq fcs* (cons v fcs*))))
 550            ((null (cddr fl))
 551             (setq fl (cadr fl))
 552             (return (setq fcs* (nconc (zff v (car fl) (cadr fl)) fcs*))))
 553            (t (setq fl (cdr fl))
 554               (setq d (zff v (caar fl) (caadr fl)))
 555               (setq v nil)
 556               (zfactsplit (car fl) (car d))
 557               (return (zfactsplit (cadr fl) (cadr d)))))))
 558
 559 (defun split2 (l)
 560   (prog (s)
 561      (setq s (nthcdr (1- (ash (length l) -1)) l))
 562      (setq dn* (copy-list (cdr s)))
 563      (rplacd s nil)
 564      (setq nn* l)))
 565
 566 (defun fsplit (l ind)
 567   (prog (nn* dn*)
 568      (cond ((null (cdr l)) (return l))
 569            ((null (cddr l))
 570             (return (list (apply #'ptimes l) l))))
 571      (split2 l)
 572      (setq nn* (fsplit nn* nil))
 573      (setq dn* (fsplit dn* nil))
 574      (return (list (cond (ind ind) (t (ptimes (car nn*) (car dn*))))
 575                    nn*
 576                    dn*))))
 577
 578 ;; this page contains routines changed for non-monic hack
 579
 580 (defun pexptmod (p n q)
 581   (prog (u x)
 582      (when (pcoefp p) (return (cexpt p n)))
 583      (setq q (cdr q) x (car p))
 584      (cond   ((oddp n)
 585               (setq p (setq u (pgcd1 (cdr p) q)))
 586               (go b))
 587              (t (setq u '(0 1))))
 588      (setq p (cdr p))
 589      a    (setq p (pgcd1 p q))
 590      b    (setq n (ash n -1))
 591      (when (zerop n) (return (cons x u)))
 592      (setq p (ptimes1 p p))
 593      (when (oddp n) (setq u (pgcd1 (ptimes1 u p) q)))
 594      (go a)))
 595
 596 (defun sqfrp (u var)
 597   (cond ((and (equal 0 (ptterm (cdr u) 0)) (equal 0 (ptterm (cdr u) 1)))
 598          nil)
 599         ((onevarp u)
 600          (setq u (pgcd u (pderivative u var)))
 601          (or (numberp u) (alg u)))
 602         (t (quick-sqfr-check u var))))
 603
 604 (declare-top (special p))
 605
 606 (defun fixvl0 (l1 l2 ov)
 607   (prog (a b c)
 608    loop (cond ((null ov) (setq subvar a subval b valist c) (return nil))
 609               ((member (car ov) l1 :test #'eq)
 610                (setq a (cons (car ov) a)
 611                      b (cons (assso (car ov) l1 l2) b)
 612                      c (cons (car b) c)))
 613               (t (setq c (cons 0 c))))
 614    (setq ov (cdr ov))
 615    (go loop)))
 616
 617 (defun assso (a l1 l2)
 618   (prog ()
 619    loop (cond ((null l1) (return nil))
 620               ((eq (car l1) a) (return (car l2))))
 621    (setq l1 (cdr l1) l2 (cdr l2))
 622    (go loop)))
 623
 624 (defun zerohk (l)
 625   (prog (ans i)
 626      (when (null l) (return nil))
 627      ag   (setq ans (car l)
 628                 i (count 0 ans))
 629      loop (setq l (cdr l))
 630      (cond ((null l) (return ans))
 631            ((> (count 0 (car l)) i) (go ag)))
 632      (go loop)))
 633
 634 (defun multfact (poly)
 635   (prog (*inl3 *i* *min* *mx* nn* *odr* lc elm listelm plim origenvar ne var valist val1
 636          ovarlist p subvar subvar1 subval1 subval dlp)
 637      ;;      (declare (special p))
 638      (setq var (car poly) elm (listovars poly)
 639            origenvar genvar
 640            genvar (intersect genvar (if algfac*
 641                                         (delete (car alpha) elm :test #'equal)
 642                                         elm))
 643            ovarlist (butlast genvar)    ;this depends on the order of the above intersection!
 644            nn* (1+ (length ovarlist)))
 645      (setq listelm 0)
 646      (setq lc (caddr poly))
 647      (setq elm 1 *i* 1 ne 1)
 648      (setq subval (reverse poly))
 649      (setq *odr*(putodr (reverse ovarlist)))
 650      (setq val1 (zerohk (nconc (degvector nil 1 lc)
 651                                (cond ((or (> (cadr subval) 0)
 652                                           (> (cadddr subval) 1))
 653                                       (degvector nil 1 (car subval)))))))
 654      (setq subval nil)
 655      (setq p poly)
 656      (when (null val1)
 657        (setq subvar1 ovarlist)
 658        (setq subval1 (polysubst (make-list (length subvar1) :initial-element 0) subvar1))
 659        (go tag))
 660      (fixvl val1 ovarlist)
 661      (fixvl1 val1 ovarlist)
 662      (when subval1 (setq subval1 (polysubst subval1 subvar1)))
 663      (setq subval (polysubst (completevector nil 0 (length subval) 1) subvar))
 664      tag  (fixvl subval1 subvar1)
 665      (setq subval1 nil subvar1 nil)
 666      (fixvl0 subvar subval (reverse ovarlist))
 667      (when algfac* (push (car alpha) genvar))
 668      (setq poly (cpber3 poly p))
 669      (setq genvar origenvar)
 670      (return poly)))
 671
 672 (defun polysubst (a b)
 673   (prog (lc *inl3 d n modulus)
 674      (when modulu* (setq modulus modulu*))
 675      (setq *inl3 t lc (caddr p) n (length a))
 676      loop (setq d (pcsubsty a b lc))
 677      (when (equal 0 d) (go inl))
 678      (let ((modulus nil))
 679        (setq d (pcsubsty a b p)))
 680      (when (sqfrp (pmod d) (car d))
 681        (setq p d) (return a))
 682      inl  (setq a (increaselist a n))
 683      (go loop)))
 684
 685 (declare-top (unspecial p))
 686
 687 (defun fixvl1 (l r)
 688   (prog nil
 689    loop (cond ((null l)
 690                (setq subval1 (nreverse subval1) subvar1 (nreverse subvar1))
 691                (return nil))
 692               ((zerop (car l))
 693                (setq subval1 (cons (car l) subval1))
 694                (setq subvar1 (cons (car r) subvar1))))
 695    (setq l (cdr l))
 696    (setq r (cdr r))
 697    (go loop)))
 698
 699 (defun fixvl (l r)
 700   (prog nil
 701    loop (cond ((null l)
 702                (setq subval (nreverse subval) subvar (nreverse subvar))
 703                (return nil))
 704               ((not (zerop (car l)))
 705                (setq subval (cons (car l) subval))
 706                (setq subvar (cons (car r) subvar))))
 707    (setq l (cdr l))
 708    (setq r (cdr r))
 709    (go loop)))
 710
 711 (defun maxcoef (p)
 712   (maxcoefficient p))
 713
 714 (defun incrlimk (p)
 715   (prog (v min-plim)
 716      (cond (modulu* (setq plim modulu* *prime modulu* limk -1) (return nil))
 717            ((null limk)(setq plim *alpha *prime *alpha limk -1)(return nil)))
 718      (setq v (butlast (pdegreevector p)))
 719      (setq v
 720            (apply
 721             '*
 722             (mapcar
 723              #'(lambda (a b)
 724                  (cond ((equal b 0) 1)
 725                        (t (max (* (simpbinocoef (list '(%binocoef)
 726                                                       a
 727                                                       (ash a -1))
 728                                                 1
 729                                                 t)
 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)))