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/division.lisp@ 68

Last change on this file since 68 was 64, checked in by Marek Rychlik, 9 years ago

* empty log message *

File size: 6.0 KB
Line 
1;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2;;
3;; An implementation of the division algorithm
4;;
5;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
6
7(defun grobner-op (ring c1 c2 m f g)
8 "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
9Assume that the leading terms will cancel."
10 #+grobner-check(funcall (ring-zerop ring)
11 (funcall (ring-sub ring)
12 (funcall (ring-mul ring) c2 (poly-lc f))
13 (funcall (ring-mul ring) c1 (poly-lc g))))
14 #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
15 ;; Note that we can drop the leading terms of f ang g
16 (poly-sub ring
17 (scalar-times-poly-1 ring c2 f)
18 (scalar-times-poly-1 ring c1 (monom-times-poly m g))))
19
20(defun poly-pseudo-divide (ring f fl)
21 "Pseudo-divide a polynomial F by the list of polynomials FL. Return
22multiple values. The first value is a list of quotients A. The second
23value is the remainder R. The third argument is a scalar coefficient
24C, such that C*F can be divided by FL within the ring of coefficients,
25which is not necessarily a field. Finally, the fourth value is an
26integer count of the number of reductions performed. The resulting
27objects satisfy the equation: C*F= sum A[i]*FL[i] + R."
28 (declare (type poly f) (list fl))
29 (do ((r (make-poly-zero))
30 (c (funcall (ring-unit ring)))
31 (a (make-list (length fl) :initial-element (make-poly-zero)))
32 (division-count 0)
33 (p f))
34 ((poly-zerop p)
35 (debug-cgb "~&~3T~d reduction~:p" division-count)
36 (when (poly-zerop r) (debug-cgb " ---> 0"))
37 (values (mapcar #'poly-nreverse a) (poly-nreverse r) c division-count))
38 (declare (fixnum division-count))
39 (do ((fl fl (rest fl)) ;scan list of divisors
40 (b a (rest b)))
41 ((cond
42 ((endp fl) ;no division occurred
43 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
44 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
45 (pop (poly-termlist p)) ;remove lt(p) from p
46 t)
47 ((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
48 (incf division-count)
49 (multiple-value-bind (gcd c1 c2)
50 (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
51 (declare (ignore gcd))
52 (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
53 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
54 (mapl #'(lambda (x)
55 (setf (car x) (scalar-times-poly ring c1 (car x))))
56 a)
57 (setf r (scalar-times-poly ring c1 r)
58 c (funcall (ring-mul ring) c c1)
59 p (grobner-op ring c2 c1 m p (car fl)))
60 (push (make-term m c2) (poly-termlist (car b))))
61 t)))))))
62
63(defun poly-exact-divide (ring f g)
64 "Divide a polynomial F by another polynomial G. Assume that exact division
65with no remainder is possible. Returns the quotient."
66 (declare (type poly f g))
67 (multiple-value-bind (quot rem coeff division-count)
68 (poly-pseudo-divide ring f (list g))
69 (declare (ignore division-count coeff)
70 (list quot)
71 (type poly rem)
72 (type fixnum division-count))
73 (unless (poly-zerop rem) (error "Exact division failed."))
74 (car quot)))
75
76
77
78;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
79;;
80;; An implementation of the normal form
81;;
82;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
83
84(defun normal-form-step (ring fl p r c division-count
85 &aux (g (find (poly-lm p) fl
86 :test #'monom-divisible-by-p
87 :key #'poly-lm)))
88 (cond
89 (g ;division possible
90 (incf division-count)
91 (multiple-value-bind (gcd cg cp)
92 (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
93 (declare (ignore gcd))
94 (let ((m (monom-div (poly-lm p) (poly-lm g))))
95 ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
96 (setf r (scalar-times-poly ring cg r)
97 c (funcall (ring-mul ring) c cg)
98 ;; p := cg*p-cp*m*g
99 p (grobner-op ring cp cg m p g))))
100 (debug-cgb "/"))
101 (t ;no division possible
102 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
103 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
104 (pop (poly-termlist p)) ;remove lt(p) from p
105 (debug-cgb "+")))
106 (values p r c division-count))
107
108;; Merge it sometime with poly-pseudo-divide
109(defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
110 ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
111 #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
112 (do ((r (make-poly-zero))
113 (c (funcall (ring-unit ring)))
114 (division-count 0))
115 ((or (poly-zerop f)
116 ;;(endp fl)
117 (and top-reduction-only (not (poly-zerop r))))
118 (progn
119 (debug-cgb "~&~3T~d reduction~:p" division-count)
120 (when (poly-zerop r)
121 (debug-cgb " ---> 0")))
122 (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
123 (values f c division-count))
124 (declare (fixnum division-count)
125 (type poly r))
126 (multiple-value-setq (f r c division-count)
127 (normal-form-step ring fl f r c division-count))))
128
129(defun buchberger-criterion (ring g)
130 "Returns T if G is a Grobner basis, by using the Buchberger
131criterion: for every two polynomials h1 and h2 in G the S-polynomial
132S(h1,h2) reduces to 0 modulo G."
133 (every
134 #'poly-zerop
135 (makelist (normal-form ring (spoly ring (elt g i) (elt g j)) g nil)
136 (i 0 (- (length g) 2))
137 (j (1+ i) (1- (length g))))))
138
139
140(defun poly-normalize (ring p &aux (c (poly-lc p)))
141 "Divide a polynomial by its leading coefficient. It assumes
142that the division is possible, which may not always be the
143case in rings which are not fields. The exact division operator
144is assumed to be provided by the RING structure of the
145COEFFICIENT-RING package."
146 (mapc #'(lambda (term)
147 (setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
148 (poly-termlist p))
149 p)
150
151(defun poly-normalize-list (ring plist)
152 "Divide every polynomial in a list PLIST by its leading coefficient. "
153 (mapcar #'(lambda (x) (poly-normalize ring x)) plist))
Note: See TracBrowser for help on using the repository browser.