source: git/Singular/LIB/ehv.lib @ a6606e6

// $Id: ehv.lib,v 1.3 2009-04-10 13:28:50 Singular Exp $
/////////////////////////////////////////////////////////////////////
// EHV.lib                                                         //
// algorithms for primary decomposition of ideals based on         //
// the algorithms of Eisenbud, Huneke, Vasconcelos                 //
// written by Kai Dehmann                                          //
//                                                                 //
/////////////////////////////////////////////////////////////////////

version="$Id: ehv.lib,v 1.3 2009-04-10 13:28:50 Singular Exp $";
category="Commutative Algebra";

info="
LIBRARY: EHV.lib      PROCEDURES FOR PRIMARY DECOMPOSITION OF IDEALS
AUTHORS: Kai Dehmann, dehmann@mathematik.uni-kl.de;

OVERVIEW:
  Algorithms for primary decomposition and radical-computation
  based on the ideas of Eisenbud, Huneke, and Vasconcelos.

PROCEDURES:
  equiMaxEHV(I);                      equidimensional part of I
  removeComponent(I,e);               intersection of the primary components
                                      of I of dimension >= e
  AssOfDim(I,e);                      an ideal such that the associated primes
                                      are exactly the associated primes of I
                                      having dimension e
  equiRadEHV(I [,Strategy]);          equidimensional radical of I
  radEHV(I [,Strategy]);              radical of I
  IntAssOfDim1(I,e);                  intersection of the associated primes of I
                                      having dimension e
  IntAssOfDim2(I,e);                  another way of computing the intersection
                                      of the associated primes of I
                                      having dimension e
  decompEHV(I);                       decomposition of a zero-dimensional
                                      radical ideal I
  AssEHV(I [,Strategy]);              associated primes of I
  minAssEHV(I [,Strategy]);           minimal associated primes of I
  localize(I,P,l);                    the contraction of the ideal generated by I
                                      in the localization w.r.t P
  componentEHV(I,P,L [,Strategy]);    a P-primary component for I
  primdecEHV(I [,Strategy]);          a minimal primary decomposition of I
  compareLists(L, K);                 procedure for comparing the output of
                                      primary decomposition algorithms (checks
                                      if the computed associated primes coincide)
";

LIB "ring.lib";
LIB "general.lib";
LIB "elim.lib";
LIB "poly.lib";
LIB "random.lib";
LIB "inout.lib";
LIB "matrix.lib";
LIB "algebra.lib";
LIB "normal.lib";


/////////////////////////////////////////////////////////////////////
//                                                                 //
//             G E N E R A L   A L G O R I T H M S                 //
//                                                                 //
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
static proc AnnExtEHV(int n,list re)
"USAGE:   AnnExtEHV(n,re); n integer, re resolution
RETURN:  ideal, the annihilator of Ext^n(R/I,R) with given
         resolution re of I"
{
  if(printlevel > 2){"Entering AnnExtEHV.";}

  if(n < 0)
    {
      ideal ann = ideal(1);
      if(printlevel > 2){"Leaving AnnExtEHV.";}
      return(ann);
    }
  int l = size(re);

  if(n < l)
    {
      matrix f = transpose(re[n+1]);
      if(n == 0)
        {
          matrix g = 0*gen(ncols(f));
        }
      else
        {
          matrix g = transpose(re[n]);
        }
      module k = syz(f);
      ideal ann = quotient1(g,k);
      if(printlevel > 2){"Leaving AnnExtEHV.";}
      return(ann);
    }

  if(n == l)
    {
      ideal ann = Ann(transpose(re[n]));
      if(printlevel > 2){"Leaving AnnExtEHV.";}
      return(ann);
    }

  ideal ann = ideal(1);
  if(printlevel > 2){"Leaving AnnExtEHV.";}
  return(ann);
}


/////////////////////////////////////////////////////////////////////
//static
proc isSubset(ideal I,ideal J)
"USAGE:   isSubset(I,J); I, J ideals
RETURN:  integer, 1 if I is a subset of J and 0 otherwise
NOTE:    if J is not a standard basis the result may be wrong
"
{
  int s = size(I);
  for(int i=1; i<=s; i++)
    {
      if(reduce(I[i],J)!=0)
        {
          return(0);
        }
    }
  return(1);
}


