(uiop:define-package :st-buchberger/src/groebner
 (:mix :cl)
 (:mix-reexport :st-buchberger/src/arithmetic)
 (:export #:s-poly #:make-index-set #:pair-member
  #:criterion #:normal-form #:groebner #:reduce-gb
  #:reduced-groebner))

(in-package :st-buchberger/src/groebner)

(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 normal-form (f G)
  (nth-value 1 (divmod f (remove-if #'ring-zero-p G))))

(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)))