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

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

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