(defun s-poly (f g)
"Returns the S-polynomial of f and g"
;; S(f, g) = x^gamma/LT(f) * f - x^gamma/LT(g) * g where x^gamma is
;; LCM(LM(f), LM(g))
(flet ((make-x-gamma (f g)
(make-instance 'term
:coefficient 1
:monomial (map 'vector #'max
(monomial (lt f)) (monomial (lt g))))))
(let ((x-gamma (make-x-gamma f g)))
(sub (mul f (div x-gamma (lt f)))
(mul g (div x-gamma (lt g)))))))
(defun make-index-set (s)
(loop :for i :below s
:nconc (loop :for j :below s :when (< i j) :collect (cons i j))))
(defun pair-member (l m B)
(or (member (cons l m) B :test #'equal)
(member (cons m l) B :test #'equal)))
(defun criterion (G B i j fi fj)
"Returns T if S_{ij} ought to be considered, NIL otherise."
(loop :with status := nil
:for fk :across G
:for k :from 0
:when (and (/= k i) (/= k j)
(divides-p (lt fk) (ring-lcm (lt fi) (lt fj)))
(not (pair-member i k B)) (not (pair-member j k B)))
:do (setf status t)
(loop-finish)
:finally (return status)))
(defun groebner (polynomials)
"Returns a Groebner basis for the ideal generated by the specified
array of polynomials."
(let* ((s (array-dimension polynomials 0))
(r s)
(G (make-array s :initial-contents polynomials
:adjustable t :fill-pointer s)))
(do ((B (make-index-set s) (rest B)))
((null B) G)
(destructuring-bind (i . j) (first B)
(let ((fi (elt G i))
(fj (elt G j)))
(when (and (not (ring-equal-p (ring-lcm (lt fi) (lt fj))
(mul (lt fi) (lt fj))))
(not (criterion G B i j fi fj)))
(let ((remainder (normal-form (s-poly fi fj) G)))
(unless (ring-zero-p remainder)
(vector-push-extend remainder G)
(setf B (append B (loop :for i :below r :collect (cons i r))))
(incf r)))))))))
(defun reduce-gb (G)
"Returns a reduced Groebner basis."
(flet ((filter (f G)
(remove-if (lambda (x) (or (ring-zero-p x) (equal f x))) G)))
(loop :with F := (map 'vector (lambda (x) (mul x (/ (lc x)))) G)
:for i :below (length F) :for fi :across F :do
(setf (elt F i)
(if (some (lambda (x) (divides-p (lt x) (lt fi)))
(filter fi F))
(make-instance 'polynomial :ring (base-ring fi)) ; zero polynomial
(normal-form fi (remove fi F))))
:finally (return (remove-if #'ring-zero-p F)))))
(defun reduced-groebner (F)
"Computes and reduces a Groebner basis of the ideal generated by F."
(reduce-gb (groebner F)))
