close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

source: branches/f4grobner/f4-matrix.lisp@ 459

Last change on this file since 459 was 160, checked in by Marek Rychlik, 9 years ago
File size: 3.7 KB
Line 
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
22(in-package :maxima)
23
24(macsyma-module f4-maxima)
25
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.
28Assume 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
38;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
39;;
40;; An implementation of the normal form
41;;
42;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
43
44(defun f4-normal-form-step (ring fl p r c division-count
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))
67
68;; Merge it sometime with poly-pseudo-divide
69(defun f4-normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
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)
87 (f4-normal-form-step ring fl f r c division-count))))
Note: See TracBrowser for help on using the repository browser.