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 1980 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module algfac) 14 15;; this is the alg factor package 16 17;;; Toplevel functions are: CPBGZASS CPTOM 18 19(load-macsyma-macros ratmac) 20 21(declare-top (special tra* trl* *xn var intbs* plim many* split* alc ind p l)) 22 23(defun ziredup (p) 24 (let ((modulus nil) (alpha nil) (minpoly* nil) (algfac* nil) 25 (gauss nil) (tellratlist nil) (many* nil) 26 (mm* 1) 27 ($gcd '$ez)) 28 (null (cddr(pfactor p))))) 29 30(defun intbasehk (p) 31 (prog (modulus) 32 (setq modulus plim) 33 (setq p (pctimes intbs* p)) 34 (setq modulus nil) 35 (return (car (ratreduce p intbs*))))) 36 37(defun findibase (p) 38 (prog (mainvar) 39 (setq mainvar (car p)) 40 (setq p (redresult p (pderivative p mainvar))) 41 (setq p (cfactorw p)) 42 (setq mainvar 1) 43 loop (when (null p) (return mainvar)) 44 (setq mainvar (* mainvar (expt (car p) (quotient (cadr p) 2)))) 45 (setq p (cddr p)) 46 (go loop))) 47 48 49(defun cpbgzass (qlist v m) 50 (prog (f y vj factors u w lc j p2 fnj fnq oldfac) 51 (cond ((equal m 1) 52 (return (list v))) 53 ((equal m (cadr v)) 54 (return (let ((var (list var 1 1))) 55 (gfsplit v))))) 56 (setq f (pmod v)) 57 (setq lc (caddr f)) 58 (setq f (monize f)) 59 (setq p2 1 60 qlist (cdr (nreverse qlist))) 61 (setq oldfac (list nil f)) 62 nextq(setq v (car qlist)) 63 (setq qlist (cdr qlist)) 64 (setq j (findses v f)) 65 (setq oldfac (nconc oldfac fnq)) 66 (setq fnq nil) 67 incrj(setq factors (nconc oldfac fnj)) 68 (setq fnj nil) 69 (setq vj (pplus v (car j)) 70 j (cdr j)) 71 tag2 (setq u (cadr factors)) 72 (setq w (pgcdu vj u)) 73 (when (or (numberp w) (and alpha (alg w))(= (cadr u) (cadr w))) 74 (go nextfac)) 75 (setq y (car (pmodquo u w))) 76 (setq fnq (cons w fnq)) 77 (setq fnj (cons y fnj)) 78 (incf p2) 79 (rplacd factors (cddr factors)) 80 (if (equal p2 m) 81 (go out) 82 (go tag1)) 83 nextfac 84 (setq factors (cdr factors)) 85 tag1 (cond ((cdr factors) 86 (go tag2)) 87 (j (go incrj)) 88 (qlist (go nextq))) 89 out (setq fnq (nconc fnq fnj (cdr oldfac))) 90 (return (cons (ptimes lc (car fnq)) (cdr fnq))))) 91 92 93;; The function PMONZ used to be defined here. It is also defined in 94;; RAT;RAT3A and BMT claims the definitions are equivalent. 95 96(defun findses (g f) 97 (prog (var tra* trl*) 98 (setq g (zassg (cdr g) (cdr f) (car g))) 99 (setq var (list (car f) 1 1)) 100 (setq f (gfsplit g)) 101 (return (mapcar #'(lambda (a) (car (last a))) f)))) 102 103(defun coefvec (p n vec) 104 (prog nil 105 loop (when (zerop n) (return vec)) 106 (decf n) 107 (push (ptterm p n) vec) 108 (go loop))) 109 110(defun zassg (g f var) 111 (prog (i mat gn ans n) 112 (setq n (car f)) 113 (setq gn g) 114 (setq i 1 115 mat (list (coefvec '(0 1) n (list 1)))) 116 (go on) 117 loop (incf i) 118 (setq gn (pgcd1 (ptimes1 gn g) f)) 119 on (setq ans (lindep mat (coefvec gn n (list (list var i 1))))) 120 (cond (ans (return ans))) 121 (go loop))) 122 123(defun divl (j a) 124 (mapcar #'(lambda (l) (car (pmodquo l a))) j)) 125 126;; (DEFUN PADDROWS (A B) (MAPCAR (FUNCTION PPLUS) A B)) 127 128(defun pdifrows (a b) 129 (mapcar #'pdifference a b)) 130 131(defun ptimesrow (var row) 132 (mapcar #'(lambda (a) (ptimes var a)) row)) 133 134(defun ddiv (j) 135 (prog (a b) 136 (setq b j) 137 ag (setq a (car b)) 138 (cond ((zerop a) 139 (setq b (cdr b)) 140 (go ag))) 141 (return (divl j a)))) 142 143(defun lindep (mat vec) 144 (prog (e d m row rowd vecd) 145 (setq m mat) 146 (cond ((equal 0. (car vec)) (setq vec (cdr vec))) 147 (t (setq vec (pdifrows (cdr vec) (ptimesrow (car vec) (cdar mat)))))) 148 loop (cond ((null (cdr m)) 149 (cond ((zerolp (cdr (reverse vec))) 150 (return (car (last vec)))) 151 (t (rplacd m (cons (ddiv vec) (cdr m))) 152 (return nil))))) 153 (setq row (cadr m)) 154 (setq rowd row vecd vec) 155 loop1(setq d (car rowd)) 156 (setq e (car vecd)) 157 (cond ((equal 0 e) 158 (cond ((equal 0 d) 159 (setq vecd (cdr vecd) rowd (cdr rowd)) 160 (go loop1)) 161 (t (setq vec (cdr vec)) (setq m (cdr m)) (go loop)))) 162 ((equal 0 d) 163 (rplacd m 164 (cons (divl vec e) (mapcar (function cdr) (cdr m)))) 165 (return nil))) 166 (setq vec (pdifrows (cdr vec) (ptimesrow e (cdr row)))) 167 (setq m (cdr m)) 168 (go loop))) 169 170(defun gfsplit (f) 171 (prog (tr fl (n 0) ans tra* (i 0) nfl) 172 (setq fl (list f) n (cadr f)) 173 loop (cond ((null fl) 174 (cond ((null nfl) 175 (cond ((= n (length ans)) 176 (setq trl* nil) 177 (return ans)) 178 (t (merror (intl:gettext "GFSPLIT: unknown error."))))) 179 (t 180 (setq fl nfl 181 nfl nil 182 i (1+ i)))))) 183 (setq f (car fl) 184 fl (cdr fl)) 185 (cond ((> i mm*) 186 (merror (intl:gettext "GFSPLIT: unknown error.")))) 187 (setq tr (tracemod0 f i)) 188 (cond ((or (pcoefp tr) (and algfac* (alg tr))) 189 (setq nfl (cons f nfl)) 190 (go loop))) 191 (setq f (cpbg0 tr f)) 192 (setq ans (nconc ans (car f))) 193 (when (null (cdr f)) (go loop)) 194 (setq nfl (nconc nfl (cdr f))) 195 (go loop))) 196 197(defun cpbg0 (tr f) 198 (prog (m f1 f2 g alc trm) 199 (setq m 0) 200 (cond ((and (not (numberp (caddr tr))) (alg (caddr tr))) 201 (setq alc (painvmod (caddr tr)) tr (ptimes alc tr))) 202 (t (setq alc 1.))) 203 bk (cond ((pcoefp f) 204 (return (cond ((and (null f1) (null f2)) 205 ;; NOTE TO TRANSLATORS: MEANING OF NEXT MESSAGE IS OBSCURE 206 (merror (intl:gettext "CPBG0: wrong trace."))) 207 (t 208 (cons f1 f2))))) 209 ((equal (cadr f) 1) 210 (return (cons (cons f f1) f2))) 211 ((equal m modulus) 212 (return (cons f1 (cons f f2))))) 213 (setq trm (pdifference tr (ptimes m alc))) 214 (setq g (pgcdu trm f)) 215 (cond ((or (numberp g) (and alpha (alg g))) 216 (incf m) 217 (go bk))) 218 (setq f (car (pmodquo f g))) 219 (cond ((equal (cadr g) 1) (setq f1 (cons g f1))) 220 (t (setq f2 (cons g f2)))) 221 (go bk))) 222 223(defun cpol2p (p var) 224 (prog((i 0) ans) 225 (setq p (nreverse p)) 226 loop (cond ((null p) (return (cons var ans))) 227 ((equal 0 (car p)) nil) 228 (t (setq ans (cons i (cons (car p) ans))))) 229 (setq p (cdr p) 230 i (1+ i)) 231 (go loop))) 232 233(defun tracemod (v) 234 (prog (ans tr qlarge term) 235 (setq ans 0 236 tr (nreverse trl*) 237 trl* nil) 238 (cond ((and (atom (caar tr)) (not (numberp (caar tr)))) 239 (setq qlarge t))) 240 loop (when (null tr) (return ans)) 241 (setq term (if qlarge 242 (car tr) 243 (cpol2p (car tr) v)) 244 tr (cdr tr)) 245 (setq ans (pplus ans term)) 246 (setq trl* (cons term trl*)) 247 (go loop))) 248 249(defun otracemod (term q m prime) 250 (prog (ans i) 251 (setq ans term 252 i 1 253 trl* (list term)) 254 loop (when (equal i m) (return ans)) 255 (setq ans (pplus ans (setq term (pexptmod term prime q)))) 256 (setq trl* (cons term trl*)) 257 (incf i) 258 (go loop))) 259 260(defun tracemod0 (q i) 261 (prog (l ans a dl) 262 (cond ((= i 0) (return (if trl* 263 (tracemod (car q)) 264 (otracemod var q mm* modulus)))) 265 (trl* (setq dl trl* 266 trl* (mapcar #'(lambda(x) 267 (cons (car x) (pgcd1 (cdr x) (cdr q)))) trl*)))) 268 (cond (tra* (go tag)) 269 (t (setq l (cdr trl*) 270 tra* (list alpha) 271 a alpha))) 272 loop (when (null l) (go tag)) 273 (setq l (cdr l) 274 a (pexpt a modulus) 275 tra* (cons a tra*)) 276 (go loop) 277 tag 278 (setq ans (tracemod1 i tra* trl*)) 279 (when dl (setq trl* dl)) 280 (return ans))) 281 282(defun tracemod1 (n a l) 283 (prog (ans) 284 (setq ans 0) 285 loop (when (null l) (return ans)) 286 (setq ans (pplus ans (ptimes (pexpt (car a) n) (car l)))) 287 (setq l (cdr l) 288 a (cdr a)) 289 (go loop))) 290 291;; The way arrays are manipulated has been changed to make this code reentrant. 292;; Previously, arrays were kept on the array properties of symbols. Now, the 293;; symbols are bound to the arrays, so they can be rebound on re-entry. 294;; The ANOTYPE, INVC, and FCTC arrays are set up in RAT;FACTOR. 295 296(declare-top (special anotype invc fctc)) 297 298(defmacro a (row col) 299 `(aref anotype ,row ,col)) 300 301(defmacro invc (ind) 302 `(aref invc ,ind)) 303 304(defmacro fctc (ind) 305 `(aref fctc ,ind)) 306 307(defun cptomexp (p m u n) 308 (prog (b1 b2 j n-1 i l) 309 (setq b1 (x**q1 (list (car u) 1 1) u m p)) 310 (cond ((equal (cdr b1) '(1 1)) 311 (setq split* t) 312 (return nil))) 313 (setq b2 b1 j 1 n-1 (1- n)) 314 (go tag1) 315 tag (incf j) 316 (when (= j n) (return nil)) 317 (setq b1 (pmodrem(ptimes b1 b2) u)) 318 tag1 (setq l (p2cpol b1 n-1) i n-1) 319 sharp2 (when (null l) (go on)) 320 (setf (a j i) (car l)) 321 (setq l (cdr l)) 322 (setq i (1- i)) 323 (go sharp2) 324 on (setf (a j j) (pdifference (a j j) 1)) 325 (go tag))) 326 327(defvar thr* 100) 328 329(defun cptom (p m u n) 330 (prog (( q (expt p m)) l s *xn (j 0) (i 0) ind n-1) 331 (declare (special q i j)) 332 (setq n-1 (1- n)) 333 (when (> q thr*) (return (cptomexp p m u n))) 334 loop (incf j) 335 (cond ((= j n) (return nil)) 336 (ind (go sa)) 337 (t 338 (setq *xn (mapcar #'pminus (p2cpol (cddr u) n-1)) 339 s (x**q (p2cpol(list var 1 1) n-1) p m) 340 ind t))) 341 (go st) 342 sa (cptimesxa s q) 343 st (cond ((and (= j 1) 344 (equal '(1 0) (last s 2)) 345 (= 1 (length (delete 0 (copy-tree s) :test #'equal)))) 346 (return (setq split* t)))) 347 (setq l s) 348 (setq i n-1) 349 sharp2 (when (null l) (go on)) 350 (setf (a j i) (car l)) 351 (setq l (cdr l)) 352 (decf i) 353 (go sharp2) 354 on (setf (a j j) (pdifference (a j j) 1)) 355 (go loop))) 356 357(defun cptimesxa (p i) 358 (prog (xn q lc) 359 ag (when (= i 0) (return p)) 360 (setq xn *xn 361 q p 362 lc (car p)) 363 loop (cond ((cdr q) 364 (rplaca q (pplus (cadr q) (ptimes lc (car xn)))) 365 (setq q (cdr q) xn (cdr xn))) 366 (t (rplaca q (ptimes lc (car xn))) 367 (decf i) 368 (go ag))) 369 (go loop))) 370 371(defun x**q (x p m) 372 (prog ((i 1) (pp 1) (d 0)) 373 (setq i 1 trl* (list x) pp 1) 374 loop (when (= i m) (return (cptimesxa x (- (* pp p) pp)))) 375 (setq d pp) 376 (cptimesxa x (- (setq pp (* pp p)) d)) 377 (setq trl* (cons (copy-tree x) trl*)) 378 (incf i) 379 (go loop))) 380 381(defun cmnull (n) 382 (prog (nullsp (sub1n (1- n)) mone (k 1) (j 0) (s 0) nullv (i 0) vj m aks) 383 (setq mone (cmod -1)) 384 sharp (when (> i sub1n) (go on)) 385 (setf (fctc i) -1) 386 (setf (invc i) -1) 387 (incf i) 388 (go sharp) 389 on (setq k 1 nullsp (list 1)) 390 n2 (when (> k sub1n) (return nullsp)) 391 (setq j 0) 392 n3a (cond ((> j sub1n) (go null)) 393 ((or (equal (a k j) 0) (> (fctc j) -1)) 394 (incf j) 395 (go n3a))) 396 (setf (invc k) j) 397 (setf (fctc j) k) 398 (setq m (a k j)) 399 (setq m (ptimes mone (painvmod m))) 400 (setq s k) 401 sharp1 (when (> s sub1n) (go on1)) 402 (setf (a s j) (ptimes m (a s j))) 403 (incf s) 404 (go sharp1) 405 on1 406 (setq s 0) 407 sharp2 (when (> s sub1n) (go nextk)) 408 (cond ((= s j) nil) 409 (t (prog (i) 410 (setq aks (a k s)) 411 (setq i k) 412 sharp3 (when (> i sub1n) (return nil)) 413 (setf (a i s) (pplus (a i s) (ptimes (a i j) aks))) 414 (incf i) 415 (go sharp3)))) 416 (incf s) 417 (go sharp2) 418 null (setq nullv nil) 419 (setq s 0) 420 sharp4 (cond ((> s sub1n) (go on4)) 421 ((= s k) (setq nullv (cons s (cons 1 nullv)))) 422 ((> (invc s) -1) 423 (setq vj (a k (invc s))) 424 (cond ((equal vj 0) nil) 425 (t (setq nullv (cons s (cons vj nullv))))))) 426 (incf s) 427 (go sharp4) 428 on4 (cond ((equal (car nullv) 0) (setq nullv (cadr nullv))) 429 ((setq nullv (cons var nullv)))) 430 (setq nullsp (cons nullv nullsp)) 431 nextk (incf k) 432 (go n2))) 433