;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; ;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik ;;; ;;; This program is free software; you can redistribute it and/or modify ;;; it under the terms of the GNU General Public License as published by ;;; the Free Software Foundation; either version 2 of the License, or ;;; (at your option) any later version. ;;; ;;; This program is distributed in the hope that it will be useful, ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;;; GNU General Public License for more details. ;;; ;;; You should have received a copy of the GNU General Public License ;;; along with this program; if not, write to the Free Software ;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. ;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defpackage "POLYNOMIAL" (:use :cl :ring :ring-and-order :monomial :term :termlist) (:export "POLY" "POLY-TERMLIST" "POLY-SUGAR" "POLY-LT" "MAKE-POLY-FROM-TERMLIST" "MAKE-POLY-ZERO" "MAKE-VARIABLE" "POLY-UNIT" "POLY-LM" "POLY-SECOND-LM" "POLY-SECOND-LT" "POLY-LC" "POLY-SECOND-LC" "POLY-ZEROP" "POLY-LENGTH" "SCALAR-TIMES-POLY" "SCALAR-TIMES-POLY-1" "MONOM-TIMES-POLY" "TERM-TIMES-POLY" "POLY-ADD" "POLY-SUB" "POLY-UMINUS" "POLY-MUL" "POLY-EXPT" "POLY-APPEND" "POLY-NREVERSE" "POLY-CONTRACT" "POLY-EXTEND" "POLY-ADD-VARIABLES" "POLY-LIST-ADD-VARIABLES" "POLY-STANDARD-EXTENSION" "SATURATION-EXTENSION" "POLYSATURATION-EXTENSION" "SATURATION-EXTENSION-1" "COERCE-COEFF" "POLY-EVAL" "SPOLY" "POLY-PRIMITIVE-PART" "POLY-CONTENT" )) (in-package :polynomial) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Polynomials ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defstruct (poly ;; ;; BOA constructor, by default constructs zero polynomial (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist)))) (:constructor make-poly-zero (&aux (termlist nil) (sugar -1))) ;; Constructor of polynomials representing a variable (:constructor make-variable (ring nvars pos &optional (power 1) &aux (termlist (list (make-term-variable ring nvars pos power))) (sugar power))) (:constructor poly-unit (ring dimension &aux (termlist (termlist-unit ring dimension)) (sugar 0)))) (termlist nil :type list) (sugar -1 :type fixnum)) ;; Leading term (defmacro poly-lt (p) `(car (poly-termlist ,p))) ;; Second term (defmacro poly-second-lt (p) `(cadar (poly-termlist ,p))) ;; Leading monomial (defun poly-lm (p) (term-monom (poly-lt p))) ;; Second monomial (defun poly-second-lm (p) (term-monom (poly-second-lt p))) ;; Leading coefficient (defun poly-lc (p) (term-coeff (poly-lt p))) ;; Second coefficient (defun poly-second-lc (p) (term-coeff (poly-second-lt p))) ;; Testing for a zero polynomial (defun poly-zerop (p) (null (poly-termlist p))) ;; The number of terms (defun poly-length (p) (length (poly-termlist p))) (defun scalar-times-poly (ring c p) (declare (type ring ring) (poly p)) (make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p))) ;; The scalar product omitting the head term (defun scalar-times-poly-1 (ring c p) (declare (type ring ring) (poly p)) (make-poly-from-termlist (scalar-times-termlist ring c (cdr (poly-termlist p))) (poly-sugar p))) (defun monom-times-poly (m p) (declare (poly p)) (make-poly-from-termlist (monom-times-termlist m (poly-termlist p)) (+ (poly-sugar p) (monom-sugar m)))) (defun term-times-poly (ring term p) (declare (type ring ring) (type term term) (type poly p)) (make-poly-from-termlist (term-times-termlist ring term (poly-termlist p)) (+ (poly-sugar p) (term-sugar term)))) (defun poly-add (ring-and-order p q) (declare (type ring-and-order ring-and-order) (type poly p q)) (make-poly-from-termlist (termlist-add ring-and-order (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q)))) (defun poly-sub (ring-and-order p q) (declare (type ring-and-order ring-and-order) (type poly p q)) (make-poly-from-termlist (termlist-sub ring-and-order (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q)))) (defun poly-uminus (ring p) (declare (type ring ring) (type poly p)) (make-poly-from-termlist (termlist-uminus ring (poly-termlist p)) (poly-sugar p))) (defun poly-mul (ring-and-order p q) (declare (type ring-and-order ring-and-order) (type poly p q)) (make-poly-from-termlist (termlist-mul ring-and-order (poly-termlist p) (poly-termlist q)) (+ (poly-sugar p) (poly-sugar q)))) (defun poly-expt (ring-and-order p n) (declare (type ring-and-order ring-and-order) (type poly p)) (make-poly-from-termlist (termlist-expt ring-and-order (poly-termlist p) n) (* n (poly-sugar p)))) (defun poly-append (&rest plist) (make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist)) (apply #'max (mapcar #'poly-sugar plist)))) (defun poly-nreverse (p) (declare (type poly p)) (setf (poly-termlist p) (nreverse (poly-termlist p))) p) (defun poly-contract (p &optional (k 1)) (declare (type poly p)) (make-poly-from-termlist (termlist-contract (poly-termlist p) k) (poly-sugar p))) (defun poly-extend (p &optional (m (make-monom :dimension 1))) (declare (type poly p)) (make-poly-from-termlist (termlist-extend (poly-termlist p) m) (+ (poly-sugar p) (monom-sugar m)))) (defun poly-add-variables (p k) (declare (type poly p)) (setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k)) p) (defun poly-list-add-variables (plist k) (mapcar #'(lambda (p) (poly-add-variables p k)) plist)) (defun poly-standard-extension (plist &aux (k (length plist))) "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]." (declare (list plist) (fixnum k)) (labels ((incf-power (g i) (dolist (x (poly-termlist g)) (incf (monom-elt (term-monom x) i))) (incf (poly-sugar g)))) (setf plist (poly-list-add-variables plist k)) (dotimes (i k plist) (incf-power (nth i plist) i)))) (defun saturation-extension (ring f plist &aux (k (length plist)) (d (monom-dimension (poly-lm (car plist))))) "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]." (setf f (poly-list-add-variables f k) plist (mapcar #'(lambda (x) (setf (poly-termlist x) (nconc (poly-termlist x) (list (make-term (make-monom :dimension d) (funcall (ring-uminus ring) (funcall (ring-unit ring))))))) x) (poly-standard-extension plist))) (append f plist)) (defun polysaturation-extension (ring f plist &aux (k (length plist)) (d (+ k (monom-dimension (poly-lm (car plist)))))) "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]." (setf f (poly-list-add-variables f k) plist (apply #'poly-append (poly-standard-extension plist)) (cdr (last (poly-termlist plist))) (list (make-term (make-monom :dimension d) (funcall (ring-uminus ring) (funcall (ring-unit ring)))))) (append f (list plist))) (defun saturation-extension-1 (ring f p) (polysaturation-extension ring f (list p))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Evaluation of polynomial (prefix) expressions ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defun coerce-coeff (ring expr vars) "Coerce an element of the coefficient ring to a constant polynomial." ;; Modular arithmetic handler by rat (make-poly-from-termlist (list (make-term (make-monom :dimension (length vars)) (funcall (ring-parse ring) expr))) 0)) (defun poly-eval (ring-and-order expr vars &optional (list-marker '[) &aux (ring (ro-ring ring-and-order))) (labels ((p-eval (arg) (poly-eval ring-and-order arg vars)) (p-eval-list (args) (mapcar #'p-eval args)) (p-add (x y) (poly-add ring-and-order x y))) (cond ((eql expr 0) (make-poly-zero)) ((member expr vars :test #'equalp) (let ((pos (position expr vars :test #'equalp))) (make-variable ring (length vars) pos))) ((atom expr) (coerce-coeff ring expr vars)) ((eq (car expr) list-marker) (cons list-marker (p-eval-list (cdr expr)))) (t (case (car expr) (+ (reduce #'p-add (p-eval-list (cdr expr)))) (- (case (length expr) (1 (make-poly-zero)) (2 (poly-uminus ring (p-eval (cadr expr)))) (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr)))) (otherwise (poly-sub ring-and-order (p-eval (cadr expr)) (reduce #'p-add (p-eval-list (cddr expr))))))) (* (if (endp (cddr expr)) ;unary (p-eval (cdr expr)) (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr))))) (expt (cond ((member (cadr expr) vars :test #'equalp) ;;Special handling of (expt var pow) (let ((pos (position (cadr expr) vars :test #'equalp))) (make-variable ring (length vars) pos (caddr expr)))) ((not (and (integerp (caddr expr)) (plusp (caddr expr)))) ;; Negative power means division in coefficient ring ;; Non-integer power means non-polynomial coefficient (coerce-coeff ring expr vars)) (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr))))) (otherwise (coerce-coeff ring expr vars))))))) (defun spoly (ring f g) "It yields the S-polynomial of polynomials F and G." (declare (type poly f g)) (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g))) (mf (monom-div lcm (poly-lm f))) (mg (monom-div lcm (poly-lm g)))) (declare (type monom mf mg)) (multiple-value-bind (c cf cg) (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g)) (declare (ignore c)) (poly-sub ring (scalar-times-poly ring cg (monom-times-poly mf f)) (scalar-times-poly ring cf (monom-times-poly mg g)))))) (defun poly-primitive-part (ring p) "Divide polynomial P with integer coefficients by gcd of its coefficients and return the result." (declare (type poly p)) (if (poly-zerop p) (values p 1) (let ((c (poly-content ring p))) (values (make-poly-from-termlist (mapcar #'(lambda (x) (make-term (term-monom x) (funcall (ring-div ring) (term-coeff x) c))) (poly-termlist p)) (poly-sugar p)) c)))) (defun poly-content (ring p) "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure to compute the greatest common divisor." (declare (type poly p)) (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p))) ;; Return the standard basis of the monomials in n variables (defun variable-basis (ring n &aux (basis (make-list n))) "Generate a list of polynomials X[i], i=0,1,...,N-1." (dotimes (i n basis) (setf (elt basis i) (make-variable ring n i)))) #| (defun poly-eval-1 (expr vars &optional (ring *ring-of-integers*) (order #'lex>) &aux (ring-and-order (make-ring-and-order :ring ring :order order)) (n (length vars)) (basis (variable-basis ring (length vars)))) "Evaluate an expression EXPR as polynomial by substituting operators + - * expt with corresponding polynomial operators and variables VARS with the corresponding polynomials in internal form. We use special versions of binary operators $poly+, $poly-, $minus-poly, $poly* and $poly-expt which work like the corresponding functions in the POLY package, but accept scalars as arguments as well. The result is a polynomial in internal form. This operation is somewhat similar to the function EXPAND in CAS." (cond ((numberp expr) (cond ((zerop expr) NIL) (t (make-poly-from-termlist (list (make-term (make-monom :dimension n) expr)))))) ((symbolp expr) (nth (position expr vars) basis)) ((consp expr) (case (car expr) (expt (if (= (length expr) 3) ($poly-expt ring-and-order (poly-eval-1 (cadr expr) vars ring order) (caddr expr) n) (error "Too many arguments to EXPT"))) (/ (if (and (= (length expr) 3) (numberp (caddr expr))) ($poly/ ring (cadr expr) (caddr expr)) (error "The second argument to / must be a number"))) (otherwise (let ((r (mapcar #'(lambda (e) (poly-eval-1 e vars ring order)) (cdr expr)))) (ecase (car expr) (+ (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) r)) (- (if (endp (cdr r)) ($minus-poly ring (car r) n) ($poly- ring-and-order (car r) (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) (cdr r)) n))) (* (reduce #'(lambda (p q) ($poly* ring-and-order p q n)) r)) ))))))) (defun poly-eval (expr vars &optional (order #'lex>) (ring *ring-of-integers*)) "Evaluate an expression EXPR, which should be a polynomial expression or a list of polynomial expressions (a list of expressions marked by prepending keyword :[ to it) given in Lisp prefix notation, in variables VARS, which should be a list of symbols. The result of the evaluation is a polynomial or a list of polynomials (marked by prepending symbol '[) in the internal alist form. This evaluator is used by the PARSE package to convert input from strings directly to internal form." (cond ((numberp expr) (unless (zerop expr) (make-poly-from-termlist (list (make-term (make-monom :dimension (length vars)) expr))))) ((or (symbolp expr) (not (eq (car expr) :[))) (poly-eval-1 expr vars ring order)) (t (cons '[ (mapcar #'(lambda (p) (poly-eval-1 p vars ring order)) (rest expr)))))) |#