[34] | 1 | ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
|
---|
| 2 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 3 | ;;;
|
---|
| 4 | ;;; Copyright (C) 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 |
|
---|
[35] | 22 | (in-package :maxima)
|
---|
| 23 |
|
---|
| 24 | (macsyma-module f4-maxima)
|
---|
| 25 |
|
---|
[160] | 26 | (defun f4-grobner-op (ring c1 c2 m f g)
|
---|
| 27 | "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
|
---|
| 28 | Assume that the leading terms will cancel."
|
---|
| 29 | #+grobner-check(funcall (ring-zerop ring)
|
---|
| 30 | (funcall (ring-sub ring)
|
---|
| 31 | (funcall (ring-mul ring) c2 (poly-lc f))
|
---|
| 32 | (funcall (ring-mul ring) c1 (poly-lc g))))
|
---|
| 33 | #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
|
---|
| 34 | (poly-sub ring
|
---|
| 35 | (scalar-times-poly ring c2 f)
|
---|
| 36 | (scalar-times-poly ring c1 (monom-times-poly m g))))
|
---|
| 37 |
|
---|
[39] | 38 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 39 | ;;
|
---|
| 40 | ;; An implementation of the normal form
|
---|
| 41 | ;;
|
---|
| 42 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
[35] | 43 |
|
---|
[40] | 44 | (defun f4-normal-form-step (ring fl p r c division-count
|
---|
[39] | 45 | &aux (g (find (poly-lm p) fl
|
---|
| 46 | :test #'monom-divisible-by-p
|
---|
| 47 | :key #'poly-lm)))
|
---|
| 48 | (cond
|
---|
| 49 | (g ;division possible
|
---|
| 50 | (incf division-count)
|
---|
| 51 | (multiple-value-bind (gcd cg cp)
|
---|
| 52 | (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
|
---|
| 53 | (declare (ignore gcd))
|
---|
| 54 | (let ((m (monom-div (poly-lm p) (poly-lm g))))
|
---|
| 55 | ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
|
---|
| 56 | (setf r (scalar-times-poly ring cg r)
|
---|
| 57 | c (funcall (ring-mul ring) c cg)
|
---|
| 58 | ;; p := cg*p-cp*m*g
|
---|
| 59 | p (grobner-op ring cp cg m p g))))
|
---|
| 60 | (debug-cgb "/"))
|
---|
| 61 | (t ;no division possible
|
---|
| 62 | (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
|
---|
| 63 | (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
|
---|
| 64 | (pop (poly-termlist p)) ;remove lt(p) from p
|
---|
| 65 | (debug-cgb "+")))
|
---|
| 66 | (values p r c division-count))
|
---|
[35] | 67 |
|
---|
[39] | 68 | ;; Merge it sometime with poly-pseudo-divide
|
---|
[40] | 69 | (defun f4-normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
|
---|
[39] | 70 | ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
|
---|
| 71 | #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
|
---|
| 72 | (do ((r (make-poly-zero))
|
---|
| 73 | (c (funcall (ring-unit ring)))
|
---|
| 74 | (division-count 0))
|
---|
| 75 | ((or (poly-zerop f)
|
---|
| 76 | ;;(endp fl)
|
---|
| 77 | (and top-reduction-only (not (poly-zerop r))))
|
---|
| 78 | (progn
|
---|
| 79 | (debug-cgb "~&~3T~d reduction~:p" division-count)
|
---|
| 80 | (when (poly-zerop r)
|
---|
| 81 | (debug-cgb " ---> 0")))
|
---|
| 82 | (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
|
---|
| 83 | (values f c division-count))
|
---|
| 84 | (declare (fixnum division-count)
|
---|
| 85 | (type poly r))
|
---|
| 86 | (multiple-value-setq (f r c division-count)
|
---|
[40] | 87 | (f4-normal-form-step ring fl f r c division-count))))
|
---|