Changeset 71
- Timestamp:
- 2015-06-05T11:32:41-07:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/f4grobner/grobner.lisp
r70 r71 134 134 135 135 136 137 138 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;139 ;;140 ;; Low-level polynomial arithmetic done on141 ;; lists of terms142 ;;143 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;144 145 (defmacro termlist-lt (p) `(car ,p))146 (defun termlist-lm (p) (term-monom (termlist-lt p)))147 (defun termlist-lc (p) (term-coeff (termlist-lt p)))148 149 (define-modify-macro scalar-mul (c) coeff-mul)150 151 (defun scalar-times-termlist (ring c p)152 "Multiply scalar C by a polynomial P. This function works153 even if there are divisors of 0."154 (mapcan155 #'(lambda (term)156 (let ((c1 (funcall (ring-mul ring) c (term-coeff term))))157 (unless (funcall (ring-zerop ring) c1)158 (list (make-term (term-monom term) c1)))))159 p))160 161 162 (defun term-mul (ring term1 term2)163 "Returns (LIST TERM) wheter TERM is the product of the terms TERM1 TERM2,164 or NIL when the product is 0. This definition takes care of divisors of 0165 in the coefficient ring."166 (let ((c (funcall (ring-mul ring) (term-coeff term1) (term-coeff term2))))167 (unless (funcall (ring-zerop ring) c)168 (list (make-term (monom-mul (term-monom term1) (term-monom term2)) c)))))169 170 (defun term-times-termlist (ring term f)171 (declare (type ring ring))172 (mapcan #'(lambda (term-f) (term-mul ring term term-f)) f))173 174 (defun termlist-times-term (ring f term)175 (mapcan #'(lambda (term-f) (term-mul ring term-f term)) f))176 177 (defun monom-times-term (m term)178 (make-term (monom-mul m (term-monom term)) (term-coeff term)))179 180 (defun monom-times-termlist (m f)181 (cond182 ((null f) nil)183 (t184 (mapcar #'(lambda (x) (monom-times-term m x)) f))))185 186 (defun termlist-uminus (ring f)187 (mapcar #'(lambda (x)188 (make-term (term-monom x) (funcall (ring-uminus ring) (term-coeff x))))189 f))190 191 (defun termlist-add (ring p q)192 (declare (type list p q))193 (do (r)194 ((cond195 ((endp p)196 (setf r (revappend r q)) t)197 ((endp q)198 (setf r (revappend r p)) t)199 (t200 (multiple-value-bind201 (lm-greater lm-equal)202 (monomial-order (termlist-lm p) (termlist-lm q))203 (cond204 (lm-equal205 (let ((s (funcall (ring-add ring) (termlist-lc p) (termlist-lc q))))206 (unless (funcall (ring-zerop ring) s) ;check for cancellation207 (setf r (cons (make-term (termlist-lm p) s) r)))208 (setf p (cdr p) q (cdr q))))209 (lm-greater210 (setf r (cons (car p) r)211 p (cdr p)))212 (t (setf r (cons (car q) r)213 q (cdr q)))))214 nil))215 r)))216 217 (defun termlist-sub (ring p q)218 (declare (type list p q))219 (do (r)220 ((cond221 ((endp p)222 (setf r (revappend r (termlist-uminus ring q)))223 t)224 ((endp q)225 (setf r (revappend r p))226 t)227 (t228 (multiple-value-bind229 (mgreater mequal)230 (monomial-order (termlist-lm p) (termlist-lm q))231 (cond232 (mequal233 (let ((s (funcall (ring-sub ring) (termlist-lc p) (termlist-lc q))))234 (unless (funcall (ring-zerop ring) s) ;check for cancellation235 (setf r (cons (make-term (termlist-lm p) s) r)))236 (setf p (cdr p) q (cdr q))))237 (mgreater238 (setf r (cons (car p) r)239 p (cdr p)))240 (t (setf r (cons (make-term (termlist-lm q) (funcall (ring-uminus ring) (termlist-lc q))) r)241 q (cdr q)))))242 nil))243 r)))244 245 ;; Multiplication of polynomials246 ;; Non-destructive version247 (defun termlist-mul (ring p q)248 (cond ((or (endp p) (endp q)) nil) ;p or q is 0 (represented by NIL)249 ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q250 ((endp (cdr p))251 (term-times-termlist ring (car p) q))252 ((endp (cdr q))253 (termlist-times-term ring p (car q)))254 (t255 (let ((head (term-mul ring (termlist-lt p) (termlist-lt q)))256 (tail (termlist-add ring (term-times-termlist ring (car p) (cdr q))257 (termlist-mul ring (cdr p) q))))258 (cond ((null head) tail)259 ((null tail) head)260 (t (nconc head tail)))))))261 262 (defun termlist-unit (ring dimension)263 (declare (fixnum dimension))264 (list (make-term (make-monom dimension :initial-element 0)265 (funcall (ring-unit ring)))))266 267 (defun termlist-expt (ring poly n &aux (dim (monom-dimension (termlist-lm poly))))268 (declare (type fixnum n dim))269 (cond270 ((minusp n) (error "termlist-expt: Negative exponent."))271 ((endp poly) (if (zerop n) (termlist-unit ring dim) nil))272 (t273 (do ((k 1 (ash k 1))274 (q poly (termlist-mul ring q q)) ;keep squaring275 (p (termlist-unit ring dim) (if (not (zerop (logand k n))) (termlist-mul ring p q) p)))276 ((> k n) p)277 (declare (fixnum k))))))278 279 280 281 282 283 284 136 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 285 137 ;; … … 289 141 (defmacro debug-cgb (&rest args) 290 142 `(when $poly_grobner_debug (format *terminal-io* ,@args))) 291 292 293 143 294 144
Note:
See TracChangeset
for help on using the changeset viewer.