path: root/src/algebra/ideal.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra ideal.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain IDEAL PolynomialIdeals}
<<domain IDEAL PolynomialIdeals>>=
)abbrev domain IDEAL PolynomialIdeals
++ Author: P. Gianni
++ Date Created: summer 1986
++ Date Last Updated: September 1996
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References: GTZ
++ Description: This domain represents polynomial ideals with coefficients in any
++ field and supports the basic ideal operations, including intersection
++ sum and quotient.
++ An ideal is represented by a list of polynomials (the generators of
++ the ideal) and a boolean that is true if the generators are a Groebner
++ basis.
++ The algorithms used are based on Groebner basis computations. The
++ ordering is determined by the datatype of the input polynomials.
++ Users may use refinements of total degree orderings.

PolynomialIdeals(F,Expon,VarSet,DPoly) : C == T
 where
   F         :  Field
   Expon     :  OrderedAbelianMonoidSup
   VarSet    :  OrderedSet
   DPoly     :  PolynomialCategory(F,Expon,VarSet)

   SUP      ==> SparseUnivariatePolynomial(DPoly)
   NNI      ==> NonNegativeInteger
   Z        ==> Integer
   P        ==> Polynomial F
   MF       ==> Matrix(F)
   ST       ==> SuchThat(List P, List Equation P)

   GenMPos  ==> Record(mval:MF,invmval:MF,genIdeal:Ideal)
   Ideal    ==> %

   C == SetCategory with

     "*"             :  (Ideal,Ideal)        -> Ideal
       ++ I*J computes the product of the ideal I and J.
     "**"            :   (Ideal,NNI)         -> Ideal
       ++ I**n computes the nth power of the ideal I.
     "+"             :  (Ideal,Ideal)        -> Ideal
       ++ I+J computes the ideal generated by the union of I and J.
     one?            :     Ideal             -> Boolean
       ++ one?(I) tests whether the ideal I is the unit ideal,
       ++ i.e. contains 1.
     zero?           :     Ideal             -> Boolean
       ++ zero?(I) tests whether the ideal I is the zero ideal
     element?        :  (DPoly,Ideal)        -> Boolean
       ++ element?(f,I) tests whether the polynomial f belongs to
       ++ the ideal I.
     in?             :  (Ideal,Ideal)        -> Boolean
       ++ in?(I,J) tests if the ideal I is contained in the ideal J.
     inRadical?      :  (DPoly,Ideal)        -> Boolean
       ++ inRadical?(f,I) tests if some power of the polynomial f
       ++ belongs to the ideal I.
     zeroDim?        : (Ideal,List VarSet)   -> Boolean
       ++ zeroDim?(I,lvar) tests if the ideal I is zero dimensional, i.e.
       ++ all its associated primes are maximal,
       ++ in the ring \spad{F[lvar]}
     zeroDim?        :     Ideal             -> Boolean
       ++ zeroDim?(I) tests if the ideal I is zero dimensional, i.e.
       ++ all its associated primes are maximal,
       ++ in the ring \spad{F[lvar]}, where lvar are the variables appearing  in I
     intersect       :  (Ideal,Ideal)        -> Ideal
       ++ intersect(I,J) computes the intersection of the ideals I and J.
     intersect       :   List(Ideal)         -> Ideal
       ++ intersect(LI) computes the intersection of the list of ideals LI.
     quotient        :  (Ideal,Ideal)        -> Ideal
       ++ quotient(I,J) computes the quotient of the ideals I and J, \spad{(I:J)}.
     quotient        :  (Ideal,DPoly)        -> Ideal
       ++ quotient(I,f) computes the quotient of the ideal I by the principal
       ++ ideal generated by the polynomial f, \spad{(I:(f))}.
     groebner        :     Ideal             -> Ideal
       ++ groebner(I) returns a set of generators of I that are a Groebner basis
       ++ for I.
     generalPosition :  (Ideal,List VarSet)      -> GenMPos
       ++ generalPosition(I,listvar) perform a random linear
       ++ transformation on the variables in listvar and returns
       ++ the transformed ideal along with the change of basis matrix.
     backOldPos      :     GenMPos           -> Ideal
       ++ backOldPos(genPos) takes the result
       ++ produced by \spadfunFrom{generalPosition}{PolynomialIdeals}
       ++ and performs the inverse transformation, returning the original ideal
       ++ \spad{backOldPos(generalPosition(I,listvar))} = I.
     dimension       : (Ideal,List VarSet)   -> Z
       ++ dimension(I,lvar) gives the dimension of the ideal I,
       ++ in the ring \spad{F[lvar]}
     dimension       :      Ideal            -> Z
       ++ dimension(I) gives the dimension of the ideal I.
       ++ in the ring \spad{F[lvar]}, where lvar are the variables appearing  in I
     leadingIdeal    :     Ideal             -> Ideal
       ++ leadingIdeal(I) is the ideal generated by the
       ++ leading terms of the elements of the ideal I.
     ideal           :   List DPoly          ->  Ideal
       ++ ideal(polyList) constructs the ideal generated by the list
       ++ of polynomials polyList.
     groebnerIdeal       :   List DPoly          ->  Ideal
       ++ groebnerIdeal(polyList) constructs the ideal generated by the list
       ++ of polynomials polyList which are assumed to be a Groebner
       ++ basis.
       ++ Note: this operation avoids a Groebner basis computation.
     groebner?           :     Ideal             ->  Boolean
       ++ groebner?(I) tests if the generators of the ideal I are a Groebner basis.
     generators      :     Ideal             ->  List DPoly
       ++ generators(I) returns a list of generators for the ideal I.
     coerce          :   List DPoly          ->  Ideal
       ++ coerce(polyList) converts the list of polynomials polyList to an ideal.

     saturate        :  (Ideal,DPoly)        -> Ideal
       ++ saturate(I,f) is the saturation of the ideal I
       ++ with respect to the multiplicative
       ++ set generated by the polynomial f.
     saturate        :(Ideal,DPoly,List VarSet)  -> Ideal
       ++ saturate(I,f,lvar) is the saturation with respect to the prime
       ++ principal ideal which is generated by f in the polynomial ring
       ++ \spad{F[lvar]}.
     if VarSet has ConvertibleTo Symbol then
       relationsIdeal  :    List DPoly         -> ST
         ++ relationsIdeal(polyList) returns the ideal of relations among the
         ++ polynomials in polyList.

   T  == add

   ---  Representation ---
     Rep := Record(idl:List DPoly,isGr:Boolean)


                ----  Local Functions  ----

     contractGrob  :    newIdeal          ->  Ideal
     npoly         :     DPoly            ->  newPoly
     oldpoly       :     newPoly          ->  Union(DPoly,"failed")
     leadterm      :   (DPoly,VarSet)     ->  DPoly
     choosel       :   (DPoly,DPoly)      ->  DPoly
     isMonic?      :   (DPoly,VarSet)     ->  Boolean
     randomat      :     List Z           ->  Record(mM:MF,imM:MF)
     monomDim      : (Ideal,List VarSet)  ->  NNI
     variables     :       Ideal          ->  List VarSet
     subset        :     List VarSet      ->  List List VarSet
     makeleast     : (List VarSet,List VarSet)  ->  List VarSet

     newExpon:  OrderedAbelianMonoidSup
     newExpon:= Product(NNI,Expon)
     newPoly := PolynomialRing(F,newExpon)

     import GaloisGroupFactorizer(SparseUnivariatePolynomial Z)
     import GroebnerPackage(F,Expon,VarSet,DPoly)
     import GroebnerPackage(F,newExpon,VarSet,newPoly)

     newIdeal ==> List(newPoly)

     npoly(f:DPoly) : newPoly ==
       f=0$DPoly => 0$newPoly
       monomial(leadingCoefficient f,makeprod(0,degree f))$newPoly +
             npoly(reductum f)

     oldpoly(q:newPoly) : Union(DPoly,"failed") ==
       q=0$newPoly => 0$DPoly
       dq:newExpon:=degree q
       n:NNI:=selectfirst (dq)
       n~=0 => "failed"
       ((g:=oldpoly reductum q) case "failed") => "failed"
       monomial(leadingCoefficient q,selectsecond dq)$DPoly + (g::DPoly)

     leadterm(f:DPoly,lvar:List VarSet) : DPoly ==
       empty?(lf:=variables f)  or lf=lvar => f
       leadterm(leadingCoefficient univariate(f,lf.first),lvar)

     choosel(f:DPoly,g:DPoly) : DPoly ==
       g=0 => f
       (f1:=f exquo g) case "failed" => f
       choosel(f1::DPoly,g)

     contractGrob(I1:newIdeal) : Ideal ==
       J1:List(newPoly):=groebner(I1)
       while (oldpoly J1.first) case "failed" repeat J1:=J1.rest
       [[(oldpoly f)::DPoly for f in J1],true]

     makeleast(fullVars: List VarSet,leastVars:List VarSet) : List VarSet ==
       n:= # leastVars
       #fullVars < n  => error "wrong vars"
       n=0 => fullVars
       append([vv for vv in fullVars| not member?(vv,leastVars)],leastVars)

     isMonic?(f:DPoly,x:VarSet) : Boolean ==
       ground? leadingCoefficient univariate(f,x)

     subset(lv : List VarSet) : List List VarSet ==
       #lv =1 => [lv,empty()]
       v:=lv.1
       ll:=subset(rest lv)
       l1:=[concat(v,set) for set in ll]
       concat(l1,ll)

     monomDim(listm:Ideal,lv:List VarSet) : NNI ==
       monvar: List List VarSet := []
       for f in generators listm repeat
         mvset := variables f
         #mvset > 1 => monvar:=concat(mvset,monvar)
         lv:=delete(lv,position(mvset.1,lv))
       empty? lv => 0
       lsubset : List List VarSet := sort(#(#1)>#(#2),subset(lv))
       for subs in lsubset repeat
         ldif:List VarSet:= lv
         for mvset in monvar while ldif ~=[] repeat
           ldif:=setDifference(mvset,subs)
         if not (empty? ldif) then  return #subs
       0

               --    Exported  Functions   ----

                 ----  is  I =  J  ?  ----
     (I:Ideal = J:Ideal) == in?(I,J) and in?(J,I)

               ----  check if f is in I  ----
     element?(f:DPoly,I:Ideal) : Boolean ==
       Id:=(groebner I).idl
       empty? Id => f = 0
       normalForm(f,Id) = 0

             ---- check if I is contained in J  ----
     in?(I:Ideal,J:Ideal):Boolean ==
       J:= groebner J
       empty?(I.idl) => true
       "and"/[element?(f,J) for f in I.idl ]


            ----  groebner base for an Ideal  ----
     groebner(I:Ideal) : Ideal  ==
       I.isGr =>
         "or"/[not zero? f for f in I.idl] => I
         [empty(),true]
       [groebner I.idl ,true]

            ----  Intersection of two ideals  ----
     intersect(I:Ideal,J:Ideal) : Ideal ==
       empty?(Id:=I.idl) => I
       empty?(Jd:=J.idl) => J
       tp:newPoly := monomial(1,makeprod(1,0$Expon))$newPoly
       tp1:newPoly:= tp-1
       contractGrob(concat([tp*npoly f for f in Id],
                     [tp1*npoly f for f in Jd]))


            ----   intersection for a list of ideals  ----

     intersect(lid:List(Ideal)) : Ideal == "intersect"/[l for l in lid]

               ----  quotient by an element  ----
     quotient(I:Ideal,f:DPoly) : Ideal ==
       --[[(g exquo f)::DPoly for g in (intersect(I,[f]::%)).idl ],true]
        import GroebnerInternalPackage(F,Expon,VarSet,DPoly)
        [minGbasis [(g exquo f)::DPoly
                 for g in (intersect(I,[f]::%)).idl ],true]

                ----  quotient of two ideals  ----
     quotient(I:Ideal,J:Ideal) : Ideal ==
       Jdl := J.idl
       empty?(Jdl) => ideal [1]
       [("intersect"/[quotient(I,f) for f in Jdl ]).idl ,true]


                ----    sum of two ideals  ----
     (I:Ideal + J:Ideal) : Ideal == [groebner(concat(I.idl ,J.idl )),true]

                ----   product of two ideals  ----
     (I:Ideal * J:Ideal):Ideal ==
       [groebner([:[f*g for f in I.idl ] for g in J.idl ]),true]

                ----  power of an ideal  ----
     (I:Ideal ** n:NNI) : Ideal ==
       n=0 => [[1$DPoly],true]
       (I * (I**(n-1):NNI))

       ----  saturation with respect to the multiplicative set f**n ----
     saturate(I:Ideal,f:DPoly) : Ideal ==
       f=0 => error "f is zero"
       tp:newPoly := (monomial(1,makeprod(1,0$Expon))$newPoly * npoly f)-1
       contractGrob(concat(tp,[npoly g for g in I.idl ]))

     ----  saturation with respect to a prime principal ideal in lvar ---
     saturate(I:Ideal,f:DPoly,lvar:List(VarSet)) : Ideal ==
       Id := I.idl
       fullVars := "setUnion"/[variables g for g in Id]
       newVars:=makeleast(fullVars,lvar)
       subVars := [monomial(1,vv,1) for vv in newVars]
       J:List DPoly:=groebner([eval(g,fullVars,subVars) for g in Id])
       ltJ:=[leadterm(g,lvar) for g in J]
       s:DPoly:=_*/[choosel(ltg,f) for ltg in ltJ]
       fullPol:=[monomial(1,vv,1) for vv in fullVars]
       [[eval(g,newVars,fullPol) for g in (saturate(J::%,s)).idl],true]

            ----  is the ideal zero dimensional?  ----
            ----      in the ring F[lvar]?        ----
     zeroDim?(I:Ideal,lvar:List VarSet) : Boolean ==
       J:=(groebner I).idl
       empty? J => false
       J = [1] => false
       n:NNI := # lvar
       #J < n => false
       for f in J while not empty?(lvar) repeat
         x:=(mainVariable f)::VarSet
         if isMonic?(f,x) then lvar:=delete(lvar,position(x,lvar))
       empty?(lvar)

            ----  is the ideal zero dimensional?  ----
     zeroDim?(I:Ideal):Boolean == zeroDim?(I,"setUnion"/[variables g for g in I.idl])

               ----  test if f is in the radical of I  ----
     inRadical?(f:DPoly,I:Ideal) : Boolean ==
       f=0$DPoly => true
       tp:newPoly :=(monomial(1,makeprod(1,0$Expon))$newPoly * npoly f)-1
       Id:=I.idl
       normalForm(1$newPoly,groebner concat(tp,[npoly g for g in Id])) = 0

              ----   dimension of an ideal  ----
              ----    in the ring F[lvar]   ----
     dimension(I:Ideal,lvar:List VarSet) : Z ==
       I:=groebner I
       empty?(I.idl) => # lvar
       element?(1,I) => -1
       truelist:="setUnion"/[variables f for f in I.idl]
       "or"/[not member?(vv,lvar) for vv in truelist] => 
          error "wrong variables"
       truelist:=setDifference(lvar,setDifference(lvar,truelist))
       ed:Z:=#lvar - #truelist
       leadid:=leadingIdeal(I)
       n1:Z:=monomDim(leadid,truelist)::Z
       ed+n1

     dimension(I:Ideal) : Z == dimension(I,"setUnion"/[variables g for g in I.idl])

     -- leading term ideal --
     leadingIdeal(I : Ideal) : Ideal ==
       Idl:= (groebner I).idl
       [[(f-reductum f) for f in Idl],true]

               ---- ideal of relations among the fi  ----
     if VarSet has ConvertibleTo Symbol then

       monompol(df:List NNI,lcf:F,lv:List VarSet) : P ==
         g:P:=lcf::P
         for dd in df for v in lv repeat
           g:= monomial(g,convert v,dd)
         g

       relationsIdeal(listf : List DPoly): ST ==
	 empty? listf  => [empty(),empty()]$ST
	 nf:=#listf
	 lvint := "setUnion"/[variables g for g in listf]
	 vl: List Symbol := [convert vv for vv in lvint]
	 nvar:List Symbol:=[new() for i in 1..nf]
	 VarSet1:=OrderedVariableList(concat(vl,nvar))
	 lv1:=[variable(vv)$VarSet1::VarSet1 for vv in nvar]
	 DirP:=DirectProduct(nf,NNI)
	 nExponent:=Product(Expon,DirP)
	 nPoly := PolynomialRing(F,nExponent)
	 gp:=GroebnerPackage(F,nExponent,VarSet1,nPoly)
	 lf:List nPoly :=[]
	 lp:List P:=[]
	 for f in listf for i in 1..  repeat
	   vec2:Vector(NNI):=new(nf,0$NNI)
	   vec2.i:=1
	   g:nPoly:=0$nPoly
	   pol:=0$P
	   while f~=0 repeat
	     df:=degree(f-reductum f,lvint)
	     lcf:=leadingCoefficient f
	     pol:=pol+monompol(df,lcf,lvint)
	     g:=g+monomial(lcf,makeprod(degree f,0))$nPoly
	     f:=reductum f
	   lp:=concat(pol,lp)
	   lf:=concat(monomial(1,makeprod(0,directProduct vec2))-g,lf)
	 npol:List P :=[v::P for v in nvar]
	 leq : List Equation P :=
	       [p = pol for p in npol for pol in reverse lp ]
	 lf:=(groebner lf)$gp
	 while lf~=[] repeat
	   q:=lf.first
	   dq:nExponent:=degree q
	   n:=selectfirst (dq)
	   if n=0 then leave "done"
	   lf:=lf.rest
	 solsn:List P:=[]
	 for q in lf repeat
	   g:Polynomial F :=0
	   while q~=0 repeat
	     dq:=degree q
	     lcq:=leadingCoefficient q
	     q:=reductum q
	     vdq:=(selectsecond dq):Vector NNI
	     g:=g+ lcq*
		_*/[p**vdq.j for p in npol for j in 1..]
	   solsn:=concat(g,solsn)
	 [solsn,leq]$ST

     coerce(Id:List DPoly) : Ideal == [Id,false]

     coerce(I:Ideal) : OutputForm ==
       Idl := I.idl
       empty? Idl => [0$DPoly] :: OutputForm
       Idl :: OutputForm

     ideal(Id:List DPoly) :Ideal  ==  [[f for f in Id|f~=0],false]

     groebnerIdeal(Id:List DPoly) : Ideal == [Id,true]

     generators(I:Ideal) : List DPoly  == I.idl

     groebner?(I:Ideal) : Boolean  ==  I.isGr

     one?(I:Ideal) : Boolean == element?(1, I)

     zero?(I:Ideal) : Boolean == empty? (groebner I).idl

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<domain IDEAL PolynomialIdeals>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}