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/pol.lisp@ 1999

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

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