/////////////////////////////////////////////////////////////////////
//                                                                 //
//         T H E   E Q U I D I M E N S I O N A L   P A R T         //
//                                                                 //
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
proc equiMaxEHV(ideal I)
"USAGE:   equiMaxEHV(I); I ideal
RETURN:  ideal, the equidimensional part of I.
NOTE:    Uses algorithm of Eisenbud, Huneke, and Vasconcelos.
EXAMPLE: example equiMaxEHV; shows an example
"
{
  if(printlevel > 2){"Entering equiMaxEHV.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  ideal J = groebner(I);
  int cod = nvars(basering)-dim(J);

  //If I is the entire ring...
  if(cod > nvars(basering))
    {
      //...then return the ideal generated by 1.
      return(ideal(1));
    }

  //Compute a resolution of I.
  if(printlevel > 2){"Computing resolution.";}
  if(homog(I)==1)
    {
      list re = sres(J,cod+1);
      re = minres(re);
    }
  else
    {
      list re = mres(I,cod+1);
    }
  if(printlevel > 2){"Finished computing resolution.";}

  //Compute the annihilator of the cod-th EXT-module.
  ideal ann = AnnExtEHV(cod,re);
  attrib(ann,"isEquidimensional",1);
  if(printlevel > 2){"Leaving equiMaxEHV.";}
  return(ann);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
  equiMaxEHV(I);
}


/////////////////////////////////////////////////////////////////////
proc removeComponent(ideal I, int e)
"USAGE:   removeComponent(I,e); I ideal, e integer
RETURN:  ideal, the intersection of the primary components
         of I of dimension >= e
EXAMPLE: example removeComponent; shows an example"
{
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }

  ideal J = groebner(I);

  //Compute a resolution of I
  if(homog(I)==1)
    {
      list re = sres(J,0);
      re = minres(re);
    }
  else
    {
      list re = mres(I,0);
    }

  int f = nvars(basering);
  int cod;
  ideal ann;
  int g = nvars(basering) - e;
  while(f > g)
    {
      ann = AnnExtEHV(f,re);
      cod = nvars(basering) - dim(groebner(ann));
      if( cod == f )
        {
          I = quotient(I,ann);
        }
      f = f-1;
    }
  return(I);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
  removeComponent(I,1);
}


/////////////////////////////////////////////////////////////////////
proc AssOfDim(ideal I, int e)
"USAGE:   AssOfDim(I,e); I ideal, e integer
RETURN:  ideal, such that the associated primes are exactly
         the associated primes of I having dimension e
EXAMPLE: example AssOfDim; shows an example"
{
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  int g = nvars(basering) - e;

  //Compute a resolution of I.
  ideal J = std(I);
  if(homog(I)==1)
    {
      list re = sres(J,g+1);
      re = minres(re);
    }
  else
    {
      list re = mres(I,g+1);
    }

  ideal ann = AnnExtEHV(g,re);
  int cod = nvars(basering) - dim(std(ann));

  //If the codimension of I_g:=Ann(Ext^g(R/I,R)) equals g...
  if(cod == g)
    {
      //...then return the equidimensional part of I_g...
      ann = equiMaxEHV(ann);
      attrib(ann,"isEquidimensional",1);
      return(ann);
    }
  //...otherwise...
  else
    {
      //...I has no associated primes of dimension e.
      return(ideal(1));
    }
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
  AssOfDim(I,1);
}

/////////////////////////////////////////////////////////////////////
//                                                                 //
//                      T H E   R A D I C A L                      //
//                                                                 //
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
static proc aJacob(ideal I, int a)
"USAGE:   aJacob(I,a); I ideal, a integer
RETURN:  ideal, the ath-Jacobian ideal of I"
{
  matrix M = jacob(I);
  int n = nvars(basering);
  if(n-a <= 0)
    {
      return(ideal(1));
    }
  if(n-a > nrows(M) or n-a > ncols(M))
    {
      return(ideal(0));
    }
  ideal J = minor(M,n-a);
  return(J);
}


