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

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

* empty log message *

File size: 20.3 KB
Line 
1;;; -*- Mode: Lisp -*-
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3;;;
4;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
5;;;
6;;; This program is free software; you can redistribute it and/or modify
7;;; it under the terms of the GNU General Public License as published by
8;;; the Free Software Foundation; either version 2 of the License, or
9;;; (at your option) any later version.
10;;;
11;;; This program is distributed in the hope that it will be useful,
12;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
13;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14;;; GNU General Public License for more details.
15;;;
16;;; You should have received a copy of the GNU General Public License
17;;; along with this program; if not, write to the Free Software
18;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19;;;
20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
21
22(defpackage "POLYNOMIAL"
23 (:use :cl :utils :ring :monom :order :term #| :infix |# )
24 (:export "POLY"
25 "POLY-TERMLIST"
26 "POLY-TERM-ORDER"
27 "CHANGE-TERM-ORDER"
28 "STANDARD-EXTENSION"
29 "STANDARD-EXTENSION-1"
30 "SATURATION-EXTENSION"
31 "ALIST->POLY")
32 (:documentation "Implements polynomials"))
33
34(in-package :polynomial)
35
36(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
37
38(defclass poly ()
39 ((termlist :initarg :termlist :accessor poly-termlist
40 :documentation "List of terms.")
41 (order :initarg :order :accessor poly-term-order
42 :documentation "Monomial/term order."))
43 (:default-initargs :termlist nil :order #'lex>)
44 (:documentation "A polynomial with a list of terms TERMLIST, ordered
45according to term order ORDER, which defaults to LEX>."))
46
47(defmethod print-object ((self poly) stream)
48 (format stream "#<POLY TERMLIST=~A ORDER=~A>"
49 (poly-termlist self)
50 (poly-term-order self)))
51
52(defgeneric change-term-order (self other)
53 (:documentation "Change term order of SELF to the term order of OTHER.")
54 (:method ((self poly) (other poly))
55 (unless (eq (poly-term-order self) (poly-term-order other))
56 (setf (poly-termlist self) (sort (poly-termlist self) (poly-term-order other))
57 (poly-term-order self) (poly-term-order other)))
58 self))
59
60(defun alist->poly (alist &aux (poly (make-instance 'poly)))
61 "It reads polynomial from an alist formatted as ( ... (exponents . coeff) ...)."
62 (dolist (x alist poly)
63 (insert-item poly (make-instance 'term :exponents (car x) :coeff (cdr x)))))
64
65
66(defmethod r-equalp ((self poly) (other poly))
67 "POLY instances are R-EQUALP if they have the same
68order and if all terms are R-EQUALP."
69 (and (every #'r-equalp (poly-termlist self) (poly-termlist other))
70 (eq (poly-term-order self) (poly-term-order other))))
71
72(defmethod insert-item ((self poly) (item term))
73 (push item (poly-termlist self))
74 self)
75
76(defmethod append-item ((self poly) (item term))
77 (setf (cdr (last (poly-termlist self))) (list item))
78 self)
79
80;; Leading term
81(defgeneric leading-term (object)
82 (:method ((self poly))
83 (car (poly-termlist self)))
84 (:documentation "The leading term of a polynomial, or NIL for zero polynomial."))
85
86;; Second term
87(defgeneric second-leading-term (object)
88 (:method ((self poly))
89 (cadar (poly-termlist self)))
90 (:documentation "The second leading term of a polynomial, or NIL for a polynomial with at most one term."))
91
92;; Leading coefficient
93(defgeneric leading-coefficient (object)
94 (:method ((self poly))
95 (r-coeff (leading-term self)))
96 (:documentation "The leading coefficient of a polynomial. It signals error for a zero polynomial."))
97
98;; Second coefficient
99(defgeneric second-leading-coefficient (object)
100 (:method ((self poly))
101 (r-coeff (second-leading-term self)))
102 (:documentation "The second leading coefficient of a polynomial. It
103 signals error for a polynomial with at most one term."))
104
105;; Testing for a zero polynomial
106(defmethod r-zerop ((self poly))
107 (null (poly-termlist self)))
108
109;; The number of terms
110(defmethod r-length ((self poly))
111 (length (poly-termlist self)))
112
113(defmethod multiply-by ((self poly) (other monom))
114 (mapc #'(lambda (term) (multiply-by term other))
115 (poly-termlist self))
116 self)
117
118(defmethod multiply-by ((self poly) (other scalar))
119 (mapc #'(lambda (term) (multiply-by term other))
120 (poly-termlist self))
121 self)
122
123
124(defmacro fast-add/subtract (p q order-fn add/subtract-fn uminus-fn)
125 "Return an expression which will efficiently adds/subtracts two
126polynomials, P and Q. The addition/subtraction of coefficients is
127performed by calling ADD/SUBTRACT-METHOD-NAME. If UMINUS-METHOD-NAME
128is supplied, it is used to negate the coefficients of Q which do not
129have a corresponding coefficient in P. The code implements an
130efficient algorithm to add two polynomials represented as sorted lists
131of terms. The code destroys both arguments, reusing the terms to build
132the result."
133 `(macrolet ((lc (x) `(r-coeff (car ,x))))
134 (do ((p ,p)
135 (q ,q)
136 r)
137 ((or (endp p) (endp q))
138 ;; NOTE: R contains the result in reverse order. Can it
139 ;; be more efficient to produce the terms in correct order?
140 (unless (endp q)
141 ;; Upon subtraction, we must change the sign of
142 ;; all coefficients in q
143 ,@(when uminus-fn
144 `((mapc #'(lambda (x) (setf x (funcall ,uminus-fn x))) q)))
145 (setf r (nreconc r q)))
146 r)
147 (multiple-value-bind
148 (greater-p equal-p)
149 (funcall ,order-fn (car p) (car q))
150 (cond
151 (greater-p
152 (rotatef (cdr p) r p)
153 )
154 (equal-p
155 (let ((s (funcall ,add/subtract-fn (lc p) (lc q))))
156 (cond
157 ((r-zerop s)
158 (setf p (cdr p))
159 )
160 (t
161 (setf (lc p) s)
162 (rotatef (cdr p) r p))))
163 (setf q (cdr q))
164 )
165 (t
166 ;;Negate the term of Q if UMINUS provided, signallig
167 ;;that we are doing subtraction
168 ,(when uminus-fn
169 `(setf (lc q) (funcall ,uminus-fn (lc q))))
170 (rotatef (cdr q) r q)))))))
171
172
173(defmacro def-add/subtract-method (add/subtract-method-name
174 uminus-method-name
175 &optional
176 (doc-string nil doc-string-supplied-p))
177 "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
178 `(defmethod ,add/subtract-method-name ((self poly) (other poly))
179 ,@(when doc-string-supplied-p `(,doc-string))
180 ;; Ensure orders are compatible
181 (change-term-order other self)
182 (setf (poly-termlist self) (fast-add/subtract
183 (poly-termlist self) (poly-termlist other)
184 (poly-term-order self)
185 #',add/subtract-method-name
186 ,(when uminus-method-name `(function ,uminus-method-name))))
187 self))
188
189(eval-when (:compile-toplevel :load-toplevel :execute)
190
191 (def-add/subtract-method add-to nil
192 "Adds to polynomial SELF another polynomial OTHER.
193This operation destructively modifies both polynomials.
194The result is stored in SELF. This implementation does
195no consing, entirely reusing the sells of SELF and OTHER.")
196
197 (def-add/subtract-method subtract-from unary-minus
198 "Subtracts from polynomial SELF another polynomial OTHER.
199This operation destructively modifies both polynomials.
200The result is stored in SELF. This implementation does
201no consing, entirely reusing the sells of SELF and OTHER.")
202
203 )
204
205
206
207(defmethod unary-minus ((self poly))
208 "Destructively modifies the coefficients of the polynomial SELF,
209by changing their sign."
210 (mapc #'unary-minus (poly-termlist self))
211 self)
212
213(defun add-termlists (p q order-fn)
214 "Destructively adds two termlists P and Q ordered according to ORDER-FN."
215 (fast-add/subtract p q order-fn #'add-to nil))
216
217(defmacro multiply-term-by-termlist-dropping-zeros (term termlist
218 &optional (reverse-arg-order-P nil))
219 "Multiplies term TERM by a list of term, TERMLIST.
220Takes into accound divisors of zero in the ring, by
221deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
222is T, change the order of arguments; this may be important
223if we extend the package to non-commutative rings."
224 `(mapcan #'(lambda (other-term)
225 (let ((prod (r*
226 ,@(cond
227 (reverse-arg-order-p
228 `(other-term ,term))
229 (t
230 `(,term other-term))))))
231 (cond
232 ((r-zerop prod) nil)
233 (t (list prod)))))
234 ,termlist))
235
236(defun multiply-termlists (p q order-fn)
237 (cond
238 ((or (endp p) (endp q))
239 ;;p or q is 0 (represented by NIL)
240 nil)
241 ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
242 ((endp (cdr p))
243 (multiply-term-by-termlist-dropping-zeros (car p) q))
244 ((endp (cdr q))
245 (multiply-term-by-termlist-dropping-zeros (car q) p t))
246 (t
247 (cons (r* (car p) (car q))
248 (add-termlists
249 (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
250 (multiply-termlists (cdr p) q order-fn)
251 order-fn)))))
252
253(defmethod multiply-by ((self poly) (other poly))
254 (change-term-order other self)
255 (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
256 (poly-termlist other)
257 (poly-term-order self)))
258 self)
259
260(defmethod r* ((poly1 poly) (poly2 poly))
261 "Non-destructively multiply POLY1 by POLY2."
262 (multiply-by (copy-instance POLY1) (copy-instance POLY2)))
263
264(defmethod left-tensor-product-by ((self poly) (other term))
265 (setf (poly-termlist self)
266 (mapcan #'(lambda (term)
267 (let ((prod (left-tensor-product-by term other)))
268 (cond
269 ((r-zerop prod) nil)
270 (t (list prod)))))
271 (poly-termlist self)))
272 self)
273
274(defmethod right-tensor-product-by ((self poly) (other term))
275 (setf (poly-termlist self)
276 (mapcan #'(lambda (term)
277 (let ((prod (right-tensor-product-by term other)))
278 (cond
279 ((r-zerop prod) nil)
280 (t (list prod)))))
281 (poly-termlist self)))
282 self)
283
284(defmethod left-tensor-product-by ((self poly) (other monom))
285 (setf (poly-termlist self)
286 (mapcan #'(lambda (term)
287 (let ((prod (left-tensor-product-by term other)))
288 (cond
289 ((r-zerop prod) nil)
290 (t (list prod)))))
291 (poly-termlist self)))
292 self)
293
294(defmethod right-tensor-product-by ((self poly) (other monom))
295 (setf (poly-termlist self)
296 (mapcan #'(lambda (term)
297 (let ((prod (right-tensor-product-by term other)))
298 (cond
299 ((r-zerop prod) nil)
300 (t (list prod)))))
301 (poly-termlist self)))
302 self)
303
304
305(defun standard-extension (plist &aux (k (length plist)) (i 0))
306 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
307is a list of polynomials. Destructively modifies PLIST elements."
308 (mapc #'(lambda (poly)
309 (left-tensor-product-by
310 poly
311 (prog1
312 (make-monom-variable k i)
313 (incf i))))
314 plist))
315
316(defmethod poly-dimension ((poly poly))
317 (cond ((r-zerop poly) -1)
318 (t (monom-dimension (leading-term poly)))))
319
320(defun standard-extension-1 (plist
321 &aux
322 (k (length plist))
323 (plist (standard-extension plist))
324 (nvars (poly-dimension (car plist))))
325 "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
326Firstly, new K variables U1, U2, ..., UK, are inserted into each
327polynomial. Subsequently, P1, P2, ..., PK are destructively modified
328tantamount to replacing PI with UI*PI-1."
329 ;; Implementation note: we use STANDARD-EXTENSION and then subtract
330 ;; 1 from each polynomial; since UI*PI has no constant term,
331 ;; we just need to append the constant term at the end
332 ;; of each termlist.
333 (flet ((subtract-1 (p)
334 (append-item p (make-instance 'term :coeff -1 :dimension (+ k nvars)))))
335 (setf plist (mapc #'subtract-1 plist)))
336 plist)
337
338#|
339
340(defun standard-sum (F plist
341 &aux
342 (k (length plist))
343 (d (+ k (monom-dimension (poly-lt (car plist)))))
344 ;; Add k variables to f
345 (f (poly-list-add-variables f k))
346 ;; Set PLIST to [U1*P1,U2*P2,...,UK*PK]
347 (plist (apply #'nconc (poly-standard-extension plist))))
348 "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
349Firstly, new K variables, U1, U2, ..., UK, are inserted into each
350polynomial. Subsequently, P1, P2, ..., PK are destructively modified
351tantamount to replacing PI with UI*PI, and the resulting polynomials
352are added. It should be noted that the term order is not modified,
353which is equivalent to using a lexicographic order on the first K
354variables."
355 (setf (cdr (last (poly-termlist plist)))
356 ;; Add -1 as the last term
357 (list (make-term :monom (make-monom :dimension d)
358 :coeff (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
359 (append f (list plist)))
360
361
362
363
364(defun saturation-extension-1 (ring f p)
365 "Calculate [F, U*P-1]. It destructively modifies F."
366 (declare (type ring ring))
367 (polysaturation-extension ring f (list p)))
368
369;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
370;;
371;; Evaluation of polynomial (prefix) expressions
372;;
373;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
374
375(defun coerce-coeff (ring expr vars)
376 "Coerce an element of the coefficient ring to a constant polynomial."
377 ;; Modular arithmetic handler by rat
378 (declare (type ring ring))
379 (make-poly-from-termlist (list (make-term :monom (make-monom :dimension (length vars))
380 :coeff (funcall (ring-parse ring) expr)))
381 0))
382
383(defun poly-eval (expr vars
384 &optional
385 (ring +ring-of-integers+)
386 (order #'lex>)
387 (list-marker :[)
388 &aux
389 (ring-and-order (make-ring-and-order :ring ring :order order)))
390 "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
391variables VARS. Return the resulting polynomial or list of
392polynomials. Standard arithmetical operators in form EXPR are
393replaced with their analogues in the ring of polynomials, and the
394resulting expression is evaluated, resulting in a polynomial or a list
395of polynomials in internal form. A similar operation in another computer
396algebra system could be called 'expand' or so."
397 (declare (type ring ring))
398 (labels ((p-eval (arg) (poly-eval arg vars ring order))
399 (p-eval-scalar (arg) (poly-eval-scalar arg))
400 (p-eval-list (args) (mapcar #'p-eval args))
401 (p-add (x y) (poly-add ring-and-order x y)))
402 (cond
403 ((null expr) (error "Empty expression"))
404 ((eql expr 0) (make-poly-zero))
405 ((member expr vars :test #'equalp)
406 (let ((pos (position expr vars :test #'equalp)))
407 (make-poly-variable ring (length vars) pos)))
408 ((atom expr)
409 (coerce-coeff ring expr vars))
410 ((eq (car expr) list-marker)
411 (cons list-marker (p-eval-list (cdr expr))))
412 (t
413 (case (car expr)
414 (+ (reduce #'p-add (p-eval-list (cdr expr))))
415 (- (case (length expr)
416 (1 (make-poly-zero))
417 (2 (poly-uminus ring (p-eval (cadr expr))))
418 (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
419 (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
420 (reduce #'p-add (p-eval-list (cddr expr)))))))
421 (*
422 (if (endp (cddr expr)) ;unary
423 (p-eval (cdr expr))
424 (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
425 (/
426 ;; A polynomial can be divided by a scalar
427 (cond
428 ((endp (cddr expr))
429 ;; A special case (/ ?), the inverse
430 (coerce-coeff ring (apply (ring-div ring) (cdr expr)) vars))
431 (t
432 (let ((num (p-eval (cadr expr)))
433 (denom-inverse (apply (ring-div ring)
434 (cons (funcall (ring-unit ring))
435 (mapcar #'p-eval-scalar (cddr expr))))))
436 (scalar-times-poly ring denom-inverse num)))))
437 (expt
438 (cond
439 ((member (cadr expr) vars :test #'equalp)
440 ;;Special handling of (expt var pow)
441 (let ((pos (position (cadr expr) vars :test #'equalp)))
442 (make-poly-variable ring (length vars) pos (caddr expr))))
443 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
444 ;; Negative power means division in coefficient ring
445 ;; Non-integer power means non-polynomial coefficient
446 (coerce-coeff ring expr vars))
447 (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
448 (otherwise
449 (coerce-coeff ring expr vars)))))))
450
451(defun poly-eval-scalar (expr
452 &optional
453 (ring +ring-of-integers+)
454 &aux
455 (order #'lex>))
456 "Evaluate a scalar expression EXPR in ring RING."
457 (declare (type ring ring))
458 (poly-lc (poly-eval expr nil ring order)))
459
460(defun spoly (ring-and-order f g
461 &aux
462 (ring (ro-ring ring-and-order)))
463 "It yields the S-polynomial of polynomials F and G."
464 (declare (type ring-and-order ring-and-order) (type poly f g))
465 (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
466 (mf (monom-div lcm (poly-lm f)))
467 (mg (monom-div lcm (poly-lm g))))
468 (declare (type monom mf mg))
469 (multiple-value-bind (c cf cg)
470 (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
471 (declare (ignore c))
472 (poly-sub
473 ring-and-order
474 (scalar-times-poly ring cg (monom-times-poly mf f))
475 (scalar-times-poly ring cf (monom-times-poly mg g))))))
476
477
478(defun poly-primitive-part (ring p)
479 "Divide polynomial P with integer coefficients by gcd of its
480coefficients and return the result."
481 (declare (type ring ring) (type poly p))
482 (if (poly-zerop p)
483 (values p 1)
484 (let ((c (poly-content ring p)))
485 (values (make-poly-from-termlist
486 (mapcar
487 #'(lambda (x)
488 (make-term :monom (term-monom x)
489 :coeff (funcall (ring-div ring) (term-coeff x) c)))
490 (poly-termlist p))
491 (poly-sugar p))
492 c))))
493
494(defun poly-content (ring p)
495 "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
496to compute the greatest common divisor."
497 (declare (type ring ring) (type poly p))
498 (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
499
500(defun read-infix-form (&key (stream t))
501 "Parser of infix expressions with integer/rational coefficients
502The parser will recognize two kinds of polynomial expressions:
503
504- polynomials in fully expanded forms with coefficients
505 written in front of symbolic expressions; constants can be optionally
506 enclosed in (); for example, the infix form
507 X^2-Y^2+(-4/3)*U^2*W^3-5
508 parses to
509 (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))
510
511- lists of polynomials; for example
512 [X-Y, X^2+3*Z]
513 parses to
514 (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
515 where the first symbol [ marks a list of polynomials.
516
517-other infix expressions, for example
518 [(X-Y)*(X+Y)/Z,(X+1)^2]
519parses to:
520 (:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
521Currently this function is implemented using M. Kantrowitz's INFIX package."
522 (read-from-string
523 (concatenate 'string
524 "#I("
525 (with-output-to-string (s)
526 (loop
527 (multiple-value-bind (line eof)
528 (read-line stream t)
529 (format s "~A" line)
530 (when eof (return)))))
531 ")")))
532
533(defun read-poly (vars &key
534 (stream t)
535 (ring +ring-of-integers+)
536 (order #'lex>))
537 "Reads an expression in prefix form from a stream STREAM.
538The expression read from the strem should represent a polynomial or a
539list of polynomials in variables VARS, over the ring RING. The
540polynomial or list of polynomials is returned, with terms in each
541polynomial ordered according to monomial order ORDER."
542 (poly-eval (read-infix-form :stream stream) vars ring order))
543
544(defun string->poly (str vars
545 &optional
546 (ring +ring-of-integers+)
547 (order #'lex>))
548 "Converts a string STR to a polynomial in variables VARS."
549 (with-input-from-string (s str)
550 (read-poly vars :stream s :ring ring :order order)))
551
552(defun poly->alist (p)
553 "Convert a polynomial P to an association list. Thus, the format of the
554returned value is ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
555MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
556corresponding coefficient in the ring."
557 (cond
558 ((poly-p p)
559 (mapcar #'term->cons (poly-termlist p)))
560 ((and (consp p) (eq (car p) :[))
561 (cons :[ (mapcar #'poly->alist (cdr p))))))
562
563(defun string->alist (str vars
564 &optional
565 (ring +ring-of-integers+)
566 (order #'lex>))
567 "Convert a string STR representing a polynomial or polynomial list to
568an association list (... (MONOM . COEFF) ...)."
569 (poly->alist (string->poly str vars ring order)))
570
571(defun poly-equal-no-sugar-p (p q)
572 "Compare polynomials for equality, ignoring sugar."
573 (declare (type poly p q))
574 (equalp (poly-termlist p) (poly-termlist q)))
575
576(defun poly-set-equal-no-sugar-p (p q)
577 "Compare polynomial sets P and Q for equality, ignoring sugar."
578 (null (set-exclusive-or p q :test #'poly-equal-no-sugar-p )))
579
580(defun poly-list-equal-no-sugar-p (p q)
581 "Compare polynomial lists P and Q for equality, ignoring sugar."
582 (every #'poly-equal-no-sugar-p p q))
583|#
Note: See TracBrowser for help on using the repository browser.