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

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

* empty log message *

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