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