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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra idecomp.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package IDECOMP IdealDecompositionPackage}
<<package IDECOMP IdealDecompositionPackage>>=
)abbrev package IDECOMP IdealDecompositionPackage
++ Author: P. Gianni
++ Date Created: summer 1986
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors: PolynomialIdeals
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This package provides functions for the primary decomposition of
++ polynomial ideals over the rational numbers. The ideals are members
++ of the \spadtype{PolynomialIdeals} domain, and the polynomial generators are
++ required to be from the \spadtype{DistributedMultivariatePolynomial} domain.

IdealDecompositionPackage(vl,nv) : C == T -- take away nv, now doesn't
                                          -- compile if it isn't there
 where
   vl      :  List Symbol
   nv      :  NonNegativeInteger
   Z      ==>  Integer  -- substitute with PFE cat
   Q      ==>  Fraction Z
   F      ==>  Fraction P
   P      ==>  Polynomial Z
   UP     ==>  SparseUnivariatePolynomial P
   Expon  ==>  DirectProduct(nv,NNI)
   OV     ==>  OrderedVariableList(vl)
   SE     ==>  Symbol
   SUP    ==>  SparseUnivariatePolynomial(DPoly)

   DPoly1 ==>  DistributedMultivariatePolynomial(vl,Q)
   DPoly  ==>  DistributedMultivariatePolynomial(vl,F)
   NNI    ==>  NonNegativeInteger

   Ideal  ==  PolynomialIdeals(Q,Expon,OV,DPoly1)
   FIdeal ==  PolynomialIdeals(F,Expon,OV,DPoly)
   Fun0   ==  Union("zeroPrimDecomp","zeroRadComp")
   GenPos ==  Record(changeval:List Z,genideal:FIdeal)

   C == with


     zeroDimPrime?       :        Ideal         -> Boolean
       ++ zeroDimPrime?(I) tests if the ideal I is a 0-dimensional prime.

     zeroDimPrimary?     :        Ideal         -> Boolean
       ++ zeroDimPrimary?(I) tests if the ideal I is 0-dimensional primary.
     prime?              :        Ideal         -> Boolean
       ++ prime?(I) tests if the ideal I is prime.
     radical             :        Ideal         -> Ideal
       ++ radical(I) returns the radical of the ideal I.
     primaryDecomp       :        Ideal         -> List(Ideal)
       ++ primaryDecomp(I) returns a list of primary ideals such that their
       ++ intersection is the ideal I.

     contract        : (Ideal,List OV   )       -> Ideal
       ++ contract(I,lvar) contracts the ideal I to the polynomial ring
       ++ \spad{F[lvar]}.

   T  == add

     import MPolyCatRationalFunctionFactorizer(Expon,OV,Z,DPoly)
     import GroebnerPackage(F,Expon,OV,DPoly)
     import GroebnerPackage(Q,Expon,OV,DPoly1)

                  ----  Local  Functions  -----
     genPosLastVar       :    (FIdeal,List OV)     -> GenPos
     zeroPrimDecomp      :    (FIdeal,List OV)     -> List(FIdeal)
     zeroRadComp         :    (FIdeal,List OV)     -> FIdeal
     zerodimcase         :    (FIdeal,List OV)     -> Boolean
     is0dimprimary       :    (FIdeal,List OV)     -> Boolean
     backGenPos          : (FIdeal,List Z,List OV) -> FIdeal
     reduceDim           : (Fun0,FIdeal,List OV)   -> List FIdeal
     findvar             :   (FIdeal,List OV)      -> OV
     testPower           :    (SUP,OV,FIdeal)      -> Boolean
     goodPower           :     (DPoly,FIdeal)  -> Record(spol:DPoly,id:FIdeal)
     pushdown            :      (DPoly,OV)        -> DPoly
     pushdterm           :     (DPoly,OV,Z)       -> DPoly
     pushup              :      (DPoly,OV)        -> DPoly
     pushuterm           :    (DPoly,SE,OV)       -> DPoly
     pushucoef           :       (UP,OV)          -> DPoly
     trueden             :        (P,SE)          -> P
     rearrange           :       (List OV)        -> List OV
     deleteunit          :      List FIdeal        -> List FIdeal
     ismonic             :      (DPoly,OV)        -> Boolean


     MPCFQF ==> MPolyCatFunctions2(OV,Expon,Expon,Q,F,DPoly1,DPoly)
     MPCFFQ ==> MPolyCatFunctions2(OV,Expon,Expon,F,Q,DPoly,DPoly1)

     convertQF(a:Q) : F == ((numer a):: F)/((denom a)::F)
     convertFQ(a:F) : Q == (ground numer a)/(ground denom a)

     internalForm(I:Ideal) : FIdeal ==
       Id:=generators I
       nId:=[map(convertQF,poly)$MPCFQF for poly in Id]
       groebner? I => groebnerIdeal nId
       ideal nId

     externalForm(I:FIdeal) : Ideal ==
       Id:=generators I
       nId:=[map(convertFQ,poly)$MPCFFQ for poly in Id]
       groebner? I => groebnerIdeal nId
       ideal nId

     lvint:=[variable(xx)::OV for xx in vl]
     nvint1:=(#lvint-1)::NNI

     deleteunit(lI: List FIdeal) : List FIdeal ==
       [I for I in lI | not element?(1$DPoly,I)]

     rearrange(vlist:List OV) :List OV ==
       vlist=[] => vlist
       sort(#1>#2,setDifference(lvint,setDifference(lvint,vlist)))

            ---- radical of a 0-dimensional ideal  ----
     zeroRadComp(I:FIdeal,truelist:List OV) : FIdeal ==
       truelist=[] => I
       Id:=generators I
       x:OV:=truelist.last
       #Id=1 =>
         f:=Id.first
	 g:= (f exquo (gcd (f,differentiate(f,x))))::DPoly
         groebnerIdeal([g])
       y:=truelist.first
       px:DPoly:=x::DPoly
       py:DPoly:=y::DPoly
       f:=Id.last
       g:= (f exquo (gcd (f,differentiate(f,x))))::DPoly
       Id:=groebner(cons(g,remove(f,Id)))
       lf:=Id.first
       pv:DPoly:=0
       pw:DPoly:=0
       while not one? degree(lf,y) repeat
         val:=random(23)$Z
         pv:=px+val*py
         pw:=px-val*py
	 Id:=groebner([(univariate(h,x)).pv for h in Id])
         lf:=Id.first
       ris:= generators(zeroRadComp(groebnerIdeal(Id.rest),truelist.rest))
       ris:=cons(lf,ris)
       if not zero? pv then
	 ris:=[(univariate(h,x)).pw for h in ris]
       groebnerIdeal(groebner ris)

          ----  find the power that stabilizes (I:s)  ----
     goodPower(s:DPoly,I:FIdeal) : Record(spol:DPoly,id:FIdeal) ==
       f:DPoly:=s
       I:=groebner I
       J:=generators(JJ:= (saturate(I,s)))
       while not in?(ideal([f*g for g in J]),I) repeat f:=s*f
       [f,JJ]

              ----  is the ideal zerodimensional?  ----
       ----     the "true variables" are  in truelist         ----
     zerodimcase(J:FIdeal,truelist:List OV) : Boolean ==
       element?(1,J) => true
       truelist=[] => true
       n:=#truelist
       Jd:=groebner generators J
       for x in truelist while Jd~=[] repeat
         f := Jd.first
         Jd:=Jd.rest
	 if ((y:=mainVariable f) case "failed") or (y::OV ~=x )
              or not (ismonic (f,x)) then return false
	 while Jd~=[] and (mainVariable Jd.first)::OV=x repeat Jd:=Jd.rest
	 if Jd=[] and position(x,truelist)<n then return false
       true

         ----  choose the variable for the reduction step  ----
                    --- J groebnerner in gen pos  ---
     findvar(J:FIdeal,truelist:List OV) : OV ==
       lmonicvar:List OV :=[]
       for f in generators J repeat
	 t:=f - reductum f
	 vt:List OV :=variables t
         if #vt=1 then lmonicvar:=setUnion(vt,lmonicvar)
       badvar:=setDifference(truelist,lmonicvar)
       badvar.first

            ---- function for the "reduction step  ----
     reduceDim(flag:Fun0,J:FIdeal,truelist:List OV) : List(FIdeal) ==
       element?(1,J) => [J]
       zerodimcase(J,truelist) =>
         (flag case "zeroPrimDecomp") => zeroPrimDecomp(J,truelist)
         (flag case "zeroRadComp") => [zeroRadComp(J,truelist)]
       x:OV:=findvar(J,truelist)
       Jnew:=[pushdown(f,x) for f in generators J]
       Jc: List FIdeal :=[]
       Jc:=reduceDim(flag,groebnerIdeal Jnew,remove(x,truelist))
       res1:=[ideal([pushup(f,x) for f in generators idp]) for idp in Jc]
       s:=pushup((_*/[leadingCoefficient f for f in Jnew])::DPoly,x)
       degree(s,x)=0 => res1
       res1:=[saturate(II,s) for II in res1]
       good:=goodPower(s,J)
       sideal := groebnerIdeal(groebner(cons(good.spol,generators J)))
       in?(good.id, sideal) => res1
       sresult:=reduceDim(flag,sideal,truelist)
       for JJ in sresult repeat
          if not(in?(good.id,JJ)) then res1:=cons(JJ,res1)
       res1

      ----  Primary Decomposition for 0-dimensional ideals  ----
     zeroPrimDecomp(I:FIdeal,truelist:List OV): List(FIdeal) ==
       truelist=[] => list I
       newJ:=genPosLastVar(I,truelist);lval:=newJ.changeval;
       J:=groebner newJ.genideal
       x:=truelist.last
       Jd:=generators J
       g:=Jd.last
       lfact:= factors factor(g)
       ris:List FIdeal:=[]
       for ef in lfact repeat
         g:DPoly:=(ef.factor)**(ef.exponent::NNI)
         J1:= groebnerIdeal(groebner cons(g,Jd))
         if not (is0dimprimary (J1,truelist)) then
                                   return zeroPrimDecomp(I,truelist)
         ris:=cons(groebner backGenPos(J1,lval,truelist),ris)
       ris

             ----  radical of an Ideal  ----
     radical(I:Ideal) : Ideal ==
       J:=groebner(internalForm I)
       truelist:=rearrange("setUnion"/[variables f for f in generators J])
       truelist=[] => externalForm J
       externalForm("intersect"/reduceDim("zeroRadComp",J,truelist))


-- the following functions are used to "push" x in the coefficient ring -

        ----  push x in the coefficient domain for a polynomial ----
     pushdown(g:DPoly,x:OV) : DPoly ==
       rf:DPoly:=0$DPoly
       i:=position(x,lvint)
       while not zero? g repeat
	 g1:=reductum g
         rf:=rf+pushdterm(g-g1,x,i)
         g := g1
       rf

      ----  push x in the coefficient domain for a term ----
     pushdterm(t:DPoly,x:OV,i:Z):DPoly ==
       n:=degree(t,x)
       xp:=convert(x)@SE
       cf:=monomial(1,xp,n)$P :: F
       newt := t exquo monomial(1,x,n)$DPoly
       cf * newt::DPoly

               ----  push back the variable  ----
     pushup(f:DPoly,x:OV) :DPoly ==
       h:=1$P
       rf:DPoly:=0$DPoly
       g := f
       xp := convert(x)@SE
       while not zero? g repeat
         h:=lcm(trueden(denom leadingCoefficient g,xp),h)
         g:=reductum g
       f:=(h::F)*f
       while not zero? f repeat
	 g:=reductum f
         rf:=rf+pushuterm(f-g,xp,x)
         f:=g
       rf

     trueden(c:P,x:SE) : P ==
       degree(c,x) = 0 => 1
       c

      ----  push x back from the coefficient domain for a term ----
     pushuterm(t:DPoly,xp:SE,x:OV):DPoly ==
       pushucoef((univariate(numer leadingCoefficient t,xp)$P), x)*
	  monomial(inv((denom leadingCoefficient t)::F),degree t)$DPoly


     pushucoef(c:UP,x:OV):DPoly ==
       c = 0 => 0
       monomial((leadingCoefficient c)::F::DPoly,x,degree c) +
		 pushucoef(reductum c,x)

           -- is the 0-dimensional ideal I primary ?  --
               ----  internal function  ----
     is0dimprimary(J:FIdeal,truelist:List OV) : Boolean ==
       element?(1,J) => true
       Jd:=generators(groebner J)
       not one?(#(factors factor Jd.last)) => return false
       i:=subtractIfCan(#truelist,1)
       (i case nothing) => return true
       k: NNI := i
       JR:=(reverse Jd);JM:=groebnerIdeal([JR.first]);JP:List(DPoly):=[]
       for f in JR.rest repeat
         if not ismonic(f,truelist.k) then
           if not inRadical?(f,JM) then return false
           JP:=cons(f,JP)
          else
           x:=truelist.k
           k := (k-1)::NNI
	   if not testPower(univariate(f,x),x,JM) then return false
           JM :=groebnerIdeal(append(cons(f,JP),generators JM))
       true

         ---- Functions for the General Position step  ----

       ----  put the ideal in general position  ----
     genPosLastVar(J:FIdeal,truelist:List OV):GenPos ==
       x := last truelist ;lv1:List OV :=remove(x,truelist)
       ranvals:List(Z):=[random(23)$Z for vv in lv1]
       val:=+/[rv*(vv::DPoly)  for vv in lv1 for rv in ranvals]
       val:=val+(x::DPoly)
       [ranvals,groebnerIdeal(groebner([(univariate(p,x)).val
                             for p in generators J]))]$GenPos


             ----  convert back the ideal  ----
     backGenPos(I:FIdeal,lval:List Z,truelist:List OV) : FIdeal ==
       lval=[] => I
       x := last truelist ;lv1:List OV:=remove(x,truelist)
       val:=-(+/[rv*(vv::DPoly) for vv in lv1 for rv in lval])
       val:=val+(x::DPoly)
       groebnerIdeal
	   (groebner([(univariate(p,x)).val for p in generators I ]))

     ismonic(f:DPoly,x:OV) : Boolean == ground? leadingCoefficient(univariate(f,x))

         ---- test if f is power of a linear mod (rad J) ----
                    ----  f is monic  ----
     testPower(uf:SUP,x:OV,J:FIdeal) : Boolean ==
       df:=degree(uf)
       trailp:DPoly := inv(df:Z ::F) *coefficient(uf,(df-1)::NNI)
       linp:SUP:=(monomial(1$DPoly,1$NNI)$SUP +
		  monomial(trailp,0$NNI)$SUP)**df
       g:DPoly:=multivariate(uf-linp,x)
       inRadical?(g,J)


                    ----  Exported Functions  ----

           -- is the 0-dimensional ideal I prime ?  --
     zeroDimPrime?(I:Ideal) : Boolean ==
       J:=groebner((genPosLastVar(internalForm I,lvint)).genideal)
       element?(1,J) => true
       n:NNI:=#vl;i:NNI:=1
       Jd:=generators J
       #Jd~=n => false
       for f in Jd repeat
         if not ismonic(f,lvint.i) then return false
	 if i<n and not one?(degree univariate(f,lvint.i)) then return false
         i:=i+1
       g:=Jd.n
       #(lfact:=factors(factor g)) >1 => false
       lfact.1.exponent =1


           -- is the 0-dimensional ideal I primary ?  --
     zeroDimPrimary?(J:Ideal):Boolean ==
       is0dimprimary(internalForm J,lvint)

             ----  Primary Decomposition of I  -----

     primaryDecomp(I:Ideal) : List(Ideal) ==
       J:=groebner(internalForm I)
       truelist:=rearrange("setUnion"/[variables f for f in generators J])
       truelist=[] => [externalForm J]
       [externalForm II for II in reduceDim("zeroPrimDecomp",J,truelist)]

          ----  contract I to the ring with lvar variables  ----
     contract(I:Ideal,lvar: List OV) : Ideal ==
       Id:= generators(groebner I)
       empty?(Id) => I
       fullVars:= "setUnion"/[variables g for g in Id]
       fullVars = lvar => I
       n:= # lvar
       #fullVars < n  => error "wrong vars"
       n=0 => I
       newVars:= append([vv for vv in fullVars
                           | not member?(vv,lvar)]$List(OV),lvar)
       subsVars := [monomial(1,vv,1)$DPoly1 for vv in newVars]
       lJ:= [eval(g,fullVars,subsVars) for g in Id]
       J := groebner(lJ)
       J=[1] => groebnerIdeal J
       J=[0] => groebnerIdeal empty()
       J:=[f for f in J| member?(mainVariable(f)::OV,newVars)]
       fullPol :=[monomial(1,vv,1)$DPoly1 for vv in fullVars]
       groebnerIdeal([eval(gg,newVars,fullPol) for gg in J])

@
\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>>

<<package IDECOMP IdealDecompositionPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}