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

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

* empty log message *
File size: 17.1 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 :ring-and-order :monom :order :term :termlist :infix)
30 (:export "POLY"
31 "POLY-TERMLIST"
32 "POLY-SUGAR"
33 "POLY-RESET-SUGAR"
34 "POLY-LT"
35 "MAKE-POLY-FROM-TERMLIST"
36 "MAKE-POLY-ZERO"
37 "MAKE-POLY-VARIABLE"
38 "POLY-UNIT"
39 "POLY-LM"
40 "POLY-SECOND-LM"
41 "POLY-SECOND-LT"
42 "POLY-LC"
43 "POLY-SECOND-LC"
44 "POLY-ZEROP"
45 "POLY-LENGTH"
46 "SCALAR-TIMES-POLY"
47 "SCALAR-TIMES-POLY-1"
48 "MONOM-TIMES-POLY"
49 "TERM-TIMES-POLY"
50 "POLY-ADD"
51 "POLY-SUB"
52 "POLY-UMINUS"
53 "POLY-MUL"
54 "POLY-EXPT"
55 "POLY-APPEND"
56 "POLY-NREVERSE"
57 "POLY-REVERSE"
58 "POLY-CONTRACT"
59 "POLY-EXTEND"
60 "POLY-ADD-VARIABLES"
61 "POLY-LIST-ADD-VARIABLES"
62 "POLY-STANDARD-EXTENSION"
63 "SATURATION-EXTENSION"
64 "POLYSATURATION-EXTENSION"
65 "SATURATION-EXTENSION-1"
66 "COERCE-COEFF"
67 "POLY-EVAL"
68 "POLY-EVAL-SCALAR"
69 "SPOLY"
70 "POLY-PRIMITIVE-PART"
71 "POLY-CONTENT"
72 "READ-INFIX-FORM"
73 "READ-POLY"
74 "STRING->POLY"
75 "POLY->ALIST"
76 "STRING->ALIST"
77 "POLY-EQUAL-NO-SUGAR-P"
78 "POLY-SET-EQUAL-NO-SUGAR-P"
79 "POLY-LIST-EQUAL-NO-SUGAR-P"
80 ))
81
82(in-package :polynomial)
83
84(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
85
86#|
87 ;;
88 ;; BOA constructor, by default constructs zero polynomial
89 (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
90 (:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
91 ;; Constructor of polynomials representing a variable
92 (:constructor make-poly-variable (ring nvars pos &optional (power 1)
93 &aux
94 (termlist (list
95 (make-term-variable ring nvars pos power)))
96 (sugar power)))
97 (:constructor poly-unit (ring dimension
98 &aux
99 (termlist (termlist-unit ring dimension))
100 (sugar 0))))
101
102|#
103
104(defclass poly ()
105 ((termlist :initarg :terms :accessor poly-termlist))
106 (:default-initargs :termlist nil))
107
108;; Leading term
109(defgeneric leading-term (object)
110 (:method ((self poly))
111 (car (poly-termlist self))))
112
113;; Second term
114(defgeneric second-leading-term (object)
115 (:method ((self poly))
116 (cadar (poly-termlist self))))
117
118;; Leading coefficient
119(defgeneric leading-coefficient (object)
120 (:method ((self poly))
121 (r-coeff (leading-term self))))
122
123;; Second coefficient
124(defgeneric second-leading-coefficient (object)
125 (:method ((self poly))
126 (term-coeff (second-leading-term self))))
127
128;; Testing for a zero polynomial
129(defun poly-zerop (p)
130 (declare (type poly p))
131 (null (poly-termlist p)))
132
133;; The number of terms
134(defun poly-length (p)
135 (declare (type poly p))
136 (length (poly-termlist p)))
137
138(defun poly-reset-sugar (p)
139 "(Re)sets the sugar of a polynomial P to the sugar of (POLY-TERMLIST P).
140Thus, the sugar is set to the maximum sugar of all monomials of P, or -1
141if P is a zero polynomial."
142 (declare (type poly p))
143 (setf (poly-sugar p) (termlist-sugar (poly-termlist p)))
144 p)
145
146(defun scalar-times-poly (ring c p)
147 "The scalar product of scalar C by a polynomial P. The sugar of the
148original polynomial becomes the sugar of the result."
149 (declare (type ring ring) (type poly p))
150 (make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p)))
151
152(defun scalar-times-poly-1 (ring c p)
153 "The scalar product of scalar C by a polynomial P, omitting the head term. The sugar of the
154original polynomial becomes the sugar of the result."
155 (declare (type ring ring) (type poly p))
156 (make-poly-from-termlist (scalar-times-termlist ring c (cdr (poly-termlist p))) (poly-sugar p)))
157
158(defun monom-times-poly (m p)
159 (declare (type monom m) (type poly p))
160 (make-poly-from-termlist
161 (monom-times-termlist m (poly-termlist p))
162 (+ (poly-sugar p) (monom-sugar m))))
163
164(defun term-times-poly (ring term p)
165 (declare (type ring ring) (type term term) (type poly p))
166 (make-poly-from-termlist
167 (term-times-termlist ring term (poly-termlist p))
168 (+ (poly-sugar p) (term-sugar term))))
169
170(defun poly-add (ring-and-order p q)
171 (declare (type ring-and-order ring-and-order) (type poly p q))
172 (make-poly-from-termlist
173 (termlist-add ring-and-order
174 (poly-termlist p)
175 (poly-termlist q))
176 (max (poly-sugar p) (poly-sugar q))))
177
178(defun poly-sub (ring-and-order p q)
179 (declare (type ring-and-order ring-and-order) (type poly p q))
180 (make-poly-from-termlist
181 (termlist-sub ring-and-order (poly-termlist p) (poly-termlist q))
182 (max (poly-sugar p) (poly-sugar q))))
183
184(defun poly-uminus (ring p)
185 (declare (type ring ring) (type poly p))
186 (make-poly-from-termlist
187 (termlist-uminus ring (poly-termlist p))
188 (poly-sugar p)))
189
190(defun poly-mul (ring-and-order p q)
191 (declare (type ring-and-order ring-and-order) (type poly p q))
192 (make-poly-from-termlist
193 (termlist-mul ring-and-order (poly-termlist p) (poly-termlist q))
194 (+ (poly-sugar p) (poly-sugar q))))
195
196(defun poly-expt (ring-and-order p n)
197 (declare (type ring-and-order ring-and-order) (type poly p))
198 (make-poly-from-termlist (termlist-expt ring-and-order (poly-termlist p) n) (* n (poly-sugar p))))
199
200(defun poly-append (&rest plist)
201 (make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist))
202 (apply #'max (mapcar #'poly-sugar plist))))
203
204(defun poly-nreverse (p)
205 "Destructively reverse the order of terms in polynomial P. Returns P"
206 (declare (type poly p))
207 (setf (poly-termlist p) (nreverse (poly-termlist p)))
208 p)
209
210(defun poly-reverse (p)
211 "Returns a copy of the polynomial P with terms in reverse order."
212 (declare (type poly p))
213 (make-poly-from-termlist (reverse (poly-termlist p))
214 (poly-sugar p)))
215
216
217(defun poly-contract (p &optional (k 1))
218 (declare (type poly p))
219 (make-poly-from-termlist (termlist-contract (poly-termlist p) k)
220 (poly-sugar p)))
221
222(defun poly-extend (p &optional (m (make-monom :dimension 1)))
223 (declare (type poly p))
224 (make-poly-from-termlist
225 (termlist-extend (poly-termlist p) m)
226 (+ (poly-sugar p) (monom-sugar m))))
227
228(defun poly-add-variables (p k)
229 (declare (type poly p))
230 (setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k))
231 p)
232
233(defun poly-list-add-variables (plist k)
234 (mapcar #'(lambda (p) (poly-add-variables p k)) plist))
235
236(defun poly-standard-extension (plist &aux (k (length plist)))
237 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
238 (declare (list plist) (fixnum k))
239 (labels ((incf-power (g i)
240 (dolist (x (poly-termlist g))
241 (incf (monom-elt (term-monom x) i)))
242 (incf (poly-sugar g))))
243 (setf plist (poly-list-add-variables plist k))
244 (dotimes (i k plist)
245 (incf-power (nth i plist) i))))
246
247(defun saturation-extension (ring f plist
248 &aux
249 (k (length plist))
250 (d (monom-dimension (poly-lm (car plist))))
251 f-x plist-x)
252 "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
253 (declare (type ring ring))
254 (setf f-x (poly-list-add-variables f k)
255 plist-x (mapcar #'(lambda (x)
256 (setf (poly-termlist x)
257 (nconc (poly-termlist x)
258 (list (make-term :monom (make-monom :dimension d)
259 :coeff (funcall (ring-uminus ring)
260 (funcall (ring-unit ring)))))))
261 x)
262 (poly-standard-extension plist)))
263 (append f-x plist-x))
264
265
266(defun polysaturation-extension (ring f plist
267 &aux
268 (k (length plist))
269 (d (+ k (monom-dimension (poly-lm (car plist)))))
270 ;; Add k variables to f
271 (f (poly-list-add-variables f k))
272 ;; Set PLIST to [U1*P1,U2*P2,...,UK*PK]
273 (plist (apply #'poly-append (poly-standard-extension plist))))
274 "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]. It destructively modifies F."
275 ;; Add -1 as the last term
276 (declare (type ring ring))
277 (setf (cdr (last (poly-termlist plist)))
278 (list (make-term :monom (make-monom :dimension d)
279 :coeff (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
280 (append f (list plist)))
281
282(defun saturation-extension-1 (ring f p)
283 "Calculate [F, U*P-1]. It destructively modifies F."
284 (declare (type ring ring))
285 (polysaturation-extension ring f (list p)))
286
287;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
288;;
289;; Evaluation of polynomial (prefix) expressions
290;;
291;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
292
293(defun coerce-coeff (ring expr vars)
294 "Coerce an element of the coefficient ring to a constant polynomial."
295 ;; Modular arithmetic handler by rat
296 (declare (type ring ring))
297 (make-poly-from-termlist (list (make-term :monom (make-monom :dimension (length vars))
298 :coeff (funcall (ring-parse ring) expr)))
299 0))
300
301(defun poly-eval (expr vars
302 &optional
303 (ring +ring-of-integers+)
304 (order #'lex>)
305 (list-marker :[)
306 &aux
307 (ring-and-order (make-ring-and-order :ring ring :order order)))
308 "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
309variables VARS. Return the resulting polynomial or list of
310polynomials. Standard arithmetical operators in form EXPR are
311replaced with their analogues in the ring of polynomials, and the
312resulting expression is evaluated, resulting in a polynomial or a list
313of polynomials in internal form. A similar operation in another computer
314algebra system could be called 'expand' or so."
315 (declare (type ring ring))
316 (labels ((p-eval (arg) (poly-eval arg vars ring order))
317 (p-eval-scalar (arg) (poly-eval-scalar arg))
318 (p-eval-list (args) (mapcar #'p-eval args))
319 (p-add (x y) (poly-add ring-and-order x y)))
320 (cond
321 ((null expr) (error "Empty expression"))
322 ((eql expr 0) (make-poly-zero))
323 ((member expr vars :test #'equalp)
324 (let ((pos (position expr vars :test #'equalp)))
325 (make-poly-variable ring (length vars) pos)))
326 ((atom expr)
327 (coerce-coeff ring expr vars))
328 ((eq (car expr) list-marker)
329 (cons list-marker (p-eval-list (cdr expr))))
330 (t
331 (case (car expr)
332 (+ (reduce #'p-add (p-eval-list (cdr expr))))
333 (- (case (length expr)
334 (1 (make-poly-zero))
335 (2 (poly-uminus ring (p-eval (cadr expr))))
336 (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
337 (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
338 (reduce #'p-add (p-eval-list (cddr expr)))))))
339 (*
340 (if (endp (cddr expr)) ;unary
341 (p-eval (cdr expr))
342 (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
343 (/
344 ;; A polynomial can be divided by a scalar
345 (cond
346 ((endp (cddr expr))
347 ;; A special case (/ ?), the inverse
348 (coerce-coeff ring (apply (ring-div ring) (cdr expr)) vars))
349 (t
350 (let ((num (p-eval (cadr expr)))
351 (denom-inverse (apply (ring-div ring)
352 (cons (funcall (ring-unit ring))
353 (mapcar #'p-eval-scalar (cddr expr))))))
354 (scalar-times-poly ring denom-inverse num)))))
355 (expt
356 (cond
357 ((member (cadr expr) vars :test #'equalp)
358 ;;Special handling of (expt var pow)
359 (let ((pos (position (cadr expr) vars :test #'equalp)))
360 (make-poly-variable ring (length vars) pos (caddr expr))))
361 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
362 ;; Negative power means division in coefficient ring
363 ;; Non-integer power means non-polynomial coefficient
364 (coerce-coeff ring expr vars))
365 (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
366 (otherwise
367 (coerce-coeff ring expr vars)))))))
368
369(defun poly-eval-scalar (expr
370 &optional
371 (ring +ring-of-integers+)
372 &aux
373 (order #'lex>))
374 "Evaluate a scalar expression EXPR in ring RING."
375 (declare (type ring ring))
376 (poly-lc (poly-eval expr nil ring order)))
377
378(defun spoly (ring-and-order f g
379 &aux
380 (ring (ro-ring ring-and-order)))
381 "It yields the S-polynomial of polynomials F and G."
382 (declare (type ring-and-order ring-and-order) (type poly f g))
383 (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
384 (mf (monom-div lcm (poly-lm f)))
385 (mg (monom-div lcm (poly-lm g))))
386 (declare (type monom mf mg))
387 (multiple-value-bind (c cf cg)
388 (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
389 (declare (ignore c))
390 (poly-sub
391 ring-and-order
392 (scalar-times-poly ring cg (monom-times-poly mf f))
393 (scalar-times-poly ring cf (monom-times-poly mg g))))))
394
395
396(defun poly-primitive-part (ring p)
397 "Divide polynomial P with integer coefficients by gcd of its
398coefficients and return the result."
399 (declare (type ring ring) (type poly p))
400 (if (poly-zerop p)
401 (values p 1)
402 (let ((c (poly-content ring p)))
403 (values (make-poly-from-termlist
404 (mapcar
405 #'(lambda (x)
406 (make-term :monom (term-monom x)
407 :coeff (funcall (ring-div ring) (term-coeff x) c)))
408 (poly-termlist p))
409 (poly-sugar p))
410 c))))
411
412(defun poly-content (ring p)
413 "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
414to compute the greatest common divisor."
415 (declare (type ring ring) (type poly p))
416 (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
417
418(defun read-infix-form (&key (stream t))
419 "Parser of infix expressions with integer/rational coefficients
420The parser will recognize two kinds of polynomial expressions:
421
422- polynomials in fully expanded forms with coefficients
423 written in front of symbolic expressions; constants can be optionally
424 enclosed in (); for example, the infix form
425 X^2-Y^2+(-4/3)*U^2*W^3-5
426 parses to
427 (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))
428
429- lists of polynomials; for example
430 [X-Y, X^2+3*Z]
431 parses to
432 (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
433 where the first symbol [ marks a list of polynomials.
434
435-other infix expressions, for example
436 [(X-Y)*(X+Y)/Z,(X+1)^2]
437parses to:
438 (:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
439Currently this function is implemented using M. Kantrowitz's INFIX package."
440 (read-from-string
441 (concatenate 'string
442 "#I("
443 (with-output-to-string (s)
444 (loop
445 (multiple-value-bind (line eof)
446 (read-line stream t)
447 (format s "~A" line)
448 (when eof (return)))))
449 ")")))
450
451(defun read-poly (vars &key
452 (stream t)
453 (ring +ring-of-integers+)
454 (order #'lex>))
455 "Reads an expression in prefix form from a stream STREAM.
456The expression read from the strem should represent a polynomial or a
457list of polynomials in variables VARS, over the ring RING. The
458polynomial or list of polynomials is returned, with terms in each
459polynomial ordered according to monomial order ORDER."
460 (poly-eval (read-infix-form :stream stream) vars ring order))
461
462(defun string->poly (str vars
463 &optional
464 (ring +ring-of-integers+)
465 (order #'lex>))
466 "Converts a string STR to a polynomial in variables VARS."
467 (with-input-from-string (s str)
468 (read-poly vars :stream s :ring ring :order order)))
469
470(defun poly->alist (p)
471 "Convert a polynomial P to an association list. Thus, the format of the
472returned value is ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
473MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
474corresponding coefficient in the ring."
475 (cond
476 ((poly-p p)
477 (mapcar #'term->cons (poly-termlist p)))
478 ((and (consp p) (eq (car p) :[))
479 (cons :[ (mapcar #'poly->alist (cdr p))))))
480
481(defun string->alist (str vars
482 &optional
483 (ring +ring-of-integers+)
484 (order #'lex>))
485 "Convert a string STR representing a polynomial or polynomial list to
486an association list (... (MONOM . COEFF) ...)."
487 (poly->alist (string->poly str vars ring order)))
488
489(defun poly-equal-no-sugar-p (p q)
490 "Compare polynomials for equality, ignoring sugar."
491 (declare (type poly p q))
492 (equalp (poly-termlist p) (poly-termlist q)))
493
494(defun poly-set-equal-no-sugar-p (p q)
495 "Compare polynomial sets P and Q for equality, ignoring sugar."
496 (null (set-exclusive-or p q :test #'poly-equal-no-sugar-p )))
497
498(defun poly-list-equal-no-sugar-p (p q)
499 "Compare polynomial lists P and Q for equality, ignoring sugar."
500 (every #'poly-equal-no-sugar-p p q))
Note: See TracBrowser for help on using the repository browser.