1 | ;;----------------------------------------------------------------
|
---|
2 | ;;; -*- Mode: Lisp -*-
|
---|
3 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
4 | ;;;
|
---|
5 | ;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
|
---|
6 | ;;;
|
---|
7 | ;;; This program is free software; you can redistribute it and/or modify
|
---|
8 | ;;; it under the terms of the GNU General Public License as published by
|
---|
9 | ;;; the Free Software Foundation; either version 2 of the License, or
|
---|
10 | ;;; (at your option) any later version.
|
---|
11 | ;;;
|
---|
12 | ;;; This program is distributed in the hope that it will be useful,
|
---|
13 | ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
14 | ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
15 | ;;; GNU General Public License for more details.
|
---|
16 | ;;;
|
---|
17 | ;;; You should have received a copy of the GNU General Public License
|
---|
18 | ;;; along with this program; if not, write to the Free Software
|
---|
19 | ;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
---|
20 | ;;;
|
---|
21 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
22 |
|
---|
23 | (defpackage "POLYNOMIAL"
|
---|
24 | (:use :cl :utils :monom :copy)
|
---|
25 | (:export "POLY"
|
---|
26 | "POLY-DIMENSION"
|
---|
27 | "POLY-TERMLIST"
|
---|
28 | "POLY-TERM-ORDER"
|
---|
29 | "POLY-INSERT-TERM"
|
---|
30 | "SCALAR-MULTIPLY-BY"
|
---|
31 | "SCALAR-DIVIDE-BY"
|
---|
32 | "LEADING-TERM"
|
---|
33 | "LEADING-MONOMIAL"
|
---|
34 | "LEADING-COEFFICIENT"
|
---|
35 | "SECOND-LEADING-TERM"
|
---|
36 | "SECOND-LEADING-MONOMIAL"
|
---|
37 | "SECOND-LEADING-COEFFICIENT"
|
---|
38 | "ADD-TO"
|
---|
39 | "ADD"
|
---|
40 | "SUBTRACT-FROM"
|
---|
41 | "SUBTRACT"
|
---|
42 | "CHANGE-TERM-ORDER"
|
---|
43 | "STANDARD-EXTENSION"
|
---|
44 | "STANDARD-EXTENSION-1"
|
---|
45 | "STANDARD-SUM"
|
---|
46 | "SATURATION-EXTENSION"
|
---|
47 | "ALIST->POLY"
|
---|
48 | "->INFIX"
|
---|
49 | "UNIVERSAL-EZGCD"
|
---|
50 | "S-POLYNOMIAL"
|
---|
51 | "POLY-CONTENT"
|
---|
52 | "POLY-PRIMITIVE-PART"
|
---|
53 | "SATURATION-EXTENSION-1"
|
---|
54 | "MAKE-POLY-VARIABLE"
|
---|
55 | "MAKE-POLY-CONSTANT"
|
---|
56 | "UNIVERSAL-EXPT"
|
---|
57 | "UNIVERSAL-EQUALP"
|
---|
58 | "POLY-LENGTH"
|
---|
59 | "POLY-P"
|
---|
60 | "+LIST-MARKER+"
|
---|
61 | "POLY-EVAL")
|
---|
62 | (:documentation "Implements polynomials. A polynomial is essentially
|
---|
63 | a mapping of monomials of the same degree to coefficients. The
|
---|
64 | momomials are ordered according to a monomial order."))
|
---|
65 |
|
---|
66 | (in-package :polynomial)
|
---|
67 |
|
---|
68 | (proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
|
---|
69 |
|
---|
70 | (defclass poly ()
|
---|
71 | ((dimension :initform nil
|
---|
72 | :initarg :dimension
|
---|
73 | :accessor poly-dimension
|
---|
74 | :documentation "Shared dimension of all terms, the number of variables")
|
---|
75 | (termlist :initform nil :initarg :termlist :accessor poly-termlist
|
---|
76 | :documentation "List of terms.")
|
---|
77 | (order :initform #'lex> :initarg :order :accessor poly-term-order
|
---|
78 | :documentation "Monomial/term order."))
|
---|
79 | (:default-initargs :dimension nil :termlist nil :order #'lex>)
|
---|
80 | (:documentation "A polynomial with a list of terms TERMLIST, ordered
|
---|
81 | according to term order ORDER, which defaults to LEX>."))
|
---|
82 |
|
---|
83 | (defmethod print-object ((self poly) stream)
|
---|
84 | (print-unreadable-object (self stream :type t :identity t)
|
---|
85 | (with-accessors ((dimension poly-dimension)
|
---|
86 | (termlist poly-termlist)
|
---|
87 | (order poly-term-order))
|
---|
88 | self
|
---|
89 | (format stream "DIMENSION=~A TERMLIST=~A ORDER=~A"
|
---|
90 | dimension termlist order))))
|
---|
91 |
|
---|
92 | (defgeneric change-term-order (self other)
|
---|
93 | (:documentation "Change term order of SELF to the term order of OTHER.")
|
---|
94 | (:method ((self poly) (other poly))
|
---|
95 | (unless (eq (poly-term-order self) (poly-term-order other))
|
---|
96 | (setf (poly-termlist self) (sort (poly-termlist self) (poly-term-order other))
|
---|
97 | (poly-term-order self) (poly-term-order other)))
|
---|
98 | self))
|
---|
99 |
|
---|
100 | (defgeneric poly-insert-term (self term)
|
---|
101 | (:documentation "Insert a term TERM into SELF before all other
|
---|
102 | terms. Order is not enforced.")
|
---|
103 | (:method ((self poly) (term term))
|
---|
104 | (cond ((null (poly-dimension self))
|
---|
105 | (setf (poly-dimension self) (monom-dimension term)))
|
---|
106 | (t (assert (= (poly-dimension self) (monom-dimension term)))))
|
---|
107 | (push term (poly-termlist self))
|
---|
108 | self))
|
---|
109 |
|
---|
110 | (defgeneric poly-append-term (self term)
|
---|
111 | (:documentation "Append a term TERM to SELF after all other terms. Order is not enforced.")
|
---|
112 | (:method ((self poly) (term term))
|
---|
113 | (cond ((null (poly-dimension self))
|
---|
114 | (setf (poly-dimension self) (monom-dimension term)))
|
---|
115 | (t (assert (= (poly-dimension self) (monom-dimension term)))))
|
---|
116 | (setf (cdr (last (poly-termlist self))) (list term))
|
---|
117 | self))
|
---|
118 |
|
---|
119 | (defun alist->poly (alist &aux (poly (make-instance 'poly)))
|
---|
120 | "It reads polynomial from an alist formatted as ( ... (exponents . coeff) ...).
|
---|
121 | It can be used to enter simple polynomials by hand, e.g the polynomial
|
---|
122 | in two variables, X and Y, given in standard notation as:
|
---|
123 |
|
---|
124 | 3*X^2*Y^3+2*Y+7
|
---|
125 |
|
---|
126 | can be entered as
|
---|
127 | (ALIST->POLY '(((2 3) . 3) ((0 1) . 2) ((0 0) . 7))).
|
---|
128 |
|
---|
129 | NOTE: The primary use is for low-level debugging of the package."
|
---|
130 | (dolist (x alist poly)
|
---|
131 | (poly-insert-term poly (make-instance 'term :exponents (car x) :coeff (cdr x)))))
|
---|
132 |
|
---|
133 | (defmethod update-instance-for-different-class :after ((old term) (new poly) &key)
|
---|
134 | "Converts OLD of class TERM to a NEW of class POLY, by making it into a 1-element TERMLIST."
|
---|
135 | (reinitialize-instance new
|
---|
136 | :dimension (monom-dimension old)
|
---|
137 | :termlist (list old)))
|
---|
138 |
|
---|
139 | (defmethod update-instance-for-different-class :after ((old monom) (new poly) &key)
|
---|
140 | "Converts OLD of class MONOM to a NEW of class POLY, by making it into a 1-element TERMLIST."
|
---|
141 | (reinitialize-instance new
|
---|
142 | :dimension (monom-dimension old)
|
---|
143 | :termlist (list (change-class old 'term))))
|
---|
144 |
|
---|
145 | (defmethod universal-equalp ((self poly) (other poly))
|
---|
146 | "Implements equality of polynomials."
|
---|
147 | (and (eql (poly-dimension self) (poly-dimension other))
|
---|
148 | (every #'universal-equalp (poly-termlist self) (poly-termlist other))
|
---|
149 | (eq (poly-term-order self) (poly-term-order other))))
|
---|
150 |
|
---|
151 | (defgeneric leading-term (object)
|
---|
152 | (:method ((self poly))
|
---|
153 | (car (poly-termlist self)))
|
---|
154 | (:documentation "The leading term of a polynomial, or NIL for zero polynomial."))
|
---|
155 |
|
---|
156 | (defgeneric second-leading-term (object)
|
---|
157 | (:method ((self poly))
|
---|
158 | (cadar (poly-termlist self)))
|
---|
159 | (:documentation "The second leading term of a polynomial, or NIL for a polynomial with at most one term."))
|
---|
160 |
|
---|
161 | (defgeneric leading-monomial (object)
|
---|
162 | (:method ((self poly))
|
---|
163 | (change-class (copy-instance (leading-term self)) 'monom))
|
---|
164 | (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
|
---|
165 |
|
---|
166 | (defgeneric second-leading-monomial (object)
|
---|
167 | (:method ((self poly))
|
---|
168 | (change-class (copy-instance (second-leading-term self)) 'monom))
|
---|
169 | (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
|
---|
170 |
|
---|
171 | (defgeneric leading-coefficient (object)
|
---|
172 | (:method ((self poly))
|
---|
173 | (term-coeff (leading-term self)))
|
---|
174 | (:documentation "The leading coefficient of a polynomial. It signals error for a zero polynomial."))
|
---|
175 |
|
---|
176 | (defgeneric second-leading-coefficient (object)
|
---|
177 | (:method ((self poly))
|
---|
178 | (term-coeff (second-leading-term self)))
|
---|
179 | (:documentation "The second leading coefficient of a polynomial. It
|
---|
180 | signals error for a polynomial with at most one term."))
|
---|
181 |
|
---|
182 | (defmethod universal-zerop ((self poly))
|
---|
183 | "Return T iff SELF is a zero polynomial."
|
---|
184 | (null (poly-termlist self)))
|
---|
185 |
|
---|
186 | (defgeneric poly-length (self)
|
---|
187 | (:documentation "Return the number of terms.")
|
---|
188 | (:method ((self poly))
|
---|
189 | (length (poly-termlist self))))
|
---|
190 |
|
---|
191 | (defgeneric scalar-multiply-by (self other)
|
---|
192 | (:documentation "Multiply vector SELF by a scalar OTHER.")
|
---|
193 | (:method ((self poly) other)
|
---|
194 | (mapc #'(lambda (term) (setf (term-coeff term) (multiply (term-coeff term) other)))
|
---|
195 | (poly-termlist self))
|
---|
196 | self))
|
---|
197 |
|
---|
198 | (defgeneric scalar-divide-by (self other)
|
---|
199 | (:documentation "Divide vector SELF by a scalar OTHER.")
|
---|
200 | (:method ((self poly) other)
|
---|
201 | (mapc #'(lambda (term) (setf (term-coeff term) (divide (term-coeff term) other)))
|
---|
202 | (poly-termlist self))
|
---|
203 | self))
|
---|
204 |
|
---|
205 | (defmethod multiply-by ((self poly) (other monom))
|
---|
206 | "Multiply a polynomial SELF by OTHER."
|
---|
207 | (mapc #'(lambda (term) (multiply-by term other))
|
---|
208 | (poly-termlist self))
|
---|
209 | self)
|
---|
210 |
|
---|
211 | (defmethod multiply-by ((self poly) (other term))
|
---|
212 | "Multiply a polynomial SELF by OTHER."
|
---|
213 | (mapc #'(lambda (term) (multiply-by term other))
|
---|
214 | (poly-termlist self))
|
---|
215 | self)
|
---|
216 |
|
---|
217 | (defmacro fast-add/subtract (p q order-fn add/subtract-fn uminus-fn)
|
---|
218 | "Return an expression which will efficiently adds/subtracts two
|
---|
219 | polynomials, P and Q. The addition/subtraction of coefficients is
|
---|
220 | performed by calling ADD/SUBTRACT-FN. If UMINUS-FN is supplied, it is
|
---|
221 | used to negate the coefficients of Q which do not have a corresponding
|
---|
222 | coefficient in P. The code implements an efficient algorithm to add
|
---|
223 | two polynomials represented as sorted lists of terms. The code
|
---|
224 | destroys both arguments, reusing the terms to build the result."
|
---|
225 | `(macrolet ((lc (x) `(term-coeff (car ,x))))
|
---|
226 | (do ((p ,p)
|
---|
227 | (q ,q)
|
---|
228 | r)
|
---|
229 | ((or (endp p) (endp q))
|
---|
230 | ;; NOTE: R contains the result in reverse order. Can it
|
---|
231 | ;; be more efficient to produce the terms in correct order?
|
---|
232 | (unless (endp q)
|
---|
233 | ;; Upon subtraction, we must change the sign of
|
---|
234 | ;; all coefficients in q
|
---|
235 | ,@(when uminus-fn
|
---|
236 | `((mapc #'(lambda (x) (setf x (funcall ,uminus-fn x))) q)))
|
---|
237 | (setf r (nreconc r q)))
|
---|
238 | (unless (endp p)
|
---|
239 | (setf r (nreconc r p)))
|
---|
240 | r)
|
---|
241 | (multiple-value-bind
|
---|
242 | (greater-p equal-p)
|
---|
243 | (funcall ,order-fn (car p) (car q))
|
---|
244 | (cond
|
---|
245 | (greater-p
|
---|
246 | (rotatef (cdr p) r p)
|
---|
247 | )
|
---|
248 | (equal-p
|
---|
249 | (let ((s (funcall ,add/subtract-fn (lc p) (lc q))))
|
---|
250 | (cond
|
---|
251 | ((universal-zerop s)
|
---|
252 | (setf p (cdr p))
|
---|
253 | )
|
---|
254 | (t
|
---|
255 | (setf (lc p) s)
|
---|
256 | (rotatef (cdr p) r p))))
|
---|
257 | (setf q (cdr q))
|
---|
258 | )
|
---|
259 | (t
|
---|
260 | ;;Negate the term of Q if UMINUS provided, signallig
|
---|
261 | ;;that we are doing subtraction
|
---|
262 | ,(when uminus-fn
|
---|
263 | `(setf (lc q) (funcall ,uminus-fn (lc q))))
|
---|
264 | (rotatef (cdr q) r q))))
|
---|
265 | ;;(format t "P:~A~%" p)
|
---|
266 | ;;(format t "Q:~A~%" q)
|
---|
267 | ;;(format t "R:~A~%" r)
|
---|
268 | )))
|
---|
269 |
|
---|
270 |
|
---|
271 |
|
---|
272 | (defgeneric add-to (self other)
|
---|
273 | (:documentation "Add OTHER to SELF.")
|
---|
274 | (:method ((self number) (other number))
|
---|
275 | (+ self other))
|
---|
276 | (:method ((self poly) (other number))
|
---|
277 | (add-to self (make-poly-constant (poly-dimension self) other)))
|
---|
278 | (:method ((self number) (other poly))
|
---|
279 | (add-to (make-poly-constant (poly-dimension other) self) other)))
|
---|
280 |
|
---|
281 |
|
---|
282 | (defgeneric subtract-from (self other)
|
---|
283 | (:documentation "Subtract OTHER from SELF.")
|
---|
284 | (:method ((self number) (other number))
|
---|
285 | (- self other))
|
---|
286 | (:method ((self poly) (other number))
|
---|
287 | (subtract-from self (make-poly-constant (poly-dimension self) other))))
|
---|
288 |
|
---|
289 |
|
---|
290 | #|
|
---|
291 | (defmacro def-add/subtract-method (add/subtract-method-name
|
---|
292 | uminus-method-name
|
---|
293 | &optional
|
---|
294 | (doc-string nil doc-string-supplied-p))
|
---|
295 | "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
|
---|
296 | `(defmethod ,add/subtract-method-name ((self poly) (other poly))
|
---|
297 | ,@(when doc-string-supplied-p `(,doc-string))
|
---|
298 | ;; Ensure orders are compatible
|
---|
299 | (change-term-order other self)
|
---|
300 | (setf (poly-termlist self) (fast-add/subtract
|
---|
301 | (poly-termlist self) (poly-termlist other)
|
---|
302 | (poly-term-order self)
|
---|
303 | #',add/subtract-method-name
|
---|
304 | ,(when uminus-method-name `(function ,uminus-method-name))))
|
---|
305 | self))
|
---|
306 |
|
---|
307 | (eval-when (:load-toplevel :execute)
|
---|
308 |
|
---|
309 | (def-add/subtract-method add-to nil
|
---|
310 | "Adds to polynomial SELF another polynomial OTHER.
|
---|
311 | This operation destructively modifies both polynomials.
|
---|
312 | The result is stored in SELF. This implementation does
|
---|
313 | no consing, entirely reusing the sells of SELF and OTHER.")
|
---|
314 |
|
---|
315 | (def-add/subtract-method subtract-from unary-minus
|
---|
316 | "Subtracts from polynomial SELF another polynomial OTHER.
|
---|
317 | This operation destructively modifies both polynomials.
|
---|
318 | The result is stored in SELF. This implementation does
|
---|
319 | no consing, entirely reusing the sells of SELF and OTHER.")
|
---|
320 | )
|
---|
321 |
|
---|
322 | |#
|
---|
323 |
|
---|
324 | (defmethod unary-minus ((self poly))
|
---|
325 | "Destructively modifies the coefficients of the polynomial SELF,
|
---|
326 | by changing their sign."
|
---|
327 | (mapc #'unary-minus (poly-termlist self))
|
---|
328 | self)
|
---|
329 |
|
---|
330 | (defun add-termlists (p q order-fn)
|
---|
331 | "Destructively adds two termlists P and Q ordered according to ORDER-FN."
|
---|
332 | (fast-add/subtract p q order-fn #'add-to nil))
|
---|
333 |
|
---|
334 | (defun subtract-termlists (p q order-fn)
|
---|
335 | "Destructively subtracts two termlists P and Q ordered according to ORDER-FN."
|
---|
336 | (fast-add/subtract p q order-fn #'subtract-from #'unary-minus))
|
---|
337 |
|
---|
338 | (defmethod add-to ((self poly) (other poly))
|
---|
339 | "Adds to polynomial SELF another polynomial OTHER.
|
---|
340 | This operation destructively modifies both polynomials.
|
---|
341 | The result is stored in SELF. This implementation does
|
---|
342 | no consing, entirely reusing the sells of SELF and OTHER."
|
---|
343 | (change-term-order other self)
|
---|
344 | (setf (poly-termlist self) (add-termlists
|
---|
345 | (poly-termlist self) (poly-termlist other)
|
---|
346 | (poly-term-order self)))
|
---|
347 | self)
|
---|
348 |
|
---|
349 |
|
---|
350 | (defmethod subtract-from ((self poly) (other poly))
|
---|
351 | "Subtracts from polynomial SELF another polynomial OTHER.
|
---|
352 | This operation destructively modifies both polynomials.
|
---|
353 | The result is stored in SELF. This implementation does
|
---|
354 | no consing, entirely reusing the sells of SELF and OTHER."
|
---|
355 | (change-term-order other self)
|
---|
356 | (setf (poly-termlist self) (subtract-termlists
|
---|
357 | (poly-termlist self) (poly-termlist other)
|
---|
358 | (poly-term-order self)))
|
---|
359 | self)
|
---|
360 |
|
---|
361 | (defmacro multiply-term-by-termlist-dropping-zeros (term termlist
|
---|
362 | &optional (reverse-arg-order-P nil))
|
---|
363 | "Multiplies term TERM by a list of term, TERMLIST.
|
---|
364 | Takes into accound divisors of zero in the ring, by
|
---|
365 | deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
|
---|
366 | is T, change the order of arguments; this may be important
|
---|
367 | if we extend the package to non-commutative rings."
|
---|
368 | `(mapcan #'(lambda (other-term)
|
---|
369 | (let ((prod (multiply
|
---|
370 | ,@(cond
|
---|
371 | (reverse-arg-order-p
|
---|
372 | `(other-term ,term))
|
---|
373 | (t
|
---|
374 | `(,term other-term))))))
|
---|
375 | (cond
|
---|
376 | ((universal-zerop prod) nil)
|
---|
377 | (t (list prod)))))
|
---|
378 | ,termlist))
|
---|
379 |
|
---|
380 | (defun multiply-termlists (p q order-fn)
|
---|
381 | "A version of polynomial multiplication, operating
|
---|
382 | directly on termlists."
|
---|
383 | (cond
|
---|
384 | ((or (endp p) (endp q))
|
---|
385 | ;;p or q is 0 (represented by NIL)
|
---|
386 | nil)
|
---|
387 | ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
|
---|
388 | ((endp (cdr p))
|
---|
389 | (multiply-term-by-termlist-dropping-zeros (car p) q))
|
---|
390 | ((endp (cdr q))
|
---|
391 | (multiply-term-by-termlist-dropping-zeros (car q) p t))
|
---|
392 | (t
|
---|
393 | (cons (multiply (car p) (car q))
|
---|
394 | (add-termlists
|
---|
395 | (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
|
---|
396 | (multiply-termlists (cdr p) q order-fn)
|
---|
397 | order-fn)))))
|
---|
398 |
|
---|
399 | (defmethod multiply-by ((self poly) (other poly))
|
---|
400 | (change-term-order other self)
|
---|
401 | (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
|
---|
402 | (poly-termlist other)
|
---|
403 | (poly-term-order self)))
|
---|
404 | self)
|
---|
405 |
|
---|
406 | (defgeneric add-2 (object1 object2)
|
---|
407 | (:documentation "Non-destructively add OBJECT1 to OBJECT2.")
|
---|
408 | (:method ((object1 t) (object2 t))
|
---|
409 | (add-to (copy-instance object1) (copy-instance object2))))
|
---|
410 |
|
---|
411 | (defun add (&rest summands)
|
---|
412 | "Non-destructively adds list SUMMANDS."
|
---|
413 | (cond ((endp summands) 0)
|
---|
414 | (t (reduce #'add-2 summands))))
|
---|
415 |
|
---|
416 | (defun subtract (minuend &rest subtrahends)
|
---|
417 | "Non-destructively subtract MINUEND and SUBTRAHENDS."
|
---|
418 | (cond ((endp subtrahends) (unary-minus minuend))
|
---|
419 | (t (subtract-from (copy-instance minuend) (reduce #'add subtrahends)))))
|
---|
420 |
|
---|
421 | (defmethod left-tensor-product-by ((self poly) (other monom))
|
---|
422 | (setf (poly-termlist self)
|
---|
423 | (mapcan #'(lambda (term)
|
---|
424 | (let ((prod (left-tensor-product-by term other)))
|
---|
425 | (cond
|
---|
426 | ((universal-zerop prod) nil)
|
---|
427 | (t (list prod)))))
|
---|
428 | (poly-termlist self)))
|
---|
429 | (incf (poly-dimension self) (monom-dimension other))
|
---|
430 | self)
|
---|
431 |
|
---|
432 | (defmethod right-tensor-product-by ((self poly) (other monom))
|
---|
433 | (setf (poly-termlist self)
|
---|
434 | (mapcan #'(lambda (term)
|
---|
435 | (let ((prod (right-tensor-product-by term other)))
|
---|
436 | (cond
|
---|
437 | ((universal-zerop prod) nil)
|
---|
438 | (t (list prod)))))
|
---|
439 | (poly-termlist self)))
|
---|
440 | (incf (poly-dimension self) (monom-dimension other))
|
---|
441 | self)
|
---|
442 |
|
---|
443 |
|
---|
444 | (defun standard-extension (plist &aux (k (length plist)) (i 0))
|
---|
445 | "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
|
---|
446 | is a list of polynomials. Destructively modifies PLIST elements."
|
---|
447 | (mapc #'(lambda (poly)
|
---|
448 | (left-tensor-product-by
|
---|
449 | poly
|
---|
450 | (prog1
|
---|
451 | (make-monom-variable k i)
|
---|
452 | (incf i))))
|
---|
453 | plist))
|
---|
454 |
|
---|
455 | (defun standard-extension-1 (plist
|
---|
456 | &aux
|
---|
457 | (plist (standard-extension plist))
|
---|
458 | (nvars (poly-dimension (car plist))))
|
---|
459 | "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
|
---|
460 | Firstly, new K variables U1, U2, ..., UK, are inserted into each
|
---|
461 | polynomial. Subsequently, P1, P2, ..., PK are destructively modified
|
---|
462 | tantamount to replacing PI with UI*PI-1. It assumes that all
|
---|
463 | polynomials have the same dimension, and only the first polynomial
|
---|
464 | is examined to determine this dimension."
|
---|
465 | ;; Implementation note: we use STANDARD-EXTENSION and then subtract
|
---|
466 | ;; 1 from each polynomial; since UI*PI has no constant term,
|
---|
467 | ;; we just need to append the constant term at the end
|
---|
468 | ;; of each termlist.
|
---|
469 | (flet ((subtract-1 (p)
|
---|
470 | (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
|
---|
471 | (setf plist (mapc #'subtract-1 plist)))
|
---|
472 | plist)
|
---|
473 |
|
---|
474 |
|
---|
475 | (defun standard-sum (plist
|
---|
476 | &aux
|
---|
477 | (plist (standard-extension plist))
|
---|
478 | (nvars (poly-dimension (car plist))))
|
---|
479 | "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
|
---|
480 | Firstly, new K variables, U1, U2, ..., UK, are inserted into each
|
---|
481 | polynomial. Subsequently, P1, P2, ..., PK are destructively modified
|
---|
482 | tantamount to replacing PI with UI*PI, and the resulting polynomials
|
---|
483 | are added. Finally, 1 is subtracted. It should be noted that the term
|
---|
484 | order is not modified, which is equivalent to using a lexicographic
|
---|
485 | order on the first K variables."
|
---|
486 | (flet ((subtract-1 (p)
|
---|
487 | (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
|
---|
488 | (subtract-1
|
---|
489 | (make-instance
|
---|
490 | 'poly
|
---|
491 | :termlist (apply #'nconc (mapcar #'poly-termlist plist))))))
|
---|
492 |
|
---|
493 | (defgeneric universal-ezgcd (x y)
|
---|
494 | (:documentation "Solves the diophantine system: X=C*X1, Y=C*X2,
|
---|
495 | C=GCD(X,Y). It returns C, X1 and Y1. The result may be obtained by
|
---|
496 | the Euclidean algorithm.")
|
---|
497 | (:method ((x integer) (y integer)
|
---|
498 | &aux (c (gcd x y)))
|
---|
499 | (values c (/ x c) (/ y c)))
|
---|
500 | )
|
---|
501 |
|
---|
502 | (defgeneric s-polynomial (object1 object2)
|
---|
503 | (:documentation "Yields the S-polynomial of OBJECT1 and OBJECT2.")
|
---|
504 | (:method ((f poly) (g poly))
|
---|
505 | (let* ((lcm (universal-lcm (leading-monomial f) (leading-monomial g)))
|
---|
506 | (mf (divide lcm (leading-monomial f)))
|
---|
507 | (mg (divide lcm (leading-monomial g))))
|
---|
508 | (multiple-value-bind (c cf cg)
|
---|
509 | (universal-ezgcd (leading-coefficient f) (leading-coefficient g))
|
---|
510 | (declare (ignore c))
|
---|
511 | (subtract
|
---|
512 | (multiply f (change-class mf 'term :coeff cg))
|
---|
513 | (multiply g (change-class mg 'term :coeff cf)))))))
|
---|
514 |
|
---|
515 | (defgeneric poly-content (object)
|
---|
516 | (:documentation "Greatest common divisor of the coefficients of the polynomial object OBJECT.")
|
---|
517 | (:method ((self poly))
|
---|
518 | (reduce #'universal-gcd
|
---|
519 | (mapcar #'term-coeff (rest (poly-termlist self)))
|
---|
520 | :initial-value (leading-coefficient self))))
|
---|
521 |
|
---|
522 | (defun poly-primitive-part (object)
|
---|
523 | "Divide polynomial OBJECT by gcd of its
|
---|
524 | coefficients. Return the resulting polynomial."
|
---|
525 | (scalar-divide-by object (poly-content object)))
|
---|
526 |
|
---|
527 | (defun poly-insert-variables (self k)
|
---|
528 | (left-tensor-product-by self (make-instance 'monom :dimension k)))
|
---|
529 |
|
---|
530 | (defun saturation-extension (f plist &aux (k (length plist)))
|
---|
531 | "Calculate [F', U1*P1-1,U2*P2-1,...,UK*PK-1], where
|
---|
532 | PLIST=[P1,P2,...,PK] and F' is F with variables U1,U2,...,UK inserted
|
---|
533 | as first K variables. It destructively modifies F and PLIST."
|
---|
534 | (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
|
---|
535 | (standard-extension-1 plist)))
|
---|
536 |
|
---|
537 | (defun polysaturation-extension (f plist &aux (k (length plist)))
|
---|
538 | "Calculate [F', U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]
|
---|
539 | and F' is F with variables U1,U2,...,UK inserted as first K
|
---|
540 | variables. It destructively modifies F and PLIST."
|
---|
541 | (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
|
---|
542 | (list (standard-sum plist))))
|
---|
543 |
|
---|
544 | (defun saturation-extension-1 (f p)
|
---|
545 | "Given family of polynomials F and a polynomial P, calculate [F',
|
---|
546 | U*P-1], where F' is F with variable inserted as the first variable. It
|
---|
547 | destructively modifies F and P."
|
---|
548 | (polysaturation-extension f (list p)))
|
---|
549 |
|
---|
550 | (defmethod multiply-by ((object1 number) (object2 poly))
|
---|
551 | (scalar-multiply-by (copy-instance object2) object1))
|
---|
552 |
|
---|
553 | (defun make-poly-variable (nvars pos &optional (power 1))
|
---|
554 | (change-class (make-monom-variable nvars pos power) 'poly))
|
---|
555 |
|
---|
556 | (defun make-poly-constant (nvars coeff)
|
---|
557 | (change-class (make-term-constant nvars coeff) 'poly))
|
---|
558 |
|
---|
559 | (defgeneric universal-expt (x y)
|
---|
560 | (:documentation "Raises X to power Y.")
|
---|
561 | (:method ((x number) (y integer)) (expt x y))
|
---|
562 | (:method ((x t) (y integer))
|
---|
563 | (declare (type fixnum y))
|
---|
564 | (cond
|
---|
565 | ((minusp y) (error "universal-expt: Negative exponent."))
|
---|
566 | ((universal-zerop x) (if (zerop y) 1))
|
---|
567 | (t
|
---|
568 | (do ((k 1 (ash k 1))
|
---|
569 | (q x (multiply q q)) ;keep squaring
|
---|
570 | (p 1 (if (not (zerop (logand k y))) (multiply p q) p)))
|
---|
571 | ((> k y) p)
|
---|
572 | (declare (fixnum k)))))))
|
---|
573 |
|
---|
574 | (defgeneric poly-p (object)
|
---|
575 | (:documentation "Checks if an object is a polynomial.")
|
---|
576 | (:method ((self poly)) t)
|
---|
577 | (:method ((self t)) nil))
|
---|
578 |
|
---|
579 | (defmethod ->infix :before ((self poly) &optional vars)
|
---|
580 | "Ensures that the number of variables in VARS maches the polynomial dimension of the
|
---|
581 | polynomial SELF."
|
---|
582 | (assert (= (length vars) (poly-dimension self))))
|
---|
583 |
|
---|
584 | (defmethod ->infix ((self poly) &optional vars)
|
---|
585 | "Converts a polynomial SELF to a sexp."
|
---|
586 | (cons '+ (mapcar #'(lambda (x) (->infix x vars))
|
---|
587 | (poly-termlist self))))
|
---|
588 |
|
---|
589 | (defparameter +list-marker+ :[
|
---|
590 | "A sexp with this head is considered a list of polynomials.")
|
---|
591 |
|
---|
592 | (defmethod ->infix ((self cons) &optional vars)
|
---|
593 | (assert (eql (car self) +list-marker+))
|
---|
594 | (cons +list-marker+ (mapcar #'(lambda (p) (->infix p vars)) (cdr self))))
|
---|
595 |
|
---|
596 |
|
---|
597 | (defun poly-eval (expr vars order)
|
---|
598 | "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
|
---|
599 | variables VARS. Return the resulting polynomial or list of
|
---|
600 | polynomials. Standard arithmetical operators in form EXPR are
|
---|
601 | replaced with their analogues in the ring of polynomials, and the
|
---|
602 | resulting expression is evaluated, resulting in a polynomial or a list
|
---|
603 | of polynomials in internal form. A similar operation in another computer
|
---|
604 | algebra system could be called 'expand' or so."
|
---|
605 | (labels ((p-eval (p) (poly-eval p vars order))
|
---|
606 | (p-eval-list (plist) (mapcar #'p-eval plist)))
|
---|
607 | (cond
|
---|
608 | ((eq expr 0)
|
---|
609 | (make-instance 'poly :dimension (length vars)))
|
---|
610 | ((member expr vars :test #'equalp)
|
---|
611 | (let ((pos (position expr vars :test #'equalp)))
|
---|
612 | (make-poly-variable (length vars) pos)))
|
---|
613 | ((atom expr)
|
---|
614 | (make-poly-constant (length vars) expr))
|
---|
615 | ((eq (car expr) +list-marker+)
|
---|
616 | (cons +list-marker+ (p-eval-list (cdr expr))))
|
---|
617 | (t
|
---|
618 | (case (car expr)
|
---|
619 | (+ (reduce #'add (p-eval-list (cdr expr))))
|
---|
620 | (- (apply #'subtract (p-eval-list (cdr expr))))
|
---|
621 | (*
|
---|
622 | (if (endp (cddr expr)) ;unary
|
---|
623 | (p-eval (cadr expr))
|
---|
624 | (reduce #'multiply (p-eval-list (cdr expr)))))
|
---|
625 | (/
|
---|
626 | ;; A polynomial can be divided by a scalar
|
---|
627 | (cond
|
---|
628 | ((endp (cddr expr))
|
---|
629 | ;; A special case (/ ?), the inverse
|
---|
630 | (divide (cadr expr)))
|
---|
631 | (t
|
---|
632 | (let ((num (p-eval (cadr expr)))
|
---|
633 | (denom-inverse (apply #'divide (mapcar #'p-eval (cddr expr)))))
|
---|
634 | (multiply denom-inverse num)))))
|
---|
635 | (expt
|
---|
636 | (cond
|
---|
637 | ((member (cadr expr) vars :test #'equalp)
|
---|
638 | ;;Special handling of (expt var pow)
|
---|
639 | (let ((pos (position (cadr expr) vars :test #'equalp)))
|
---|
640 | (make-poly-variable (length vars) pos (caddr expr))))
|
---|
641 | ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
|
---|
642 | ;; Negative power means division in coefficient ring
|
---|
643 | ;; Non-integer power means non-polynomial coefficient
|
---|
644 | expr)
|
---|
645 | (t (universal-expt (p-eval (cadr expr)) (caddr expr)))))
|
---|
646 | (otherwise
|
---|
647 | (error "Cannot evaluate as polynomial: ~A" expr)))))))
|
---|