;;; -*-  Mode: Lisp -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>		 
;;;  		       								 
;;;  This program is free software; you can redistribute it and/or modify	 
;;;  it under the terms of the GNU General Public License as published by	 
;;;  the Free Software Foundation; either version 2 of the License, or		 
;;;  (at your option) any later version.					 
;;; 		       								 
;;;  This program is distributed in the hope that it will be useful,		 
;;;  but WITHOUT ANY WARRANTY; without even the implied warranty of		 
;;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the		 
;;;  GNU General Public License for more details.				 
;;; 		       								 
;;;  You should have received a copy of the GNU General Public License		 
;;;  along with this program; if not, write to the Free Software 		 
;;;  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.	 
;;;										 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Operations in ideal theory
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defpackage "IDEAL"
  (:use :cl :ring :monomial :order :term :polynomial :division :grobner-wrap :ring-and-order)
  (:export "POLY-DEPENDS-P"
	   "RING-INTERSECTION" 
	   "ELIMINATION-IDEAL"
	   "COLON-IDEAL"
	   "COLON-IDEAL-1"
	   "IDEAL-INTERSECTION"
	   "POLY-LCM"
	   "GROBNER-EQUAL"
	   "GROBNER-SUBSETP"
	   "GROBNER-MEMBER"
	   "IDEAL-SATURATION-1"
	   "IDEAL-SATURATION"
	   "IDEAL-POLYSATURATION-1"
	   "IDEAL-POLYSATURATION"
	   ))

(in-package :ideal)

;; Does the term depend on variable K?
(defun term-depends-p (term k)
  "Return T if the term TERM depends on variable number K."
  (monom-depends-p (term-monom term) k))

