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