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

Last change on this file since 4227 was 4219, checked in by Marek Rychlik, 8 years ago

* empty log message *

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