- Timestamp:
- 2016-06-14T19:35:34-07:00 (9 years ago)
- Location:
- branches/f4grobner
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/f4grobner/5am-buchberger.lisp
r4440 r4463 107 107 (with-fixture buchberger-advanced-context () 108 108 (is-true (grobner-test gb fl)) 109 (is ( every #'universal-equalp(buchberger fl) gb))109 (is (grobner-equal (buchberger fl) gb)) 110 110 ;;(is (every #'universal-equalp (parallel-buchberger fl) gb)) 111 111 ) -
branches/f4grobner/5am-symbolic-poly.lisp
r4325 r4463 70 70 (is (universal-equalp (change-class p 'symbolic-poly :vars '(x)) p-symbolic ))) 71 71 (with-fixture sym-poly-context () 72 (signals 73 (error "Number of variables does not equal dimension.") 74 (universal-equalp (change-class p 'symbolic-poly :vars '(x y)) p-symbolic )))) 72 (universal-equalp (change-class p 'symbolic-poly :vars '(x y)) p-symbolic ))) 75 73 76 74 ;; The following example uses the Grobner basis obtained with Maxima -
branches/f4grobner/buchberger.lisp
r4462 r4463 26 26 :ring 27 27 ) 28 (:export "BUCHBERGER" "PARALLEL-BUCHBERGER" )28 (:export "BUCHBERGER" "PARALLEL-BUCHBERGER" "GROBNER-MEMBER" "GROBNER-SUBSETP" "GROBNER-EQUAL") 29 29 (:documentation "Buchberger Algorithm Implementation.")) 30 30 … … 49 49 (grobner-test (subseq f 0 start) (subseq f 0 start))) 50 50 ;;Initialize critical pairs 51 (let ((b (make-critical-pair-queue * normal-strategy* f start))51 (let ((b (make-critical-pair-queue *min-total-degree-strategy* f start)) 52 52 (b-done (make-hash-table :test #'equal))) 53 53 (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) 86 86 t))))) 87 88 89 (defun grobner-member (p g) 90 "Returns T if a polynomial P belongs to the ideal generated by the 91 polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise." 92 (universal-zerop (normal-form p g nil))) 93 94 95 (defun grobner-subsetp (g1 g2) 96 "Returns T if a list of polynomials G1 generates 97 an ideal contained in the ideal generated by a polynomial list G2, 98 both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise." 99 (every #'(lambda (p) (grobner-member p g2)) g1)) 100 101 102 (defun grobner-equal (g1 g2) 103 "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases, 104 generate the same ideal, and NIL otherwise." 105 (and (grobner-subsetp g1 g2) (grobner-subsetp g2 g1))) -
branches/f4grobner/division.lisp
r4434 r4463 43 43 44 44 (defun grobner-op (c1 c2 m f g) 45 "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial. 46 Assume that the leading terms will cancel." 47 (declare (type monom m) 48 (type poly f g)) 49 (assert (universal-zerop 50 (subtract 51 (multiply c2 (leading-coefficient f)) 52 (multiply c1 (leading-coefficient g))))) 53 (assert (universal-equalp (leading-monomial f) (multiply m (leading-monomial g)))) 54 ;; Note that below we can drop the leading terms of f ang g for the 55 ;; purpose of polynomial arithmetic. 56 ;; 57 ;; TODO: Make sure that the sugar calculation is correct if leading 58 ;; terms are dropped. 59 (subtract 60 (multiply f c2) 61 (multiply g m c1))) 45 "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial." 46 (declare (type monom m) (type poly f g)) 47 (subtract (multiply f c2) (multiply g m c1))) 62 48 63 49 (defun check-loop-invariant (c f a fl r p &aux (p-zero (make-zero-for f))) … … 74 60 |# 75 61 (let* ((prod (inner-product a fl add multiply p-zero)) 76 (succeeded-p (universal-zerop (subtract (multiply f c) (add prod rp)))))62 (succeeded-p (universal-zerop (subtract (multiply f c) (add prod (make-instance 'poly :termlist (reverse r)) p))))) 77 63 (unless succeeded-p 78 64 (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)) 93 79 ;; Loop invariant: c*f=sum ai*fi+r+p, where p must eventually become 0 94 (do ((r (make-zero-for f))80 (do ((r nil) 95 81 (c (make-unit-for (leading-coefficient f))) 96 82 (a (make-list (length fl) :initial-element (make-zero-for f))) … … 100 86 #+grobner-check(check-loop-invariant c f a fl r p) 101 87 (debug-cgb "~&~3T~d reduction~:p" division-count) 102 (when ( universal-zeropr) (debug-cgb " ---> 0"))103 (values a rc division-count))88 (when (null r) (debug-cgb " ---> 0")) 89 (values a (make-instance 'poly :termlist r) c division-count)) 104 90 (declare (fixnum division-count)) 105 91 ;; Check the loop invariant here … … 108 94 (b a (rest b))) 109 95 ((cond 110 ((endp fl) ;no division occurred 111 (setf r (add-to r (leading-term p)) ;move lt(p) to remainder 112 p (subtract-from p (leading-term p))) ;remove lt(p) from p 96 ((endp fl) ;no division occurred 97 (push (poly-remove-term p) r) ;move lt(p) to remainder 113 98 t) 114 99 ((divides-p (leading-monomial (car fl)) (leading-monomial p)) ;division occurred … … 122 107 (setf (car x) (multiply-by (car x) c1))) 123 108 a) 124 (setf r (m ultiply-by r c1)109 (setf r (mapc #'multiply-by r c1) 125 110 c (multiply-by c c1) 126 111 p (grobner-op c2 c1 m p (car fl))) … … 154 139 :test #'divisible-by-p 155 140 :key #'leading-monomial))) 141 ;; NOTE: Currently R is a list of terms of the remainder 156 142 (cond 157 143 (g ;division possible … … 162 148 (let ((m (divide (leading-monomial p) (leading-monomial g)))) 163 149 ;; Multiply the equation c*f=sum ai*fi+r+p by cg. 164 (setf r (m ultiply-by r cg)150 (setf r (mapc #'(lambda (trm) (multiply-by trm cg)) r) 165 151 c (multiply-by c cg) 166 152 ;; p := cg*p-cp*m*g 167 153 p (grobner-op cp cg m p g)))) 168 154 (debug-cgb "/")) 169 (t ;no division possible 170 (setf r (add-to r (leading-term p))) ;move lt(p) to remainder 171 (setf p (subtract-from p (leading-term p))) ;move lt(p) to remainder 155 (t ;no division possible 156 (setf r (push (poly-remove-term p) r)) ;move lt(p) to remainder 172 157 (debug-cgb "+"))) 173 158 (values p r c division-count)) … … 187 172 ;; unit in the coefficient field. 188 173 (return-from normal-form (values f nil 0))) 189 (do ((r (make-zero-for f))174 (do ((r nil) 190 175 (c (make-unit-for (leading-coefficient f))) 191 176 (division-count 0)) 192 177 ((or (universal-zerop f) 193 178 ;;(endp fl) 194 (and top-reduction-only (not ( universal-zeropr))))179 (and top-reduction-only (not (null r)))) 195 180 (progn 196 181 (debug-cgb "~&~3T~D reduction~:P" division-count) 197 (when ( universal-zeropr)182 (when (null r) 198 183 (debug-cgb " ---> 0"))) 199 (setf f (add-to f r))184 (setf (poly-termlist f) (nreconc r (poly-termlist f))) 200 185 (values f c division-count)) 201 (declare (fixnum division-count) 202 (type poly r)) 186 (declare (fixnum division-count)) 203 187 (multiple-value-setq (f r c division-count) 204 188 (normal-form-step fl f r c division-count)))) -
branches/f4grobner/fast-add.lisp
r4450 r4463 12 12 result-last (cdr result-last))))) 13 13 (loop 14 (format t "RESULT: ~S~%" (mapcar #'->list result)) 15 (format t "RESULT-LAST: ~S~%" (mapcar #'->list result-last)) 16 (format t "P: ~S~%" (mapcar #'->list p)) 17 (format t "Q: ~S~%" (mapcar #'->list q)) 14 ;;TODO: Eliminate this debugging code 15 #+nil 16 (progn 17 (format t "~V@{~C~:*~}~%" 80 #\-) 18 (format t "RESULT: ~S~%" (mapcar #'->list result)) 19 (format t "RESULT-LAST: ~S~%" (mapcar #'->list result-last)) 20 (format t "P: ~S~%" (mapcar #'->list p)) 21 (format t "Q: ~S~%" (mapcar #'->list q))) 18 22 (cond 19 23 ((endp p) (unless (endp q) (add-to-result q)) (return result)) -
branches/f4grobner/monom.lisp
r4446 r4463 568 568 copy)) 569 569 570 #| 571 (defmethod multiply-by ((self term) (other number)) 570 (defmethod multiply-by ((self term) (other ring)) 572 571 (reinitialize-instance self :coeff (multiply-by (term-coeff self) other))) 573 572 574 (defmethod divide-by ((self term) (other number))573 (defmethod divide-by ((self term) (other ring)) 575 574 (reinitialize-instance self :coeff (divide-by (term-coeff self) other))) 576 |# 575 577 576 578 577 (defmethod unary-inverse :after ((self term)) -
branches/f4grobner/polynomial-eval.lisp
r4386 r4463 38 38 (cond 39 39 ((eq expr 0) 40 (make-instance 'poly :dimension (length vars)))40 (make-instance 'poly)) 41 41 ((member expr vars :test #'equalp) 42 42 (let ((pos (position expr vars :test #'equalp))) -
branches/f4grobner/polynomial.lisp
r4457 r4463 738 738 "Ensures that the number of variables in VARS maches the polynomial dimension of the 739 739 polynomial SELF." 740 (let ((dimension (poly-dimension self))) 741 (assert (= (length vars) dimension) 742 nil 743 "Number of variables ~S does not match the dimension ~S" 744 vars dimension))) 740 (unless (endp (poly-termlist self)) 741 (let ((dimension (poly-dimension self))) 742 (assert (= (length vars) dimension) 743 nil 744 "Number of variables ~S does not match the dimension ~S" 745 vars dimension)))) 745 746 746 747 (defmethod ->sexp ((self poly) &optional vars) -
branches/f4grobner/symbolic-polynomial.lisp
r4442 r4463 38 38 (defmethod print-object ((self symbolic-poly) stream) 39 39 (print-unreadable-object (self stream :type t :identity t) 40 (with-accessors ((dimension poly-dimension) 41 (termlist poly-termlist) 40 (with-accessors ((termlist poly-termlist) 42 41 (order poly-term-order) 43 42 (vars symbolic-poly-vars)) 44 43 self 45 (format stream " DIMENSION=~ATERMLIST=~A ORDER=~A VARS=~A"46 dimensiontermlist order vars))))44 (format stream "TERMLIST=~A ORDER=~A VARS=~A" 45 termlist order vars)))) 47 46 48 47 … … 59 58 (defmethod universal-equalp ((self symbol) (other symbol)) 60 59 (eq self other)) 61 62 (defmethod update-instance-for-different-class :after ((old poly) (new symbolic-poly) &key)63 "After adding variables to NEW, we need to make sure that the number64 of variables given by POLY-DIMENSION is consistent with VARS."65 (assert (= (length (symbolic-poly-vars new)) (poly-dimension new))))66 67 60 68 61 #| … … 168 161 (format nil "[~{~a~^, ~}]" (mapcar #'(lambda (p) (poly->string p)) (cdr self)))))) 169 162 (:method ((self poly) &optional (vars nil)) 170 ;; Ensure that the number of variables matches the dimension171 (assert (= (length vars) (poly-dimension self)))172 163 (infix-print-to-string (->sexp self vars))) 173 164 (:method ((self symbolic-poly) &optional (vars (symbolic-poly-vars self))) -
branches/f4grobner/test0.lisp
r4442 r4463 4 4 5 5 (setf p (ADD-TO 6 (make-instance 'poly : dimension 1 :termlist (list (make-instance 'term :exponents #(0) :coeff 5)) :order #'lex>)7 (make-instance 'poly : dimension 1 :termlist (list (make-instance 'term :exponents #(1) :coeff 1)) :order #'lex>)))6 (make-instance 'poly :termlist (list (make-instance 'term :exponents #(0) :coeff 5)) :order #'lex>) 7 (make-instance 'poly :termlist (list (make-instance 'term :exponents #(1) :coeff 1)) :order #'lex>))) 8 8 9 9 (setf r (POLYNOMIAL::ADD-TERMLISTS (list (make-instance 'term :exponents #(0) :coeff 5)) -
branches/f4grobner/test8.lisp
r4447 r4463 1 1 (in-package :polynomial) 2 2 3 (proclaim '(special p0 q0 p q p+q-good p+q-risky)) 3 4 4 (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)))))5 (setf q0 (alist->poly ( nreverse '(((2 0 0) . 2) ((0 1 0) . -1) ((0 0 0) . 2)))))5 (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))))) 6 (setf q0 (alist->poly (reverse '(((2 0 0) . 2) ((0 1 0) . -1) ((0 0 0) . 2))))) 6 7 7 8 (setf p (poly-termlist p0)) 8 9 (setf q (poly-termlist q0)) 9 10 10 (setf p+q-good (mapcar #'->list (slow-add p q#'lex> #'add-to)))11 (setf p+q-good (mapcar #'->list (slow-add (mapcar #'copy:copy-instance p) (mapcar #'copy:copy-instance q) #'lex> #'add-to))) 11 12 12 (setf p+q-risky (mapcar #'->list (fast-and-risky-add p q#'lex> #'add-to)))13 (setf p+q-risky (mapcar #'->list (fast-and-risky-add (mapcar #'copy:copy-instance p) (mapcar #'copy:copy-instance q) #'lex> #'add-to))) 13 14 14 ;;(assert (equal p+q-good p+q-risky) nil "Two methods gave different results.")15 (assert (equal p+q-good p+q-risky) nil "Two methods gave different results.")
Note:
See TracChangeset
for help on using the changeset viewer.