src/ezgcd.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 ezgcd)
  14
  15 (declare-top (special genvar lcprod svals svars oldsvars oldsvals
  16                       valflag $gcd pl0 d0 degd0 var
  17                       many* tempprime ovarlist valist modulus
  18                       zl *prime plim nn* ne nn*-1 dlp
  19                       ez1skip svalsl nsvals $algebraic
  20                       lc1 oldlc limk *alpha $ratfac))
  21
  22 (load-macsyma-macros ratmac)
  23
  24 (defun ezgcd2 (f g)
  25   (prog (allvars)
  26      (setq allvars (union* (listovars f) (listovars g)))
  27      (cond ((> 2 (length allvars))
  28             (setq allvars (newgcd f g modulus))
  29             (cond ((cdr allvars) (return allvars))
  30                   (t (return (list (setq allvars (car allvars))
  31                                    (pquotient f allvars)
  32                                    (pquotient g allvars)))))))
  33      (setq allvars (sort allvars #'pointergp))
  34      (return (ezgcd (list f g) allvars modulus))))
  35
  36 (defun newgcdcall (p q)
  37   (car (newgcd p q modulus)))
  38
  39 (defun gcdl (pl)
  40   (do ((d (car pl) (pgcd d (car l)))
  41        (l (cdr pl) (cdr l)))
  42       ((or (null l) (eql d 1)) d)))
  43
  44 (defun newgcdl (pl)
  45   (let (($gcd '$mod))
  46     (gcdl pl)))
  47
  48 (defun oldgcdl (elt pl)
  49   (let (($gcd '$red))
  50     (gcdl (cons elt pl))))
  51
  52 (defun oldgcdcall (pfl)
  53   (let ((a (oldgcdl (car pfl) (cdr pfl))))
  54     (cons a (mapcar #'(lambda (h) (pquotient h a)) pfl))))
  55
  56 (defun non0rand (modulus)
  57   (do ((r))
  58       ((not (zerop (setq r (cmod (random 1000))))) r)))
  59
  60 (defun getgoodvals (varl lcp)
  61   (mapcar #'(lambda (v) (do ((val 0 (non0rand tempprime)) (temp))
  62                             ((not (pzerop (setq temp (pcsubsty val v lcp))))
  63                              (setq lcp temp) val)))
  64           varl))
  65
  66 (defun evmap (vals pl)
  67   (prog (pl0 d0)
  68      (cond ((equal nsvals (length svalsl)) (return nil)))
  69      (cond (valflag (go newvals)))
  70      (setq vals (getgoodvals svars lcprod))
  71      again (cond ((member vals svalsl :test #'equal)
  72                  (setq vals (rand (length svars) tempprime))
  73                  (go again)))
  74      (setq valflag t svalsl (cons vals svalsl))
  75      (go end)
  76      newvals
  77      (setq pl0 (rand (length svars) tempprime))
  78      (cond ((member pl0 svalsl :test #'equal) (go newvals))
  79            (t (setq vals pl0 svalsl (cons vals svalsl))))
  80      (cond ((equal 0 (pcsub lcprod vals svars))
  81             (cond ((equal nsvals (length svalsl))
  82                    (return nil))
  83                   (t (go newvals)))))
  84      ;;      END  (GETD0 PL(SETQ VALS(SUBST 1. 0.  VALS))) WHAT WAS SUBST FOR?
  85      end  (getd0 pl vals)
  86      (return (list vals pl0 d0))))
  87
  88 (defun degodr (a b)
  89   (cond ((numberp a) nil)
  90         ((numberp b) t)
  91         (t (> (cadr a) (cadr b)))))
  92
  93 (defun evtildegless (pl)
  94   (prog (evout npl0 nd0 ndeg)
  95    again (setq evout (evmap svals pl))
  96    (cond (evout (setq npl0 (cadr evout) nd0 (caddr evout)))
  97          (t (return nil)))
  98    (cond ((numberp nd0) (setq ndeg 0)) (t (setq ndeg (cadr nd0))))
  99    (when (or (= degd0 ndeg) (> ndeg degd0)) (go again))
 100    (return (setq degd0 ndeg pl0 npl0 d0 nd0 svals (car evout)))))
 101
 102 (defun ptimesmerge (pl1 pl2)
 103   (cond (pl1 (cons (ptimes (car pl1) (car pl2))
 104                    (ptimesmerge (cdr pl1) (cdr pl2))))
 105         (t nil)))
 106
 107 (defun ez1call (builder factrs lc1 valist ovarlist)
 108   (prog (*prime plim nn* ne nn*-1 zl zfactr oldlc lcd0
 109          dlp limk genvar mult subval subvar)
 110      (declare (special subval subvar))
 111      (setq oldlc (caddr builder))
 112      (cond ((not (equal 1 lc1))
 113             (setq builder (ptimes builder lc1))))
 114      (setq genvar (append ovarlist (list (car builder))))
 115      (cond (modulus
 116             (setq *prime modulus plim modulus limk -1)
 117             (go mod))
 118            (t (setq *prime (max (norm builder)
 119                                 (maxcoefficient (car factrs))
 120                                 (maxcoefficient (cadr factrs))))))
 121      (cond ((> *prime *alpha)
 122             (prog (newmodulus)
 123                (setq newmodulus (* *alpha *alpha)
 124                      limk 0)
 125                again(cond ((> newmodulus *prime)
 126                            (setq *prime *alpha plim newmodulus))
 127                           (t (setq limk (1+ limk) newmodulus (* newmodulus newmodulus))
 128                              (go again)))))
 129            (t (setq limk -1 *prime *alpha plim *alpha)))
 130      mod  (setq nn* (1+ (setq ne (setq nn*-1 (length ovarlist)))))
 131      (setq zl (completevector nil 1 nn* 0))
 132      (fixvl valist ovarlist)
 133      (cond ((equal 1 lc1)
 134             (setq modulus plim builder (newrep builder))
 135             (setq dlp  (loop for x in (cdr (oddelm builder))
 136                           maximize (multideg x)))
 137             (setq zfactr (z1 builder (car factrs)(cadr factrs)))
 138             (setq zfactr (restorelc zfactr (caddr builder)))
 139             (return (oldrep(cadr zfactr)))))
 140      (setq modulus plim lcd0 (caddar factrs))
 141      (setq mult (ctimes (pcsub lc1 svals svars)
 142                         (crecip lcd0)))
 143      (setq factrs (list (ptimes mult (car factrs))
 144                         (ptimes lcd0 (cadr factrs))))
 145      (setq builder (newrep builder))
 146      (setq dlp (loop for x in (cdr (oddelm builder))
 147                   maximize (multideg x)))
 148      (setq zfactr (z1 builder (car factrs) (cadr factrs)))
 149      (setq zfactr (pmod (oldrep (car zfactr))))
 150      (return (cadr (let ((modulus nil))
 151                      (fastcont zfactr))))))
 152
 153 (defun getd0 (tpl tvals)
 154   (prog (c)
 155      (setq d0 (pcsub (car tpl) tvals svars)
 156            pl0 (list d0) tpl (cdr tpl))
 157      loop (cond ((null tpl) (return d0)))
 158      (setq c (pcsub (car tpl) tvals svars)
 159            d0 (newgcdcall d0 c))
 160      (cond ((numberp d0) (return (setq d0 1))))
 161      (setq pl0 (append pl0 (list c)) tpl (cdr tpl))
 162      (go loop)))
 163
 164 (defun numberinlistp (l)
 165   (do ((l l (cdr l))) ((null l))
 166     (and (numberp (car l)) (return (car l)))))
 167
 168 (defun ezgcd (pfl vl modulus)
 169   (prog (svars svals valflag tempprime pfcontl contgcd contcofactl
 170          pl nsvars nsvals svalsl lcprod gcdlcs lcpl evmapout
 171          pl0 d0 d degd0 degd0n d0n pl0n temp tryagain cofact0
 172          pcofactl ith builder var termcont tcontl $algebraic)
 173      (cond ((setq temp (numberinlistp pfl))
 174             (cond ((or (member 1 pfl) (member -1 pfl))
 175                    (return (cons 1 pfl))))
 176             (setq temp (oldgcdl temp (remove temp pfl :test #'equal))
 177                   pl (mapcar #'(lambda(h) (pquotient h temp)) pfl))
 178             (return (cons temp pl))))
 179      (setq svars (cdr vl) var (car vl))
 180      (cond (svars (setq many* t))
 181            (t (return (cons (setq d (newgcdl pfl))
 182                             (mapcar #'(lambda(h) (pquotient h d)) pfl)))))
 183      (cond (modulus (setq tempprime modulus))
 184            (t (setq tempprime 13.)))
 185      (setq tcontl (mapcar #'ptermcont pfl)
 186            pfl (mapcar #'cadr tcontl)
 187            tcontl (mapcar #'car tcontl))
 188      (setq termcont (oldgcdcall tcontl)
 189            tcontl (cdr termcont)
 190            termcont (car termcont))
 191      (cond ((setq temp (numberinlistp pfl))
 192             (setq d (oldgcdl temp (remove temp pfl :test #'equal))
 193                   pcofactl
 194                   (mapcar #'(lambda(h) (pquotient h d)) pfl))
 195             (setq contgcd termcont contcofactl tcontl)
 196             (go out)))
 197      (setq pfcontl (mapcar #'(lambda(h) (if (eq var (car h))
 198                                             (fastcont h)
 199                                             (list h 1)))
 200                            pfl))
 201      (setq pfl (mapcar #'cadr pfcontl)
 202            pfcontl (mapcar #'car pfcontl))
 203      (setq contgcd (ezgcd pfcontl svars modulus) pfcontl nil
 204            contcofactl (ptimesmerge tcontl (cdr contgcd))
 205            contgcd (ptimes termcont (car contgcd)))
 206      (cond ((numberinlistp pfl)
 207             (setq d 1 pcofactl pfl) (go out)))
 208      (setq temp (listovarsl pfl))
 209      (cond ((setq temp (intersection svars temp :test #'equal)) nil)
 210            (t (setq d (newgcdl pfl)) (go end)))
 211      (setq pl (bbsort pfl #'degodr))
 212      (setq nsvars (length svars))
 213      (do ((i nsvars (1- i)))
 214          ((zerop i))
 215        (push 0 svals))
 216      (setq lcprod 1 svalsl (list svals)
 217            nsvals (expt tempprime (length svars)))
 218      (do ((l (mapcar #'caddr pl) (cdr l)))
 219          ((null l))
 220        (setq lcprod (ptimes lcprod (car l))))
 221      (cond ((equal 0 (pcsub lcprod svals svars))
 222             (setq evmapout (evmap svals pl))
 223             (cond (evmapout (setq svals (car evmapout)
 224                                   pl0 (cadr evmapout)
 225                                   d0 (caddr evmapout)))
 226                   (t (desetq (d . pcofactl) (oldgcdcall pfl))
 227                      (go out))))
 228            (t (setq valflag t) (getd0 pl svals)))
 229      (cond ((numberp d0) (setq degd0 0))
 230            (t (setq degd0 (cadr d0))))
 231      testd0
 232      (cond ((equal 1 d0) (setq d 1)
 233             (setq d 1 pcofactl pfl) (go out)))
 234      (cond (degd0n (go testcofact)))
 235      anothersvals
 236      (setq evmapout (evmap svals pl))
 237      (cond (evmapout (setq pl0n (cadr evmapout)
 238                            d0n (caddr evmapout)
 239                            evmapout (car evmapout)))
 240            (t (desetq (d . pcofactl) (oldgcdcall pfl))
 241               (go out)))
 242      (cond ((numberp d0n) (setq degd0n 0))
 243            (t (setq degd0n (cadr d0n))))
 244      (cond ((> degd0 degd0n)
 245             (setq degd0 degd0n pl0 pl0n d0 d0n svals evmapout)
 246             (go anothersvals)))
 247      (cond ((equal degd0 degd0n) (go testd0)) (t (go anothersvals)))
 248      testcofact
 249      (cond ((equal degd0 (cadar pl0)) nil) (t (go testgcd)))
 250      (setq d (car pl) temp pfl pcofactl nil)
 251      loop (cond (temp (setq d0n (eztestdivide (car temp) d)))
 252                 (t (setq ez1skip t) (go out)))
 253      (cond (d0n (setq pcofactl (append pcofactl (list d0n)))
 254                 (setq temp (cdr temp)) (go loop))
 255            (t (cond ((evtildegless pl)
 256                      (setq degd0n nil) (go testd0))
 257                     (t (desetq (d . pcofactl) (oldgcdcall pfl))
 258                        (go out)))))
 259      testgcd
 260      (setq ith 1. temp pl0)
 261      next (cond (temp nil)
 262                 (t (cond (tryagain  (setq d (nonsqfrcase pfl vl)
 263                                           pcofactl (cdr d)
 264                                           d (car d))
 265                                     (go out))
 266                          (t (setq degd0 degd0n pl0 pl0n d0 d0n
 267                                   degd0n nil pl0n nil d0n nil
 268                                   svals evmapout tryagain t)
 269                             (go testgcd)))))
 270      (setq cofact0 (pquotient (car temp) d0))
 271      (cond ((numberp (newgcdcall d0 cofact0))
 272             (setq builder (ith pl ith))
 273             (cond ((intersection (listovars builder) svars :test #'equal)
 274                    (go callez1)))))
 275      (setq temp (cdr temp) ith (1+ ith)) (go next)
 276      callez1
 277      (setq lcpl (mapcar #'caddr pl)
 278            gcdlcs (car (ezgcd lcpl svars modulus))
 279            lcpl nil)
 280      (setq d (ez1call builder
 281                       (list d0 cofact0) gcdlcs
 282                       (reverse svals)
 283                       (reverse svars)))
 284      (setq modulus nil)
 285      end  (setq pcofactl nil temp pfl)
 286      (cond ((pminusp d) (setq d (pminus d))))
 287      loop1(cond (temp (setq cofact0 (eztestdivide (car temp) d)))
 288                 (t (setq ez1skip nil) (go out)))
 289      (cond (cofact0 (setq pcofactl (append pcofactl (list cofact0)))
 290                     (setq temp (cdr temp)) (go loop1))
 291            (t (cond ((evtildegless pl)
 292                      (setq degd0n nil) (go testd0))
 293                     (t (desetq (d . pcofactl) (oldgcdcall pfl))
 294                        (go out)))))
 295      out  (setq oldsvars svars oldsvals svals)
 296      (return (cons (ptimes contgcd d)
 297                    (ptimesmerge contcofactl pcofactl)))))
 298
 299 (defun listovarsl (plist)
 300   (prog (allvarsl allvars)
 301      (setq allvarsl (mapcar #'listovars plist))
 302      (setq allvars (car allvarsl))
 303      (do ((l (cdr allvarsl) (cdr l)))
 304          ((null l))
 305        (setq allvars (union* allvars (car l))))
 306      (return allvars)))
 307
 308 (defmfun $ezgcd (&rest args)
 309   (prog (pfl allvars presult flag genvar denom pfl2)
 310      ;;need if genvar doesn't shrink
 311      (when (null args)
 312        (wna-err '$ezgcd))
 313      (when (some #'$ratp args)
 314        (setq flag t))
 315      (setq pfl (mapcar #'(lambda (h) (cdr (ratf h))) args))
 316      (setq pfl2 (list 1))
 317      (do ((lcm (cdar pfl))
 318           (l (cdr pfl) (cdr l))
 319           (cof1)
 320           (cof2))
 321          ((null l) (setq denom lcm))
 322        (desetq (lcm cof1 cof2) (plcmcofacts lcm (cdar l)))
 323        (unless (equal cof1 1)
 324          (mapcar #'(lambda (x) (ptimes x cof1)) pfl2))
 325        (push cof2 pfl2))
 326      (setq pfl (mapcar #'car pfl))
 327      (setq allvars (sort (listovarsl pfl) #'pointergp))
 328      (setq presult (if $ratfac
 329                        (let (($gcd '$ez))
 330                          (facmgcd pfl))
 331                        (ezgcd pfl allvars modulus)))
 332      (setq presult (cons (cons (car presult) denom)
 333                          (if (equal denom 1)
 334                              (cdr presult)
 335                              (mapcar #'ptimes (cdr presult) pfl2))))
 336      (setq presult (cons '(mlist)
 337                          (cons (rdis* (car presult))
 338                                (mapcar #'pdis* (cdr presult)))))
 339      (return (if flag presult ($totaldisrep presult)))))
 340
 341 (defun insrt (nth elt l)
 342   (if (eql nth 1)
 343       (cons elt l)
 344       (cons (car l) (insrt (1- nth) elt (cdr l)))))
 345
 346 (defun nonsqfrcase(pl vl)
 347   (prog (d f ptr)
 348      (do ((dl pl (cdr dl))
 349           (pt 1 (1+ pt)))
 350          ((intersection (cdr vl) (listovars (car dl)) :test #'equal)
 351           (setq f (car dl) ptr pt)))
 352      (setq d (ezgcd (list f (pderivative f (car f))) vl modulus)
 353            pl (ezgcd (cons (cadr d) (remove f pl :test #'equal)) vl modulus)
 354            pl (cons (car pl) (cons (car d) (cddr pl)))
 355            d (car pl))
 356      loop  (setq pl (ezgcd pl vl modulus))
 357      (cond ((equal 1 (car pl))
 358             (return (cons d (insrt ptr (pquotient f d) (cdddr pl))))))
 359      (setq d (ptimes (car pl) d)
 360            pl (cons (car pl) (cddr pl)))
 361      (go loop)))
 362
 363 (defun eztestdivide (x y)
 364   (when (or (pcoefp x) (pcoefp y)
 365             (ignore-rat-err (pquotient (car (last x)) (car (last y)))))
 366     (ignore-rat-err (pquotient x y))))
 367
 368 (defun noterms (p)
 369   (cond ((pcoefp p) 1)
 370         (t (do ((nt (noterms (caddr p)) (+ nt (noterms (cadr p))))
 371                 (p (cdddr p) (cddr p)))
 372                ((null p) nt)))))
 373
 374 (defun fastcont (p)
 375   (prog (oldgenvar var tppl tcontl tcont coefvarl temp small1 small2 ans minus?)
 376      (cond ((univar (cdr p)) (return (oldcontent p)))
 377            (t (setq oldgenvar genvar)
 378               (setq var (car p))
 379               ;; intersect must maintain order of genvar list
 380               (setq genvar (remove var (intersect (cdr genvar) (listovars p))
 381                                    :test #'equal))))
 382      (cond ((pminusp p) (setq p (pminus p) minus? t)))
 383      (setq tppl (oddelm (cddr p)))
 384      (cond ((null (cdr tppl))
 385             (setq tcont 1)
 386             (setq ans (car tppl))
 387             (go out)))
 388      (setq tcontl (mapcar #'pmindegvec tppl))
 389      (setq tppl (mapcar #'(lambda(x y) (pquotient x (degvecdisrep y)))
 390                         tppl tcontl))
 391      (setq tcont (car tcontl))
 392      (do ((l (cdr tcontl) (cdr l))) ((null l))
 393        (setq tcont (mapcar #'(lambda (x y) (min x y))
 394                            tcont (car l))))
 395      (setq tcontl nil)
 396      (setq tcont (degvecdisrep tcont))
 397      (setq genvar oldgenvar)
 398      (cond ((setq temp (numberinlistp tppl))
 399             (cond ((or (member 1 tppl) (member -1 tppl))
 400                    (setq ans 1))
 401                   (t (setq ans (oldgcdl temp (delete temp tppl :test #'equal)))))
 402             (go out)))
 403      (cond ((> 4 (length tppl))
 404             (setq tppl (bbsort tppl #'(lambda(a b) (> (length a) (length b)))))
 405             (go skip)))
 406      (setq coefvarl (mapcar #'listovars tppl))
 407      (setq temp (car coefvarl))
 408      (setq coefvarl (cdr coefvarl))
 409      loop (cond ((null coefvarl) nil)
 410                 (t (cond ((null (setq temp (intersection temp (car coefvarl) :test #'equal)))
 411                           (setq ans 1) (go out))
 412                          (t (setq coefvarl (cdr coefvarl)) (go loop)))))
 413      (setq temp (mapcar #'noterms tppl))
 414      (setq tppl (mapcar #'(lambda (x y) (list x y))
 415                         temp tppl))
 416      (setq tppl (bbsort tppl #'(lambda(x y) (> (car x) (car y)))))
 417      (setq tppl (mapcar #'cadr tppl))
 418      skip (setq small1 (car tppl))
 419      (setq small2 (cadr tppl))
 420      (setq ans (pgcd small1 small2))
 421      (cond ((eql 1 ans) (go out))
 422            ((eql -1 ans) (setq ans 1) (go out)))
 423      (cond ((cddr tppl) (setq ans (cons ans (cddr tppl))))
 424            (t (go out)))
 425      (setq temp (sort (listovarsl ans) #'pointergp))
 426      (setq ans (car (ezgcd ans temp modulus)))
 427      out (setq tcont (ptimes tcont ans))
 428      (setq p (pquotient p tcont))
 429      (cond (minus? (setq tcont (pminus tcont))))
 430      (return (list tcont p))))
 431
 432 (declare-top (unspecial lcprod svals svars oldsvars oldsvals
 433                         valflag pl0 d0 degd0 var many* tempprime ovarlist valist
 434                         zl *prime plim nn* ne nn*-1 dlp
 435                         ez1skip svalsl nsvals lc1 oldlc limk))