[3400] | 1 | ;;----------------------------------------------------------------
|
---|
| 2 | ;; File: polynomial.lisp
|
---|
| 3 | ;;----------------------------------------------------------------
|
---|
| 4 | ;;
|
---|
| 5 | ;; Author: Marek Rychlik (rychlik@u.arizona.edu)
|
---|
| 6 | ;; Date: Thu Aug 27 09:41:24 2015
|
---|
| 7 | ;; Copying: (C) Marek Rychlik, 2010. All rights reserved.
|
---|
| 8 | ;;
|
---|
| 9 | ;;----------------------------------------------------------------
|
---|
[1201] | 10 | ;;; -*- Mode: Lisp -*-
|
---|
[77] | 11 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 12 | ;;;
|
---|
| 13 | ;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
|
---|
| 14 | ;;;
|
---|
| 15 | ;;; This program is free software; you can redistribute it and/or modify
|
---|
| 16 | ;;; it under the terms of the GNU General Public License as published by
|
---|
| 17 | ;;; the Free Software Foundation; either version 2 of the License, or
|
---|
| 18 | ;;; (at your option) any later version.
|
---|
| 19 | ;;;
|
---|
| 20 | ;;; This program is distributed in the hope that it will be useful,
|
---|
| 21 | ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 22 | ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 23 | ;;; GNU General Public License for more details.
|
---|
| 24 | ;;;
|
---|
| 25 | ;;; You should have received a copy of the GNU General Public License
|
---|
| 26 | ;;; along with this program; if not, write to the Free Software
|
---|
| 27 | ;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
---|
| 28 | ;;;
|
---|
| 29 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 30 |
|
---|
[431] | 31 | (defpackage "POLYNOMIAL"
|
---|
[3643] | 32 | (:use :cl :utils :monom :copy)
|
---|
[2596] | 33 | (:export "POLY"
|
---|
[3270] | 34 | "POLY-DIMENSION"
|
---|
[2596] | 35 | "POLY-TERMLIST"
|
---|
[3016] | 36 | "POLY-TERM-ORDER"
|
---|
[3509] | 37 | "POLY-INSERT-TERM"
|
---|
[3690] | 38 | "SCALAR-MULTIPLY-BY"
|
---|
| 39 | "SCALAR-DIVIDE-BY"
|
---|
[3642] | 40 | "LEADING-TERM"
|
---|
[3657] | 41 | "LEADING-MONOMIAL"
|
---|
[3642] | 42 | "LEADING-COEFFICIENT"
|
---|
[3657] | 43 | "SECOND-LEADING-TERM"
|
---|
| 44 | "SECOND-LEADING-MONOMIAL"
|
---|
| 45 | "SECOND-LEADING-COEFFICIENT"
|
---|
[3642] | 46 | "ADD-TO"
|
---|
[3646] | 47 | "ADD"
|
---|
[3642] | 48 | "SUBTRACT-FROM"
|
---|
[3646] | 49 | "SUBTRACT"
|
---|
[3071] | 50 | "CHANGE-TERM-ORDER"
|
---|
[3099] | 51 | "STANDARD-EXTENSION"
|
---|
[3101] | 52 | "STANDARD-EXTENSION-1"
|
---|
[3109] | 53 | "STANDARD-SUM"
|
---|
[3094] | 54 | "SATURATION-EXTENSION"
|
---|
[3655] | 55 | "ALIST->POLY"
|
---|
| 56 | "UNIVERSAL-EZGCD"
|
---|
[3678] | 57 | "S-POLYNOMIAL"
|
---|
[3686] | 58 | "POLY-CONTENT"
|
---|
| 59 | "POLY-PRIMITIVE-PART")
|
---|
[3489] | 60 | (:documentation "Implements polynomials. A polynomial is essentially
|
---|
| 61 | a mapping of monomials of the same degree to coefficients. The
|
---|
| 62 | momomials are ordered according to a monomial order."))
|
---|
[143] | 63 |
|
---|
[431] | 64 | (in-package :polynomial)
|
---|
| 65 |
|
---|
[1927] | 66 | (proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
|
---|
[52] | 67 |
|
---|
[2442] | 68 | (defclass poly ()
|
---|
[3253] | 69 | ((dimension :initform nil
|
---|
[3250] | 70 | :initarg :dimension
|
---|
| 71 | :accessor poly-dimension
|
---|
[3242] | 72 | :documentation "Shared dimension of all terms, the number of variables")
|
---|
[3250] | 73 | (termlist :initform nil :initarg :termlist :accessor poly-termlist
|
---|
[3619] | 74 | :documentation "List of terms.")
|
---|
[3250] | 75 | (order :initform #'lex> :initarg :order :accessor poly-term-order
|
---|
[2697] | 76 | :documentation "Monomial/term order."))
|
---|
[3262] | 77 | (:default-initargs :dimension nil :termlist nil :order #'lex>)
|
---|
[2695] | 78 | (:documentation "A polynomial with a list of terms TERMLIST, ordered
|
---|
[2696] | 79 | according to term order ORDER, which defaults to LEX>."))
|
---|
[2442] | 80 |
|
---|
[2471] | 81 | (defmethod print-object ((self poly) stream)
|
---|
[3241] | 82 | (print-unreadable-object (self stream :type t :identity t)
|
---|
[3243] | 83 | (with-accessors ((dimension poly-dimension)
|
---|
| 84 | (termlist poly-termlist)
|
---|
| 85 | (order poly-term-order))
|
---|
[3237] | 86 | self
|
---|
[3244] | 87 | (format stream "DIMENSION=~A TERMLIST=~A ORDER=~A"
|
---|
| 88 | dimension termlist order))))
|
---|
[2469] | 89 |
|
---|
[3015] | 90 | (defgeneric change-term-order (self other)
|
---|
[3012] | 91 | (:documentation "Change term order of SELF to the term order of OTHER.")
|
---|
[3010] | 92 | (:method ((self poly) (other poly))
|
---|
| 93 | (unless (eq (poly-term-order self) (poly-term-order other))
|
---|
[3620] | 94 | (setf (poly-termlist self) (sort (poly-termlist self) (poly-term-order other))
|
---|
[3010] | 95 | (poly-term-order self) (poly-term-order other)))
|
---|
[3012] | 96 | self))
|
---|
[3010] | 97 |
|
---|
[3621] | 98 | (defgeneric poly-insert-term (self term)
|
---|
[3622] | 99 | (:documentation "Insert a term TERM into SELF before all other
|
---|
[3621] | 100 | terms. Order is not enforced.")
|
---|
| 101 | (:method ((self poly) (term term))
|
---|
[3510] | 102 | (cond ((null (poly-dimension self))
|
---|
[3621] | 103 | (setf (poly-dimension self) (monom-dimension term)))
|
---|
| 104 | (t (assert (= (poly-dimension self) (monom-dimension term)))))
|
---|
| 105 | (push term (poly-termlist self))
|
---|
[3510] | 106 | self))
|
---|
| 107 |
|
---|
[3622] | 108 | (defgeneric poly-append-term (self term)
|
---|
| 109 | (:documentation "Append a term TERM to SELF after all other terms. Order is not enforced.")
|
---|
| 110 | (:method ((self poly) (term term))
|
---|
[3510] | 111 | (cond ((null (poly-dimension self))
|
---|
[3622] | 112 | (setf (poly-dimension self) (monom-dimension term)))
|
---|
| 113 | (t (assert (= (poly-dimension self) (monom-dimension term)))))
|
---|
| 114 | (setf (cdr (last (poly-termlist self))) (list term))
|
---|
[3510] | 115 | self))
|
---|
| 116 |
|
---|
[3095] | 117 | (defun alist->poly (alist &aux (poly (make-instance 'poly)))
|
---|
[3126] | 118 | "It reads polynomial from an alist formatted as ( ... (exponents . coeff) ...).
|
---|
| 119 | It can be used to enter simple polynomials by hand, e.g the polynomial
|
---|
| 120 | in two variables, X and Y, given in standard notation as:
|
---|
| 121 |
|
---|
| 122 | 3*X^2*Y^3+2*Y+7
|
---|
| 123 |
|
---|
| 124 | can be entered as
|
---|
| 125 | (ALIST->POLY '(((2 3) . 3) ((0 1) . 2) ((0 0) . 7))).
|
---|
| 126 |
|
---|
| 127 | NOTE: The primary use is for low-level debugging of the package."
|
---|
[3099] | 128 | (dolist (x alist poly)
|
---|
[3623] | 129 | (poly-insert-term poly (make-instance 'term :exponents (car x) :coeff (cdr x)))))
|
---|
[3092] | 130 |
|
---|
[3401] | 131 | (defmethod update-instance-for-different-class :after ((old monom) (new poly) &key)
|
---|
| 132 | "Converts OLD of class MONOM to a NEW of class POLY, by making it into a 1-element TERMLIST."
|
---|
| 133 | (reinitialize-instance new
|
---|
| 134 | :dimension (monom-dimension old)
|
---|
[3511] | 135 | :termlist (list (cons old 1))))
|
---|
[3403] | 136 |
|
---|
[3624] | 137 | (defmethod universal-equalp ((self poly) (other poly))
|
---|
| 138 | "Implements equality of polynomials."
|
---|
| 139 | (and (eql (poly-dimension self) (poly-dimension other))
|
---|
| 140 | (every #'universal-equalp (poly-termlist self) (poly-termlist other))
|
---|
| 141 | (eq (poly-term-order self) (poly-term-order other))))
|
---|
[2650] | 142 |
|
---|
[3624] | 143 | (defgeneric leading-term (object)
|
---|
[2442] | 144 | (:method ((self poly))
|
---|
[2525] | 145 | (car (poly-termlist self)))
|
---|
| 146 | (:documentation "The leading term of a polynomial, or NIL for zero polynomial."))
|
---|
[52] | 147 |
|
---|
[3625] | 148 | (defgeneric second-leading-term (object)
|
---|
[2442] | 149 | (:method ((self poly))
|
---|
[2525] | 150 | (cadar (poly-termlist self)))
|
---|
| 151 | (:documentation "The second leading term of a polynomial, or NIL for a polynomial with at most one term."))
|
---|
[52] | 152 |
|
---|
[3656] | 153 | (defgeneric leading-monomial (object)
|
---|
| 154 | (:method ((self poly))
|
---|
| 155 | (change-class (copy-instance (leading-term self)) 'monom))
|
---|
| 156 | (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
|
---|
| 157 |
|
---|
| 158 | (defgeneric second-leading-monomial (object)
|
---|
| 159 | (:method ((self poly))
|
---|
| 160 | (change-class (copy-instance (second-leading-term self)) 'monom))
|
---|
| 161 | (:documentation "The leading monomial of a polynomial, or NIL for zero polynomial."))
|
---|
| 162 |
|
---|
[3625] | 163 | (defgeneric leading-coefficient (object)
|
---|
[2442] | 164 | (:method ((self poly))
|
---|
[3642] | 165 | (term-coeff (leading-term self)))
|
---|
[2545] | 166 | (:documentation "The leading coefficient of a polynomial. It signals error for a zero polynomial."))
|
---|
[52] | 167 |
|
---|
[2442] | 168 | (defgeneric second-leading-coefficient (object)
|
---|
| 169 | (:method ((self poly))
|
---|
[3645] | 170 | (term-coeff (second-leading-term self)))
|
---|
[2906] | 171 | (:documentation "The second leading coefficient of a polynomial. It
|
---|
| 172 | signals error for a polynomial with at most one term."))
|
---|
[52] | 173 |
|
---|
[3629] | 174 | (defmethod universal-zerop ((self poly))
|
---|
| 175 | "Return T iff SELF is a zero polynomial."
|
---|
[3639] | 176 | (null (poly-termlist self)))
|
---|
[52] | 177 |
|
---|
[3518] | 178 | (defgeneric poly-length (self)
|
---|
[3630] | 179 | (:documentation "Return the number of terms.")
|
---|
[3518] | 180 | (:method ((self poly))
|
---|
| 181 | (length (poly-termlist self))))
|
---|
[52] | 182 |
|
---|
[3689] | 183 | (defgeneric scalar-multiply-by (self other)
|
---|
| 184 | (:documentation "Multiply vector SELF by a scalar OTHER.")
|
---|
| 185 | (:method ((self poly) other)
|
---|
| 186 | (mapc #'(lambda (term) (setf (term-coeff term) (multiply (term-coeff term) other)))
|
---|
| 187 | (poly-termlist self))
|
---|
| 188 | self))
|
---|
| 189 |
|
---|
| 190 | (defgeneric scalar-divide-by (self other)
|
---|
| 191 | (:documentation "Divide vector SELF by a scalar OTHER.")
|
---|
| 192 | (:method ((self poly) other)
|
---|
| 193 | (mapc #'(lambda (term) (setf (term-coeff term) (divide (term-coeff term) other)))
|
---|
| 194 | (poly-termlist self))
|
---|
| 195 | self))
|
---|
| 196 |
|
---|
[3663] | 197 | (defmethod multiply-by ((self poly) (other monom))
|
---|
[3630] | 198 | "Multiply a polynomial SELF by OTHER."
|
---|
| 199 | (mapc #'(lambda (term) (multiply-by term other))
|
---|
| 200 | (poly-termlist self))
|
---|
| 201 | self)
|
---|
[2469] | 202 |
|
---|
[3672] | 203 | (defmethod multiply-by ((self poly) (other term))
|
---|
| 204 | "Multiply a polynomial SELF by OTHER."
|
---|
| 205 | (mapc #'(lambda (term) (multiply-by term other))
|
---|
| 206 | (poly-termlist self))
|
---|
| 207 | self)
|
---|
| 208 |
|
---|
[2761] | 209 | (defmacro fast-add/subtract (p q order-fn add/subtract-fn uminus-fn)
|
---|
[2755] | 210 | "Return an expression which will efficiently adds/subtracts two
|
---|
| 211 | polynomials, P and Q. The addition/subtraction of coefficients is
|
---|
| 212 | performed by calling ADD/SUBTRACT-METHOD-NAME. If UMINUS-METHOD-NAME
|
---|
| 213 | is supplied, it is used to negate the coefficients of Q which do not
|
---|
[2756] | 214 | have a corresponding coefficient in P. The code implements an
|
---|
| 215 | efficient algorithm to add two polynomials represented as sorted lists
|
---|
| 216 | of terms. The code destroys both arguments, reusing the terms to build
|
---|
| 217 | the result."
|
---|
[3631] | 218 | `(macrolet ((lc (x) `(term-coeff (car ,x))))
|
---|
[2742] | 219 | (do ((p ,p)
|
---|
| 220 | (q ,q)
|
---|
| 221 | r)
|
---|
| 222 | ((or (endp p) (endp q))
|
---|
| 223 | ;; NOTE: R contains the result in reverse order. Can it
|
---|
| 224 | ;; be more efficient to produce the terms in correct order?
|
---|
[2774] | 225 | (unless (endp q)
|
---|
[2776] | 226 | ;; Upon subtraction, we must change the sign of
|
---|
| 227 | ;; all coefficients in q
|
---|
[2774] | 228 | ,@(when uminus-fn
|
---|
[2775] | 229 | `((mapc #'(lambda (x) (setf x (funcall ,uminus-fn x))) q)))
|
---|
[2774] | 230 | (setf r (nreconc r q)))
|
---|
[2742] | 231 | r)
|
---|
| 232 | (multiple-value-bind
|
---|
| 233 | (greater-p equal-p)
|
---|
[3632] | 234 | (funcall ,order-fn (car p) (car q))
|
---|
[2742] | 235 | (cond
|
---|
| 236 | (greater-p
|
---|
| 237 | (rotatef (cdr p) r p)
|
---|
| 238 | )
|
---|
| 239 | (equal-p
|
---|
[2766] | 240 | (let ((s (funcall ,add/subtract-fn (lc p) (lc q))))
|
---|
[2742] | 241 | (cond
|
---|
[3640] | 242 | ((universal-zerop s)
|
---|
[2742] | 243 | (setf p (cdr p))
|
---|
| 244 | )
|
---|
| 245 | (t
|
---|
| 246 | (setf (lc p) s)
|
---|
| 247 | (rotatef (cdr p) r p))))
|
---|
| 248 | (setf q (cdr q))
|
---|
| 249 | )
|
---|
| 250 | (t
|
---|
[2743] | 251 | ;;Negate the term of Q if UMINUS provided, signallig
|
---|
| 252 | ;;that we are doing subtraction
|
---|
[2908] | 253 | ,(when uminus-fn
|
---|
| 254 | `(setf (lc q) (funcall ,uminus-fn (lc q))))
|
---|
[2743] | 255 | (rotatef (cdr q) r q)))))))
|
---|
[2585] | 256 |
|
---|
[2655] | 257 |
|
---|
[3647] | 258 | (defgeneric add-to (self other)
|
---|
| 259 | (:documentation "Add OTHER to SELF.")
|
---|
| 260 | (:method ((self number) (other number))
|
---|
| 261 | (+ self other)))
|
---|
| 262 |
|
---|
| 263 | (defgeneric subtract-from (self other)
|
---|
[3648] | 264 | (:documentation "Subtract OTHER from SELF.")
|
---|
| 265 | (:method ((self number) (other number))
|
---|
| 266 | (- self other)))
|
---|
[3647] | 267 |
|
---|
[2763] | 268 | (defmacro def-add/subtract-method (add/subtract-method-name
|
---|
[2752] | 269 | uminus-method-name
|
---|
| 270 | &optional
|
---|
[2913] | 271 | (doc-string nil doc-string-supplied-p))
|
---|
[3647] | 272 | "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
|
---|
[2749] | 273 | `(defmethod ,add/subtract-method-name ((self poly) (other poly))
|
---|
[2615] | 274 | ,@(when doc-string-supplied-p `(,doc-string))
|
---|
[2769] | 275 | ;; Ensure orders are compatible
|
---|
[3015] | 276 | (change-term-order other self)
|
---|
[2772] | 277 | (setf (poly-termlist self) (fast-add/subtract
|
---|
| 278 | (poly-termlist self) (poly-termlist other)
|
---|
| 279 | (poly-term-order self)
|
---|
| 280 | #',add/subtract-method-name
|
---|
| 281 | ,(when uminus-method-name `(function ,uminus-method-name))))
|
---|
[2609] | 282 | self))
|
---|
[2487] | 283 |
|
---|
[2916] | 284 | (eval-when (:compile-toplevel :load-toplevel :execute)
|
---|
[2777] | 285 |
|
---|
[3639] | 286 | (def-add/subtract-method add-to nil
|
---|
[2777] | 287 | "Adds to polynomial SELF another polynomial OTHER.
|
---|
[2610] | 288 | This operation destructively modifies both polynomials.
|
---|
| 289 | The result is stored in SELF. This implementation does
|
---|
[2752] | 290 | no consing, entirely reusing the sells of SELF and OTHER.")
|
---|
[2609] | 291 |
|
---|
[3639] | 292 | (def-add/subtract-method subtract-from unary-minus
|
---|
[2753] | 293 | "Subtracts from polynomial SELF another polynomial OTHER.
|
---|
[2610] | 294 | This operation destructively modifies both polynomials.
|
---|
| 295 | The result is stored in SELF. This implementation does
|
---|
[2752] | 296 | no consing, entirely reusing the sells of SELF and OTHER.")
|
---|
[2916] | 297 | )
|
---|
[2777] | 298 |
|
---|
[2691] | 299 | (defmethod unary-minus ((self poly))
|
---|
[2694] | 300 | "Destructively modifies the coefficients of the polynomial SELF,
|
---|
| 301 | by changing their sign."
|
---|
[2692] | 302 | (mapc #'unary-minus (poly-termlist self))
|
---|
[2683] | 303 | self)
|
---|
[52] | 304 |
|
---|
[2795] | 305 | (defun add-termlists (p q order-fn)
|
---|
[2794] | 306 | "Destructively adds two termlists P and Q ordered according to ORDER-FN."
|
---|
[3639] | 307 | (fast-add/subtract p q order-fn #'add-to nil))
|
---|
[2794] | 308 |
|
---|
[2800] | 309 | (defmacro multiply-term-by-termlist-dropping-zeros (term termlist
|
---|
[2927] | 310 | &optional (reverse-arg-order-P nil))
|
---|
[2799] | 311 | "Multiplies term TERM by a list of term, TERMLIST.
|
---|
[2792] | 312 | Takes into accound divisors of zero in the ring, by
|
---|
[2927] | 313 | deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
|
---|
[2928] | 314 | is T, change the order of arguments; this may be important
|
---|
[2927] | 315 | if we extend the package to non-commutative rings."
|
---|
[2800] | 316 | `(mapcan #'(lambda (other-term)
|
---|
[3633] | 317 | (let ((prod (multiply
|
---|
[2923] | 318 | ,@(cond
|
---|
[2930] | 319 | (reverse-arg-order-p
|
---|
[2925] | 320 | `(other-term ,term))
|
---|
| 321 | (t
|
---|
| 322 | `(,term other-term))))))
|
---|
[2800] | 323 | (cond
|
---|
[3633] | 324 | ((universal-zerop prod) nil)
|
---|
[2800] | 325 | (t (list prod)))))
|
---|
| 326 | ,termlist))
|
---|
[2790] | 327 |
|
---|
[2796] | 328 | (defun multiply-termlists (p q order-fn)
|
---|
[3127] | 329 | "A version of polynomial multiplication, operating
|
---|
| 330 | directly on termlists."
|
---|
[2787] | 331 | (cond
|
---|
[2917] | 332 | ((or (endp p) (endp q))
|
---|
| 333 | ;;p or q is 0 (represented by NIL)
|
---|
| 334 | nil)
|
---|
[2789] | 335 | ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
|
---|
[2787] | 336 | ((endp (cdr p))
|
---|
[2918] | 337 | (multiply-term-by-termlist-dropping-zeros (car p) q))
|
---|
| 338 | ((endp (cdr q))
|
---|
[2919] | 339 | (multiply-term-by-termlist-dropping-zeros (car q) p t))
|
---|
| 340 | (t
|
---|
[3633] | 341 | (cons (multiply (car p) (car q))
|
---|
[2949] | 342 | (add-termlists
|
---|
| 343 | (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
|
---|
| 344 | (multiply-termlists (cdr p) q order-fn)
|
---|
| 345 | order-fn)))))
|
---|
[2793] | 346 |
|
---|
[2803] | 347 | (defmethod multiply-by ((self poly) (other poly))
|
---|
[3014] | 348 | (change-term-order other self)
|
---|
[2803] | 349 | (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
|
---|
| 350 | (poly-termlist other)
|
---|
| 351 | (poly-term-order self)))
|
---|
| 352 | self)
|
---|
| 353 |
|
---|
[3634] | 354 | (defun add (object1 object2)
|
---|
[3374] | 355 | "Non-destructively add POLY1 by POLY2."
|
---|
[3634] | 356 | (add-to (copy-instance object1) (copy-instance object2)))
|
---|
[3374] | 357 |
|
---|
[3634] | 358 | (defun subtract (minuend &rest subtrahends)
|
---|
[3427] | 359 | "Non-destructively subtract MINUEND and SUBTRAHENDS."
|
---|
[3634] | 360 | (subtract-from (copy-instance minuend) (reduce #'add subtrahends)))
|
---|
[3374] | 361 |
|
---|
[3062] | 362 | (defmethod left-tensor-product-by ((self poly) (other monom))
|
---|
| 363 | (setf (poly-termlist self)
|
---|
| 364 | (mapcan #'(lambda (term)
|
---|
| 365 | (let ((prod (left-tensor-product-by term other)))
|
---|
| 366 | (cond
|
---|
[3640] | 367 | ((universal-zerop prod) nil)
|
---|
[3062] | 368 | (t (list prod)))))
|
---|
| 369 | (poly-termlist self)))
|
---|
[3249] | 370 | (incf (poly-dimension self) (monom-dimension other))
|
---|
[3062] | 371 | self)
|
---|
[3044] | 372 |
|
---|
[3062] | 373 | (defmethod right-tensor-product-by ((self poly) (other monom))
|
---|
| 374 | (setf (poly-termlist self)
|
---|
| 375 | (mapcan #'(lambda (term)
|
---|
| 376 | (let ((prod (right-tensor-product-by term other)))
|
---|
| 377 | (cond
|
---|
[3640] | 378 | ((universal-zerop prod) nil)
|
---|
[3062] | 379 | (t (list prod)))))
|
---|
| 380 | (poly-termlist self)))
|
---|
[3249] | 381 | (incf (poly-dimension self) (monom-dimension other))
|
---|
[3062] | 382 | self)
|
---|
| 383 |
|
---|
| 384 |
|
---|
[3084] | 385 | (defun standard-extension (plist &aux (k (length plist)) (i 0))
|
---|
[2716] | 386 | "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
|
---|
[3060] | 387 | is a list of polynomials. Destructively modifies PLIST elements."
|
---|
[3061] | 388 | (mapc #'(lambda (poly)
|
---|
[3085] | 389 | (left-tensor-product-by
|
---|
| 390 | poly
|
---|
| 391 | (prog1
|
---|
| 392 | (make-monom-variable k i)
|
---|
| 393 | (incf i))))
|
---|
[3061] | 394 | plist))
|
---|
[52] | 395 |
|
---|
[3087] | 396 | (defun standard-extension-1 (plist
|
---|
| 397 | &aux
|
---|
[3096] | 398 | (plist (standard-extension plist))
|
---|
[3087] | 399 | (nvars (poly-dimension (car plist))))
|
---|
[3081] | 400 | "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
|
---|
[3087] | 401 | Firstly, new K variables U1, U2, ..., UK, are inserted into each
|
---|
| 402 | polynomial. Subsequently, P1, P2, ..., PK are destructively modified
|
---|
[3105] | 403 | tantamount to replacing PI with UI*PI-1. It assumes that all
|
---|
[3106] | 404 | polynomials have the same dimension, and only the first polynomial
|
---|
| 405 | is examined to determine this dimension."
|
---|
[3089] | 406 | ;; Implementation note: we use STANDARD-EXTENSION and then subtract
|
---|
| 407 | ;; 1 from each polynomial; since UI*PI has no constant term,
|
---|
| 408 | ;; we just need to append the constant term at the end
|
---|
| 409 | ;; of each termlist.
|
---|
[3064] | 410 | (flet ((subtract-1 (p)
|
---|
[3641] | 411 | (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
|
---|
[3083] | 412 | (setf plist (mapc #'subtract-1 plist)))
|
---|
[3077] | 413 | plist)
|
---|
[52] | 414 |
|
---|
| 415 |
|
---|
[3107] | 416 | (defun standard-sum (plist
|
---|
| 417 | &aux
|
---|
| 418 | (plist (standard-extension plist))
|
---|
| 419 | (nvars (poly-dimension (car plist))))
|
---|
[3087] | 420 | "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
|
---|
| 421 | Firstly, new K variables, U1, U2, ..., UK, are inserted into each
|
---|
| 422 | polynomial. Subsequently, P1, P2, ..., PK are destructively modified
|
---|
| 423 | tantamount to replacing PI with UI*PI, and the resulting polynomials
|
---|
[3117] | 424 | are added. Finally, 1 is subtracted. It should be noted that the term
|
---|
| 425 | order is not modified, which is equivalent to using a lexicographic
|
---|
| 426 | order on the first K variables."
|
---|
[3107] | 427 | (flet ((subtract-1 (p)
|
---|
[3641] | 428 | (poly-append-term p (make-instance 'term :dimension nvars :coeff -1))))
|
---|
[3108] | 429 | (subtract-1
|
---|
| 430 | (make-instance
|
---|
| 431 | 'poly
|
---|
[3115] | 432 | :termlist (apply #'nconc (mapcar #'poly-termlist plist))))))
|
---|
[52] | 433 |
|
---|
[3653] | 434 | (defgeneric universal-ezgcd (x y)
|
---|
| 435 | (:documentation "Solves the diophantine system: X=C*X1, Y=C*X2,
|
---|
| 436 | C=GCD(X,Y). It returns C, X1 and Y1. The result may be obtained by
|
---|
| 437 | the Euclidean algorithm.")
|
---|
| 438 | (:method ((x integer) (y integer)
|
---|
| 439 | &aux (c (gcd x y)))
|
---|
| 440 | (values c (/ x c) (/ y c)))
|
---|
| 441 | )
|
---|
| 442 |
|
---|
[3655] | 443 | (defgeneric s-polynomial (object1 object2)
|
---|
[3651] | 444 | (:documentation "Yields the S-polynomial of OBJECT1 and OBJECT2.")
|
---|
| 445 | (:method ((f poly) (g poly))
|
---|
| 446 | (let* ((lcm (universal-lcm (leading-monomial f) (leading-monomial g)))
|
---|
| 447 | (mf (divide lcm (leading-monomial f)))
|
---|
| 448 | (mg (divide lcm (leading-monomial g))))
|
---|
| 449 | (multiple-value-bind (c cf cg)
|
---|
[3652] | 450 | (universal-ezgcd (leading-coefficient f) (leading-coefficient g))
|
---|
[3651] | 451 | (declare (ignore c))
|
---|
| 452 | (subtract
|
---|
[3673] | 453 | (multiply f (change-class mf 'term :coeff cg))
|
---|
| 454 | (multiply g (change-class mg 'term :coeff cf)))))))
|
---|
[3651] | 455 |
|
---|
[3676] | 456 | (defgeneric poly-content (object)
|
---|
| 457 | (:documentation "Greatest common divisor of the coefficients of the polynomial object OBJECT.")
|
---|
[3677] | 458 | (:method ((self poly))
|
---|
| 459 | (reduce #'universal-gcd
|
---|
[3679] | 460 | (mapcar #'term-coeff (rest (poly-termlist self)))
|
---|
| 461 | :initial-value (leading-coefficient self))))
|
---|
[3676] | 462 |
|
---|
[3684] | 463 | (defun poly-primitive-part (object)
|
---|
[3685] | 464 | "Divide polynomial OBJECT by gcd of its
|
---|
[3684] | 465 | coefficients. Return the resulting polynomial."
|
---|
[3688] | 466 | (scalar-divide-by object (poly-content object)))
|
---|
[3682] | 467 |
|
---|
[3122] | 468 | #|
|
---|
| 469 |
|
---|
[1477] | 470 | (defun saturation-extension-1 (ring f p)
|
---|
[1497] | 471 | "Calculate [F, U*P-1]. It destructively modifies F."
|
---|
[1908] | 472 | (declare (type ring ring))
|
---|
[1477] | 473 | (polysaturation-extension ring f (list p)))
|
---|
[53] | 474 |
|
---|
[55] | 475 |
|
---|
[1066] | 476 |
|
---|
[2456] | 477 | |#
|
---|