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)) 436