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(macsyma-module ufact) 13 14(load-macsyma-macros ratmac rzmac) 15 16;; Dense Polynomial Representation 17 18(defun dprep (p) 19 (do ((n (car p)) 20 (e (car p) (f1- e)) 21 (l)) 22 ((< e 0) (cons n (nreverse l))) 23 (cond ((equal e (car p)) 24 (push (cadr p) l) 25 (setq p (cddr p))) 26 (t (push 0 l))))) 27 28(defun dpdisrep (l) 29 (cond ((zerop (car l)) (cadr l)) 30 ((do ((l (nreverse (cdr l)) (cdr l)) 31 (n 0 (f1+ n)) 32 (ll)) 33 ((null l) ll) 34 (or (= (car l) 0) 35 (setq ll (cons n (cons (car l) ll)))))))) 36 37;; not currently called 38;;(DEFUN PGCDU* (P Q) 39;; (COND ((OR (PCOEFP P) (PCOEFP Q)) 1) 40;; ((NULL MODULUS) 41;; (merror "Illegal CALL TO PGCDU")) 42;; ((> (CADR P) (CADR Q)) 43;; (PSIMP (CAR P) (DPDISREP (DPGCD (DPREP (CDR P)) (DPREP (CDR Q)))))) 44;; ((PSIMP (CAR P) (DPDISREP (DPGCD (DPREP (CDR Q)) (DPREP (CDR P)))))))) 45;; 46;;(DEFUN PMODSQFRU (P) 47;; (DO ((DPL (DPSQFR (DPREP (CDR P))) (CDR DPL)) 48;; (PL NIL (CONS (PSIMP (CAR P) (DPDISREP (CDAR DPL))) (CONS (CAAR DPL) PL)))) 49;; ((NULL DPL) PL))) 50 51(defun dpgcd (p q) 52 (if (< (car p) (car q)) (rotatef p q)) 53 (do ((p (copy-list p) q) 54 (q (copy-list q) (dpremquo p q nil))) 55 ((= (car q) 0) 56 (if (= (cadr q) 0) p '(0 1))))) 57 58(defun dpdif (p q) 59 (cond ((> (car p) (car q)) 60 (do ((i (car p) (f1- i)) 61 (pl (cdr p) (cdr pl)) 62 (l nil (cons (car pl) l))) 63 ((= i (car q)) (dpdif1 pl (cdr q) l)) )) 64 ((< (car p) (car q)) 65 (do ((i (car q) (f1- i)) 66 (ql (cdr q) (cdr ql)) 67 (l nil (cons (cminus (car ql)) l))) 68 ((= i (car p)) (dpdif1 (cdr p) ql l)))) 69 (t (dpdif1 (cdr p) (cdr q) nil)))) 70 71(defun dpdif1 (p1 q1 l) 72 (do ((pl p1 (cdr pl)) 73 (ql q1 (cdr ql)) 74 (ll l (cons (cdifference (car pl) (car ql)) ll))) 75 ((null pl) (dpsimp (nreverse ll))))) 76 77(defun dpsimp (pl) (setq pl (ufact-strip-zeroes pl)) 78 (cond ((null pl) '(0 0)) 79 (t (cons (f1- (length pl)) pl)))) 80 81(defun dpderiv (p) 82 (cond ((= 0 (car p)) '(0 0)) 83 (t (do ((l (cdr p) (cdr l)) 84 (i (car p) (f1- i)) 85 (dp nil (cons (ctimes i (car l)) dp))) 86 ((= i 0) (cons (f1- (car p)) (nreverse dp))))))) 87 88(defun dpsqfr (q) ;ASSUMES MOD > DEGREE 89 (do ((c q (dpmodquo c p)) 90 (d (dpderiv q) (dpmodquo d p)) 91 (i 0 (f1+ i)) 92 (p) 93 (pl)) 94 ((= 0 (car c)) pl) 95 (cond (p (setq d (dpdif d (dpderiv c)) 96 p (dpgcd c d)) 97 (and (> (car p) 0) 98 (setq pl (cons (cons i p) pl)))) 99 (t (setq p (dpgcd c d)) 100 (cond ((= (car p) 0) (return (ncons (cons 1 c))))))))) 101 102 103 104(defun dpmodrem (p q) 105 (cond ((< (car p) (car q)) p) 106 ((= (car q) 0) '(0 0)) 107 ((dpremquo (copy-list p) (copy-list q) nil)))) 108 109(defun dpmodquo (p q) 110 (cond ((< (car p) (car q)) '(0 0)) 111 ((= (car q) 0) 112 (cond ((equal (cadr q) 1) p) 113 (t (cons (car p) 114 (mapcar #'(lambda (c) (cquotient c (cadr q))) (cdr p)) 115 )))) 116 ((dpremquo (copy-list p) (copy-list q) t)))) 117 118;; If FLAG is T, return quotient. Otherwise return remainder. 119 120(defun dpremquo (p q flag) 121 (prog (lp lq l alpha) 122 (cond ((= (cadr q) 1) 123 (setq alpha 1)) 124 (t (setq alpha (crecip (cadr q))) 125 (do ((l (cddr q) (cdr l))) 126 ((null l) 127 (rplaca (cdr q) 1)) 128 (rplaca l (ctimes (car l) alpha))))) 129 a (and flag (setq l (cons (ctimes (cadr p) alpha) l))) 130 (setq lp (cddr p) lq (cddr q)) 131 b (rplaca lp (cdifference (car lp) (ctimes (car lq) (cadr p)))) 132 (cond ((null (setq lq (cdr lq))) 133 (do ((e (f1- (car p)) (f1- e)) 134 (pp (cddr p) (cdr pp))) 135 ((null pp) (setq p '(0 0))) 136 (cond ((signp e (car pp)) 137 (and flag (not (< e (car q))) 138 (setq l (cons 0 l)))) 139 ((return (setq p (cons e pp)))))) 140 (cond ((< (car p) (car q)) 141 (return (cond (flag (dpsimp (nreverse l)));GET EXP? 142 (p)))) 143 ((go a)))) 144 (t (setq lp (cdr lp)) 145 (go b))))) 146 147(defun ufact-strip-zeroes (l) 148 (do ((l l (cdr l))) 149 ((null (pzerop (car l))) l))) 150 151(defun cpres1 (a b) 152 (prog (res (v 0) a3) (declare (fixnum v)) 153 (setq a (dprep a) b (dprep b)) 154 (setq res 1) 155 again (setq a3 (dpmodrem a b)) 156 (setq v (boole boole-xor v (logand 1 (car a) (car b) ))) 157 (setq res (ctimes res (cexpt (cadr b) 158 (f- (car a) (car a3))))) 159 (cond ((= 0 (car a3)) 160 (setq res (ctimes res (cexpt (cadr a3) (car b)))) 161 (return (cond ((oddp v) (cminus res)) 162 (t res))) )) 163 (setq a b) 164 (setq b a3) 165 (go again) )) 166