close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

source: branches/f4grobner/polynomial.lisp@ 4451

Last change on this file since 4451 was 4451, checked in by Marek Rychlik, 8 years ago

* empty log message *

File size: 25.9 KB
Line 
1;;----------------------------------------------------------------
2;;; -*- Mode: Lisp -*-
3;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4;;;
5;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
6;;;
7;;; This program is free software; you can redistribute it and/or modify
8;;; it under the terms of the GNU General Public License as published by
9;;; the Free Software Foundation; either version 2 of the License, or
10;;; (at your option) any later version.
11;;;
12;;; This program is distributed in the hope that it will be useful,
13;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
14;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15;;; GNU General Public License for more details.
16;;;
17;;; You should have received a copy of the GNU General Public License
18;;; along with this program; if not, write to the Free Software
19;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20;;;
21;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
22
23(defpackage "POLYNOMIAL"
24 (:use :cl :utils :monom :copy :ring)
25 (:export "POLY"
26 "POLY-DIMENSION"
27 "POLY-TERMLIST"
28 "POLY-TERM-ORDER"
29 "POLY-INSERT-TERM"
30 "SCALAR-MULTIPLY-BY"
31 "SCALAR-DIVIDE-BY"
32 "LEADING-TERM"
33 "LEADING-MONOMIAL"
34 "LEADING-COEFFICIENT"
35 "SECOND-LEADING-TERM"
36 "SECOND-LEADING-MONOMIAL"
37 "SECOND-LEADING-COEFFICIENT"
38 "ADD-TO"
39 "ADD"
40 "SUBTRACT-FROM"
41 "SUBTRACT"
42 "CHANGE-TERM-ORDER"
43 "STANDARD-EXTENSION"
44 "STANDARD-EXTENSION-1"
45 "STANDARD-SUM"
46 "SATURATION-EXTENSION"
47 "ALIST->POLY"
48 "POLY->ALIST"
49 "->INFIX"
50 "UNIVERSAL-EZGCD"
51 "S-POLYNOMIAL"
52 "POLY-CONTENT"
53 "POLY-PRIMITIVE-PART"
54 "SATURATION-EXTENSION-1"
55 "MAKE-POLY-VARIABLE"
56 "MAKE-POLY-CONSTANT"
57 "MAKE-ZERO-FOR"
58 "MAKE-UNIT-FOR"
59 "UNIVERSAL-EXPT"
60 "UNIVERSAL-EQUALP"
61 "UNIVERSAL-ZEROP"
62 "POLY-LENGTH"
63 "POLY-REVERSE"
64 "POLY-P"
65 "+LIST-MARKER+"
66 "POLY-EVAL"
67 "*COEFFICIENT-CLASS*")
68 (:documentation "Implements polynomials. A polynomial is essentially
69a mapping of monomials of the same degree to coefficients. The
70momomials are ordered according to a monomial order."))
71
72(in-package :polynomial)
73
74(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
75
76(defclass poly (ring)
77 ((dimension :initform nil
78 :initarg :dimension
79 :accessor poly-dimension
80 :documentation "Shared dimension of all terms, the number of variables")
81 (termlist :initform nil :initarg :termlist :accessor poly-termlist
82 :documentation "List of terms.")
83 (order :initform #'lex> :initarg :order :accessor poly-term-order
84 :documentation "Monomial/term order."))
85 (:default-initargs :dimension nil :termlist nil :order #'lex>)
86 (:documentation "A polynomial with a list of terms TERMLIST, ordered
87according to term order ORDER, which defaults to LEX>."))
88
89(defmethod print-object ((self poly) stream)
90 (print-unreadable-object (self stream :type t :identity t)
91 (with-accessors ((dimension poly-dimension)
92 (termlist poly-termlist)
93 (order poly-term-order))
94 self
95 (format stream "DIMENSION=~A TERMLIST=~A ORDER=~A"
96 dimension termlist order))))
97
98(defmethod copy-instance :around ((object poly) &rest initargs &key &allow-other-keys)
99 "Returns a deep copy of the polynomial POLY, by copying the TERMLIST and its terms."
100 (declare (ignore object initargs))
101 (let ((copy (call-next-method)))
102 (with-slots (termlist)
103 copy
104 (setf termlist (mapcar #'copy-instance termlist)))
105 copy))
106
107
108(defgeneric change-term-order (self other)
109 (:documentation "Change term order of SELF to the term order of OTHER.")
110 (:method ((self poly) (other poly))
111 (unless (eq (poly-term-order self) (poly-term-order other))
112 (setf (poly-termlist self) (sort (poly-termlist self) (poly-term-order other))
113 (poly-term-order self) (poly-term-order other)))
114 self))
115
116(defgeneric poly-insert-term (self term)
117 (:documentation "Insert a term TERM into SELF before all other
118terms. Order is not enforced.")
119 (:method ((self poly) (term term))
120 (cond ((null (poly-dimension self))
121 (setf (poly-dimension self) (monom-dimension term)))
122 (t (assert (= (poly-dimension self) (monom-dimension term)))))
123 (push term (poly-termlist self))
124 self))
125
126(defgeneric poly-append-term (self term)
127 (:documentation "Append a term TERM to SELF after all other terms. Order is not enforced.")
128 (:method ((self poly) (term term))
129 (cond ((null (poly-dimension self))
130 (setf (poly-dimension self) (monom-dimension term)))
131 (t (assert (= (poly-dimension self) (monom-dimension term)))))
132 (setf (cdr (last (poly-termlist self))) (list term))
133 self))
134
135(defun alist->poly (alist &aux (poly (make-instance 'poly)))
136 "It reads polynomial from an alist formatted as ( ... (exponents . coeff) ...).
137It can be used to enter simple polynomials by hand, e.g the polynomial
138in two variables, X and Y, given in standard notation as:
139
140 3*X^2*Y^3+2*Y+7
141
142can be entered as
143(ALIST->POLY '(((0 0) . 7) ((0 1) . 2) ((2 3) . 3) )). NOTE: the
144terms are entered in the increasing order.
145
146NOTE: The primary use is for low-level debugging of the package."
147 (dolist (x alist poly)
148 (poly-insert-term poly (make-instance 'term :exponents (car x) :coeff (cdr x)))))
149
150(defun poly->alist (p)
151 "Convert a polynomial P to an association list. Thus, the format of the
152returned value is ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
153MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
154corresponding coefficient in the ring."
155 (cond
156 ((poly-p p)
157 (mapcar #'->list (poly-termlist p)))
158 ((and (consp p) (eq (car p) :[))
159 (cons :[ (mapcar #'poly->alist (cdr p))))))
160
161
162
163(defmethod update-instance-for-different-class :after ((old term) (new poly) &key)
164 "Converts OLD of class TERM to a NEW of class POLY, by making it into a 1-element TERMLIST."
165 (reinitialize-instance new
166 :dimension (monom-dimension old)
167 :termlist (list old)))
168
169(defmethod update-instance-for-different-class :after ((old monom) (new poly) &key)
170 "Converts OLD of class MONOM to a NEW of class POLY, by making it into a 1-element TERMLIST."
171 (reinitialize-instance new
172 :dimension (monom-dimension old)
173 :termlist (list (change-class old 'term))))
174
175(defmethod universal-equalp ((self poly) (other poly))
176 "Implements equality of polynomials."
177 (and (eql (poly-dimension self) (poly-dimension other))
178 (every #'universal-equalp (poly-termlist self) (poly-termlist other))
179 (eq (poly-term-order self) (poly-term-order other))))
180
181(defgeneric leading-term (object)
182 (:method ((self poly))
183 (car (poly-termlist self)))
184 (:documentation "The leading term of a polynomial, or NIL for zero polynomial."))
185
186(defgeneric second-leading-term (object)
187 (:method ((self poly))
188 (cadar (poly-termlist self)))
189 (:documentation "The second leading term of a polynomial, or NIL for a polynomial with at most one term."))
190
191(defgeneric leading-monomial (object)
192 (:method ((self poly))
193 (change-class (copy-instance (leading-term self)) 'monom))
194 (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
195
196(defgeneric second-leading-monomial (object)
197 (:method ((self poly))
198 (change-class (copy-instance (second-leading-term self)) 'monom))
199 (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
200
201(defgeneric leading-coefficient (object)
202 (:method ((self poly))
203 (term-coeff (leading-term self)))
204 (:documentation "The leading coefficient of a polynomial. It signals error for a zero polynomial."))
205
206(defgeneric second-leading-coefficient (object)
207 (:method ((self poly))
208 (term-coeff (second-leading-term self)))
209 (:documentation "The second leading coefficient of a polynomial. It
210 signals error for a polynomial with at most one term."))
211
212(defmethod universal-zerop ((self poly))
213 "Return T iff SELF is a zero polynomial."
214 (null (poly-termlist self)))
215
216(defgeneric poly-length (self)
217 (:documentation "Return the number of terms.")
218 (:method ((self poly))
219 (length (poly-termlist self))))
220
221(defgeneric scalar-multiply-by (self other)
222 (:documentation "Multiply vector SELF by a scalar OTHER.")
223 (:method ((self poly) other)
224 (mapc #'(lambda (term) (setf (term-coeff term) (multiply-by (term-coeff term) other)))
225 (poly-termlist self))
226 self))
227
228(defgeneric scalar-divide-by (self other)
229 (:documentation "Divide vector SELF by a scalar OTHER.")
230 (:method ((self poly) other)
231 (mapc #'(lambda (term) (setf (term-coeff term) (divide-by (term-coeff term) other)))
232 (poly-termlist self))
233 self))
234
235(defmethod unary-inverse :before ((self poly))
236 "Checks invertibility of a polynomial SELF. To be invertable, the
237polynomial must be an invertible, constant polynomial."
238 (with-slots (termlist)
239 self
240 (assert (and (= (length termlist) 1) (zerop (total-degree (car termlist))))
241 nil
242 "To be invertible, the polynomial must have 1 term of total degree 0.")))
243
244(defmethod unary-inverse ((self poly))
245 "Returns the unary inverse of a polynomial SELF."
246 (with-slots (termlist)
247 self
248 (setf (car termlist) (unary-inverse (car termlist)))
249 self))
250
251(defmethod multiply-by ((self poly) (other monom))
252 "Multiply a polynomial SELF by OTHER."
253 (mapc #'(lambda (term) (multiply-by term other))
254 (poly-termlist self))
255 self)
256
257(defmethod multiply-by ((self poly) (other term))
258 "Multiply a polynomial SELF by OTHER."
259 (mapc #'(lambda (term) (multiply-by term other))
260 (poly-termlist self))
261 self)
262
263#|
264(defmacro fast-add/subtract (p q order-fn add/subtract-fn uminus-fn)
265 "Return an expression which will efficiently adds/subtracts two
266polynomials, P and Q. The addition/subtraction of coefficients is
267performed by calling ADD/SUBTRACT-FN. If UMINUS-FN is supplied, it is
268used to negate the coefficients of Q which do not have a corresponding
269coefficient in P. The code implements an efficient algorithm to add
270two polynomials represented as sorted lists of terms. The code
271destroys both arguments, reusing the terms to build the result."
272 `(macrolet ((lc (x) `(term-coeff (car ,x))))
273 (do ((p ,p)
274 (q ,q)
275 r)
276 ((or (endp p) (endp q))
277 ;; NOTE: R contains the result in reverse order. Can it
278 ;; be more efficient to produce the terms in correct order?
279 (unless (endp q)
280 ;; Upon subtraction, we must change the sign of
281 ;; all coefficients in q
282 ,@(when uminus-fn
283 `((mapc #'(lambda (x) (setf x (funcall ,uminus-fn x))) q)))
284 (setf r (nreconc r q)))
285 (unless (endp p)
286 (setf r (nreconc r p)))
287 r)
288 (multiple-value-bind
289 (greater-p equal-p)
290 (funcall ,order-fn (car p) (car q))
291 (cond
292 (greater-p
293 (rotatef (cdr p) r p)
294 )
295 (equal-p
296 (let ((s (funcall ,add/subtract-fn (lc p) (lc q))))
297 (cond
298 ((universal-zerop s)
299 (setf p (cdr p))
300 )
301 (t
302 (setf (lc p) s)
303 (rotatef (cdr p) r p))))
304 (setf q (cdr q))
305 )
306 (t
307 ;;Negate the term of Q if UMINUS provided, signallig
308 ;;that we are doing subtraction
309 ,(when uminus-fn
310 `(setf (lc q) (funcall ,uminus-fn (lc q))))
311 (rotatef (cdr q) r q))))
312 ;;(format t "P:~A~%" p)
313 ;;(format t "Q:~A~%" q)
314 ;;(format t "R:~A~%" r)
315 )))
316|#
317
318
319
320
321#|
322(defun fast-add (p q order-fn add-fn)
323 "Add two polynomials, P and Q, represented as lists of terms.
324The operation is destructive to both polynomials, as the terms
325of both lists are combined into the result. The operation does not
326create any new instance of TERM."
327 (macrolet ((lc (x) `(term-coeff (car ,x))))
328 (do (r)
329 ((or (endp p) (endp q))
330 ;; NOTE: R contains the result in reverse order. Can it
331 ;; be more efficient to produce the terms in correct order?
332 (unless (endp q)
333 (setf r (nreconc r q)))
334 (unless (endp p)
335 (setf r (nreconc r p)))
336 r)
337 (multiple-value-bind
338 (greater-p equal-p)
339 (funcall order-fn (car p) (car q))
340 (cond
341 (greater-p
342 (rotatef (cdr p) r p)
343 )
344 (equal-p
345 (let ((s (funcall add-fn (lc p) (lc q))))
346 (cond
347 ((universal-zerop s)
348 (setf p (cdr p))
349 )
350 (t
351 (setf (lc p) s)
352 (rotatef (cdr p) r p))))
353 (setf q (cdr q))
354 )
355 (t
356 (rotatef (cdr q) r q)))))))
357|#
358
359;; Getter/setter of leading coefficient
360(defun lc (x) (term-coeff (car x)))
361(defun (setf lc) (new-value x) (setf (term-coeff (car x)) new-value))
362
363
364(defun slow-add (p q order-fn add-fn)
365 (cond
366 ((endp p) q)
367 ((endp q) p)
368 (t
369 (multiple-value-bind
370 (greater-p equal-p)
371 (funcall order-fn (car p) (car q))
372 (cond
373 (greater-p ; (> (car p) (car q))
374 (cons (car p) (slow-add (cdr p) q order-fn add-fn))
375 )
376 (equal-p ; (= (car p)) (car q))
377 (let ((s (funcall add-fn (lc p) (lc q))))
378 (cond
379 ((universal-zerop s)
380 (slow-add (cdr p) (cdr q) order-fn add-fn))
381 (t
382 ;; Adjust the lc of p
383 (setf (lc p) s)
384 (cons (car p) (slow-add (cdr p) (cdr q) order-fn add-fn))
385 ))))
386 (t ;(< (car p) (car q))
387 (cons (car q) (slow-add p (cdr q) order-fn add-fn))
388 ))))))
389
390
391(defun fast-and-risky-add (p q order-fn add-fn &aux result result-last)
392 (when (and p q (eq p q)) (warn "FAST-AND-RISKY-ADD: ~S is EQ to ~S" p q))
393 (flet ((add-to-result (x)
394 (assert (consp x))
395 (setf (cdr x) nil)
396 (if (endp result)
397 (setf result x
398 result-last x)
399 (setf (cdr result-last) x
400 result-last (cdr result-last)))))
401 (loop
402 (cond
403 ((endp p) (unless (endp q) (add-to-result q)) (return result))
404 ((endp q) (unless (endp p) (add-to-result p)) (return result))
405 (t
406 (multiple-value-bind
407 (greater-p equal-p)
408 (funcall order-fn (car p) (car q))
409 (cond
410 (greater-p ; (> (car p) (car q))
411 (let ((tmp (cdr p)))
412 (add-to-result p)
413 (setf p tmp)))
414 (equal-p ; (= (car p)) (car q))
415 (let ((s (funcall add-fn (lc p) (lc q))))
416 (cond
417 ((universal-zerop s)
418 ;; Terms cancel, discard both
419 (setf p (cdr p)
420 q (cdr q)))
421 (t
422 ;; Terms do not cancel, store the
423 ;; sum of coefficients in (lc p)
424 (setf (lc p) s)
425 (let ((tmp (cdr p)))
426 (add-to-result p)
427 (setf p tmp
428 q (cdr q)))))))
429 (t ;(< (car p) (car q))
430 (let ((tmp (cdr q)))
431 (add-to-result q)
432 (setf q tmp))
433 ))))))))
434
435(defun fast-add (p q order-fn add-fn)
436 "This version calls SLOW-ADD and is bullet-proof."
437 ;;(slow-add p q order-fn add-fn)
438 (fast-and-risky-add p q order-fn add-fn)
439 )
440
441#|
442;; NOTE: The stuff below works, but may not be worth the trouble.
443
444(defmacro def-add/subtract-method (add/subtract-method-name
445 uminus-method-name
446 &optional
447 (doc-string nil doc-string-supplied-p))
448 "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
449 `(defmethod ,add/subtract-method-name ((self poly) (other poly))
450 ,@(when doc-string-supplied-p `(,doc-string))
451 ;; Ensure orders are compatible
452 (change-term-order other self)
453 (setf (poly-termlist self) (fast-add/subtract
454 (poly-termlist self) (poly-termlist other)
455 (poly-term-order self)
456 #',add/subtract-method-name
457 ,(when uminus-method-name `(function ,uminus-method-name))))
458 self))
459
460(eval-when (:load-toplevel :execute)
461
462 (def-add/subtract-method add-to nil
463 "Adds to polynomial SELF another polynomial OTHER.
464This operation destructively modifies both polynomials.
465The result is stored in SELF. This implementation does
466no consing, entirely reusing the sells of SELF and OTHER.")
467
468 (def-add/subtract-method subtract-from unary-minus
469 "Subtracts from polynomial SELF another polynomial OTHER.
470This operation destructively modifies both polynomials.
471The result is stored in SELF. This implementation does
472no consing, entirely reusing the sells of SELF and OTHER.")
473 )
474
475|#
476
477(defmethod unary-minus ((self poly))
478 "Destructively modifies the coefficients of the polynomial SELF,
479by changing their sign."
480 (mapc #'unary-minus (poly-termlist self))
481 self)
482
483(defun add-termlists (p q order-fn)
484 "Destructively adds two termlists P and Q ordered according to ORDER-FN."
485 (fast-add p q order-fn #'add-to))
486
487(defun subtract-termlists (p q order-fn)
488 "Destructively subtracts two termlists P and Q ordered according to ORDER-FN."
489 (setf q (mapc #'unary-minus q))
490 (add-termlists p q order-fn))
491
492(defmethod add-to ((self poly) (other poly) &aux (other-copy (copy-instance other)))
493 "Adds to polynomial SELF another polynomial OTHER.
494This operation destructively modifies both polynomials.
495The result is stored in SELF. This implementation does
496no consing, entirely reusing the sells of SELF and OTHER."
497 (change-term-order other-copy self)
498 (setf (poly-termlist self) (add-termlists
499 (poly-termlist self) (poly-termlist other-copy)
500 (poly-term-order self)))
501 self)
502
503
504(defmethod subtract-from ((self poly) (other poly))
505 "Subtracts from polynomial SELF another polynomial OTHER.
506This operation destructively modifies both polynomials.
507The result is stored in SELF. This implementation does
508no consing, entirely reusing the sells of SELF and OTHER."
509 (change-term-order other self)
510 (setf (poly-termlist self) (subtract-termlists
511 (poly-termlist self) (poly-termlist other)
512 (poly-term-order self)))
513 self)
514
515
516(defmethod add-to ((self poly) (other term) &aux (other-copy (copy-instance other)))
517 "Adds to a polynomial SELF a term OTHER. The term OTHER is not
518modified."
519 (add-to self (change-class other-copy 'poly)))
520
521(defmethod subtract-from ((self poly) (other term) &aux (other-copy (copy-instance other)))
522 "Subtracts from a polynomial SELF a term OTHER. The term OTHER is not
523modified."
524 (subtract-from self (change-class other-copy 'poly)))
525
526
527(defmacro multiply-term-by-termlist-dropping-zeros (term termlist
528 &optional (reverse-arg-order-P nil))
529 "Multiplies term TERM by a list of term, TERMLIST.
530Takes into accound divisors of zero in the ring, by
531deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
532is T, change the order of arguments; this may be important
533if we extend the package to non-commutative rings."
534 `(mapcan #'(lambda (other-term)
535 (let ((prod (multiply
536 ,@(cond
537 (reverse-arg-order-p
538 `(other-term ,term))
539 (t
540 `(,term other-term))))))
541 (cond
542 ((universal-zerop prod) nil)
543 (t (list prod)))))
544 ,termlist))
545
546(defun multiply-termlists (p q order-fn)
547 "A version of polynomial multiplication, operating
548directly on termlists."
549 (cond
550 ((or (endp p) (endp q))
551 ;;p or q is 0 (represented by NIL)
552 nil)
553 ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
554 ((endp (cdr p))
555 (multiply-term-by-termlist-dropping-zeros (car p) q))
556 ((endp (cdr q))
557 (multiply-term-by-termlist-dropping-zeros (car q) p t))
558 (t
559 (cons (multiply (car p) (car q))
560 (add-termlists
561 (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
562 (multiply-termlists (cdr p) q order-fn)
563 order-fn)))))
564
565(defmethod multiply-by ((self poly) (other poly) &aux (other-copy (copy-instance other)))
566 (change-term-order other-copy self)
567 (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
568 (poly-termlist other-copy)
569 (poly-term-order self)))
570 self)
571
572(defmethod left-tensor-product-by ((self poly) (other monom))
573 (setf (poly-termlist self)
574 (mapcan #'(lambda (term)
575 (let ((prod (left-tensor-product-by term other)))
576 (cond
577 ((universal-zerop prod) nil)
578 (t (list prod)))))
579 (poly-termlist self)))
580 (incf (poly-dimension self) (monom-dimension other))
581 self)
582
583(defmethod right-tensor-product-by ((self poly) (other monom))
584 (setf (poly-termlist self)
585 (mapcan #'(lambda (term)
586 (let ((prod (right-tensor-product-by term other)))
587 (cond
588 ((universal-zerop prod) nil)
589 (t (list prod)))))
590 (poly-termlist self)))
591 (incf (poly-dimension self) (monom-dimension other))
592 self)
593
594
595(defun standard-extension (plist &aux (k (length plist)) (i 0))
596 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
597is a list of polynomials. Destructively modifies PLIST elements."
598 (mapc #'(lambda (poly)
599 (left-tensor-product-by
600 poly
601 (prog1
602 (make-monom-variable k i)
603 (incf i))))
604 plist))
605
606(defun standard-extension-1 (plist
607 &aux
608 (plist (standard-extension plist))
609 (nvars (poly-dimension (car plist))))
610 "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
611Firstly, new K variables U1, U2, ..., UK, are inserted into each
612polynomial. Subsequently, P1, P2, ..., PK are destructively modified
613tantamount to replacing PI with UI*PI-1. It assumes that all
614polynomials have the same dimension, and only the first polynomial
615is examined to determine this dimension."
616 ;; Implementation note: we use STANDARD-EXTENSION and then subtract
617 ;; 1 from each polynomial; since UI*PI has no constant term,
618 ;; we just need to append the constant term at the end
619 ;; of each termlist.
620 (flet ((subtract-1 (p)
621 (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
622 (setf plist (mapc #'subtract-1 plist)))
623 plist)
624
625
626(defun standard-sum (plist
627 &aux
628 (plist (standard-extension plist))
629 (nvars (poly-dimension (car plist))))
630 "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
631Firstly, new K variables, U1, U2, ..., UK, are inserted into each
632polynomial. Subsequently, P1, P2, ..., PK are destructively modified
633tantamount to replacing PI with UI*PI, and the resulting polynomials
634are added. Finally, 1 is subtracted. It should be noted that the term
635order is not modified, which is equivalent to using a lexicographic
636order on the first K variables."
637 (flet ((subtract-1 (p)
638 (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
639 (subtract-1
640 (make-instance
641 'poly
642 :termlist (apply #'nconc (mapcar #'poly-termlist plist))))))
643
644(defgeneric s-polynomial (object1 object2)
645 (:documentation "Yields the S-polynomial of OBJECT1 and OBJECT2.")
646 (:method ((f poly) (g poly))
647 (let* ((lcm (universal-lcm (leading-monomial f) (leading-monomial g)))
648 (mf (divide lcm (leading-monomial f)))
649 (mg (divide lcm (leading-monomial g))))
650 (multiple-value-bind (c cf cg)
651 (universal-ezgcd (leading-coefficient f) (leading-coefficient g))
652 (declare (ignore c))
653 (subtract
654 (multiply f mf cg)
655 (multiply g mg cf))))))
656
657(defgeneric poly-content (object)
658 (:documentation "Greatest common divisor of the coefficients of the polynomial object OBJECT.")
659 (:method ((self poly))
660 (reduce #'universal-gcd
661 (mapcar #'term-coeff (rest (poly-termlist self)))
662 :initial-value (leading-coefficient self))))
663
664(defun poly-primitive-part (self)
665 "Divide polynomial SELF by gcd of its
666coefficients. Return the resulting polynomial."
667 (scalar-divide-by self (poly-content self)))
668
669(defun poly-insert-variables (self k)
670 (left-tensor-product-by self (make-instance 'monom :dimension k)))
671
672(defun saturation-extension (f plist &aux (k (length plist)))
673 "Calculate [F', U1*P1-1,U2*P2-1,...,UK*PK-1], where
674PLIST=[P1,P2,...,PK] and F' is F with variables U1,U2,...,UK inserted
675as first K variables. It destructively modifies F and PLIST."
676 (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
677 (standard-extension-1 plist)))
678
679(defun polysaturation-extension (f plist &aux (k (length plist)))
680 "Calculate [F', U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]
681and F' is F with variables U1,U2,...,UK inserted as first K
682variables. It destructively modifies F and PLIST."
683 (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
684 (list (standard-sum plist))))
685
686(defun saturation-extension-1 (f p)
687 "Given family of polynomials F and a polynomial P, calculate [F',
688U*P-1], where F' is F with variable inserted as the first variable. It
689destructively modifies F and P."
690 (polysaturation-extension f (list p)))
691
692(defmethod multiply-by ((self poly) (other ring))
693 (scalar-multiply-by self other))
694
695(defun make-poly-variable (nvars pos &optional (power 1))
696 (change-class (make-monom-variable nvars pos power) 'poly))
697
698(defun make-poly-constant (nvars coeff)
699 (change-class (make-term-constant nvars coeff) 'poly))
700
701(defgeneric universal-expt (x y)
702 (:documentation "Raises X to power Y.")
703 (:method ((x number) (y integer)) (expt x y))
704 (:method ((x t) (y integer))
705 (declare (type fixnum y))
706 (cond
707 ((minusp y) (error "universal-expt: Negative exponent."))
708 ((universal-zerop x) (if (zerop y) 1))
709 (t
710 (do ((k 1 (ash k 1))
711 (q x (multiply q q)) ;keep squaring
712 (p (make-unit-for x) (if (not (zerop (logand k y))) (multiply p q) p)))
713 ((> k y) p)
714 (declare (fixnum k)))))))
715
716(defgeneric poly-p (object)
717 (:documentation "Checks if an object is a polynomial.")
718 (:method ((self poly)) t)
719 (:method ((self t)) nil))
720
721(defmethod ->sexp :before ((self poly) &optional vars)
722 "Ensures that the number of variables in VARS maches the polynomial dimension of the
723polynomial SELF."
724 (with-slots (dimension)
725 self
726 (assert (= (length vars) dimension)
727 nil
728 "Number of variables ~S does not match the dimension ~S"
729 vars dimension)))
730
731(defmethod ->sexp ((self poly) &optional vars)
732 "Converts a polynomial SELF to a sexp."
733 (let ((m (mapcar #'(lambda (trm) (->sexp trm vars))
734 (poly-termlist self))))
735 (cond ((endp m) 0)
736 ((endp (cdr m)) (car m))
737 (t (cons '+ m)))))
738
739(defconstant +list-marker+ :[
740 "A sexp with this head is considered a list of polynomials.")
741
742(defmethod ->sexp ((self cons) &optional vars)
743 (assert (eql (car self) +list-marker+))
744 (cons +list-marker+ (mapcar #'(lambda (p) (->sexp p vars)) (cdr self))))
745
746(defmethod make-zero-for ((self poly))
747 (make-instance 'poly :dimension (poly-dimension self)))
748
749(defmethod make-unit-for ((self poly))
750 (make-poly-constant (poly-dimension self) 1))
751
752(defgeneric poly-reverse (self)
753 (:documentation "Reverse the order of terms in a polynomial SELF.")
754 (:method ((self poly))
755 (with-slots (termlist)
756 self
757 (setf termlist (nreverse termlist)))
758 self))
759
760
761
Note: See TracBrowser for help on using the repository browser.