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

Last change on this file since 2485 was 2483, checked in by Marek Rychlik, 10 years ago

* empty log message *
File size: 13.4 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;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
23;;
24;; Polynomials
25;;
26;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
27
28(defpackage "POLYNOMIAL"
29 (:use :cl :ring :monom :order :term #| :infix |# )
30 (:export "POLY"
31 ))
32
33(in-package :polynomial)
34
35(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
36
37#|
38 ;;
39 ;; BOA constructor, by default constructs zero polynomial
40 (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
41 (:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
42 ;; Constructor of polynomials representing a variable
43 (:constructor make-poly-variable (ring nvars pos &optional (power 1)
44 &aux
45 (termlist (list
46 (make-term-variable ring nvars pos power)))
47 (sugar power)))
48 (:constructor poly-unit (ring dimension
49 &aux
50 (termlist (termlist-unit ring dimension))
51 (sugar 0))))
52
53|#
54
55(defclass poly ()
56 ((termlist :initarg :termlist :accessor poly-termlist))
57 (:default-initargs :termlist nil))
58
59(defmethod print-object ((self poly) stream)
60 (format stream "#<POLY TERMLIST=~A >" (poly-termlist self)))
61
62(defgeneric insert-item (object item)
63 (:method ((self poly) (item term))
64 (push item (poly-termlist self))
65 self))
66
67
68(defgeneric append-item (object item)
69 (:method ((self poly) (item term))
70 (setf (cdr (last (poly-termlist self))) (list item))
71 self))
72
73;; Leading term
74(defgeneric leading-term (object)
75 (:method ((self poly))
76 (car (poly-termlist self))))
77
78;; Second term
79(defgeneric second-leading-term (object)
80 (:method ((self poly))
81 (cadar (poly-termlist self))))
82
83;; Leading coefficient
84(defgeneric leading-coefficient (object)
85 (:method ((self poly))
86 (r-coeff (leading-term self))))
87
88;; Second coefficient
89(defgeneric second-leading-coefficient (object)
90 (:method ((self poly))
91 (r-coeff (second-leading-term self))))
92
93;; Testing for a zero polynomial
94(defmethod r-zerop ((self poly))
95 (null (poly-termlist self)))
96
97;; The number of terms
98(defmethod r-length ((self poly))
99 (length (poly-termlist self)))
100
101(defmethod multiply-by ((self poly) (other scalar))
102 (mapc #'(lambda (term) (multiply-by term other)) (poly-termlist self))
103 self)
104
105(defmethod multiply-by ((self poly) (other monom))
106 (mapc #'(lambda (term) (multiply-by term monom)) (poly-termlist self))
107 self)
108
109
110(defgeneric add-to (self other)
111 (:method ((self poly) (other poly))))
112
113(defgeneric subtract-from (self other)
114 (:method ((self poly) (other poly))))
115
116(defmethod unary-uminus (self))
117
118(defun poly-standard-extension (plist &aux (k (length plist)))
119 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
120 (declare (list plist) (fixnum k))
121 (labels ((incf-power (g i)
122 (dolist (x (poly-termlist g))
123 (incf (monom-elt (term-monom x) i)))
124 (incf (poly-sugar g))))
125 (setf plist (poly-list-add-variables plist k))
126 (dotimes (i k plist)
127 (incf-power (nth i plist) i))))
128
129(defun saturation-extension (ring f plist
130 &aux
131 (k (length plist))
132 (d (monom-dimension (poly-lm (car plist))))
133 f-x plist-x)
134 "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
135 (declare (type ring ring))
136 (setf f-x (poly-list-add-variables f k)
137 plist-x (mapcar #'(lambda (x)
138 (setf (poly-termlist x)
139 (nconc (poly-termlist x)
140 (list (make-term :monom (make-monom :dimension d)
141 :coeff (funcall (ring-uminus ring)
142 (funcall (ring-unit ring)))))))
143 x)
144 (poly-standard-extension plist)))
145 (append f-x plist-x))
146
147
148(defun polysaturation-extension (ring f plist
149 &aux
150 (k (length plist))
151 (d (+ k (monom-dimension (poly-lm (car plist)))))
152 ;; Add k variables to f
153 (f (poly-list-add-variables f k))
154 ;; Set PLIST to [U1*P1,U2*P2,...,UK*PK]
155 (plist (apply #'poly-append (poly-standard-extension plist))))
156 "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]. It destructively modifies F."
157 ;; Add -1 as the last term
158 (declare (type ring ring))
159 (setf (cdr (last (poly-termlist plist)))
160 (list (make-term :monom (make-monom :dimension d)
161 :coeff (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
162 (append f (list plist)))
163
164(defun saturation-extension-1 (ring f p)
165 "Calculate [F, U*P-1]. It destructively modifies F."
166 (declare (type ring ring))
167 (polysaturation-extension ring f (list p)))
168
169;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
170;;
171;; Evaluation of polynomial (prefix) expressions
172;;
173;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
174
175(defun coerce-coeff (ring expr vars)
176 "Coerce an element of the coefficient ring to a constant polynomial."
177 ;; Modular arithmetic handler by rat
178 (declare (type ring ring))
179 (make-poly-from-termlist (list (make-term :monom (make-monom :dimension (length vars))
180 :coeff (funcall (ring-parse ring) expr)))
181 0))
182
183(defun poly-eval (expr vars
184 &optional
185 (ring +ring-of-integers+)
186 (order #'lex>)
187 (list-marker :[)
188 &aux
189 (ring-and-order (make-ring-and-order :ring ring :order order)))
190 "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
191variables VARS. Return the resulting polynomial or list of
192polynomials. Standard arithmetical operators in form EXPR are
193replaced with their analogues in the ring of polynomials, and the
194resulting expression is evaluated, resulting in a polynomial or a list
195of polynomials in internal form. A similar operation in another computer
196algebra system could be called 'expand' or so."
197 (declare (type ring ring))
198 (labels ((p-eval (arg) (poly-eval arg vars ring order))
199 (p-eval-scalar (arg) (poly-eval-scalar arg))
200 (p-eval-list (args) (mapcar #'p-eval args))
201 (p-add (x y) (poly-add ring-and-order x y)))
202 (cond
203 ((null expr) (error "Empty expression"))
204 ((eql expr 0) (make-poly-zero))
205 ((member expr vars :test #'equalp)
206 (let ((pos (position expr vars :test #'equalp)))
207 (make-poly-variable ring (length vars) pos)))
208 ((atom expr)
209 (coerce-coeff ring expr vars))
210 ((eq (car expr) list-marker)
211 (cons list-marker (p-eval-list (cdr expr))))
212 (t
213 (case (car expr)
214 (+ (reduce #'p-add (p-eval-list (cdr expr))))
215 (- (case (length expr)
216 (1 (make-poly-zero))
217 (2 (poly-uminus ring (p-eval (cadr expr))))
218 (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
219 (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
220 (reduce #'p-add (p-eval-list (cddr expr)))))))
221 (*
222 (if (endp (cddr expr)) ;unary
223 (p-eval (cdr expr))
224 (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
225 (/
226 ;; A polynomial can be divided by a scalar
227 (cond
228 ((endp (cddr expr))
229 ;; A special case (/ ?), the inverse
230 (coerce-coeff ring (apply (ring-div ring) (cdr expr)) vars))
231 (t
232 (let ((num (p-eval (cadr expr)))
233 (denom-inverse (apply (ring-div ring)
234 (cons (funcall (ring-unit ring))
235 (mapcar #'p-eval-scalar (cddr expr))))))
236 (scalar-times-poly ring denom-inverse num)))))
237 (expt
238 (cond
239 ((member (cadr expr) vars :test #'equalp)
240 ;;Special handling of (expt var pow)
241 (let ((pos (position (cadr expr) vars :test #'equalp)))
242 (make-poly-variable ring (length vars) pos (caddr expr))))
243 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
244 ;; Negative power means division in coefficient ring
245 ;; Non-integer power means non-polynomial coefficient
246 (coerce-coeff ring expr vars))
247 (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
248 (otherwise
249 (coerce-coeff ring expr vars)))))))
250
251(defun poly-eval-scalar (expr
252 &optional
253 (ring +ring-of-integers+)
254 &aux
255 (order #'lex>))
256 "Evaluate a scalar expression EXPR in ring RING."
257 (declare (type ring ring))
258 (poly-lc (poly-eval expr nil ring order)))
259
260(defun spoly (ring-and-order f g
261 &aux
262 (ring (ro-ring ring-and-order)))
263 "It yields the S-polynomial of polynomials F and G."
264 (declare (type ring-and-order ring-and-order) (type poly f g))
265 (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
266 (mf (monom-div lcm (poly-lm f)))
267 (mg (monom-div lcm (poly-lm g))))
268 (declare (type monom mf mg))
269 (multiple-value-bind (c cf cg)
270 (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
271 (declare (ignore c))
272 (poly-sub
273 ring-and-order
274 (scalar-times-poly ring cg (monom-times-poly mf f))
275 (scalar-times-poly ring cf (monom-times-poly mg g))))))
276
277
278(defun poly-primitive-part (ring p)
279 "Divide polynomial P with integer coefficients by gcd of its
280coefficients and return the result."
281 (declare (type ring ring) (type poly p))
282 (if (poly-zerop p)
283 (values p 1)
284 (let ((c (poly-content ring p)))
285 (values (make-poly-from-termlist
286 (mapcar
287 #'(lambda (x)
288 (make-term :monom (term-monom x)
289 :coeff (funcall (ring-div ring) (term-coeff x) c)))
290 (poly-termlist p))
291 (poly-sugar p))
292 c))))
293
294(defun poly-content (ring p)
295 "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
296to compute the greatest common divisor."
297 (declare (type ring ring) (type poly p))
298 (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
299
300(defun read-infix-form (&key (stream t))
301 "Parser of infix expressions with integer/rational coefficients
302The parser will recognize two kinds of polynomial expressions:
303
304- polynomials in fully expanded forms with coefficients
305 written in front of symbolic expressions; constants can be optionally
306 enclosed in (); for example, the infix form
307 X^2-Y^2+(-4/3)*U^2*W^3-5
308 parses to
309 (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))
310
311- lists of polynomials; for example
312 [X-Y, X^2+3*Z]
313 parses to
314 (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
315 where the first symbol [ marks a list of polynomials.
316
317-other infix expressions, for example
318 [(X-Y)*(X+Y)/Z,(X+1)^2]
319parses to:
320 (:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
321Currently this function is implemented using M. Kantrowitz's INFIX package."
322 (read-from-string
323 (concatenate 'string
324 "#I("
325 (with-output-to-string (s)
326 (loop
327 (multiple-value-bind (line eof)
328 (read-line stream t)
329 (format s "~A" line)
330 (when eof (return)))))
331 ")")))
332
333(defun read-poly (vars &key
334 (stream t)
335 (ring +ring-of-integers+)
336 (order #'lex>))
337 "Reads an expression in prefix form from a stream STREAM.
338The expression read from the strem should represent a polynomial or a
339list of polynomials in variables VARS, over the ring RING. The
340polynomial or list of polynomials is returned, with terms in each
341polynomial ordered according to monomial order ORDER."
342 (poly-eval (read-infix-form :stream stream) vars ring order))
343
344(defun string->poly (str vars
345 &optional
346 (ring +ring-of-integers+)
347 (order #'lex>))
348 "Converts a string STR to a polynomial in variables VARS."
349 (with-input-from-string (s str)
350 (read-poly vars :stream s :ring ring :order order)))
351
352(defun poly->alist (p)
353 "Convert a polynomial P to an association list. Thus, the format of the
354returned value is ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
355MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
356corresponding coefficient in the ring."
357 (cond
358 ((poly-p p)
359 (mapcar #'term->cons (poly-termlist p)))
360 ((and (consp p) (eq (car p) :[))
361 (cons :[ (mapcar #'poly->alist (cdr p))))))
362
363(defun string->alist (str vars
364 &optional
365 (ring +ring-of-integers+)
366 (order #'lex>))
367 "Convert a string STR representing a polynomial or polynomial list to
368an association list (... (MONOM . COEFF) ...)."
369 (poly->alist (string->poly str vars ring order)))
370
371(defun poly-equal-no-sugar-p (p q)
372 "Compare polynomials for equality, ignoring sugar."
373 (declare (type poly p q))
374 (equalp (poly-termlist p) (poly-termlist q)))
375
376(defun poly-set-equal-no-sugar-p (p q)
377 "Compare polynomial sets P and Q for equality, ignoring sugar."
378 (null (set-exclusive-or p q :test #'poly-equal-no-sugar-p )))
379
380(defun poly-list-equal-no-sugar-p (p q)
381 "Compare polynomial lists P and Q for equality, ignoring sugar."
382 (every #'poly-equal-no-sugar-p p q))
383|#
Note: See TracBrowser for help on using the repository browser.