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

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

* empty log message *
File size: 9.5 KB
RevLine 
[77]1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
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
[431]23(defpackage "POLYNOMIAL"
24 (:use :cl :monomial :term :termlist))
[143]25
[431]26(in-package :polynomial)
27
[52]28;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
29;;
30;; Polynomials
31;;
32;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
33
34(defstruct (poly
35 ;;
36 ;; BOA constructor, by default constructs zero polynomial
37 (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
38 (:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
39 ;; Constructor of polynomials representing a variable
40 (:constructor make-variable (ring nvars pos &optional (power 1)
[53]41 &aux
42 (termlist (list
43 (make-term-variable ring nvars pos power)))
44 (sugar power)))
45 (:constructor poly-unit (ring dimension
46 &aux
47 (termlist (termlist-unit ring dimension))
48 (sugar 0))))
[52]49 (termlist nil :type list)
50 (sugar -1 :type fixnum))
51
52;; Leading term
53(defmacro poly-lt (p) `(car (poly-termlist ,p)))
54
55;; Second term
56(defmacro poly-second-lt (p) `(cadar (poly-termlist ,p)))
57
58;; Leading monomial
59(defun poly-lm (p) (term-monom (poly-lt p)))
60
61;; Second monomial
62(defun poly-second-lm (p) (term-monom (poly-second-lt p)))
63
64;; Leading coefficient
65(defun poly-lc (p) (term-coeff (poly-lt p)))
66
67;; Second coefficient
68(defun poly-second-lc (p) (term-coeff (poly-second-lt p)))
69
70;; Testing for a zero polynomial
71(defun poly-zerop (p) (null (poly-termlist p)))
72
73;; The number of terms
74(defun poly-length (p) (length (poly-termlist p)))
75
76(defun scalar-times-poly (ring c p)
77 (make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p)))
78
79;; The scalar product omitting the head term
80(defun scalar-times-poly-1 (ring c p)
81 (make-poly-from-termlist (scalar-times-termlist ring c (cdr (poly-termlist p))) (poly-sugar p)))
[53]82
[52]83(defun monom-times-poly (m p)
84 (make-poly-from-termlist (monom-times-termlist m (poly-termlist p)) (+ (poly-sugar p) (monom-sugar m))))
85
86(defun term-times-poly (ring term p)
87 (make-poly-from-termlist (term-times-termlist ring term (poly-termlist p)) (+ (poly-sugar p) (term-sugar term))))
88
89(defun poly-add (ring p q)
90 (make-poly-from-termlist (termlist-add ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
91
92(defun poly-sub (ring p q)
93 (make-poly-from-termlist (termlist-sub ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
94
95(defun poly-uminus (ring p)
96 (make-poly-from-termlist (termlist-uminus ring (poly-termlist p)) (poly-sugar p)))
97
98(defun poly-mul (ring p q)
99 (make-poly-from-termlist (termlist-mul ring (poly-termlist p) (poly-termlist q)) (+ (poly-sugar p) (poly-sugar q))))
100
101(defun poly-expt (ring p n)
102 (make-poly-from-termlist (termlist-expt ring (poly-termlist p) n) (* n (poly-sugar p))))
103
104(defun poly-append (&rest plist)
105 (make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist))
[53]106 (apply #'max (mapcar #'poly-sugar plist))))
[52]107
108(defun poly-nreverse (p)
109 (setf (poly-termlist p) (nreverse (poly-termlist p)))
110 p)
111
112(defun poly-contract (p &optional (k 1))
113 (make-poly-from-termlist (termlist-contract (poly-termlist p) k)
[53]114 (poly-sugar p)))
[52]115
116(defun poly-extend (p &optional (m (make-monom 1 :initial-element 0)))
117 (make-poly-from-termlist
118 (termlist-extend (poly-termlist p) m)
119 (+ (poly-sugar p) (monom-sugar m))))
120
121(defun poly-add-variables (p k)
122 (setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k))
123 p)
124
125(defun poly-list-add-variables (plist k)
126 (mapcar #'(lambda (p) (poly-add-variables p k)) plist))
127
128(defun poly-standard-extension (plist &aux (k (length plist)))
129 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
130 (declare (list plist) (fixnum k))
131 (labels ((incf-power (g i)
132 (dolist (x (poly-termlist g))
133 (incf (monom-elt (term-monom x) i)))
134 (incf (poly-sugar g))))
135 (setf plist (poly-list-add-variables plist k))
136 (dotimes (i k plist)
137 (incf-power (nth i plist) i))))
138
139(defun saturation-extension (ring f plist &aux (k (length plist)) (d (monom-dimension (poly-lm (car plist)))))
140 "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
141 (setf f (poly-list-add-variables f k)
142 plist (mapcar #'(lambda (x)
143 (setf (poly-termlist x) (nconc (poly-termlist x)
144 (list (make-term (make-monom d :initial-element 0)
145 (funcall (ring-uminus ring) (funcall (ring-unit ring)))))))
146 x)
147 (poly-standard-extension plist)))
148 (append f plist))
149
150
151(defun polysaturation-extension (ring f plist &aux (k (length plist))
[53]152 (d (+ k (length (poly-lm (car plist))))))
[52]153 "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]."
154 (setf f (poly-list-add-variables f k)
155 plist (apply #'poly-append (poly-standard-extension plist))
156 (cdr (last (poly-termlist plist))) (list (make-term (make-monom d :initial-element 0)
157 (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
158 (append f (list plist)))
159
160(defun saturation-extension-1 (ring f p) (polysaturation-extension ring f (list p)))
[53]161
162
163
164;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
165;;
166;; Evaluation of polynomial (prefix) expressions
167;;
168;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
169
170(defun coerce-coeff (ring expr vars)
171 "Coerce an element of the coefficient ring to a constant polynomial."
172 ;; Modular arithmetic handler by rat
173 (make-poly-from-termlist (list (make-term (make-monom (length vars) :initial-element 0)
174 (funcall (ring-parse ring) expr)))
175 0))
176
177(defun poly-eval (ring expr vars &optional (list-marker '[))
178 (labels ((p-eval (arg) (poly-eval ring arg vars))
179 (p-eval-list (args) (mapcar #'p-eval args))
180 (p-add (x y) (poly-add ring x y)))
181 (cond
182 ((eql expr 0) (make-poly-zero))
183 ((member expr vars :test #'equalp)
184 (let ((pos (position expr vars :test #'equalp)))
185 (make-variable ring (length vars) pos)))
186 ((atom expr)
187 (coerce-coeff ring expr vars))
188 ((eq (car expr) list-marker)
189 (cons list-marker (p-eval-list (cdr expr))))
190 (t
191 (case (car expr)
192 (+ (reduce #'p-add (p-eval-list (cdr expr))))
193 (- (case (length expr)
194 (1 (make-poly-zero))
195 (2 (poly-uminus ring (p-eval (cadr expr))))
196 (3 (poly-sub ring (p-eval (cadr expr)) (p-eval (caddr expr))))
197 (otherwise (poly-sub ring (p-eval (cadr expr))
198 (reduce #'p-add (p-eval-list (cddr expr)))))))
199 (*
200 (if (endp (cddr expr)) ;unary
201 (p-eval (cdr expr))
202 (reduce #'(lambda (p q) (poly-mul ring p q)) (p-eval-list (cdr expr)))))
203 (expt
204 (cond
205 ((member (cadr expr) vars :test #'equalp)
206 ;;Special handling of (expt var pow)
207 (let ((pos (position (cadr expr) vars :test #'equalp)))
208 (make-variable ring (length vars) pos (caddr expr))))
209 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
210 ;; Negative power means division in coefficient ring
211 ;; Non-integer power means non-polynomial coefficient
212 (coerce-coeff ring expr vars))
213 (t (poly-expt ring (p-eval (cadr expr)) (caddr expr)))))
214 (otherwise
215 (coerce-coeff ring expr vars)))))))
[55]216
217(defun spoly (ring f g)
218 "It yields the S-polynomial of polynomials F and G."
219 (declare (type poly f g))
220 (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
221 (mf (monom-div lcm (poly-lm f)))
222 (mg (monom-div lcm (poly-lm g))))
223 (declare (type monom mf mg))
224 (multiple-value-bind (c cf cg)
225 (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
226 (declare (ignore c))
227 (poly-sub
228 ring
229 (scalar-times-poly ring cg (monom-times-poly mf f))
[53]230 (scalar-times-poly ring cf (monom-times-poly mg g))))))
231
[55]232
233(defun poly-primitive-part (ring p)
234 "Divide polynomial P with integer coefficients by gcd of its
235coefficients and return the result."
236 (declare (type poly p))
237 (if (poly-zerop p)
238 (values p 1)
239 (let ((c (poly-content ring p)))
240 (values (make-poly-from-termlist (mapcar
241 #'(lambda (x)
242 (make-term (term-monom x)
243 (funcall (ring-div ring) (term-coeff x) c)))
244 (poly-termlist p))
245 (poly-sugar p))
246 c))))
247
248(defun poly-content (ring p)
249 "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
250to compute the greatest common divisor."
251 (declare (type poly p))
[57]252 (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
253
Note: See TracBrowser for help on using the repository browser.