| 1 | ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
|
---|
| 2 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 3 | ;;;
|
---|
| 4 | ;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
|
---|
| 5 | ;;;
|
---|
| 6 | ;;; This program is free software; you can redistribute it and/or modify
|
---|
| 7 | ;;; it under the terms of the GNU General Public License as published by
|
---|
| 8 | ;;; the Free Software Foundation; either version 2 of the License, or
|
---|
| 9 | ;;; (at your option) any later version.
|
---|
| 10 | ;;;
|
---|
| 11 | ;;; This program is distributed in the hope that it will be useful,
|
---|
| 12 | ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 13 | ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 14 | ;;; GNU General Public License for more details.
|
---|
| 15 | ;;;
|
---|
| 16 | ;;; You should have received a copy of the GNU General Public License
|
---|
| 17 | ;;; along with this program; if not, write to the Free Software
|
---|
| 18 | ;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
---|
| 19 | ;;;
|
---|
| 20 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 21 |
|
---|
| 22 | (in-package :ngrobner)
|
---|
| 23 |
|
---|
| 24 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 25 | ;;
|
---|
| 26 | ;; Low-level polynomial arithmetic done on
|
---|
| 27 | ;; lists of terms
|
---|
| 28 | ;;
|
---|
| 29 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 30 |
|
---|
| 31 | (defmacro termlist-lt (p) `(car ,p))
|
---|
| 32 | (defun termlist-lm (p) (term-monom (termlist-lt p)))
|
---|
| 33 | (defun termlist-lc (p) (term-coeff (termlist-lt p)))
|
---|
| 34 |
|
---|
| 35 | (define-modify-macro scalar-mul (c) coeff-mul)
|
---|
| 36 |
|
---|
| 37 | (defun scalar-times-termlist (ring c p)
|
---|
| 38 | "Multiply scalar C by a polynomial P. This function works
|
---|
| 39 | even if there are divisors of 0."
|
---|
| 40 | (mapcan
|
---|
| 41 | #'(lambda (term)
|
---|
| 42 | (let ((c1 (funcall (ring-mul ring) c (term-coeff term))))
|
---|
| 43 | (unless (funcall (ring-zerop ring) c1)
|
---|
| 44 | (list (make-term (term-monom term) c1)))))
|
---|
| 45 | p))
|
---|
| 46 |
|
---|
| 47 |
|
---|
| 48 | (defun term-mul-lst (ring term1 term2)
|
---|
| 49 | "Returns (LIST TERM) where TERM is the product of the terms TERM1 TERM2,
|
---|
| 50 | or NIL when the product is 0. This definition takes care of divisors of 0
|
---|
| 51 | in the coefficient ring."
|
---|
| 52 | (let ((c (funcall (ring-mul ring) (term-coeff term1) (term-coeff term2))))
|
---|
| 53 | (unless (funcall (ring-zerop ring) c)
|
---|
| 54 | (list (make-term (monom-mul (term-monom term1) (term-monom term2)) c)))))
|
---|
| 55 |
|
---|
| 56 | (defun term-times-termlist (ring term f)
|
---|
| 57 | (declare (type ring ring))
|
---|
| 58 | (mapcan #'(lambda (term-f) (term-mul-lst ring term term-f)) f))
|
---|
| 59 |
|
---|
| 60 | (defun termlist-times-term (ring f term)
|
---|
| 61 | (mapcan #'(lambda (term-f) (term-mul-lst ring term-f term)) f))
|
---|
| 62 |
|
---|
| 63 | (defun monom-times-term (m term)
|
---|
| 64 | (make-term (monom-mul m (term-monom term)) (term-coeff term)))
|
---|
| 65 |
|
---|
| 66 | (defun monom-times-termlist (m f)
|
---|
| 67 | (cond
|
---|
| 68 | ((null f) nil)
|
---|
| 69 | (t
|
---|
| 70 | (mapcar #'(lambda (x) (monom-times-term m x)) f))))
|
---|
| 71 |
|
---|
| 72 | (defun termlist-uminus (ring f)
|
---|
| 73 | (mapcar #'(lambda (x)
|
---|
| 74 | (make-term (term-monom x) (funcall (ring-uminus ring) (term-coeff x))))
|
---|
| 75 | f))
|
---|
| 76 |
|
---|
| 77 | (defun termlist-add (ring p q)
|
---|
| 78 | (declare (type list p q))
|
---|
| 79 | (do (r)
|
---|
| 80 | ((cond
|
---|
| 81 | ((endp p)
|
---|
| 82 | (setf r (revappend r q)) t)
|
---|
| 83 | ((endp q)
|
---|
| 84 | (setf r (revappend r p)) t)
|
---|
| 85 | (t
|
---|
| 86 | (multiple-value-bind
|
---|
| 87 | (lm-greater lm-equal)
|
---|
| 88 | (monomial-order (termlist-lm p) (termlist-lm q))
|
---|
| 89 | (cond
|
---|
| 90 | (lm-equal
|
---|
| 91 | (let ((s (funcall (ring-add ring) (termlist-lc p) (termlist-lc q))))
|
---|
| 92 | (unless (funcall (ring-zerop ring) s) ;check for cancellation
|
---|
| 93 | (setf r (cons (make-term (termlist-lm p) s) r)))
|
---|
| 94 | (setf p (cdr p) q (cdr q))))
|
---|
| 95 | (lm-greater
|
---|
| 96 | (setf r (cons (car p) r)
|
---|
| 97 | p (cdr p)))
|
---|
| 98 | (t (setf r (cons (car q) r)
|
---|
| 99 | q (cdr q)))))
|
---|
| 100 | nil))
|
---|
| 101 | r)))
|
---|
| 102 |
|
---|
| 103 | (defun termlist-sub (ring p q)
|
---|
| 104 | (declare (type list p q))
|
---|
| 105 | (do (r)
|
---|
| 106 | ((cond
|
---|
| 107 | ((endp p)
|
---|
| 108 | (setf r (revappend r (termlist-uminus ring q)))
|
---|
| 109 | t)
|
---|
| 110 | ((endp q)
|
---|
| 111 | (setf r (revappend r p))
|
---|
| 112 | t)
|
---|
| 113 | (t
|
---|
| 114 | (multiple-value-bind
|
---|
| 115 | (mgreater mequal)
|
---|
| 116 | (monomial-order (termlist-lm p) (termlist-lm q))
|
---|
| 117 | (cond
|
---|
| 118 | (mequal
|
---|
| 119 | (let ((s (funcall (ring-sub ring) (termlist-lc p) (termlist-lc q))))
|
---|
| 120 | (unless (funcall (ring-zerop ring) s) ;check for cancellation
|
---|
| 121 | (setf r (cons (make-term (termlist-lm p) s) r)))
|
---|
| 122 | (setf p (cdr p) q (cdr q))))
|
---|
| 123 | (mgreater
|
---|
| 124 | (setf r (cons (car p) r)
|
---|
| 125 | p (cdr p)))
|
---|
| 126 | (t (setf r (cons (make-term (termlist-lm q) (funcall (ring-uminus ring) (termlist-lc q))) r)
|
---|
| 127 | q (cdr q)))))
|
---|
| 128 | nil))
|
---|
| 129 | r)))
|
---|
| 130 |
|
---|
| 131 | ;; Multiplication of polynomials
|
---|
| 132 | ;; Non-destructive version
|
---|
| 133 | (defun termlist-mul (ring p q)
|
---|
| 134 | (cond ((or (endp p) (endp q)) nil) ;p or q is 0 (represented by NIL)
|
---|
| 135 | ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
|
---|
| 136 | ((endp (cdr p))
|
---|
| 137 | (term-times-termlist ring (car p) q))
|
---|
| 138 | ((endp (cdr q))
|
---|
| 139 | (termlist-times-term ring p (car q)))
|
---|
| 140 | (t
|
---|
| 141 | (let ((head (term-mul-lst ring (termlist-lt p) (termlist-lt q)))
|
---|
| 142 | (tail (termlist-add ring (term-times-termlist ring (car p) (cdr q))
|
---|
| 143 | (termlist-mul ring (cdr p) q))))
|
---|
| 144 | (cond ((null head) tail)
|
---|
| 145 | ((null tail) head)
|
---|
| 146 | (t (nconc head tail)))))))
|
---|
| 147 |
|
---|
| 148 | (defun termlist-unit (ring dimension)
|
---|
| 149 | (declare (fixnum dimension))
|
---|
| 150 | (list (make-term (make-monom dimension :initial-element 0)
|
---|
| 151 | (funcall (ring-unit ring)))))
|
---|
| 152 |
|
---|
| 153 | (defun termlist-expt (ring poly n &aux (dim (monom-dimension (termlist-lm poly))))
|
---|
| 154 | (declare (type fixnum n dim))
|
---|
| 155 | (cond
|
---|
| 156 | ((minusp n) (error "termlist-expt: Negative exponent."))
|
---|
| 157 | ((endp poly) (if (zerop n) (termlist-unit ring dim) nil))
|
---|
| 158 | (t
|
---|
| 159 | (do ((k 1 (ash k 1))
|
---|
| 160 | (q poly (termlist-mul ring q q)) ;keep squaring
|
---|
| 161 | (p (termlist-unit ring dim) (if (not (zerop (logand k n))) (termlist-mul ring p q) p)))
|
---|
| 162 | ((> k n) p)
|
---|
| 163 | (declare (fixnum k))))))
|
---|