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@ 3729

Last change on this file since 3729 was 3721, checked in by Marek Rychlik, 9 years ago

* empty log message *

File size: 18.4 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)
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 "UNIVERSAL-EZGCD"
49 "S-POLYNOMIAL"
50 "POLY-CONTENT"
51 "POLY-PRIMITIVE-PART"
52 "SATURATION-EXTENSION-1"
53 "UNIVERSAL-EXPT")
54 (:documentation "Implements polynomials. A polynomial is essentially
55a mapping of monomials of the same degree to coefficients. The
56momomials are ordered according to a monomial order."))
57
58(in-package :polynomial)
59
60(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
61
62(defclass poly ()
63 ((dimension :initform nil
64 :initarg :dimension
65 :accessor poly-dimension
66 :documentation "Shared dimension of all terms, the number of variables")
67 (termlist :initform nil :initarg :termlist :accessor poly-termlist
68 :documentation "List of terms.")
69 (order :initform #'lex> :initarg :order :accessor poly-term-order
70 :documentation "Monomial/term order."))
71 (:default-initargs :dimension nil :termlist nil :order #'lex>)
72 (:documentation "A polynomial with a list of terms TERMLIST, ordered
73according to term order ORDER, which defaults to LEX>."))
74
75(defmethod print-object ((self poly) stream)
76 (print-unreadable-object (self stream :type t :identity t)
77 (with-accessors ((dimension poly-dimension)
78 (termlist poly-termlist)
79 (order poly-term-order))
80 self
81 (format stream "DIMENSION=~A TERMLIST=~A ORDER=~A"
82 dimension termlist order))))
83
84(defgeneric change-term-order (self other)
85 (:documentation "Change term order of SELF to the term order of OTHER.")
86 (:method ((self poly) (other poly))
87 (unless (eq (poly-term-order self) (poly-term-order other))
88 (setf (poly-termlist self) (sort (poly-termlist self) (poly-term-order other))
89 (poly-term-order self) (poly-term-order other)))
90 self))
91
92(defgeneric poly-insert-term (self term)
93 (:documentation "Insert a term TERM into SELF before all other
94 terms. Order is not enforced.")
95 (:method ((self poly) (term term))
96 (cond ((null (poly-dimension self))
97 (setf (poly-dimension self) (monom-dimension term)))
98 (t (assert (= (poly-dimension self) (monom-dimension term)))))
99 (push term (poly-termlist self))
100 self))
101
102(defgeneric poly-append-term (self term)
103 (:documentation "Append a term TERM to SELF after all other terms. Order is not enforced.")
104 (:method ((self poly) (term term))
105 (cond ((null (poly-dimension self))
106 (setf (poly-dimension self) (monom-dimension term)))
107 (t (assert (= (poly-dimension self) (monom-dimension term)))))
108 (setf (cdr (last (poly-termlist self))) (list term))
109 self))
110
111(defun alist->poly (alist &aux (poly (make-instance 'poly)))
112 "It reads polynomial from an alist formatted as ( ... (exponents . coeff) ...).
113It can be used to enter simple polynomials by hand, e.g the polynomial
114in two variables, X and Y, given in standard notation as:
115
116 3*X^2*Y^3+2*Y+7
117
118can be entered as
119(ALIST->POLY '(((2 3) . 3) ((0 1) . 2) ((0 0) . 7))).
120
121NOTE: The primary use is for low-level debugging of the package."
122 (dolist (x alist poly)
123 (poly-insert-term poly (make-instance 'term :exponents (car x) :coeff (cdr x)))))
124
125(defmethod update-instance-for-different-class :after ((old monom) (new poly) &key)
126 "Converts OLD of class MONOM to a NEW of class POLY, by making it into a 1-element TERMLIST."
127 (reinitialize-instance new
128 :dimension (monom-dimension old)
129 :termlist (list (cons old 1))))
130
131(defmethod universal-equalp ((self poly) (other poly))
132 "Implements equality of polynomials."
133 (and (eql (poly-dimension self) (poly-dimension other))
134 (every #'universal-equalp (poly-termlist self) (poly-termlist other))
135 (eq (poly-term-order self) (poly-term-order other))))
136
137(defgeneric leading-term (object)
138 (:method ((self poly))
139 (car (poly-termlist self)))
140 (:documentation "The leading term of a polynomial, or NIL for zero polynomial."))
141
142(defgeneric second-leading-term (object)
143 (:method ((self poly))
144 (cadar (poly-termlist self)))
145 (:documentation "The second leading term of a polynomial, or NIL for a polynomial with at most one term."))
146
147(defgeneric leading-monomial (object)
148 (:method ((self poly))
149 (change-class (copy-instance (leading-term self)) 'monom))
150 (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
151
152(defgeneric second-leading-monomial (object)
153 (:method ((self poly))
154 (change-class (copy-instance (second-leading-term self)) 'monom))
155 (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
156
157(defgeneric leading-coefficient (object)
158 (:method ((self poly))
159 (term-coeff (leading-term self)))
160 (:documentation "The leading coefficient of a polynomial. It signals error for a zero polynomial."))
161
162(defgeneric second-leading-coefficient (object)
163 (:method ((self poly))
164 (term-coeff (second-leading-term self)))
165 (:documentation "The second leading coefficient of a polynomial. It
166 signals error for a polynomial with at most one term."))
167
168(defmethod universal-zerop ((self poly))
169 "Return T iff SELF is a zero polynomial."
170 (null (poly-termlist self)))
171
172(defgeneric poly-length (self)
173 (:documentation "Return the number of terms.")
174 (:method ((self poly))
175 (length (poly-termlist self))))
176
177(defgeneric scalar-multiply-by (self other)
178 (:documentation "Multiply vector SELF by a scalar OTHER.")
179 (:method ((self poly) other)
180 (mapc #'(lambda (term) (setf (term-coeff term) (multiply (term-coeff term) other)))
181 (poly-termlist self))
182 self))
183
184(defgeneric scalar-divide-by (self other)
185 (:documentation "Divide vector SELF by a scalar OTHER.")
186 (:method ((self poly) other)
187 (mapc #'(lambda (term) (setf (term-coeff term) (divide (term-coeff term) other)))
188 (poly-termlist self))
189 self))
190
191(defmethod multiply-by ((self poly) (other monom))
192 "Multiply a polynomial SELF by OTHER."
193 (mapc #'(lambda (term) (multiply-by term other))
194 (poly-termlist self))
195 self)
196
197(defmethod multiply-by ((self poly) (other term))
198 "Multiply a polynomial SELF by OTHER."
199 (mapc #'(lambda (term) (multiply-by term other))
200 (poly-termlist self))
201 self)
202
203(defmacro fast-add/subtract (p q order-fn add/subtract-fn uminus-fn)
204 "Return an expression which will efficiently adds/subtracts two
205polynomials, P and Q. The addition/subtraction of coefficients is
206performed by calling ADD/SUBTRACT-METHOD-NAME. If UMINUS-METHOD-NAME
207is supplied, it is used to negate the coefficients of Q which do not
208have a corresponding coefficient in P. The code implements an
209efficient algorithm to add two polynomials represented as sorted lists
210of terms. The code destroys both arguments, reusing the terms to build
211the result."
212 `(macrolet ((lc (x) `(term-coeff (car ,x))))
213 (do ((p ,p)
214 (q ,q)
215 r)
216 ((or (endp p) (endp q))
217 ;; NOTE: R contains the result in reverse order. Can it
218 ;; be more efficient to produce the terms in correct order?
219 (unless (endp q)
220 ;; Upon subtraction, we must change the sign of
221 ;; all coefficients in q
222 ,@(when uminus-fn
223 `((mapc #'(lambda (x) (setf x (funcall ,uminus-fn x))) q)))
224 (setf r (nreconc r q)))
225 r)
226 (multiple-value-bind
227 (greater-p equal-p)
228 (funcall ,order-fn (car p) (car q))
229 (cond
230 (greater-p
231 (rotatef (cdr p) r p)
232 )
233 (equal-p
234 (let ((s (funcall ,add/subtract-fn (lc p) (lc q))))
235 (cond
236 ((universal-zerop s)
237 (setf p (cdr p))
238 )
239 (t
240 (setf (lc p) s)
241 (rotatef (cdr p) r p))))
242 (setf q (cdr q))
243 )
244 (t
245 ;;Negate the term of Q if UMINUS provided, signallig
246 ;;that we are doing subtraction
247 ,(when uminus-fn
248 `(setf (lc q) (funcall ,uminus-fn (lc q))))
249 (rotatef (cdr q) r q)))))))
250
251
252(defgeneric add-to (self other)
253 (:documentation "Add OTHER to SELF.")
254 (:method ((self number) (other number))
255 (+ self other)))
256
257(defgeneric subtract-from (self other)
258 (:documentation "Subtract OTHER from SELF.")
259 (:method ((self number) (other number))
260 (- self other)))
261
262(defmacro def-add/subtract-method (add/subtract-method-name
263 uminus-method-name
264 &optional
265 (doc-string nil doc-string-supplied-p))
266 "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
267 `(defmethod ,add/subtract-method-name ((self poly) (other poly))
268 ,@(when doc-string-supplied-p `(,doc-string))
269 ;; Ensure orders are compatible
270 (change-term-order other self)
271 (setf (poly-termlist self) (fast-add/subtract
272 (poly-termlist self) (poly-termlist other)
273 (poly-term-order self)
274 #',add/subtract-method-name
275 ,(when uminus-method-name `(function ,uminus-method-name))))
276 self))
277
278(eval-when (:compile-toplevel :load-toplevel :execute)
279
280 (def-add/subtract-method add-to nil
281 "Adds to polynomial SELF another polynomial OTHER.
282This operation destructively modifies both polynomials.
283The result is stored in SELF. This implementation does
284no consing, entirely reusing the sells of SELF and OTHER.")
285
286 (def-add/subtract-method subtract-from unary-minus
287 "Subtracts from polynomial SELF another polynomial OTHER.
288This operation destructively modifies both polynomials.
289The result is stored in SELF. This implementation does
290no consing, entirely reusing the sells of SELF and OTHER.")
291 )
292
293(defmethod unary-minus ((self poly))
294 "Destructively modifies the coefficients of the polynomial SELF,
295by changing their sign."
296 (mapc #'unary-minus (poly-termlist self))
297 self)
298
299(defun add-termlists (p q order-fn)
300 "Destructively adds two termlists P and Q ordered according to ORDER-FN."
301 (fast-add/subtract p q order-fn #'add-to nil))
302
303(defmacro multiply-term-by-termlist-dropping-zeros (term termlist
304 &optional (reverse-arg-order-P nil))
305 "Multiplies term TERM by a list of term, TERMLIST.
306Takes into accound divisors of zero in the ring, by
307deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
308is T, change the order of arguments; this may be important
309if we extend the package to non-commutative rings."
310 `(mapcan #'(lambda (other-term)
311 (let ((prod (multiply
312 ,@(cond
313 (reverse-arg-order-p
314 `(other-term ,term))
315 (t
316 `(,term other-term))))))
317 (cond
318 ((universal-zerop prod) nil)
319 (t (list prod)))))
320 ,termlist))
321
322(defun multiply-termlists (p q order-fn)
323 "A version of polynomial multiplication, operating
324directly on termlists."
325 (cond
326 ((or (endp p) (endp q))
327 ;;p or q is 0 (represented by NIL)
328 nil)
329 ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
330 ((endp (cdr p))
331 (multiply-term-by-termlist-dropping-zeros (car p) q))
332 ((endp (cdr q))
333 (multiply-term-by-termlist-dropping-zeros (car q) p t))
334 (t
335 (cons (multiply (car p) (car q))
336 (add-termlists
337 (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
338 (multiply-termlists (cdr p) q order-fn)
339 order-fn)))))
340
341(defmethod multiply-by ((self poly) (other poly))
342 (change-term-order other self)
343 (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
344 (poly-termlist other)
345 (poly-term-order self)))
346 self)
347
348(defun add (object1 object2)
349 "Non-destructively add POLY1 by POLY2."
350 (add-to (copy-instance object1) (copy-instance object2)))
351
352(defun subtract (minuend &rest subtrahends)
353 "Non-destructively subtract MINUEND and SUBTRAHENDS."
354 (subtract-from (copy-instance minuend) (reduce #'add subtrahends)))
355
356(defmethod left-tensor-product-by ((self poly) (other monom))
357 (setf (poly-termlist self)
358 (mapcan #'(lambda (term)
359 (let ((prod (left-tensor-product-by term other)))
360 (cond
361 ((universal-zerop prod) nil)
362 (t (list prod)))))
363 (poly-termlist self)))
364 (incf (poly-dimension self) (monom-dimension other))
365 self)
366
367(defmethod right-tensor-product-by ((self poly) (other monom))
368 (setf (poly-termlist self)
369 (mapcan #'(lambda (term)
370 (let ((prod (right-tensor-product-by term other)))
371 (cond
372 ((universal-zerop prod) nil)
373 (t (list prod)))))
374 (poly-termlist self)))
375 (incf (poly-dimension self) (monom-dimension other))
376 self)
377
378
379(defun standard-extension (plist &aux (k (length plist)) (i 0))
380 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
381is a list of polynomials. Destructively modifies PLIST elements."
382 (mapc #'(lambda (poly)
383 (left-tensor-product-by
384 poly
385 (prog1
386 (make-monom-variable k i)
387 (incf i))))
388 plist))
389
390(defun standard-extension-1 (plist
391 &aux
392 (plist (standard-extension plist))
393 (nvars (poly-dimension (car plist))))
394 "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
395Firstly, new K variables U1, U2, ..., UK, are inserted into each
396polynomial. Subsequently, P1, P2, ..., PK are destructively modified
397tantamount to replacing PI with UI*PI-1. It assumes that all
398polynomials have the same dimension, and only the first polynomial
399is examined to determine this dimension."
400 ;; Implementation note: we use STANDARD-EXTENSION and then subtract
401 ;; 1 from each polynomial; since UI*PI has no constant term,
402 ;; we just need to append the constant term at the end
403 ;; of each termlist.
404 (flet ((subtract-1 (p)
405 (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
406 (setf plist (mapc #'subtract-1 plist)))
407 plist)
408
409
410(defun standard-sum (plist
411 &aux
412 (plist (standard-extension plist))
413 (nvars (poly-dimension (car plist))))
414 "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
415Firstly, new K variables, U1, U2, ..., UK, are inserted into each
416polynomial. Subsequently, P1, P2, ..., PK are destructively modified
417tantamount to replacing PI with UI*PI, and the resulting polynomials
418are added. Finally, 1 is subtracted. It should be noted that the term
419order is not modified, which is equivalent to using a lexicographic
420order on the first K variables."
421 (flet ((subtract-1 (p)
422 (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
423 (subtract-1
424 (make-instance
425 'poly
426 :termlist (apply #'nconc (mapcar #'poly-termlist plist))))))
427
428(defgeneric universal-ezgcd (x y)
429 (:documentation "Solves the diophantine system: X=C*X1, Y=C*X2,
430C=GCD(X,Y). It returns C, X1 and Y1. The result may be obtained by
431the Euclidean algorithm.")
432 (:method ((x integer) (y integer)
433 &aux (c (gcd x y)))
434 (values c (/ x c) (/ y c)))
435 )
436
437(defgeneric s-polynomial (object1 object2)
438 (:documentation "Yields the S-polynomial of OBJECT1 and OBJECT2.")
439 (:method ((f poly) (g poly))
440 (let* ((lcm (universal-lcm (leading-monomial f) (leading-monomial g)))
441 (mf (divide lcm (leading-monomial f)))
442 (mg (divide lcm (leading-monomial g))))
443 (multiple-value-bind (c cf cg)
444 (universal-ezgcd (leading-coefficient f) (leading-coefficient g))
445 (declare (ignore c))
446 (subtract
447 (multiply f (change-class mf 'term :coeff cg))
448 (multiply g (change-class mg 'term :coeff cf)))))))
449
450(defgeneric poly-content (object)
451 (:documentation "Greatest common divisor of the coefficients of the polynomial object OBJECT.")
452 (:method ((self poly))
453 (reduce #'universal-gcd
454 (mapcar #'term-coeff (rest (poly-termlist self)))
455 :initial-value (leading-coefficient self))))
456
457(defun poly-primitive-part (object)
458 "Divide polynomial OBJECT by gcd of its
459coefficients. Return the resulting polynomial."
460 (scalar-divide-by object (poly-content object)))
461
462(defun poly-insert-variables (self k)
463 (left-tensor-product-by self (make-instance 'monom :dimension k)))
464
465(defun saturation-extension (f plist &aux (k (length plist)))
466 "Calculate [F', U1*P1-1,U2*P2-1,...,UK*PK-1], where
467PLIST=[P1,P2,...,PK] and F' is F with variables U1,U2,...,UK inserted
468as first K variables. It destructively modifies F and PLIST."
469 (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
470 (standard-extension-1 plist)))
471
472(defun polysaturation-extension (f plist &aux (k (length plist)))
473 "Calculate [F', U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]
474and F' is F with variables U1,U2,...,UK inserted as first K
475variables. It destructively modifies F and PLIST."
476 (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
477 (list (standard-sum plist))))
478
479(defun saturation-extension-1 (f p)
480 "Given family of polynomials F and a polynomial P, calculate [F',
481U*P-1], where F' is F with variable inserted as the first variable. It
482destructively modifies F and P."
483 (polysaturation-extension f (list p)))
484
485(defmethod multiply-by ((object1 number) (object2 poly))
486 (scalar-multiply-by (copy-instance object2) object1))
487
488(defgeneric universal-expt (x y)
489 (:documentation "Raises X to power Y.")
490 (:method ((x number) (y integer)) (expt x y))
491 (:method ((x t) (y integer))
492 (declare (type fixnum y))
493 (cond
494 ((minusp y) (error "universal-expt: Negative exponent."))
495 ((universal-zerop x) (if (zerop y) 1))
496 (t
497 (do ((k 1 (ash k 1))
498 (q x (multiply q q)) ;keep squaring
499 (p 1 (if (not (zerop (logand k y))) (multiply p q) p)))
500 ((> k y) p)
501 (declare (fixnum k)))))))
Note: See TracBrowser for help on using the repository browser.