source: branches/f4grobner/ideal.lisp@ 4061

;;; -*-  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 :monom :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 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=(P[0],P[1],...,P[J-1]) is a Grobner basis and it calculates the

intersection of Id({P[0],P[1],...,P[J-1]}) with the ring

R[X[K],...,X[N-1]], 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))

  "Given a list of polynomials FLIST, and an integer K, tt finds and

returns the Groebner basis the elimination ideal of Id({FLIST})

obtained by eliminating the first K variables. Optional argument

TOP-REDUCTION-ONLY indicates whether to fully reduce or only

top-reduce. Optional argument START, defaulting to 0, is used to

indicate that the first START elements of F form a Groebner basis."

  (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

                       :monom (make-monom :dimension (monom-dimension (poly-lm (find-if-not #'poly-zerop f))))

                       :coeff (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 :monom m 

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