Changeset 4463 for branches/f4grobner

Timestamp:

2016-06-14T19:35:34-07:00 (9 years ago)

Author:

Marek Rychlik

Message:

 

Location:

branches/f4grobner

Files:

• 5am-buchberger.lisp (modified) (1 diff)
• 5am-symbolic-poly.lisp (modified) (1 diff)
• buchberger.lisp (modified) (3 diffs)
• division.lisp (modified) (9 diffs)
• fast-add.lisp (modified) (1 diff)
• monom.lisp (modified) (1 diff)
• polynomial-eval.lisp (modified) (1 diff)
• polynomial.lisp (modified) (1 diff)
• symbolic-polynomial.lisp (modified) (3 diffs)
• test0.lisp (modified) (1 diff)
• test8.lisp (modified) (1 diff)

branches/f4grobner/5am-buchberger.lisp

@@ -107 +107 @@
   (with-fixture buchberger-advanced-context ()
     (is-true (grobner-test gb fl))
-    (is (every #'universal-equalp (buchberger fl) gb))
+    (is (grobner-equal (buchberger fl) gb))
     ;;(is (every #'universal-equalp (parallel-buchberger fl) gb))
     )

branches/f4grobner/5am-symbolic-poly.lisp

@@ -70 +70 @@
     (is (universal-equalp (change-class p 'symbolic-poly :vars '(x)) p-symbolic )))
   (with-fixture sym-poly-context ()
-    (signals
-        (error "Number of variables does not equal dimension.")
-      (universal-equalp (change-class p 'symbolic-poly :vars '(x y)) p-symbolic ))))
+    (universal-equalp (change-class p 'symbolic-poly :vars '(x y)) p-symbolic )))
 
 ;; The following example uses the Grobner basis obtained with Maxima

branches/f4grobner/buchberger.lisp

@@ -26 +26 @@
         :ring
         )
-  (:export "BUCHBERGER" "PARALLEL-BUCHBERGER")
+  (:export "BUCHBERGER" "PARALLEL-BUCHBERGER" "GROBNER-MEMBER" "GROBNER-SUBSETP" "GROBNER-EQUAL")
   (:documentation "Buchberger Algorithm Implementation."))
 
@@ -49 +49 @@
                     (grobner-test (subseq f 0 start) (subseq f 0 start)))
   ;;Initialize critical pairs