;; Does the polynomial P depend on variable K?
(defun poly-depends-p (p k)
  "Return T if the term polynomial P depends on variable number K."
  (some #'(lambda (term) (term-depends-p term k)) (poly-termlist p)))

(defun ring-intersection (plist k)
  "This function assumes that polynomial list PLIST is a Grobner basis
and it calculates the intersection with the ring R[X[K],...,X[N]], i.e.
it discards polynomials which depend on variables X[0], X[1], ..., X[K-1]."
  (dotimes (i k plist)
    (setf plist
	  (remove-if #'(lambda (p)
			 (poly-depends-p p i))
		     plist))))

(defun elimination-ideal (ring-and-order flist k &optional (top-reduction-only $poly_top_reduction_only) (start 0))
  (ring-intersection (reduced-grobner ring-and-order flist start top-reduction-only) k))

(defun colon-ideal (ring-and-order f g 
		    &optional 
		      (top-reduction-only $poly_top_reduction_only)
		    &aux
		      (ring (ro-ring ring-and-order)))
  "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G),
where F and G are two lists of polynomials. The colon ideal I:J is
defined as the set of polynomials H such that for all polynomials W in
J the polynomial W*H belongs to I."
  (declare (type ring-and-order ring-and-order))
  (cond
    ((endp g)
     ;;Id(G) consists of 0 only so W*0=0 belongs to Id(F)
     (if (every #'poly-zerop f)
         (error "First ideal must be non-zero.")
         (list (make-poly-from-termlist
                (list (make-term
                       (make-monom :dimension (monom-dimension (poly-lm (find-if-not #'poly-zerop f))))
                       (funcall (ring-unit ring))))))))
    ((endp (cdr g))
     (colon-ideal-1 ring-and-order f (car g) top-reduction-only))
    (t
     (ideal-intersection ring-and-order
                         (colon-ideal-1 ring-and-order f (car g) top-reduction-only)
                         (colon-ideal ring-and-order f (rest g) top-reduction-only)
                         top-reduction-only))))

(defun colon-ideal-1 (ring-and-order f g 
		      &optional 
			(top-reduction-only $poly_top_reduction_only))
  "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
F is a list of polynomials and G is a polynomial."
  (declare (type ring-and-order ring-and-order))
  (mapcar #'(lambda (x) 
	      (poly-exact-divide ring-and-order x g)) 
	  (ideal-intersection ring-and-order f (list g) top-reduction-only)))

(defun ideal-intersection (ring-and-order f g 
			   &optional 
			     (top-reduction-only $poly_top_reduction_only)
			     (ring (ro-ring ring-and-order)))
  (declare (type ring-and-order ring-and-order))
  (mapcar #'poly-contract
	  (ring-intersection
	   (reduced-grobner
	    ring-and-order
	    (append (mapcar #'(lambda (p) (poly-extend p (make-monom :dimension 1 :initial-exponent 1))) f)
		    (mapcar #'(lambda (p)
				(poly-append (poly-extend (poly-uminus ring p)
							  (make-monom :dimension 1 :initial-exponent 1))
					     (poly-extend p)))
			    g))
	    0
	    top-reduction-only)
	   1)))

(defun poly-lcm (ring-and-order f g &aux (ring (ro-ring ring-and-order)))
  "Return LCM (least common multiple) of two polynomials F and G.
The polynomials must be ordered according to monomial order PRED
and their coefficients must be compatible with the RING structure
defined in the COEFFICIENT-RING package."
  (cond
    ((poly-zerop f) f)
    ((poly-zerop g) g)
    ((and (endp (cdr (poly-termlist f))) (endp (cdr (poly-termlist g))))
     (let ((m (monom-lcm (poly-lm f) (poly-lm g))))
       (make-poly-from-termlist (list (make-term m (funcall (ring-lcm ring) (poly-lc f) (poly-lc g)))))))
    (t
     (multiple-value-bind (f f-cont)
	 (poly-primitive-part ring f)
       (multiple-value-bind (g g-cont)
	   (poly-primitive-part ring g)
	 (scalar-times-poly
	  ring
	  (funcall (ring-lcm ring) f-cont g-cont)
	  (poly-primitive-part ring (car (ideal-intersection ring-and-order (list f) (list g) nil)))))))))

;; Do two Grobner bases yield the same ideal?
(defun grobner-equal (ring-and-order g1 g2)
  "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
generate  the same ideal, and NIL otherwise."
  (declare (type ring-and-order ring-and-order))
  (and (grobner-subsetp ring-and-order g1 g2) (grobner-subsetp ring-and-order g2 g1)))

(defun grobner-subsetp (ring-and-order 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."
  (declare (type ring-and-order ring-and-order))
  (every #'(lambda (p) (grobner-member ring-and-order p g2)) g1))

(defun grobner-member (ring-and-order 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."
  (declare (type ring-and-order ring-and-order))
  (poly-zerop (normal-form ring-and-order p g nil)))

;; Calculate F : p^inf
(defun ideal-saturation-1 (ring-and-order f p 
			   &optional 
			     (start 0)
			     (top-reduction-only $poly_top_reduction_only)
			   &aux
			     (ring (ro-ring ring-and-order)))
  "Returns the reduced Grobner basis of the saturation of the ideal
generated by a polynomial list F in the ideal generated by a single
polynomial P. The saturation ideal is defined as the set of
polynomials H such for some natural number n (* (EXPT P N) H) is in
the ideal spanned by F. Geometrically, over an algebraically closed
field, this is the set of polynomials in the ideal generated by F
which do not identically vanish on the variety of P."
  (declare (type ring-and-order ring-and-order))
  (mapcar
   #'poly-contract
   (ring-intersection
    (reduced-grobner
     ring-and-order
     (saturation-extension-1 ring f p)
     start top-reduction-only)
    1)))


;; Calculate F : p1^inf : p2^inf : ... : ps^inf
(defun ideal-polysaturation-1 (ring-and-order f plist 
			       &optional 
				 (start 0)
				 (top-reduction-only $poly_top_reduction_only))
  "Returns the reduced Grobner basis of the ideal obtained by a
sequence of successive saturations in the polynomials
of the polynomial list PLIST of the ideal generated by the
polynomial list F."
  (cond
    ((endp plist) 
     (reduced-grobner ring-and-order f start top-reduction-only))
    (t (let ((g (ideal-saturation-1 ring-and-order f (car plist) start top-reduction-only)))
	 (ideal-polysaturation-1 ring-and-order g (rest plist) (length g) top-reduction-only)))))

(defun ideal-saturation (ring-and-order f g
			 &optional 
			   (start 0)
			   (top-reduction-only $poly_top_reduction_only)
			 &aux
			   (k (length g))
			   (ring (ro-ring ring-and-order)))
  "Returns the reduced Grobner basis of the saturation of the ideal
generated by a polynomial list F in the ideal generated a polynomial
list G. The saturation ideal is defined as the set of polynomials H
such for some natural number n and some P in the ideal generated by G
the polynomial P**N * H is in the ideal spanned by F.  Geometrically,
over an algebraically closed field, this is the set of polynomials in
the ideal generated by F which do not identically vanish on the
variety of G."
  (declare (type ring-and-order ring-and-order))
  (mapcar
   #'(lambda (q) (poly-contract q k))
   (ring-intersection
    (reduced-grobner ring-and-order
		     (polysaturation-extension ring f g)
		     start
		     top-reduction-only)
    k)))

(defun ideal-polysaturation (ring-and-order f ideal-list
			     &optional 
			       (start 0)
			       (top-reduction-only $poly_top_reduction_only))
  "Returns the reduced Grobner basis of the ideal obtained by a
successive applications of IDEAL-SATURATION to F and lists of
polynomials in the list IDEAL-LIST."
  (declare (type ring-and-order ring-and-order))
  (cond
    ((endp ideal-list) f)
    (t (let ((h (ideal-saturation ring-and-order f (car ideal-list) start top-reduction-only)))
	 (ideal-polysaturation ring-and-order h (rest ideal-list) (length h) top-reduction-only)))))