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

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

* empty log message *

File size: 16.5 KB
Line 
1;;----------------------------------------------------------------
2;; File: polynomial.lisp
3;;----------------------------------------------------------------
4;;
5;; Author: Marek Rychlik (rychlik@u.arizona.edu)
6;; Date: Thu Aug 27 09:41:24 2015
7;; Copying: (C) Marek Rychlik, 2010. All rights reserved.
8;;
9;;----------------------------------------------------------------
10;;; -*- Mode: Lisp -*-
11;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
12;;;
13;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
14;;;
15;;; This program is free software; you can redistribute it and/or modify
16;;; it under the terms of the GNU General Public License as published by
17;;; the Free Software Foundation; either version 2 of the License, or
18;;; (at your option) any later version.
19;;;
20;;; This program is distributed in the hope that it will be useful,
21;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
22;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23;;; GNU General Public License for more details.
24;;;
25;;; You should have received a copy of the GNU General Public License
26;;; along with this program; if not, write to the Free Software
27;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
28;;;
29;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
30
31(defpackage "POLYNOMIAL"
32 (:use :cl :utils :monom)
33 (:export "POLY"
34 "POLY-DIMENSION"
35 "POLY-TERMLIST"
36 "POLY-TERM-ORDER"
37 "POLY-INSERT-TERM"
38 "CHANGE-TERM-ORDER"
39 "STANDARD-EXTENSION"
40 "STANDARD-EXTENSION-1"
41 "STANDARD-SUM"
42 "SATURATION-EXTENSION"
43 "ALIST->POLY")
44 (:documentation "Implements polynomials. A polynomial is essentially
45a mapping of monomials of the same degree to coefficients. The
46momomials are ordered according to a monomial order."))
47
48(in-package :polynomial)
49
50(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
51
52(defclass poly ()
53 ((dimension :initform nil
54 :initarg :dimension
55 :accessor poly-dimension
56 :documentation "Shared dimension of all terms, the number of variables")
57 (termlist :initform nil :initarg :termlist :accessor poly-termlist
58 :documentation "List of terms. This is an association
59list mapping monomials to coefficients, ordered by this polynomial's
60monomial order.")
61 (order :initform #'lex> :initarg :order :accessor poly-term-order
62 :documentation "Monomial/term order."))
63 (:default-initargs :dimension nil :termlist nil :order #'lex>)
64 (:documentation "A polynomial with a list of terms TERMLIST, ordered
65according to term order ORDER, which defaults to LEX>."))
66
67(defmethod print-object ((self poly) stream)
68 (print-unreadable-object (self stream :type t :identity t)
69 (with-accessors ((dimension poly-dimension)
70 (termlist poly-termlist)
71 (order poly-term-order))
72 self
73 (format stream "DIMENSION=~A TERMLIST=~A ORDER=~A"
74 dimension termlist order))))
75
76(defgeneric change-term-order (self other)
77 (:documentation "Change term order of SELF to the term order of OTHER.")
78 (:method ((self poly) (other poly))
79 (unless (eq (poly-term-order self) (poly-term-order other))
80 (setf (poly-termlist self) (sort (poly-termlist self) (poly-term-order other))
81 (poly-term-order self) (poly-term-order other)))
82 self))
83
84(defgeneric poly-insert-term (self monom coeff)
85 (:documentation "Insert a term with monomial MONOM and coefficient COEFF
86before all other terms. Order is not enforced.")
87 (:method ((self poly) (monom monom) coeff)
88 (cond ((null (poly-dimension self))
89 (setf (poly-dimension self) (monom-dimension monom)))
90 (t (assert (= (poly-dimension self) (monom-dimension monom)))))
91 (push (cons monom coeff) (poly-termlist self))
92 self))
93
94(defgeneric poly-append-term (self monom coeff)
95 (:documentation "Append a term with monomial MONOM and coefficient COEFF
96after all other terms. Order is not enforced.")
97 (:method ((self poly) (monom monom) coeff)
98 (cond ((null (poly-dimension self))
99 (setf (poly-dimension self) (monom-dimension monom)))
100 (t (assert (= (poly-dimension self) (monom-dimension monom)))))
101 (setf (cdr (last (poly-termlist self))) (list (cons monom coeff)))
102 self))
103
104(defun alist->poly (alist &aux (poly (make-instance 'poly)))
105 "It reads polynomial from an alist formatted as ( ... (exponents . coeff) ...).
106It can be used to enter simple polynomials by hand, e.g the polynomial
107in two variables, X and Y, given in standard notation as:
108
109 3*X^2*Y^3+2*Y+7
110
111can be entered as
112(ALIST->POLY '(((2 3) . 3) ((0 1) . 2) ((0 0) . 7))).
113
114NOTE: The primary use is for low-level debugging of the package."
115 (dolist (x alist poly)
116 (poly-insert-term poly (make-instance 'monom :exponents (car x)) (cdr x))))
117
118(defmethod update-instance-for-different-class :after ((old monom) (new poly) &key)
119 "Converts OLD of class MONOM to a NEW of class POLY, by making it into a 1-element TERMLIST."
120 (reinitialize-instance new
121 :dimension (monom-dimension old)
122 :termlist (list (cons old 1))))
123
124(defgeneric poly-equalp (self other)
125 (:documentation "Implements equality of polynomials.")
126 (:method ((self poly) (other poly))
127 (and (eql (poly-dimension self) (poly-dimension other))
128 (every #'r-equalp (poly-termlist self) (poly-termlist other))
129 (eq (poly-term-order self) (poly-term-order other)))))
130
131;; Leading term
132(defgeneric leading-term (object)
133 (:method ((self poly))
134 (car (poly-termlist self)))
135 (:documentation "The leading term of a polynomial, or NIL for zero polynomial."))
136
137;; Second term
138(defgeneric second-leading-term (object)
139 (:method ((self poly))
140 (cadar (poly-termlist self)))
141 (:documentation "The second leading term of a polynomial, or NIL for a polynomial with at most one term."))
142
143;; Leading coefficient
144(defgeneric leading-coefficient (object)
145 (:method ((self poly))
146 (scalar-coeff (leading-term self)))
147 (:documentation "The leading coefficient of a polynomial. It signals error for a zero polynomial."))
148
149;; Second coefficient
150(defgeneric second-leading-coefficient (object)
151 (:method ((self poly))
152 (scalar-coeff (second-leading-term self)))
153 (:documentation "The second leading coefficient of a polynomial. It
154 signals error for a polynomial with at most one term."))
155
156;; Testing for a zero polynomial
157(defmethod r-zerop ((self poly))
158 (null (poly-termlist self)))
159
160;; The number of terms
161(defmethod r-length ((self poly))
162 (length (poly-termlist self)))
163
164(defmethod multiply-by ((self poly) (other monom))
165 (mapc #'(lambda (term) (multiply-by term other))
166 (poly-termlist self))
167 self)
168
169(defmethod multiply-by ((self poly) other)
170 (mapc #'(lambda (term) (multiply-by term other))
171 (poly-termlist self))
172 self)
173
174(defmacro fast-add/subtract (p q order-fn add/subtract-fn uminus-fn)
175 "Return an expression which will efficiently adds/subtracts two
176polynomials, P and Q. The addition/subtraction of coefficients is
177performed by calling ADD/SUBTRACT-METHOD-NAME. If UMINUS-METHOD-NAME
178is supplied, it is used to negate the coefficients of Q which do not
179have a corresponding coefficient in P. The code implements an
180efficient algorithm to add two polynomials represented as sorted lists
181of terms. The code destroys both arguments, reusing the terms to build
182the result."
183 `(macrolet ((lc (x) `(scalar-coeff (car ,x))))
184 (do ((p ,p)
185 (q ,q)
186 r)
187 ((or (endp p) (endp q))
188 ;; NOTE: R contains the result in reverse order. Can it
189 ;; be more efficient to produce the terms in correct order?
190 (unless (endp q)
191 ;; Upon subtraction, we must change the sign of
192 ;; all coefficients in q
193 ,@(when uminus-fn
194 `((mapc #'(lambda (x) (setf x (funcall ,uminus-fn x))) q)))
195 (setf r (nreconc r q)))
196 r)
197 (multiple-value-bind
198 (greater-p equal-p)
199 (funcall ,order-fn (car p) (car q))
200 (cond
201 (greater-p
202 (rotatef (cdr p) r p)
203 )
204 (equal-p
205 (let ((s (funcall ,add/subtract-fn (lc p) (lc q))))
206 (cond
207 ((r-zerop s)
208 (setf p (cdr p))
209 )
210 (t
211 (setf (lc p) s)
212 (rotatef (cdr p) r p))))
213 (setf q (cdr q))
214 )
215 (t
216 ;;Negate the term of Q if UMINUS provided, signallig
217 ;;that we are doing subtraction
218 ,(when uminus-fn
219 `(setf (lc q) (funcall ,uminus-fn (lc q))))
220 (rotatef (cdr q) r q)))))))
221
222
223(defmacro def-add/subtract-method (add/subtract-method-name
224 uminus-method-name
225 &optional
226 (doc-string nil doc-string-supplied-p))
227 "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
228 `(defmethod ,add/subtract-method-name ((self poly) (other poly))
229 ,@(when doc-string-supplied-p `(,doc-string))
230 ;; Ensure orders are compatible
231 (change-term-order other self)
232 (setf (poly-termlist self) (fast-add/subtract
233 (poly-termlist self) (poly-termlist other)
234 (poly-term-order self)
235 #',add/subtract-method-name
236 ,(when uminus-method-name `(function ,uminus-method-name))))
237 self))
238
239(eval-when (:compile-toplevel :load-toplevel :execute)
240
241 (def-add/subtract-method add-to nil
242 "Adds to polynomial SELF another polynomial OTHER.
243This operation destructively modifies both polynomials.
244The result is stored in SELF. This implementation does
245no consing, entirely reusing the sells of SELF and OTHER.")
246
247 (def-add/subtract-method subtract-from unary-minus
248 "Subtracts from polynomial SELF another polynomial OTHER.
249This operation destructively modifies both polynomials.
250The result is stored in SELF. This implementation does
251no consing, entirely reusing the sells of SELF and OTHER.")
252 )
253
254(defmethod unary-minus ((self poly))
255 "Destructively modifies the coefficients of the polynomial SELF,
256by changing their sign."
257 (mapc #'unary-minus (poly-termlist self))
258 self)
259
260(defun add-termlists (p q order-fn)
261 "Destructively adds two termlists P and Q ordered according to ORDER-FN."
262 (fast-add/subtract p q order-fn #'add-to nil))
263
264(defmacro multiply-term-by-termlist-dropping-zeros (term termlist
265 &optional (reverse-arg-order-P nil))
266 "Multiplies term TERM by a list of term, TERMLIST.
267Takes into accound divisors of zero in the ring, by
268deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
269is T, change the order of arguments; this may be important
270if we extend the package to non-commutative rings."
271 `(mapcan #'(lambda (other-term)
272 (let ((prod (r*
273 ,@(cond
274 (reverse-arg-order-p
275 `(other-term ,term))
276 (t
277 `(,term other-term))))))
278 (cond
279 ((r-zerop prod) nil)
280 (t (list prod)))))
281 ,termlist))
282
283(defun multiply-termlists (p q order-fn)
284 "A version of polynomial multiplication, operating
285directly on termlists."
286 (cond
287 ((or (endp p) (endp q))
288 ;;p or q is 0 (represented by NIL)
289 nil)
290 ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
291 ((endp (cdr p))
292 (multiply-term-by-termlist-dropping-zeros (car p) q))
293 ((endp (cdr q))
294 (multiply-term-by-termlist-dropping-zeros (car q) p t))
295 (t
296 (cons (r* (car p) (car q))
297 (add-termlists
298 (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
299 (multiply-termlists (cdr p) q order-fn)
300 order-fn)))))
301
302(defmethod multiply-by ((self poly) (other poly))
303 (change-term-order other self)
304 (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
305 (poly-termlist other)
306 (poly-term-order self)))
307 self)
308
309(defmethod r+ ((poly1 poly) poly2)
310 "Non-destructively add POLY1 by POLY2."
311 (add-to (copy-instance POLY1) (change-class (copy-instance POLY2) 'poly)))
312
313(defmethod r- ((minuend poly) &rest subtrahends)
314 "Non-destructively subtract MINUEND and SUBTRAHENDS."
315 (subtract-from (copy-instance minuend)
316 (change-class (reduce #'r+ subtrahends) 'poly)))
317
318(defmethod r+ ((poly1 monom) poly2)
319 "Non-destructively add POLY1 by POLY2."
320 (add-to (change-class (copy-instance poly1) 'poly)
321 (change-class (copy-instance poly2) 'poly)))
322
323(defmethod r- ((minuend monom) &rest subtrahends)
324 "Non-destructively subtract MINUEND and SUBTRAHENDS."
325 (subtract-from (change-class (copy-instance minuend) 'poly)
326 (change-class (reduce #'r+ subtrahends) 'poly)))
327
328(defmethod r* ((poly1 poly) (poly2 poly))
329 "Non-destructively multiply POLY1 by POLY2."
330 (multiply-by (copy-instance poly1) (copy-instance poly2)))
331
332(defmethod left-tensor-product-by ((self poly) (other monom))
333 (setf (poly-termlist self)
334 (mapcan #'(lambda (term)
335 (let ((prod (left-tensor-product-by term other)))
336 (cond
337 ((r-zerop prod) nil)
338 (t (list prod)))))
339 (poly-termlist self)))
340 (incf (poly-dimension self) (monom-dimension other))
341 self)
342
343(defmethod right-tensor-product-by ((self poly) (other monom))
344 (setf (poly-termlist self)
345 (mapcan #'(lambda (term)
346 (let ((prod (right-tensor-product-by term other)))
347 (cond
348 ((r-zerop prod) nil)
349 (t (list prod)))))
350 (poly-termlist self)))
351 (incf (poly-dimension self) (monom-dimension other))
352 self)
353
354
355(defun standard-extension (plist &aux (k (length plist)) (i 0))
356 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
357is a list of polynomials. Destructively modifies PLIST elements."
358 (mapc #'(lambda (poly)
359 (left-tensor-product-by
360 poly
361 (prog1
362 (make-monom-variable k i)
363 (incf i))))
364 plist))
365
366(defun standard-extension-1 (plist
367 &aux
368 (plist (standard-extension plist))
369 (nvars (poly-dimension (car plist))))
370 "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
371Firstly, new K variables U1, U2, ..., UK, are inserted into each
372polynomial. Subsequently, P1, P2, ..., PK are destructively modified
373tantamount to replacing PI with UI*PI-1. It assumes that all
374polynomials have the same dimension, and only the first polynomial
375is examined to determine this dimension."
376 ;; Implementation note: we use STANDARD-EXTENSION and then subtract
377 ;; 1 from each polynomial; since UI*PI has no constant term,
378 ;; we just need to append the constant term at the end
379 ;; of each termlist.
380 (flet ((subtract-1 (p)
381 (poly-append-term p (make-instance 'monom :dimension nvars) -1)))
382 (setf plist (mapc #'subtract-1 plist)))
383 plist)
384
385
386(defun standard-sum (plist
387 &aux
388 (plist (standard-extension plist))
389 (nvars (poly-dimension (car plist))))
390 "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
391Firstly, new K variables, U1, U2, ..., UK, are inserted into each
392polynomial. Subsequently, P1, P2, ..., PK are destructively modified
393tantamount to replacing PI with UI*PI, and the resulting polynomials
394are added. Finally, 1 is subtracted. It should be noted that the term
395order is not modified, which is equivalent to using a lexicographic
396order on the first K variables."
397 (flet ((subtract-1 (p)
398 (poly-append-term p (make-instance 'monom :dimension nvars) -1)))
399 (subtract-1
400 (make-instance
401 'poly
402 :termlist (apply #'nconc (mapcar #'poly-termlist plist))))))
403
404#|
405
406(defun saturation-extension-1 (ring f p)
407 "Calculate [F, U*P-1]. It destructively modifies F."
408 (declare (type ring ring))
409 (polysaturation-extension ring f (list p)))
410
411
412
413
414(defun spoly (ring-and-order f g
415 &aux
416 (ring (ro-ring ring-and-order)))
417 "It yields the S-polynomial of polynomials F and G."
418 (declare (type ring-and-order ring-and-order) (type poly f g))
419 (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
420 (mf (monom-div lcm (poly-lm f)))
421 (mg (monom-div lcm (poly-lm g))))
422 (declare (type monom mf mg))
423 (multiple-value-bind (c cf cg)
424 (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
425 (declare (ignore c))
426 (poly-sub
427 ring-and-order
428 (scalar-times-poly ring cg (monom-times-poly mf f))
429 (scalar-times-poly ring cf (monom-times-poly mg g))))))
430
431
432(defun poly-primitive-part (ring p)
433 "Divide polynomial P with integer coefficients by gcd of its
434coefficients and return the result."
435 (declare (type ring ring) (type poly p))
436 (if (poly-zerop p)
437 (values p 1)
438 (let ((c (poly-content ring p)))
439 (values (make-poly-from-termlist
440 (mapcar
441 #'(lambda (x)
442 (make-term :monom (term-monom x)
443 :coeff (funcall (ring-div ring) (term-coeff x) c)))
444 (poly-termlist p))
445 (poly-sugar p))
446 c))))
447
448(defun poly-content (ring p)
449 "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
450to compute the greatest common divisor."
451 (declare (type ring ring) (type poly p))
452 (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
453
454|#
Note: See TracBrowser for help on using the repository browser.