/////////////////////////////////////////////////////////////////////
proc equiRadEHV(ideal I, list #)
"USAGE:   equiRadEHV(I [,Strategy]); I ideal, Strategy list
RETURN:  ideal, the equidimensional radical of I,
         i.e. the intersection of the minimal associated primes of I
         having the same dimension as I
NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos,
         Works only in characteristic 0 or p large.
         The (optional) second argument determines the strategy used:
         Strategy[1]        > strategy for the equidimensional part
                     = 0    : uses equiMaxEHV
                     = 1    : uses equidimMax
         Strategy[2]        > strategy for the radical
                     = 0    : combination of strategy 1 and 2
                     = 1    : computation of the radical just with the
                              help of regular sequences
                     = 2    : does not try to find a regular sequence
         Strategy[3]        > strategy for the computation of ideal quotients
                     = n    : uses quot(.,.,n) for the ideal quotient computations
         If no second argument is given then Strategy=(0,0,0) is used.
EXAMPLE: example equiRadEHV; shows an example"
{
  if(printlevel > 2){"Entering equiRadEHV.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  if((char(basering)<100)&&(char(basering)!=0))
    {
      "WARNING: The characteristic is too small, the result may be wrong";
    }

  //Define the Strategy to be used.
  if(size(#) > 0)
    {
      if(#[1]!=1)
        {
          int equStr = 0;
        }
      else
        {
          int equStr = 1;
        }
      if(size(#) > 1)
        {
          if(#[2]!=1 and #[2]!=2)
            {
              int strategy = 0;
            }
          else
            {
              int strategy = #[2];
            }
          if(size(#) > 2)
            {
              int quoStr = #[3];
            }
          else
            {
              int quoStr = 0;
            }
        }
      else
        {
          int strategy = 0;
          int quoStr = 0;
        }
    }
  else
    {
      int equStr = 0;
      int strategy = 0;
      int quoStr = 0;
    }

  ideal J,I0,radI0,L,radI1,I2,radI2;
  int l,n;
  intvec op = option(get);
  matrix M;

  option(redSB);
  list m = mstd(I);
  option(set,op);

  int d = dim(m[1]);
  if(d==-1)
    {
      return(ideal(1));
    }

  if(strategy != 2)
    {
      /////////////////////////////////////////////
      // Computing the equidimensional radical   //
      // via regular sequenves                   //
      /////////////////////////////////////////////

      if(printlevel > 2){"Trying to find a regular sequence.";}
      int cod = nvars(basering)-d;

      //Complete intersection case:
      if(cod==size(m[2]))
        {
          J = aJacob(m[2],d);
          if(printlevel > 2){"Leaving equiRadEHV.";}
          return(quot1(m[2],J,quoStr));
        }

      //First codim elements of I are a complete intersection:
      for(l=1; l<=cod; l++)
        {
          I0[l] = m[2][l];
        }
      n = dim(groebner(I0))+cod-nvars(basering);

      //Last codim elements of I are a complete intersection:
      if(n!=0)
              {
          for(l=1; l<=cod; l++)
            {
              I0[l] = m[2][size(m[2])-l+1];
            }
          n = dim(groebner(I0))+cod-nvars(basering);
        }

      //Taking a generic linear combination of the input:
      if(n!=0)
              {
          M = transpose(sparsetriag(size(m[2]),cod,95,1));
          I0 = ideal(M*transpose(m[2]));
          n = dim(groebner(I0))+cod-nvars(basering);
        }

      //Taking a more generic linear combination of the input:
      if(n!=0)
        {
          while(strategy == 1 and n!=0)
            {
              M = transpose(sparsetriag(size(m[2]),cod,0,100));
              I0 = ideal(M*transpose(m[2]));
              n = dim(groebner(I0))+cod-nvars(basering);
            }
        }

      if(n==0)
        {
          J = aJacob(I0,d);
          if(printlevel > 2){"1st quotient.";}
          radI0 = quot1(I0,J,quoStr);
          if(printlevel > 2){"2nd quotient.";}
          L = quot1(radI0,m[2],quoStr);
          if(printlevel > 2){"3rd quotient.";}
          radI1 = quot1(radI0,L,quoStr);
          attrib(radI1,"isEquidimensional",1);
          attrib(radI1,"isRadical",1);
          if(printlevel > 2){"Leaving equiRadEHV.";}
          return(radI1);
        }
    }

  ////////////////////////////////////////////////////
  // Computing the equidimensional radical directly //
  ////////////////////////////////////////////////////

  if(printlevel > 2){"Computing the equidimensional radical directly";}

  //Compute the equidimensional part depending on the chosen strategy
  if(equStr == 0)
    {
      I = equiMaxEHV(I);
    }
  if(equStr == 1)
    {
      I = equidimMax(I);
    }
  int a = nvars(basering)-1;

  while(a > d)
    {
      if(printlevel > 2){"While-Loop: "+string(a);}
      J = aJacob(I,a);
      while(dim(groebner(J+I))==d)
        {
          if(printlevel > 2){"Quotient-Computation.";}
          I = quot1(I,J,quoStr);
          if(printlevel > 2){"Computing the a-th Jacobian";}
            J = aJacob(I,a);
        }
      a = a-1;
    }
  if(printlevel > 2){"We left While-Loop.";}
  if(printlevel > 2){"Computing the a-th Jacobian";}
  J = aJacob(I,d);
  if(printlevel > 2){"Quotient-Computation.";}
  I = quot1(I,J,quoStr);
  attrib(I,"isEquidimensional",1);
  attrib(I,"isRadical",1);
  if(printlevel > 2){"Leaving equiRadEHV.";}
  return(I);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  poly  p = z2+1;
  poly  q = z3+2;
  ideal i = p*q^2,y-z2;
  ideal pr= equiRadEHV(i);
  pr;
}

/////////////////////////////////////////////////////////////////////
proc radEHV(ideal I, list #)
"USAGE:   radEHV(I [,Strategy]); ideal I, Strategy list
RETURN:  ideal, the radical of I
NOTE:    uses the algorithm of Eisenbud/Huneke/Vasconcelos
         Works only in characteristic 0 or p large.
         The (optional) second argument determines the strategy used:
         Strategy[1]        > strategy for the equidimensional part
                     = 0    : uses equiMaxEHV
                     = 1    : uses equidimMax
         Strategy[2]        > strategy for the radical
                     = 0    : combination of strategy 1 and 2
                     = 1    : computation of the radical just with the
                              help of regular sequences
                     = 2    : does not try to find a regular sequence
         Strategy[3]        > strategy for the computation of ideal quotients
                     = n    : uses quot(.,.,n) for the ideal quotient computations
         If no second argument is given then Strategy=(0,0,0) is used.
EXAMPLE: example radEHV; shows an example"
{
  if(printlevel > 2){"Entering radEHV.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }

  //Compute the equidimensional radical J of I.
  ideal J = equiRadEHV(I,#);

  //If I is the entire ring...
  if(deg(J[1]) <= 0)
    {
      //...then return the ideal generated by 1...
      return(ideal(1));
    }

  //...else remove the maximal dimensional components and
  //compute the radical K of the lower dimensional components ...
  ideal K = radEHV(sat(I,J)[1],#);

  //..and intersect it with J.
  K = intersect(J,K);
  attrib(K,"isRadical",1);
  return(K);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  poly  p = z2+1;
  poly  q = z3+2;
  ideal i = p*q^2,y-z2;
  ideal pr= radical(i);
  pr;
}

/////////////////////////////////////////////////////////////////////
proc IntAssOfDim1(ideal I, int e)
"USAGE:   IntAssOfDim1(I,e); I idea, e integer
RETURN:  ideal, the intersection of the associated primes of I having dimension e
EXAMPLE: example IntAssOfDim1; shows an example"
{
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  int g = nvars(basering) - e;

  //Compute a resolution of I.
  ideal J = groebner(I);
  if(homog(I)==1)
    {
      list re = sres(J,g+1);
      re = minres(re);
    }
  else
    {
      list re = mres(I,g+1);
    }

  ideal ann = AnnExtEHV(g,re);
  int cod = nvars(basering) - dim(groebner(ann));
  //If the codimension of I_g:=Ann(Ext^g(R/I,I)) equals g...
  if(cod == g)
    {
      //...then return the equidimensional radical of I_g...
      ann = equiRadEHV(ann);
      attrib(ann,"isEquidimensional",1);
      attrib(ann,"isRadical",1);
      return(ann);
    }
  else
    {
      //...else I has no associated primes of dimension e.
      return(ideal(1));
    }
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
  IntAssOfDim1(I,1);
}

/////////////////////////////////////////////////////////////////////
proc IntAssOfDim2(ideal I, int e)
"USAGE:   IntAssOfDim2(I,e); I ideal, e integer
RETURN:  ideal, the intersection of the associated primes of I having dimension e
EXAMPLE: example IntAssOfDim2; shows an example"
{
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  ideal I1 = removeComponent(I,e);
  ideal I2 = removeComponent(I,e+1);
  ideal b = quotient(I1,I2);
  b = equiRadEHV(b);
  return(b);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
  IntAssOfDim2(I,1);
}


/////////////////////////////////////////////////////////////////////
//                                                                 //
//   Z E R O - D I M E N S I O N A L   D E C O M P O S I T I O N   //
//                                                                 //
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
proc decompEHV(ideal I)
"USAGE:   decompEHV(I); I zero-dimensional radical ideal
RETURN:  list, the associated primes of I
EXAMPLE: example decompEHV; shows an example"
{
  if(printlevel > 2){"Entering decompEHV.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }

  int e,m,vd,nfact;
  list l,k,L;
  poly f,h;
  ideal J;
  def base = basering;

  if( attrib(I,"isSB")!=1 )
    {
      I = groebner(I);
    }

  while(1)
    {
      //Choose a random polynomial f from R.
      e = random(0,100);
      f = sparsepoly(e);

      //Check if f lies not in I.
      if(reduce(f,I)!=0)
        {
          J = quotient1(I,f);

          //If f is a zerodivisor modulo I...
          if( isSubset(J,I) == 0 )
            {
              //...then use recursion...
              if(printlevel > 2){"We found a zerodivisor -- recursion";}
              l = decompEHV(J) + decompEHV(I+f);
              return(l);
            }
          if(printlevel > 2){"We found a non-zero-divisor.";}

          //...else compute the vectorspace dimension vd of I and...
          vd = vdim(I);

          //...compute m minimal such that 1,f,f^2,...,f^m are linearly dependent.
          qring Q = I;
          poly g = fetch(base,f);
          k = algDependent(g);
          def R = k[2];
          setring R;
          if(size(ker)!=1)
            {
              //Calculate a generator for ker.
              ker = mstd(ker)[2];
            }
          poly g = ker[1];
          m = deg(g);

          //If m and vd coincide...
          if(m==vd)
            {
              if(printlevel > 2){"We have a good candidate.";}
              //...then factorize g.
              if(printlevel > 2){"Factorizing.";}
              L = factorize(g,2);    //returns non-constant factors and multiplicities
              nfact = size(L[1]);

              //If g is irreducible...
              if(nfact==1 and L[2][1]==1)
                {
                  if(printlevel > 2){"The element is irreducible.";}
                  setring base;
                  //..then I is a maximal ideal...
                  l[1] = I;
                  kill R,Q;
                  return(l);
                }
              //...else...
              else
                {

                  if(printlevel > 2){"The element is not irreducible -- recursion.";}
                  //...take a non-trivial factor g1 of g ...
                  poly g1 = L[1][1];
                  //..and insert f in g1...
                  execute("ring newR = (" + charstr(R) + "),(y(1)),(" + ordstr(R) + ");");
                  poly g2 = imap(R,g1);
                  setring base;
                  h = fetch(newR,g2);
                  h = subst(h,var(1),f);
                  //...and use recursion...
                  kill R,Q,newR;
                  l = l + decompEHV(quotient1(I,h));
                  l = l + decompEHV(I+h);
                  return(l);
                }
            }
          setring base;
          kill R,Q;
          }
    }
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring r = 32003,(x,y),dp;
  ideal i = x2+y2-10,x2+xy+2y2-16;
  decompEHV(i);
}


/////////////////////////////////////////////////////////////////////
//                                                                 //
//           A S S O C I A T E D   P R I M E S                     //
//                                                                 //
/////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////
static proc idempotent(ideal I)
"USAGE:   idempotent(I); ideal I (weighted) homogeneous radical ideal,
         I intersected K[x(1),...,x(k)] zero-dimensional
         where deg(x(i))=0 for all i <= k and deg(x(i))>0 for all i>k.
RETURN:  a list of ideals I(1),...,I(t) such that
         K[x]/I = K[x]/I(1) x ... x K[x]/I(t)"
{
  if(printlevel > 2){"Entering idempotent.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  int n = nvars(basering);
  poly f,g;
  string j;
  int i,k,splits;
  list l;

  //Collect all variables of degree 0.
  for(i=1; i<=n; i++)
    {
      if(deg(var(i)) > 0)
        {
          f = f*var(i);
        }
      else
        {
          splits = 1;
          j = j +  string(var(i)) + "," ;
        }
    }
  //If there are no variables of degree 0
  //then there are no idempotents and we are done...
  if(splits == 0)
    {
      l[1]=I;
      return(l);
    }

  //...else compute J = I intersected K[x(1),...,x(k)]...
  ideal J = eliminate(I,f);
  def base = basering;
  j = j[1,size(j)-1];
  j = "ring @r = (" + charstr(basering) + "),(" + j + "),(" + ordstr(basering) + ");";
  execute(j);
  ideal J = imap(base,J);

  //...and compute the associated primes of the zero-dimensional ideal J.
  list L = decompEHV(J);
  int s = size(L);
  ideal K,Z;
  poly g;
  //For each associated prime ideal P_i of J...
  for(i=1; i<=s; i++)
    {
      K = ideal(1);
      //...comnpute the intersection K of the other associated prime ideals...
      for(k=1; k<=s; k++)
        {
          if(i!=k)
            {
              K = intersect(K,L[k]);
            }
        }

      //...and find an element that lies in K but not in P_i...
      g = randomid(K,1)[1];
      Z = L[i];
      Z = std(Z);
      while(reduce(g,Z)==0)
        {
          g = randomid(K,1)[1];
        }
      setring base;
      g = imap(@r,g);
      //...and compute the corresponding ideal I(i)
      l[i] = quotient(I,g);
      setring @r;

    }
  setring base;
  return(l);
}


/////////////////////////////////////////////////////////////////////
static proc equiAssEHV(ideal I)
"USAGE:   equiAssEHV(I); I equidimensional, radical, and homogeneous ideal
RETURN:  a list, the associated prime ideals of I"
{
  if(printlevel > 2){"Entering equiAssEHV.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  list L;
  def base = basering;
  int n = nvars(basering);
  int i,j;

  //Compute the normalization of I.
  if(printlevel > 2){"Entering Normalization.";}
  list norOut = normal(I, "noDeco");
  list K = norOut[1];
  if(printlevel > 2){"Leaving Normalisation.";}

  //The normalization algorithm returns k factors.
  int k = size(K);
  if(printlevel > 1){"Normalization algorithm splits ideal in " + string(k) + " factors.";}

  //Examine each factor of the normalization.
  def P;
  for(i=1; i<=k; i++)
    {
      P = K[i];
      setring P;

      //Use procedure idempotent to split the i-th factor of
      //the normalization in a product of integral domains.
      if(printlevel > 1){"Examining " + string(i) +". factor.";}
      list l = idempotent(norid);
      int leng = size(l);
      if(printlevel > 1){"Idempotent algorithm splits factor " + string(i) + " in " + string(leng) + " factors.";}

      //Intersect the minimal primes corresponding to the
      //integral domains obtained from idempotent with the groundring,
      //i.e. compute the preimages w.r.t. the corresponding normalization map
      ideal J;
      for(j=1; j<=leng; j++)
        {
          J = l[j];
          setring base;
          L[size(L)+1] = preimage(P,normap,J);
          setring P;
        }
      kill l, leng, J;
      setring base;
    }

  return(L);
}


/////////////////////////////////////////////////////////////////////
proc AssEHV(ideal I, list #)
"USAGE:  AssEHV(I [,Strategy]); I Ideal, Strategy list
RETURN:  a list, the associated prime ideals of I
NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos.
         The (optional) second argument determines the strategy used:
         Strategy[1]        > strategy for the equidimensional part
                     = 0    : uses equiMaxEHV
                     = 1    : uses equidimMax
         Strategy[2]        > strategy for the equidimensional radical
                     = 0    : uses equiRadEHV
                     = 1    : uses equiRadical
         Strategy[3]        > strategy for equiRadEHV
                     = 0    : combination of strategy 1 and 2
                     = 1    : computation of the radical just with the
                              help of regular sequences
                     = 2    : does not try to find a regular sequence
         Strategy[4]        > strategy for the computation of ideal quotients
                     = n    : uses quot(.,.,n) for the ideal quotient computations
         If no second argument is given, Strategy=(0,0,0,0) is used.
EXAMPLE: example AssEHV; shows an example"
{
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  if(printlevel > 2){"Entering AssEHV";}

  //Specify the strategy to be used.
  if(size(#)==0)
    {
      # = 0,0,0,0;
    }
  if(size(#)==1)
    {
      # = #[1],0,0,0;
    }
  if(size(#)==2)
    {
      # = #[1],#[2],0,0;
    }
  if(size(#)==3)
    {
      # = #[1],#[2],#[3],0;
    }

  list L;
  ideal K;
  def base = basering;
  int m,j;
  ideal J = groebner(I);
  int n = nvars(basering);
  int d = dim(J);

  //Compute a resolution of I.
  if(printlevel > 2){"Computing resolution.";}
  if(homog(I)==1)
    {
      list re = sres(J,0);
      re = minres(re);
    }
  else
    {
      list re = mres(I,0);
    }
  ideal ann;
  int cod;

  //For 0<=i<= dim(I) compute the intersection of the i-dimensional associated primes of I
  for(int i=0; i<=d; i++)
    {
      if(printlevel > 1){"Are there components of dimension " + string(i) + "?";}
      ann = AnnExtEHV(n-i,re);
      cod = n - dim(groebner(ann));

      //If there are associated primes of dimension i...
      if(cod == n-i)
        {
          if(printlevel > 1){"Yes. There are components of dimension " + string(i) + ".";}
          //...then compute the intersection K of all associated primes of I of dimension i
          if(#[2]==0)
            {
              if(size(#) > 3)
                {
                  K = equiRadEHV(ann,#[1],#[3],#[4]);
                }
              if(size(#) > 2)
                {
                  K = equiRadEHV(ann,#[1],#[3]);
                }
              else
                {
                  K = equiRadEHV(ann,#[1]);
                }
            }
          if(#[2]==1)
            {
              K = equiRadical(ann);
            }
          attrib(K,"isEquidimensional",1);
          attrib(K,"isRadical",1);

          //If K is already homogeneous then use equiAssEHV to recover the associated primes of K...
           if(homog(K)==1)
            {
              if(printlevel > 2){"Input already homogeneous.";}
              L = L + equiAssEHV(K);
            }

          //...else...
          else
            {
              //...homogenize K w.r.t. t,...
              if(printlevel > 2){"Input not homogeneous; must homogenize.";}
              changeord("homoR","dp");
              ideal homoJ = fetch(base,K);
              homoJ = groebner(homoJ);
              execute("ring newR = (" + charstr(base) + "),(x(1..n),t),dp;");
              ideal homoK = fetch(homoR,homoJ);
              homoK = homog(homoK,t);
              attrib(homoK,"isEquidimensional",1);
              attrib(homoK,"isRadical",1);

              //...compute the associated primes of the homogenization using equiAssEHV,...
              list l = equiAssEHV(homoK);

              //...and set t=1 in the generators of the associated primes just computed.
              ideal Z;
              for(j=1; j<=size(l); j++)
                {
                  Z = subst(l[j],t,1);
                  setring base;
                  L[size(L)+1] = fetch(newR,Z);
                  setring newR;
                }
              setring base;
              kill homoR;
              kill newR;
            }
        }
      else
        {
          if(printlevel > 1){"No. There are no components of dimension " + string(i) + ".";}
        }
    }
  return(L);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  poly  p = z2+1;
  poly  q = z3+2;
  ideal i = p*q^2,y-z2;
  list pr = AssEHV(i);
  pr;
}

/////////////////////////////////////////////////////////////////////
proc minAssEHV(ideal I, list #)
"USAGE:  minAssEHV(I [,Strategy]); I ideal, Strategy list
RETURN:  a list, the minimal associated prime ideals of I
NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos.
         The (optional) second argument determines the strategy used:
         Strategy[1]        > strategy for the equidimensional part
                     = 0    : uses equiMaxEHV
                     = 1    : uses equidimMax
         Strategy[2]        > strategy for the equidimensional radical
                     = 0    : uses equiRadEHV, resp. radicalEHV
                     = 1    : uses equiRadical, resp. radical
         Strategy[3]        > strategy for equiRadEHV
                     = 0    : combination of strategy 1 and 2
                     = 1    : computation of the radical just with the
                              help of regular sequences
                     = 2    : does not try to find a regular sequence
         Strategy[4]        > strategy for the computation of ideal quotients
                     = n    : uses quot(.,.,n) for the ideal quotient computations
         If no second argument is given, Strategy=(0,0,0,0) is used.
EXAMPLE: example minAssEHV; shows an example"
{
  if(ord_test(basering)!=1)
    {ERROR("// Not implemented for this ordering, please change to global ordering.");}

  //Specify the strategy to be used.
  if(size(#)==0)
    {
      # = 0,0,0,0;
    }
  if(size(#)==1)
    {
      # = #[1],0,0,0;
    }
  if(size(#)==2)
    {
      # = #[1],#[2],0,0;
    }
  if(size(#)==3)
    {
      # = #[1],#[2],#[3],0;
    }

  //Compute the radical of I.
  if(#[2]==0)
    {
      I = radEHV(I,#);
    }
  if(#[2]==1)
    {
      I = radical(I);
    }

  //Compute the associated primes of the radical.
  return(AssEHV(I,#));
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  poly  p = z2+1;
  poly  q = z3+2;
  ideal i = p*q^2,y-z2;
  list pr = minAssEHV(i);
  pr;
}


/////////////////////////////////////////////////////////////////////
//                                                                 //
//        P R I M A R Y   D E C O M P O S I T I O N                //
//                                                                 //
/////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////
proc localize(ideal I, ideal P, list l)
"USAGE:   localize(I,P,l); I ideal, P an associated prime ideal of I,
         l list of all associated primes of I
RETURN:  ideal,  the contraction of the ideal generated by I
         in the localization w.r.t P
EXAMPLE: example localize; shows an example"
{
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }

  ideal Intersection = ideal(1);
  int s = size(l);
  if(attrib(P,"isSB")!=1)
    {
      P = groebner(P);
    }

  //Compute the intersection of all associated primes of I that are not contained in P.
  for(int i=1; i<=s; i++)
    {
      if(isSubset(l[i],P)!=1)
        {
          Intersection = intersect(Intersection,l[i]);
        }
    }
  Intersection = groebner(Intersection);

  //If the intersection is the entire ring...
  if(reduce(1,Intersection)==0)
    {
      //...then return I...
      return(I);
    }
  //...else try to find an element f that lies in the intersection but outside P...
  poly f = 0;
  while(reduce(f,P) == 0)
    {
      f = randomid(Intersection,1)[1];
    }

  //...and saturate I w.r.t. f.
  I = sat(I,f)[1];
  return(I);
}

/////////////////////////////////////////////////////////////////////
proc componentEHV(ideal I, ideal P, list L, list #)
"USAGE:   componentEHV(I,P,L [,Strategy]);
         I ideal, P associated prime of I,
         L list of all associated primes of I, Strategy list
RETURN:  ideal, a P-primary component Q for I
NOTE:    The (optional) second argument determines the strategy used:
         Strategy[1]        > strategy for equidimensional part
                     = 0    : uses equiMaxEHV
                     = 1    : uses equidimMax
         If no second argument is given then Strategy=0 is used.
EXAMPLE: example componentEHV; shows an example"
{
  if(printlevel > 2){"Entering componentEHV.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }

  //If no strategy is specified use standard strategy.
  if(size(#)==0)
    {
      # = 0;
    }

  ideal T = P;
  ideal Q;

  //Compute the localization of I at P...
  ideal IP = groebner(localize(I,P,L));

  //...and compute the saturation of the localization w.r.t. P.
  ideal IP2 = sat(IP,P)[1];

  //As long as we have not found a primary component...
  int isPrimaryComponent = 0;
  while(isPrimaryComponent!=1)
    {
      //...compute the equidimensional part Q of I+P^n...
      if(#[1]==0)
        {
          Q = equiMaxEHV(I+T);
        }
      if(#[1]==1)
        {
          Q = equidimMax(I+T);
        }
      //...and check if it is a primary component for P.
      if(isSubset(intersect(IP2,Q),IP)==1)
        {
          isPrimaryComponent = 1;
        }
      else
        {
          T = T*P;
        }
    }
  if(printlevel > 2){"Leaving componentEHV.";}
  return(Q);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  poly  p = z2+1;
  poly  q = z3+2;
  ideal i = p*q^2,y-z2;
  list pr = AssEHV(i);
  componentEHV(i,pr[1],pr);
}


/////////////////////////////////////////////////////////////////////
proc primdecEHV(ideal I, list #)
"USAGE:   primdecEHV(I [,Strategy]); I ideal, Strategy list
RETURN:  a list pr of primary ideals and their associated primes:
                pr[i][1]   the i-th primary component,
                pr[i][2]   the i-th prime component.
NOTE:    Algorithm of Eisenbud/Huneke/Vasconcelos.
         The (optional) second argument determines the strategy used:
         Strategy[1]        > strategy for equidimensional part
                     = 0    : uses equiMaxEHV
                     = 1    : uses equidimMax
         Strategy[2]        > strategy for equidimensional radical
                     = 0    : uses equiRadEHV, resp. radicalEHV
                     = 1    : uses equiRadical, resp. radical
         Strategy[3]        > strategy for equiRadEHV
                     = 0    : combination of strategy 1 and 2
                     = 1    : computation of the radical just with the
                              help of regular sequences
                     = 2    : does not try to find a regular sequence
         Strategy[4]        > strategy for the computation of ideal quotients
                     = n    : uses quot(.,.,n) for the ideal quotient computations
         If no second argument is given then Strategy=(0,0,0,0) is used.
EXAMPLE: example primdecEHV; shows an example"
{
  if(printlevel > 2){"Entering primdecEHV.";}
  if(ord_test(basering)!=1)
    {
      ERROR("// Not implemented for this ordering, please change to global ordering.");
    }
  list L,K;

  //Specify the strategy to be used.
  if(size(#)==0)
    {
      # = 0,0,0,0;
    }
  if(size(#)==1)
    {
      # = #[1],0,0,0;
    }
  if(size(#)==2)
    {
      # = #[1],#[2],0,0;
    }
  if(size(#)==3)
    {
      # = #[1],#[2],#[3],0;
    }

  //Compute the associated primes of I...
  L = AssEHV(I,#);
  if(printlevel > 0){"We have " + string(size(L)) + " prime components.";}

  //...and compute for each associated prime of I a corresponding primary component.
  int l = size(L);
  for(int i=1; i<=l; i++)
    {
      K[i] = list();
      K[i][2] = L[i];
      K[i][1] = componentEHV(I,L[i],L,#[1]);
    }
  if(printlevel > 2){"Leaving primdecEHV.";}
  return(K);

}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y,z),dp;
  poly  p = z2+1;
  poly  q = z3+2;
  ideal i = p*q^2,y-z2;
  list pr = primdecEHV(i);
  pr;
}


/////////////////////////////////////////////////////////////////////
//                                                                 //
//        A L G O R I T H M S   F O R   T E S T I N G              //
//                                                                 //
/////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////
proc compareLists(list L, list K)
"USAGE:   checkLists(L,K); L,K list of ideals
RETURN:  integer, 1 if the lists are the same up to ordering and 0 otherwise
EXAMPLE: example checkLists; shows an example"
{
  int s1 = size(L);
  int s2 = size(K);
  if(s1!=s2)
    {
      return(0);
    }
  list L1, K1;
  int i,j,t;
  list N;
  for(i=1; i<=s1; i++)
    {
      L1[i]=std(L[i][2]);
      K1[i]=std(K[i][2]);
    }
  for(i=1; i<=s1; i++)
    {
      for(j=1; j<=s1; j++)
        {
          if(isSubset(L1[i],K1[j]))
            {
              if(isSubset(K1[j],L1[i]))
                {
                  for(t=1; t<=size(N); t++)
                    {
                      if(N[t]==j)
                        {
                          return(0);
                        }
                    }
                  N[size(N)+1]=j;
                }

            }
        }
    }
  return(1);
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring  r = 0,(x,y),dp;
  ideal i = x2,xy;
  list L1 = primdecGTZ(i);
  list L2 = primdecEHV(i);
  compareLists(L1,L2);
}