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

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

* empty log message *

File size: 25.1 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 :ring)
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(defmacro def-add/subtract-method (add/subtract-method-name
310 uminus-method-name
311 &optional
312 (doc-string nil doc-string-supplied-p))
313 "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
314 `(defmethod ,add/subtract-method-name ((self poly) (other poly))
315 ,@(when doc-string-supplied-p `(,doc-string))
316 ;; Ensure orders are compatible
317 (change-term-order other self)
318 (setf (poly-termlist self) (fast-add/subtract
319 (poly-termlist self) (poly-termlist other)
320 (poly-term-order self)
321 #',add/subtract-method-name
322 ,(when uminus-method-name `(function ,uminus-method-name))))
323 self))
324
325(eval-when (:load-toplevel :execute)
326
327 (def-add/subtract-method add-to nil
328 "Adds to polynomial SELF another polynomial OTHER.
329This operation destructively modifies both polynomials.
330The result is stored in SELF. This implementation does
331no consing, entirely reusing the sells of SELF and OTHER.")
332
333 (def-add/subtract-method subtract-from unary-minus
334 "Subtracts from polynomial SELF another polynomial OTHER.
335This operation destructively modifies both polynomials.
336The result is stored in SELF. This implementation does
337no consing, entirely reusing the sells of SELF and OTHER.")
338 )
339
340|#
341
342(defmethod unary-minus ((self poly))
343 "Destructively modifies the coefficients of the polynomial SELF,
344by changing their sign."
345 (mapc #'unary-minus (poly-termlist self))
346 self)
347
348(defun add-termlists (p q order-fn)
349 "Destructively adds two termlists P and Q ordered according to ORDER-FN."
350 (fast-add/subtract p q order-fn #'add-to nil))
351
352(defun subtract-termlists (p q order-fn)
353 "Destructively subtracts two termlists P and Q ordered according to ORDER-FN."
354 (fast-add/subtract p q order-fn #'subtract-from #'unary-minus))
355
356(defmethod add-to ((self poly) (other poly) &aux (other-copy (copy-instance other)))
357 "Adds to polynomial SELF another polynomial OTHER.
358This operation destructively modifies both polynomials.
359The result is stored in SELF. This implementation does
360no consing, entirely reusing the sells of SELF and OTHER."
361 (change-term-order other-copy self)
362 (setf (poly-termlist self) (add-termlists
363 (poly-termlist self) (poly-termlist other-copy)
364 (poly-term-order self)))
365 self)
366
367
368(defmethod subtract-from ((self poly) (other poly) &aux (other-copy (copy-instance other)))
369 "Subtracts from polynomial SELF another polynomial OTHER.
370This operation destructively modifies both polynomials.
371The result is stored in SELF. This implementation does
372no consing, entirely reusing the sells of SELF and OTHER."
373 (change-term-order other-copy self)
374 (setf (poly-termlist self) (subtract-termlists
375 (poly-termlist self) (poly-termlist other-copy)
376 (poly-term-order self)))
377 self)
378
379
380(defmethod add-to ((self poly) (other term) &aux (other-copy (copy-instance other)))
381 "Adds to a polynomial SELF a term OTHER. The term OTHER is not
382modified."
383 (add-to self (change-class other-copy 'poly)))
384
385(defmethod subtract-from ((self poly) (other term) &aux (other-copy (copy-instance other)))
386 "Subtracts from a polynomial SELF a term OTHER. The term OTHER is not
387modified."
388 (subtract-from self (change-class other-copy 'poly)))
389
390
391(defmacro multiply-term-by-termlist-dropping-zeros (term termlist
392 &optional (reverse-arg-order-P nil))
393 "Multiplies term TERM by a list of term, TERMLIST.
394Takes into accound divisors of zero in the ring, by
395deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
396is T, change the order of arguments; this may be important
397if we extend the package to non-commutative rings."
398 `(mapcan #'(lambda (other-term)
399 (let ((prod (multiply
400 ,@(cond
401 (reverse-arg-order-p
402 `(other-term ,term))
403 (t
404 `(,term other-term))))))
405 (cond
406 ((universal-zerop prod) nil)
407 (t (list prod)))))
408 ,termlist))
409
410(defun multiply-termlists (p q order-fn)
411 "A version of polynomial multiplication, operating
412directly on termlists."
413 (cond
414 ((or (endp p) (endp q))
415 ;;p or q is 0 (represented by NIL)
416 nil)
417 ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
418 ((endp (cdr p))
419 (multiply-term-by-termlist-dropping-zeros (car p) q))
420 ((endp (cdr q))
421 (multiply-term-by-termlist-dropping-zeros (car q) p t))
422 (t
423 (cons (multiply (car p) (car q))
424 (add-termlists
425 (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
426 (multiply-termlists (cdr p) q order-fn)
427 order-fn)))))
428
429(defmethod multiply-by ((self poly) (other poly))
430 (change-term-order other self)
431 (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
432 (poly-termlist other)
433 (poly-term-order self)))
434 self)
435
436(defun add (summand &rest more-summands)
437 "Successively Adds to SUMMAND the elements of MORE-SUMMANDS."
438 (reduce #'add-to more-summands :initial-value summand))
439
440(defun subtract (minuend &rest subtrahends)
441 "Non-destructively subtract MINUEND and SUBTRAHENDS."
442 (cond ((endp subtrahends) (unary-minus minuend))
443 (t (subtract-from (copy-instance minuend) (reduce #'add subtrahends)))))
444
445(defmethod left-tensor-product-by ((self poly) (other monom))
446 (setf (poly-termlist self)
447 (mapcan #'(lambda (term)
448 (let ((prod (left-tensor-product-by term other)))
449 (cond
450 ((universal-zerop prod) nil)
451 (t (list prod)))))
452 (poly-termlist self)))
453 (incf (poly-dimension self) (monom-dimension other))
454 self)
455
456(defmethod right-tensor-product-by ((self poly) (other monom))
457 (setf (poly-termlist self)
458 (mapcan #'(lambda (term)
459 (let ((prod (right-tensor-product-by term other)))
460 (cond
461 ((universal-zerop prod) nil)
462 (t (list prod)))))
463 (poly-termlist self)))
464 (incf (poly-dimension self) (monom-dimension other))
465 self)
466
467
468(defun standard-extension (plist &aux (k (length plist)) (i 0))
469 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
470is a list of polynomials. Destructively modifies PLIST elements."
471 (mapc #'(lambda (poly)
472 (left-tensor-product-by
473 poly
474 (prog1
475 (make-monom-variable k i)
476 (incf i))))
477 plist))
478
479(defun standard-extension-1 (plist
480 &aux
481 (plist (standard-extension plist))
482 (nvars (poly-dimension (car plist))))
483 "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
484Firstly, new K variables U1, U2, ..., UK, are inserted into each
485polynomial. Subsequently, P1, P2, ..., PK are destructively modified
486tantamount to replacing PI with UI*PI-1. It assumes that all
487polynomials have the same dimension, and only the first polynomial
488is examined to determine this dimension."
489 ;; Implementation note: we use STANDARD-EXTENSION and then subtract
490 ;; 1 from each polynomial; since UI*PI has no constant term,
491 ;; we just need to append the constant term at the end
492 ;; of each termlist.
493 (flet ((subtract-1 (p)
494 (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
495 (setf plist (mapc #'subtract-1 plist)))
496 plist)
497
498
499(defun standard-sum (plist
500 &aux
501 (plist (standard-extension plist))
502 (nvars (poly-dimension (car plist))))
503 "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
504Firstly, new K variables, U1, U2, ..., UK, are inserted into each
505polynomial. Subsequently, P1, P2, ..., PK are destructively modified
506tantamount to replacing PI with UI*PI, and the resulting polynomials
507are added. Finally, 1 is subtracted. It should be noted that the term
508order is not modified, which is equivalent to using a lexicographic
509order on the first K variables."
510 (flet ((subtract-1 (p)
511 (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
512 (subtract-1
513 (make-instance
514 'poly
515 :termlist (apply #'nconc (mapcar #'poly-termlist plist))))))
516
517(defgeneric s-polynomial (object1 object2)
518 (:documentation "Yields the S-polynomial of OBJECT1 and OBJECT2.")
519 (:method ((f poly) (g poly))
520 (let* ((lcm (universal-lcm (leading-monomial f) (leading-monomial g)))
521 (mf (divide lcm (leading-monomial f)))
522 (mg (divide lcm (leading-monomial g))))
523 (multiple-value-bind (c cf cg)
524 (universal-ezgcd (leading-coefficient f) (leading-coefficient g))
525 (declare (ignore c))
526 (subtract
527 (multiply f (change-class mf 'term :coeff cg))
528 (multiply g (change-class mg 'term :coeff cf)))))))
529
530(defgeneric poly-content (object)
531 (:documentation "Greatest common divisor of the coefficients of the polynomial object OBJECT.")
532 (:method ((self poly))
533 (reduce #'universal-gcd
534 (mapcar #'term-coeff (rest (poly-termlist self)))
535 :initial-value (leading-coefficient self))))
536
537(defun poly-primitive-part (object)
538 "Divide polynomial OBJECT by gcd of its
539coefficients. Return the resulting polynomial."
540 (scalar-divide-by object (poly-content object)))
541
542(defun poly-insert-variables (self k)
543 (left-tensor-product-by self (make-instance 'monom :dimension k)))
544
545(defun saturation-extension (f plist &aux (k (length plist)))
546 "Calculate [F', U1*P1-1,U2*P2-1,...,UK*PK-1], where
547PLIST=[P1,P2,...,PK] and F' is F with variables U1,U2,...,UK inserted
548as first K variables. It destructively modifies F and PLIST."
549 (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
550 (standard-extension-1 plist)))
551
552(defun polysaturation-extension (f plist &aux (k (length plist)))
553 "Calculate [F', U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]
554and F' is F with variables U1,U2,...,UK inserted as first K
555variables. It destructively modifies F and PLIST."
556 (nconc (mapc #'(lambda (x) (poly-insert-variables x k)) f)
557 (list (standard-sum plist))))
558
559(defun saturation-extension-1 (f p)
560 "Given family of polynomials F and a polynomial P, calculate [F',
561U*P-1], where F' is F with variable inserted as the first variable. It
562destructively modifies F and P."
563 (polysaturation-extension f (list p)))
564
565;; (defmethod multiply-by ((object1 number) (object2 poly))
566;; (scalar-multiply-by (copy-instance object2) object1))
567
568(defmethod multiply-by ((object1 poly) (object2 number))
569 (scalar-multiply-by (copy-instance object1) object2))
570
571(defun make-poly-variable (nvars pos &optional (power 1))
572 (change-class (make-monom-variable nvars pos power) 'poly))
573
574(defun make-poly-constant (nvars coeff)
575 (change-class (make-term-constant nvars coeff) 'poly))
576
577(defgeneric universal-expt (x y)
578 (:documentation "Raises X to power Y.")
579 (:method ((x number) (y integer)) (expt x y))
580 (:method ((x t) (y integer))
581 (declare (type fixnum y))
582 (cond
583 ((minusp y) (error "universal-expt: Negative exponent."))
584 ((universal-zerop x) (if (zerop y) 1))
585 (t
586 (do ((k 1 (ash k 1))
587 (q x (multiply q q)) ;keep squaring
588 (p (make-unit-for x) (if (not (zerop (logand k y))) (multiply p q) p)))
589 ((> k y) p)
590 (declare (fixnum k)))))))
591
592(defgeneric poly-p (object)
593 (:documentation "Checks if an object is a polynomial.")
594 (:method ((self poly)) t)
595 (:method ((self t)) nil))
596
597(defmethod ->sexp :before ((self poly) &optional vars)
598 "Ensures that the number of variables in VARS maches the polynomial dimension of the
599polynomial SELF."
600 (with-slots (dimension)
601 self
602 (assert (= (length vars) dimension)
603 nil
604 "Number of variables ~S does not match the dimension ~S"
605 vars dimension)))
606
607(defmethod ->sexp ((self poly) &optional vars)
608 "Converts a polynomial SELF to a sexp."
609 (let ((m (mapcar #'(lambda (x) (->sexp x vars))
610 (poly-termlist self))))
611 (cond ((endp m) 0)
612 ((endp (cdr m)) (car m))
613 (t (cons '+ m)))))
614
615(defparameter +list-marker+ :[
616 "A sexp with this head is considered a list of polynomials.")
617
618(defmethod ->sexp ((self cons) &optional vars)
619 (assert (eql (car self) +list-marker+))
620 (cons +list-marker+ (mapcar #'(lambda (p) (->sexp p vars)) (cdr self))))
621
622
623(defun poly-eval (expr vars order)
624 "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
625variables VARS. Return the resulting polynomial or list of
626polynomials. Standard arithmetical operators in form EXPR are
627replaced with their analogues in the ring of polynomials, and the
628resulting expression is evaluated, resulting in a polynomial or a list
629of polynomials in internal form. A similar operation in another computer
630algebra system could be called 'expand' or so."
631 (labels ((p-eval (p) (poly-eval p vars order))
632 (p-eval-list (plist) (mapcar #'p-eval plist)))
633 (cond
634 ((eq expr 0)
635 (make-instance 'poly :dimension (length vars)))
636 ((member expr vars :test #'equalp)
637 (let ((pos (position expr vars :test #'equalp)))
638 (make-poly-variable (length vars) pos)))
639 ((integerp expr)
640 (make-poly-constant (length vars) (make-instance 'integer-ring :value expr)))
641 ((eq (car expr) +list-marker+)
642 (cons +list-marker+ (p-eval-list (cdr expr))))
643 (t
644 (case (car expr)
645 (+ (reduce #'add (p-eval-list (cdr expr))))
646 (- (apply #'subtract (p-eval-list (cdr expr))))
647 (*
648 (if (endp (cddr expr)) ;unary
649 (p-eval (cadr expr))
650 (apply #'multiply (p-eval-list (cdr expr)))))
651 (/
652 ;; A polynomial can be divided by a scalar
653 (cond
654 ((endp (cddr expr))
655 ;; A special case (/ ?), the inverse
656 (divide (cadr expr)))
657 (t
658 (let ((num (p-eval (cadr expr)))
659 (denom-inverse (apply #'divide (mapcar #'p-eval (cddr expr)))))
660 (multiply denom-inverse num)))))
661 (expt
662 (cond
663 ((member (cadr expr) vars :test #'equalp)
664 ;;Special handling of (expt var pow)
665 (let ((pos (position (cadr expr) vars :test #'equalp)))
666 (make-poly-variable (length vars) pos (caddr expr))))
667 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
668 ;; Negative power means division in coefficient ring
669 ;; Non-integer power means non-polynomial coefficient
670 expr)
671 (t (universal-expt (p-eval (cadr expr)) (caddr expr)))))
672 (otherwise
673 (error "Cannot evaluate as polynomial: ~A" expr)))))))
674
675(defmethod make-zero-for ((self poly))
676 (make-instance 'poly :dimension (poly-dimension self)))
677
678(defmethod make-unit-for ((self poly))
679 (make-poly-constant (poly-dimension self) 1))
680
681(defgeneric poly-reverse (self)
682 (:documentation "Reverse the order of terms in a polynomial SELF.")
683 (:method ((self poly))
684 (with-slots (termlist)
685 self
686 (setf termlist (nreverse termlist)))
687 self))
688
689
690
Note: See TracBrowser for help on using the repository browser.