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

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

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