Changeset 39

Timestamp:

2015-06-02T00:11:37-07:00 (10 years ago)

Author:

Marek Rychlik

Message:

* empty log message *

File:

• branches/f4grobner/f4-matrix.lisp (modified) (1 diff)

branches/f4grobner/f4-matrix.lisp

@@ -24 +24 @@
 (macsyma-module f4-maxima)
 
-(defun poly-monomials (p) 
-  (mapcar #'term-monom (poly-termlist p)))
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;
+;; An implementation of the normal form
+;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
-(defun make-f4-matrix (f fl &aux ml))
-  
-  
+(defun normal-form-step (ring fl p r c division-count
+                         &aux (g (find (poly-lm p) fl
+                                       :test #'monom-divisible-by-p
+                                       :key #'poly-lm)))
+  (cond
+   (g                                   ;division possible
+    (incf division-count)
+    (multiple-value-bind (gcd cg cp)
+        (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
+      (declare (ignore gcd))
+      (let ((m (monom-div (poly-lm p) (poly-lm g))))
+        ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
+        (setf r (scalar-times-poly ring cg r)
+              c (funcall (ring-mul ring) c cg)
+              ;; p := cg*p-cp*m*g
+              p (grobner-op ring cp cg m p g))))
+    (debug-cgb "/"))
+   (t                                                   ;no division possible
+    (push (poly-lt p) (poly-termlist r))                ;move lt(p) to remainder
+    (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
+    (pop (poly-termlist p))                             ;remove lt(p) from p
+    (debug-cgb "+")))
+  (values p r c division-count))
 
+;; Merge it sometime with poly-pseudo-divide
+(defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
+  ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
+  #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
+  (do ((r (make-poly-zero))
+       (c (funcall (ring-unit ring)))
+       (division-count 0))
+      ((or (poly-zerop f)
+           ;;(endp fl)
+           (and top-reduction-only (not (poly-zerop r))))
+       (progn
+         (debug-cgb "~&~3T~d reduction~:p" division-count)
+         (when (poly-zerop r)
+           (debug-cgb " ---> 0")))
+       (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
+       (values f c division-count))
+    (declare (fixnum division-count)
+             (type poly r))
+    (multiple-value-setq (f r c division-count)
+      (normal-form-step ring fl f r c division-count))))