-  (let ((b (make-critical-pair-queue *normal-strategy* f start))
+  (let ((b (make-critical-pair-queue *min-total-degree-strategy* f start))
         (b-done (make-hash-table :test #'equal)))
     (declare (type critical-pair-queue b) (type hash-table b-done))
@@ -85 +85 @@
         (setf (gethash (list (critical-pair-first pair) (critical-pair-second pair)) b-done)
               t)))))
+
+
+(defun grobner-member (p g)
+  "Returns T if a polynomial P belongs to the ideal generated by the
+polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
+  (universal-zerop (normal-form p g nil)))
+
+
+(defun grobner-subsetp (g1 g2)
+  "Returns T if a list of polynomials G1 generates
+an ideal contained in the ideal generated by a polynomial list G2,
+both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
+  (every #'(lambda (p) (grobner-member p g2)) g1))
+
+
+(defun grobner-equal (g1 g2)
+  "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
+generate  the same ideal, and NIL otherwise."
+  (and (grobner-subsetp g1 g2) (grobner-subsetp g2 g1)))

branches/f4grobner/division.lisp

@@ -43 +43 @@
 
 (defun grobner-op (c1 c2 m f g)
-  "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
-Assume that the leading terms will cancel."
-  (declare (type monom m)
-           (type poly f g))
-  (assert (universal-zerop
-           (subtract
-            (multiply c2 (leading-coefficient f))
-            (multiply c1 (leading-coefficient g)))))
-  (assert (universal-equalp (leading-monomial f) (multiply m (leading-monomial g))))
-  ;; Note that below we can drop the leading terms of f ang g for the
-  ;; purpose of polynomial arithmetic.  
-  ;;
-  ;; TODO: Make sure that the sugar calculation is correct if leading
-  ;; terms are dropped.
-  (subtract
-   (multiply f c2)
-   (multiply g m c1)))
+  "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial."
+  (declare (type monom m) (type poly f g))
+  (subtract (multiply f c2) (multiply g m c1)))
 
 (defun check-loop-invariant (c f a fl r p &aux (p-zero (make-zero-for f)))
@@ -74 +60 @@
   |#
   (let* ((prod (inner-product a fl add multiply p-zero))
-         (succeeded-p (universal-zerop (subtract (multiply f c) (add prod r p)))))
+         (succeeded-p (universal-zerop (subtract (multiply f c) (add prod (make-instance 'poly :termlist (reverse r)) p)))))
     (unless succeeded-p 
       (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%" 
@@ -92 +78 @@
   (declare (type poly f) (list fl))
   ;; Loop invariant: c*f=sum ai*fi+r+p, where p must eventually become 0
-  (do ((r (make-zero-for f))
+  (do ((r nil)
        (c (make-unit-for (leading-coefficient f)))
        (a (make-list (length fl) :initial-element (make-zero-for f)))
@@ -100 +86 @@
        #+grobner-check(check-loop-invariant c f a fl r p)
        (debug-cgb "~&~3T~d reduction~:p" division-count)
-       (when (universal-zerop r) (debug-cgb " ---> 0"))
-       (values a r c division-count))
+       (when (null r) (debug-cgb " ---> 0"))
+       (values a (make-instance 'poly :termlist r) c division-count))
     (declare (fixnum division-count))
     ;; Check the loop invariant here
@@ -108 +94 @@
          (b a (rest b)))
         ((cond
-           ((endp fl)                           ;no division occurred
-            (setf r (add-to r (leading-term p)) ;move lt(p) to remainder
-                  p (subtract-from p (leading-term p))) ;remove lt(p) from p
+           ((endp fl)                     ;no division occurred
+            (push (poly-remove-term p) r) ;move lt(p) to remainder
             t)
            ((divides-p (leading-monomial (car fl)) (leading-monomial p)) ;division occurred
@@ -122 +107 @@
                           (setf (car x) (multiply-by (car x) c1)))
                       a)
-                (setf r (multiply-by r c1)
+                (setf r (mapc #'multiply-by r c1)
                       c (multiply-by c c1)
                       p (grobner-op c2 c1 m p (car fl)))
@@ -154 +139 @@
                                     :test #'divisible-by-p
                                     :key #'leading-monomial)))
+  ;; NOTE: Currently R is a list of terms of the remainder
   (cond
    (g                                   ;division possible
@@ -162 +148 @@
       (let ((m (divide (leading-monomial p) (leading-monomial g))))
         ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
-        (setf r (multiply-by r cg)
+        (setf r (mapc #'(lambda (trm) (multiply-by trm cg)) r)
               c (multiply-by c cg)
               ;; p := cg*p-cp*m*g
               p (grobner-op cp cg m p g))))
     (debug-cgb "/"))
-   (t                                                   ;no division possible
-    (setf r (add-to r (leading-term p))) ;move lt(p) to remainder
-    (setf p (subtract-from p (leading-term p))) ;move lt(p) to remainder
+   (t                                      ;no division possible
+    (setf r (push (poly-remove-term p) r))         ;move lt(p) to remainder
     (debug-cgb "+")))
   (values p r c division-count))
@@ -187 +172 @@
     ;; unit in the coefficient field.
     (return-from normal-form (values f nil 0)))
-  (do ((r (make-zero-for f))
+  (do ((r nil)
        (c (make-unit-for (leading-coefficient f)))
        (division-count 0))
       ((or (universal-zerop f)
            ;;(endp fl)
-           (and top-reduction-only (not (universal-zerop r))))
+           (and top-reduction-only (not (null r))))
        (progn
          (debug-cgb "~&~3T~D reduction~:P" division-count)
-         (when (universal-zerop r)
+         (when (null r)
            (debug-cgb " ---> 0")))
-       (setf f (add-to f r))
+       (setf (poly-termlist f) (nreconc r (poly-termlist f)))
        (values f c division-count))
-    (declare (fixnum division-count)
-             (type poly r))
+    (declare (fixnum division-count))
     (multiple-value-setq (f r c division-count)
       (normal-form-step fl f r c division-count))))

branches/f4grobner/fast-add.lisp

@@ -12 +12 @@
                      result-last (cdr result-last)))))
     (loop 
-       (format t "RESULT:      ~S~%" (mapcar #'->list result))
-       (format t "RESULT-LAST: ~S~%" (mapcar #'->list result-last))
-       (format t "P:           ~S~%" (mapcar #'->list p))
-       (format t "Q:           ~S~%" (mapcar #'->list q))
+       ;;TODO: Eliminate this debugging code
+       #+nil
+       (progn
+         (format t "~V@{~C~:*~}~%" 80 #\-)
+         (format t "RESULT:      ~S~%" (mapcar #'->list result))
+         (format t "RESULT-LAST: ~S~%" (mapcar #'->list result-last))
+         (format t "P:           ~S~%" (mapcar #'->list p))
+         (format t "Q:           ~S~%" (mapcar #'->list q)))
        (cond
          ((endp p) (unless (endp q) (add-to-result q)) (return result))

branches/f4grobner/monom.lisp

@@ -568 +568 @@
       copy))
 
-#|
-(defmethod multiply-by ((self term) (other number))
+(defmethod multiply-by ((self term) (other ring))
   (reinitialize-instance self :coeff (multiply-by (term-coeff self) other)))
 
-(defmethod divide-by ((self term) (other number))
+(defmethod divide-by ((self term) (other ring))
   (reinitialize-instance self :coeff (divide-by (term-coeff self) other)))
-|#
+
 
 (defmethod unary-inverse :after ((self term))

branches/f4grobner/polynomial-eval.lisp

@@ -38 +38 @@
     (cond 
       ((eq expr 0) 
-       (make-instance 'poly :dimension (length vars)))
+       (make-instance 'poly))
       ((member expr vars :test #'equalp)
        (let ((pos (position expr vars :test #'equalp)))

branches/f4grobner/polynomial.lisp

@@ -738 +738 @@
   "Ensures that the number of variables in VARS maches the polynomial dimension of the
 polynomial SELF."
-  (let ((dimension (poly-dimension self)))
-    (assert (= (length vars) dimension)
-            nil
-            "Number of variables ~S does not match the dimension ~S"
-            vars dimension)))
+  (unless (endp (poly-termlist self))
+    (let ((dimension (poly-dimension self)))
+      (assert (= (length vars) dimension)
+              nil
+              "Number of variables ~S does not match the dimension ~S"
+              vars dimension))))
 
 (defmethod ->sexp ((self poly) &optional vars)

branches/f4grobner/symbolic-polynomial.lisp

@@ -38 +38 @@
 (defmethod print-object ((self symbolic-poly) stream)
   (print-unreadable-object (self stream :type t :identity t)
-    (with-accessors ((dimension poly-dimension)
-                     (termlist poly-termlist) 
+    (with-accessors ((termlist poly-termlist) 
                      (order poly-term-order)
                      (vars symbolic-poly-vars))
         self
-      (format stream "DIMENSION=~A TERMLIST=~A ORDER=~A VARS=~A" 
-              dimension termlist order vars))))
+      (format stream "TERMLIST=~A ORDER=~A VARS=~A" 
+              termlist order vars))))
 
 
@@ -59 +58 @@
 (defmethod universal-equalp ((self symbol) (other symbol))
   (eq self other))
-
-(defmethod update-instance-for-different-class :after ((old poly) (new  symbolic-poly) &key)
-  "After adding variables to NEW, we need to make sure that the number
-of variables given by POLY-DIMENSION is consistent with VARS."
-  (assert (= (length (symbolic-poly-vars new)) (poly-dimension new))))
-
 
 #|
@@ -168 +161 @@
            (format nil "[~{~a~^, ~}]" (mapcar #'(lambda (p) (poly->string p)) (cdr self))))))
   (:method ((self poly) &optional (vars nil))
-    ;; Ensure that the number of variables matches the dimension
-    (assert (= (length vars) (poly-dimension self)))
     (infix-print-to-string (->sexp self vars)))
   (:method ((self symbolic-poly) &optional (vars (symbolic-poly-vars self)))

branches/f4grobner/test0.lisp

@@ -4 +4 @@
 
 (setf p (ADD-TO
-         (make-instance 'poly :dimension 1 :termlist (list (make-instance 'term :exponents #(0) :coeff 5)) :order #'lex>)
-         (make-instance 'poly :dimension 1 :termlist (list (make-instance 'term :exponents #(1) :coeff 1)) :order #'lex>)))
+         (make-instance 'poly :termlist (list (make-instance 'term :exponents #(0) :coeff 5)) :order #'lex>)
+         (make-instance 'poly :termlist (list (make-instance 'term :exponents #(1) :coeff 1)) :order #'lex>)))
 
 (setf r (POLYNOMIAL::ADD-TERMLISTS (list (make-instance 'term :exponents #(0) :coeff 5))

branches/f4grobner/test8.lisp

@@ -1 +1 @@
 (in-package :polynomial)
+
 (proclaim '(special p0 q0 p q p+q-good p+q-risky))
 
-(setf p0 (alist->poly (nreverse '(((2 0 0) . 2) ((1 1 0) . 1) ((1 0 0) . 1) ((0 1 0) . -1) ((0 0 0) . 2)))))
-(setf q0 (alist->poly (nreverse '(((2 0 0) . 2) ((0 1 0) . -1) ((0 0 0) . 2)))))
+(setf p0 (alist->poly (reverse '(((2 0 0) . 2) ((1 1 0) . 1) ((1 0 0) . 1) ((0 1 0) . -1) ((0 0 0) . 2)))))
+(setf q0 (alist->poly (reverse '(((2 0 0) . 2) ((0 1 0) . -1) ((0 0 0) . 2)))))
 
 (setf p (poly-termlist p0))
 (setf q (poly-termlist q0))
 
-(setf p+q-good (mapcar #'->list (slow-add p q #'lex> #'add-to)))
+(setf p+q-good (mapcar #'->list (slow-add (mapcar #'copy:copy-instance p) (mapcar #'copy:copy-instance q) #'lex> #'add-to)))
 
-(setf p+q-risky (mapcar #'->list (fast-and-risky-add p q #'lex> #'add-to)))
+(setf p+q-risky (mapcar #'->list (fast-and-risky-add (mapcar #'copy:copy-instance p) (mapcar #'copy:copy-instance q) #'lex> #'add-to)))
 
-;;(assert (equal p+q-good p+q-risky) nil "Two methods gave different results.")
+(assert (equal p+q-good p+q-risky) nil "Two methods gave different results.")