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