source: git/Singular/LIB/monomial.lib @ 0dd77c2

// last modified:  10.06.2010
//////////////////////////////////////////////////////////////////////////////
version = "$Id: monomial.lib$";
category = "Commutative algebra";
info = "
LIBRARY: monomial.lib   Primary and irreducible decompositions of monomial
                        ideals
AUTHORS: I.Bermejo,           ibermejo@ull.es
@*       E.Garcia-Llorente,   evgarcia@ull.es
@*       Ph.Gimenez,          pgimenez@agt.uva.es
OVERVIEW:
 A library for computing a primary and the irreducible decompositions of a
 monomial ideal using several methods.
 In this library we also take advantage of the fact that the ideal is
 monomial to make some computations that are Grobner free in this case
 (radical, intersection, quotient...).
PROCEDURES:
 isMonomial(id);       checks whether an ideal id is monomial
 minbaseMon(id);       computes the minimal monomial generating set of a
                       monomial ideal id
 gcdMon(f,g);          computes the gcd of two monomials f, g
 lcmMon(f,g);          computes the lcm of two monomials f, g
 membershipMon(f,id);  checks whether a polynomial f belongs to a monomial
                       ideal id
 intersectMon(id1,id2);intersection of monomial ideals id1 and id2
 quotientMon(id1,id2); quotient ideal id1:id2
 radicalMon(id);       computes the radical of a monomial ideal id
 isprimeMon(id);       checks whether a monomial ideal id is prime
 isprimaryMon(id);     checks whether a monomial ideal id is primary
 isirreducibleMon(id); checks whether a monomial ideal id is irreducible
 isartinianMon(id);    checks whether a monomial ideal id is artininan
 isgenericMon(id);     checks whether a monomial ideal id is generic
 dimMon(id);           dimension of a monomial ideal id
 irreddecMon(id,..);   computes the irreducible decomposition of a monomial
                       ideal id
 primdecMon(id,..);    computes a minimal primary decomposition of a monomial
                       ideal id
";
LIB "poly.lib";  // Para "maxdeg1" en "isprimeMon"
//---------------------------------------------------------------------------        
//-----------------------   INTERNOS    -------------------------------------
//---------------------------------------------------------------------------        
/////////////////////////////////////////////////////////////////////////////
//
static proc checkIdeal (ideal I)
"
USAGE:    checkIdeal (I); I ideal.
RETURN:   1, if ideal is generated by monomials; 0, otherwise.
"
// Aqui NO estoy quitando el caso de que el ideal sea el trivial.
{
  int i,n;
  n = ncols(I);
  for (i = 1 ; i <= n ; i ++)
    {
      if ( size(I[i]) > 1 )
        {
          return (0);
        }
    }
  return (1);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc quotientIdealMon (ideal I,poly f)
"
USAGE:    quotientIdealMon(I,f); I ideal, f polynomial.
RETURN:   an ideal, the quotient ideal I:(f).
ASSUME:   I is an ideal generated by a list of monomials and f is a monomial
          of the basering.
"
{
  // Variables
  int i,j;
  poly g,generator;
  intvec v;
  ideal J;
  J = 0;

  int sizI = size(I);
  for (i = 1 ; i <= sizI ; i++)
    {
      g = gcdMon(I[i],f);
      // Cociente de dos monomios: restamos los exponentes, y en el
      // denominador va el mcd
      v = leadexp(I[i]) - leadexp(g);
      generator = monomial (v);
      if (membershipMon(generator,J) == 0)
        {
          J=J,generator;
        }
    }
  // minimal monomial basis
  return ( minbase(J) );
}
/////////////////////////////////////////////////////////////////////////////
//
static proc soporte (poly f)
"
USAGE:    soporte(f); f polynomial.
RETURN:   0, if the monomial f is product of more than one variable;
          otherwise, an integer j, 1<=j<=n, if the monomial f is a power of
          x(j).
ASSUME:   f is a monomial of the basering K[x(1)..x(n)].
"
{
  def R = basering;
  // Variables
  int nvar,i,cont,sop;
  intvec expf;
  nvar = nvars(R);
  expf = leadexp(f);
  cont = 0;
  // cont va a contar el numero de componentes del vector no nulas.
  // En sop guardamos el subindice de la componente no nula.
  for (i = 1 ; i <= nvar ; i++)
    {
      if (expf[i] > 0)
        {
          cont ++;
          sop = i;
          // Si cont > 1 ==> aparece mas de una variable, devolvemos 0
          if (cont > 1)
            {
              return (0);
            }
        }
    }
  return(sop);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc irredAux (ideal I)
"
USAGE:    irredAux (I); I ideal.
RETURN:   1, if I is irreducible; otherwise, an intvec whose fist entry is
          the position of a generator which is the product of more than one
          variable, the next entries are the indexes of those variables.
ASSUME:   I is a monomial ideal of the basering K[x(1)..x(n)] and it is
          generated by its minimal monomial generators.
NOTE:     This procedure is a modification of isirreducibleMon to give
          more information when ideal is not irreducible.
"
{
  def R = basering;
  // Variables
  int sizI,i,nvar,j,sum;
  intvec w,exp;
  sizI = size(I);
  nvar = nvars(R);
  for (i = 1 ; i <= sizI ; i++)
    {
      sum = 0;
      exp = leadexp(I[i]);
      for (j = 1 ; j <= nvar ; j++)
        {
          // Al menos tenemos una variable en cada generador, luego
          // entramos minimo 1 vez, luego sum >= 1.
          if (exp[j] <> 0)
            {
              sum = sum + 1;
              w[sum] = j;
            }
        }
      // Si hay mas de una variable la suma sera mayor que 1; y ya
      // sabemos que I no es irreducible.
      if (sum <> 1)
        {
          return(i,w);
        }
    }
  return(1);
}
//////////////////////////////////////////////////////////////////////
//
static proc contents (ideal I,ideal J)
"
USAGE:    contents (I,J); I,J ideals.
RETURN:   1, if I is contained in J; 0, otherwise.
ASSUME:   I,J are monomial ideals of the basering.
"
{
  // Variables
  poly f;
  int i,resp;
  int n = size(I);
  // Desde que haya un generador que no pertenzca al ideal, ya no se da
  // el contenido y terminamos.
  for (i = 1 ; i <= n ; i++)
    {
      resp = membershipMon(I[i],J);
      if (resp == 0)
        {
          return(0);
        }
    }
  return(1);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc equal (ideal I,ideal J)
"
USAGE:    equal (I,J); I,J ideals.
RETURN:   1, if I and J are the same ideal; 0, otherwise.
ASSUME:   I,J are monomial ideals of the basering and are defined by their
          minimal monomial generators.
"
{
  // Variables
  int n,i,j;
  intvec resps;
  // Si no tienen el mismo numero de generadores, no pueden ser iguales; ya
  // que vienen dados por el sistema minimal de generadores.
  if (size(I) <> size(J))
    {
      return(0);
    }
  n = size(I);
  // Como ambos ideales vienen dados por la base minimal, no vamos a
  // tener problemas con que comparemos uno de I con otro de J, pues
  // no puede haber generadores iguales en el mismo ideal.
  // Si los ordenamos, se puede comparar uno a uno
  I = sort(I)[1];
  J = sort(J)[1];
  for (i = 1 ; i <= n ; i++)
    {
      if (I[i] <> J[i])
        {
          return(0);
        }
    }
  return(1);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc radicalAux (ideal I)
"
USAGE:    radicalAux (I); I ideal.
RETURN:   an ideal, the radical ideal of I
ASSUME:   I is an irreducible monomial ideal of the basering given by its
          minimal monomial generators.
"
{
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Variables
  int i,cont;
  intvec exp;
  ideal rad;
  ideal I = fetch(R,I);
  // Como en cada generador aparece solo una variable,  y ademas la
  // la misma variable no va a aparecer dos veces, es suficiente
  // con sumar los exponentes de todos los generadores para saber que
  // variables aparecen.
  int n = size(I);
  for (i = 1 ; i <= n ; i++)
    {
      exp = exp + leadexp (I[i]);
    }
  cont = 1;
  for (i = 1 ; i <= nvar ; i++)
    {
      if (exp[i] <> 0)
        {
          rad[cont] = x(i);
          cont ++;
        }
    }
  setring R;
  ideal rad = fetch(R1,rad);
  kill R1;
  return (rad);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc primAux (ideal I)
"
USAGE:    primAux (I); I ideal.
RETURN:   1, if I is primary; otherwise, an intvec, whose first element is
          0, the second is the index of one variable such that a power of it
          does not appear as a generator of I, the rest of the elements are
          the situation in the ideal of that elements of I which
          are product of more than one variable.
ASSUME:   I is a monomial ideal of the basering K[x(1)..x(n)].
NOTE:     This procedure detects if the ideal is primary, when the
          ideal is not primary, it gives some additional information.
"
{
  def R = basering;
  // Variables
  int control,nvar,i,sub_in,l,j;
  intvec v,w,exp_gen;
  // El ideal ya entra generado por el sistema minimal
  nvar = nvars(R);
  int sizI = size(I);
  v[nvar] = 0;
  int cont = 1;
  // v = 1 en la posicion en el ideal de variables sueltas, que son
  // las que no hay que tocar, 0 en las demas. w = posiciones de los
  // generadores de I que hay que comprobar.
  for (i = 1 ; i <= sizI ; i++ )
    {
      sub_in = soporte(I[i]);
      if ( sub_in <> 0)
        {
          v[sub_in] = 1;
        }
      else
        {
          w[cont] = i;
          cont ++;
        }
    }
  l = size(w);
  // No hay ningun generador que tenga productos de variables, luego
  //  este ideal ya es primario.
  if (l == 1 && w[1] == 0)
    {
      return (1);
    }
  for (i = 1 ; i <= l ; i++)
    {
      exp_gen = leadexp(I[w[i]]);
      // Ahora hay que ver que valor tiene el exponente de los
      // generadores oportunos en la posicion que es cero dentro del
      // vector v.
      for (j = 1 ; j <= nvar ; j++)
        {
          if (v[j] == 0)
            {
              if (exp_gen[j] <> 0)
                {
                  return (0,j,w);
                }
            }
        }
    }
  // Si hemos llegado hasta aqui hemos recorrido todo el ideal y por tanto
  // es primario.
  return (1);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc maxExp (ideal I,intvec v)
"
USAGE:    maxExp (I,v); I ideal, v integer vector.
RETURN:   an integer, the greatest power of a variable in the minimal
          monomial set of generators of I.
ASSUME:   I is a monomial ideal of the basering, v=primAux(I) and the
          variable considered is v[2].
          If the ideal I is primary, it returns 0.
NOTE:     The elements of the vector shows what variable and what
          generators we must consider to look for the greatest power
          of this variable.
"
{
  // Variables
  int n,i,max;
  intvec exp;
  // Ponemos el tama?o de v menos 2 porque en el vector v a partir de
  // la tercera componente es donde tenemos la posicion de los
  // generadores que tenemos que estudiar.
  n = size(v)-2;
  // Buscamos el maximo de la variable que no aparece "sola" en los
  // generadores del ideal (donde nos indica v).
  max = 0;
  for (i = 1 ; i <= n ; i++)
    {
      exp = leadexp (I[v[i+2]]);
      if (exp[v[2]] > max)
        {
          max = exp[v[2]];
        }
    }
  return (max);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc irredundant (list l)
"
USAGE:    irredundant (l); l, list.
RETURN:   a list such that the intersection of the elements in that list has
          no redundant component.
ASSUME:   elements of l are monomial ideals of the basering.
"
{
  // Variables
  int i,j,resp;
  ideal J;
  // Recalculamos el tamano de l cuando modificamos la lista (sizl)
  int sizl = size(l);
  for (i = 1 ; i <= sizl ; i++)
    {
      J = 1;
      for (j = 1 ; j <= sizl ; j++)
        {
          // Hacemos la interseccion de todos los ideales menos uno y
          // luego se estudia el contenido.
          if (j <> i)
            {
              J = intersect (J,l[j]);
            }
        }
      J = minbase(J);
      resp = contents(J,l[i]);
      if (resp == 1)
        {
          l = delete (l,i);
          i--;
          sizl = size(l);
        }
    }
  return (l);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc alexDif (intvec v,ideal I)
"
USAGE:    alexDif (v,I); v, intvec; I, ideal.
RETURN:   a list, irreducible monomial ideals whose intersection is the
          Alexander dual of I with respect to v.
ASSUME:   I is a monomial ideal of the basering K[x(1),...,x(n)] given by
          its minimal monomial generators and v is an integer vector with
          n entries s.t.monomial(v) is a multiple of all minimal monomial
          generators of I.
"
{
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Variables
  int i,j;
  intvec exp_I,exp;
  list l;
  ideal J;
  ideal I = fetch (R,I);
  int sizI = size(I);
  // Vamos a tener tantas componentes como generadores minimales tiene el
  // ideal, por eso el bucle es de 1 a size(I).
 for (i = 1 ; i <= sizI ; i++)
   {
     J = 0;
     exp_I = leadexp (I[i]);
     for (j = 1 ; j <= nvar ; j++)
       {
         if (exp_I[j] <> 0)
           {
             exp[j] = v[j] + 1 - exp_I[j];
             J = J, x(j)^exp[j];
           }
       }
     // Tenemos siempre un cero por la inicializacion de J, entonces
     // lo quitamos.
     J = simplify (J,2);
     l = insert (l,J);
   }
 setring R;
 list l = fetch (R1,l);
 kill R1;
 return (l);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc irredPrimary (list l1)
"
USAGE:    irredPrimary (l1); l1, list of ideals.
RETURN:   a list, primary monomial ideals whose intersection is an
          irredundant primary decomposition.
ASSUME:   list l1 is the list of the irredundant irreducible components of a
          monomial ideal I of the basering.
"
{
  // Variables
  int i,sizl1,sizl2,j;
  ideal J,K;
  list l2,l3;
//----- irredundant primary decomposition
  sizl1 = size(l1);
  for (i = 1 ; i <= sizl1 ; i++)
    {
      l2[i] = radicalAux (l1[i]);
    }
  sizl2 = size(l2);
  int sizl2i, sizl2j;
  // Looking for irreducible components whose radicals are equal.
  // l1 = irreducible components list
  // l2 = radical of irreducible components list
  // l3 = primary components list
  for (i = 1 ; i <= sizl1 ; i++)
    {
      J = l2[i];
      sizl2i = size(l2[i]);
      K = l1[i];
      for (j = i+1 ; j <= sizl2 ; j++)
        {
          sizl2j = size(l2[j]);
          if (sizl2i == sizl2j)
            {
              if (equal (J,l2[j]) == 1)
                {
                  K = minbase(intersect (K,l1[j]));
                  l1 = delete (l1,j);
                  sizl1 = size(l1);
                  l2 = delete (l2,j);
                  sizl2 = size(l2);
                  j--;
                }
            }
        }
      l3 = insert (l3,K);
    }
  return (l3);
}
/////////////////////////////////////////////////////////////////////////////
//
static proc isMinimal (ideal I)        
"
USAGE:    isMinimal (I); I ideal.
RETURN:   1, if the generators of I are the minimal ones;
          0 & minimal generators of I, otherwise.
ASSUME:   I is an ideal of the basering generated by monomials.
"
{
  // VARIABLES
  int i;
  ideal J;
 // Quitamos los ceros del sistema de generadores.
  I = simplify(I,2);
  int resp = 1;
  int sizI = size(I);
  // Cambiamos el tamano de I cuando eliminamos generadores
  for (i = 1 ; i <= sizI ; i++)
    {
      if (sizI <> 1)
        {
          if (i == 1)
            {
              J = I[2..sizI];
            }
          else
            {
              if (i > 1 && i < sizI)
                {
                  J = I[1..i-1], I[i+1..sizI];
                }
              else
                {
                  J = I[1..sizI-1];
                }
            }
          // Si quitamos el generador del lugar "i", luego el que
          // ocupa ahora el lugar "i" es el "i+1", de ahi que restemos
          // 1 al i para volver al for de manera que recorramos los
          // generadores como debemos.                
          if (membershipMon(I[i],J) == 1)
            {
              resp = 0;
              I = J;
              i--;
              sizI = size(I);
            }
        }
    }
  if (resp == 1)
    {
      return (1);
    }
  else
    {
      return (0,I);
    }
}
/////////////////////////////////////////////////////////////////////////////
//        
static proc isMonomialGB (ideal I)        
"
USAGE:    isMonomialGB (I); I ideal.
RETURN:   a list, 1 & the minimal generators of I, if I is a monomial ideal;
          0, otherwise.
ASSUME:   I is an ideal of the basering which is not generated by
          monomials.
NOTE:     this procedure is NOT Grobner free and should be used only if the
          ideal has non-monomial generators (use first checkIdeal)
"
{
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
        
  // Variables
  int resp,m,k;
  // Queremos la base de Grobner reducida, para uncidad.
  ideal I = fetch(R,I);
  option(redSB);        
  // Si el ideal es cero, no es monomial.
  if ( size(I) == 0)
    {
      return(0);
    }
  // Base de Grobner
  I = std(I);
  // Una vez que tenemos la GB, no es mas que comprobar que el ideal
  // esta generado por monomios.
  resp = checkIdeal(I);
  if (resp == 0)
    {
      setring R;
      kill R1;
      return (0);
    }
  else
    {
      setring R;
      I = fetch (R1,I);
      kill R1;
      return (1,I);
    }
}
/////////////////////////////////////////////////////////////////////////////
//
// Comparing irreducible decompsitions
// WARNING: this is not a test, when the answer is 1 and the decompositions
//          may not coincide but it is fast and easy and when the answer is
//          0 the decomposition do not coincide.
//
proc areEqual(list l1,list l2)
{
  def R = basering;
  l1 = fetch(R,l1);
  l2 = fetch(R,l2);
  int i,j,sizIdeal;
  poly generator;
  ideal l1Ideal,l2Ideal;
  int sizl1 = size(l1);
  for (i = 1 ; i <= sizl1 ; i ++)
    {
      sizIdeal = size(l1[i]);
      generator = 1;
      for (j = 1 ; j <= sizIdeal ; j ++)
        {
          generator = generator*l1[i][j];
        }
      l1Ideal[i] = generator;
    }
  int sizl2 = size(l2);
  for (i = 1 ; i <= sizl2 ; i ++)
    {
      sizIdeal = size(l2[i]);
      generator = 1;
      for (j = 1 ; j <= sizIdeal ; j ++)
        {
          generator = generator*l2[i][j];
        }
      l2Ideal[i] = generator;
    }
  return (equal(l1Ideal,l2Ideal));
}
/////////////////////////////////////////////////////////////////////////////
//-------------------------------------------------------------------------//
//-----------------------   EXTERNOS   ------------------------------------//
//-------------------------------------------------------------------------//
/////////////////////////////////////////////////////////////////////////////
//        
proc isMonomial (ideal I)        
"USAGE:    isMonomial (I); I ideal.
RETURN:   1, if I is monomial ideal; 0, otherwise.
ASSUME:   I is an ideal of the basering.
EXAMPLE:  example isMonomial; shows some examples.
"
{
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);        
  // Variables
  int resp,m,k;
  // Queremos la base de Grobner reducida, para uncidad.
  ideal I = fetch(R,I);
  option(redSB);        
  // Si el ideal es cero, no es monomial.
  if ( size(I) == 0)
    {
      return(0);
    }
  // Si ya viene dado por sistema de generadores monomiales, devolvemos 1
  if (checkIdeal (I) == 1)
    {
      return(1);
    }
  // Hallamos GB
  I = std(I);
  // Una vez que tenemos la GB, no es mas que comprobar si el ideal
  // esta generado por monomios.
  resp = checkIdeal(I);
  // Volvemos a dejar el comando "std" devolviendo una GB estandar.
  option(noredSB);
  setring R;
  kill R1;
  return(resp);
}
example
{ "EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^3*x*y, w^2*x^2*y^2*z^2 - y^3*z+x^4*z^4*t^4, w*x*y*z*t - w*x^6*y^5*z^4, x^2*z^4*t^3 , w^6*y^4*z^2 + x^2*y^2*z^2;
 isMonomial(I);
 ideal J = w^3*x*y + x^3*y^9*t^3, w^2*x^2*y^2*z^2 - y^3*z, w*x*y*z*t - w*x^6*y^5*z^4, x^2*z^4*t^3 + y^4*z^4*t^4, w^6*y^4*z^2 + x^2*y^2*z^2;
 isMonomial(J);
}
/////////////////////////////////////////////////////////////////////////////
//
proc minbaseMon (ideal I)
"USAGE:    minbaseMon (I); I ideal.
RETURN:   an ideal, the minimal monomial generators of I.
          (-1 if the generators of I are not monomials)
ASSUME:   I is an  ideal generated by a list of monomials of the basering.
EXAMPLE:  example minbaseMon; shows an example.
"
{
  // VARIABLES
  int i;
  ideal J;
  // Si no esta generado por monomios este metodo no vale
  int control = checkIdeal(I);
  if (control == 0)
    {
      return (-1);
    }
  // Quitamos los ceros del sistema de generadores.
  I = simplify(I,2);
  int sizI = size(I);
  for (i = 1 ; i <= sizI ; i++)
    {
      if (sizI > 1)
        {
          if (i == 1)
            {
              J = I[2..sizI];
            }
          else
            {
              if (i > 1 && i < sizI)
                {
                  J = I[1..i-1], I[i+1..sizI];
                }
              else
                {
                  if (i == sizI)
                    {
                      J = I[1..sizI-1];
                    }
                }
            }
          // Si quitamos el generador del lugar "i", luego el que
          // ocupa ahora el lugar "i" es el "i+1", de ahi que restemos
          // 1 al i para volver al for de manera que recorramos los
          // generadores como debemos.        
          if (membershipMon(I[i],J) == 1)
            {
              I = J;
              i--;
              sizI = size(I);
            }
        }
    }
  return (I);
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^3*x*y, w^2*x^2*y^2*z^2, y^3*z,x^4*z^4*t^4, w*x*y*z*t,w*x^6*y^5*z^4, x^2*z^4*t^3 , w^6*y^4*z^2,x^2*y^2*z^2;
 minbaseMon(I);
}
/////////////////////////////////////////////////////////////////////////////
//
proc gcdMon (poly f,poly g)
"USAGE:    gcdMon (f,g); f,g polynomials.
RETURN:   a monomial, the greatest common divisor of f and g.
ASSUME:   f and g are monomials of the basering.
EXAMPLE:  example gcdMon; shows an example.
"
{
  if (size(f) <> 1 or size(g) <> 1)
  {
    ERROR ("the input must be 2 monomials.");
  }
  return(gcd(f,g));

//  // Variables
//  int k;
//  intvec exp;
//  int nvar = nvars(basering);
//  // intput: monomials
//
//  intvec expf = leadexp(f);
//  intvec expg = leadexp(g);
//  // Nos quedamos con el menor exponente de cada variable.
//  for (k = 1 ; k <= nvar ; k++)
//    {
//      if (expf[k] <= expg[k])
//        {
//          exp[k] = expf[k];
//        }
//      else
//        {
//          exp[k] = expg[k];
//        }
//    }
//  // Devolvemos el monomio correspondiente al exponente que hemos
//  // calculado.
//  return(monomial(exp));
}
example
{"EXAMPLE:"; echo = 2;
  ring R = 0,(x,y,z,t),dp;
  poly f = x^3*z^5*t^2;
  poly g = y^6*z^3*t^3;
  gcdMon(f,g);
}
/////////////////////////////////////////////////////////////////////////////
//
proc lcmMon (poly f,poly g)
"USAGE:    lcmMon (f,g); f,g polynomials.
RETURN:   a monomial,the least common multiple of f and g.
ASSUME:   f,g are monomials of the basering.
EXAMPLE:  example lcmMon; shows an example.
"
{
  // Hay que verificar que son monomios
  if (size(f) <> 1 or size(g) <> 1)
  {
    ERROR ("the input must be 2 monomials.");
  }
  // Variables.
  int k;
  intvec exp;
  int nvar = nvars(basering);

  // No tenemos mas que tomar el mayor exponente.
  intvec expf = leadexp (f);
  intvec expg = leadexp (g);
        
  for (k = 1 ; k <= nvar ; k ++)
  {
    if (expf[k] <= expg[k])
    {
      exp[k] = expg[k];
    }
    else
    {
       exp[k] = expf[k];
    }
  }
  // Transformamos el vector de exponentes al monomio correspondiente.
  return(monomial (exp));
}
example
{"EXAMPLE:"; echo = 2;
  ring R = 0,(x,y,z,t),dp;
  poly f = x^3*z^5*t^2;
  poly g = y^6*z^3*t^3;
  lcmMon(f,g);
}
//////////////////////////////////////////////////////////////////////
//
proc membershipMon(poly f,ideal I)        
"USAGE:    membershipMon(f,I); f polynomial, I ideal.
RETURN:   1, if f lies in I; 0 otherwise.
          (-1 if I and f are nonzero and I is not a monomial ideal)
ASSUME:   I is a monomial ideal of the basering.
EXAMPLE:  example membershipMon; shows some examples
"
{
  def R = basering;
  // Variables
  int i,j,resp,k,control;
  intvec restf;
  // Si el ideal es cero no es monomial, pero la pertenencia no se da si
  // el polinomio es no nulo
  if ( size(I) == 0 && size(f) > 0)
  {
    return (0);
  }
 // El cero esta en todos los ideales.
  if (f == 0)
  {
    return (1);
  }
  // COMPROBACIONES
  control = checkIdeal(I);
  if (control == 0)
  {
    list isMon = isMonomialGB (I);
    if (isMon[1] == 0)
    {
      ERROR ("the ideal is not monomial.");
    }
    else
    {
      // Sabemos que I ya viene dado por el sistema minimal de
      // generadores, aunque aqui no sea necesario.
      I = isMon[2];
    }
  }
  // Quitamos ceros.
  I = simplify(I,2);
  int n = size(I);
  int m = size(f);
  int nvar = nvars(R);
  for (i=1 ; i<=m ; i++)
  {
    resp = 0;
    for (j=1 ; j<=n  ;j++)
    {
      // Vamos termino a termino viendo si son divididos por algun
      // generador. Trabajamos con los exponentes, pues es mas
      // simple.
      restf = leadexp(f) - leadexp(I[j]);
      for (k=1 ; k<=nvar; k++)
      {
        // Si hay una componente negativa es que este generador
        // no divide. Queremos entonces ir al siguiente
        // generador.
        if (restf[k] < 0)
        {
          break;
        }
      }
      // Si no ha habido componente nula, hemos encontrado un
      // divisor para el actual termino de f. En esta situacion
      // queremos pasar a otro termino de f, no seguir con los
      // otros generadores.
      if (k == nvar+1)
      {
        resp = 1;
        break;
      }
    }
    // Si hemos encontrado para el anterior termino, voy al
    // siguiente; en caso contrario salimos, pues desde que un
    // termino no sea dividido, f no esta en I.
    if (resp == 1)
    {
      f = f - lead(f);
    }
    else
    {
      break;
    }
  }
  // Si hemos recorrido todo el bucle, f esta en I.
  if (i == m+1)
  {
    return (1);
  }
  return (0);
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I =  w*x, x^2, y*z*t, y^5*t;
 poly f =  3*x^2*y + 6*t^5*z*y^6 - 4*x^2 + 8*w*x^5*y^6 - 10*y^10*t^10;
 membershipMon(f,I);
 poly g = 3*w^2*t^3 - 4*y^3*z*t^3 - 2*x^2*y^5*t + 4*x*y^3;
 membershipMon(g,I);
}
//////////////////////////////////////////////////////////////////////
//
proc intersectMon (ideal I,ideal J)        
"USAGE:    intersectMon (I,J); I,J ideals.
RETURN:   an ideal, the intersection of I and J.
          (it returns -1 if I or J are not monomial ideals)
ASSUME:   I,J are monomial ideals of the basering.
NOTE:     the minimal monomial generating set is returned.
EXAMPLE:  example intersectMon; shows some examples
"
{
  // Variables
  ideal K;
  int i,j,control;
  list isMon;
  // El ideal trivial no es monomial.
  if ( size(I) == 0 || size(J) == 0)
    {
      ERROR("one of the ideals is zero, hence not monomial.");
    }
 // COMPROBACIONES
  control = checkIdeal(I);
  if (control == 0)
    {
      isMon = isMonomialGB(I);
      if (isMon[1] == 0)
        {
          ERROR ("the first ideal is not monomial.");
        }
      else
        {
          // Sabemos que I ya viene dado por el sistema minimal de
          // generadores, aunque aqui no sea necesario.
          I = isMon[2];
        }
    }
  else
    {
      // Los generadores son monomiales, hallamos el sistema minimal
      I = minbase(I);
    }
  control = checkIdeal(J);
  if (control == 0)
    {
      isMon = isMonomialGB (J);
      if (isMon[1] == 0)
        {
          ERROR ("the second ideal is not monomial.");
        }
      else
        {
          // Sabemos que J ya viene dado por el sistema minimal de
          // generadores, aunque aqui no sea necesario.
          J = isMon[2];
        }
    }
  else
    {
      // Los generadores son monomiales,hallamos la base minimal
      J = minbase(J);
    }
  // Hemos asegurado que los ideales sean monomiales.
  // Quitamos ceros de la base para no alargar calculos.
  int n = size(I);
  int m = size(J);
  // Hallamos el m.c.m.de cada par de generadores de uno y otro ideal
  // y los a?adimos al ideal interseccion.
  for (i=1 ; i<=n ; i++)
    {
      for (j=1 ; j<=m ; j++)
        {
          K = K , lcmMon (I[i],J[j]);
        }
    }
  // Devolvemos el ideal ya dado por la base minimal porque sabemos
  // que este procedimiento genera muchos generadores redundantes.
  return(minbase(K));
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
 ideal J = w*x, x^2, y*z*t, y^5*t;
 intersectMon (I,J);
}
//////////////////////////////////////////////////////////////////////
//
proc quotientMon (ideal I,ideal J)        
"USAGE:    quotientMon (I,J); I,J ideals.
RETURN:   an ideal, the quotient I:J.
          (returns -1 if I or J is not monomial)
ASSUME:   I,J are monomial ideals of the basering.
NOTE:     the minimal monomial generating set is returned.
EXAMPLE:  example quotientMon; shows an example.
"
{
  // Variables
  int i,control,n;
  poly f;
  list isMon;
  //COMPROBACIONES
  if (size(I) == 0 || size(J) == 0)
    {
      ERROR("one of the ideals is zero, hence not monomial.");
    }
  control = checkIdeal(I);
  if (control == 0)
    {
      isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the first ideal is not monomial.");
        }
      else
        {
          // Sabemos que I ya viene dado por el sistema minimal de
          // generadores, aunque aqui no sea necesario.
          I = isMon[2];
        }
    }
  else
    {
      // Los generadores son monomiales,hallamos sistema minimal
      I = minbase(I);
    }
  control = checkIdeal(J);
  if (control == 0)
    {
      isMon = isMonomialGB (J);
      if (isMon[1] == 0)
        {
          ERROR ("the second ideal is not monomial.");
          return (-1);
        }
      else
        {
          // Sabemos que J ya viene dado por el sistema minimal de
          // generadores, aunque aqui no sea necesario.
          J = isMon[2];
        }
    }
  else
    {
      // Los generadores son monomiales, hallamos el sistema minimal
      J = minbase(J);
    }
  // Tenemos los ideales dados por su sistema minimal (aunque no es necesario, ahorra
  // calculos. En K vamos a tener la interseccion de los ideales.
  ideal K = 1;
  // ------ Empezamos a hacer el cociente.
  // Los ideales estan dados por su sistema minimal, con lo que no aparecen ceros.
  // Luego podemos usar size(J).
  n = size(J);
  for (i = 1 ; i <= n ; i++)
    {
      K = intersect (K ,quotientIdealMon(I,J[i]));
    }
  // Aqui tambien surgen muchos generadores que no forman parte de la
  // base minimal del ideal obtenido.
  return(minbase(K));
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
 ideal J = w*x, x^2, y*z*t, y^5*t;
 quotientMon (I,J);
}
//////////////////////////////////////////////////////////////////////
//
proc radicalMon(ideal I)
"USAGE:    radicalMon(I); I ideal
RETURN:   an ideal, the radical ideal of the ideal I.
          (returns -1 if I is not a monomial ideal)
ASSUME:   I is a monomial ideal of the basering.
NOTE:     the minimal monomial generating set is returned.
EXAMPLE:  example radicalMon; shows an example.
"
{
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Variables
  int i,m,j,control;
  poly f;
  intvec v;
  ideal rad_I;
  ideal I = fetch(R,I);
  // COMPROBACIONES
  control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores monomiales minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos sistema minimal
      I = minbase(I);
    }
  // Ya tenemos el ideal generado por la BASE MINIMAL MONOMIAL
  m = size(I);
  // Solo hay que poner exponente 1 a todas las variables que tengan
  // exponente >1.
  for (i = 1 ; i <= m ; i++)
    {
      v = leadexp(I[i]);
      f = 1;
      for (j = 1 ; j <= nvar ; j++)
        {
          if (v[j] <> 0)
            {
              f = f*x(j);
            }
        }
      rad_I = rad_I,f;
    }
  setring R;        
  // Hay que devolver el ideal dado por la base minimal, pues se
  // producen muchos generadores redundantes.
  ideal rad_I = fetch(R1,rad_I);
  kill R1;
  return( minbase (rad_I));
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
 radicalMon(I);
}
//////////////////////////////////////////////////////////////////////
//
proc isprimeMon (ideal I)        
"USAGE:    isprimeMon (I); I ideal
RETURN:   1, if I is prime; 0, otherwise.
          (returns -1 if I is not a monomial ideal)
ASSUME:   I is a monomial ideal of the basering.
EXAMPLE:  example isprimeMon; shows some example.
"
{
  // Variables
  int control,i,j,suma;
  intvec expin;
 // COMPROBACIONES
  control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos el sistema minimal
      I = minbase(I);
    }
  // Ya tenemos el ideal generado por la BASE MINIMAL
  if (maxdeg1(I) == 1)
    {
      return (1);
    }
  return (0);                
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w,y,t;
 isprimeMon (I);
 ideal J = w,y,t,x*z;
 isprimeMon (J);
}
//////////////////////////////////////////////////////////////////////
//
proc isprimaryMon (ideal I)        
"USAGE:    isprimaryMon (I); I ideal
RETURN:   1, if I is primary; 0, otherwise.
          (returns -1 if I is not a monomial ideal)
ASSUME:   I is a monomial ideal of the basering.
EXAMPLE:  example isprimaryMon; shows some examples
"
{
  def R = basering;
  // Variables
  int nvar,control,m,l,sub_in,i,j,k;
  intvec v,w,exp_gen;
  // COMPROBACIONES
  control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos el sistema minimal
      I = minbase(I);
    }
  // Ya tenemos el ideal generado por la BASE MINIMAL
  // Usamos la funcion "soporte" que hemos creado, para saber que
  // variables aparecen en la base de generadores como producto de si
  // mismas y tambien cuales son los generadores del ideal donde
  // tenemos que comprobar si aparecen tales variables.
  nvar = nvars(R);
  m = size(I);
  // Inicializo la ultima componente del vector para que contenga a
  // todas las variables (el subindice de la variable es el lugar
  // que ocupa como componente de v).
  v[nvar] = 0;
  k = 1;
  // v = 1 en variables solas y 0 en el resto.
  // w = lugar de los generadores de I que son producto de mas de una variable.
  for (i = 1 ; i <= m ; i++)
    {
      sub_in = soporte(I[i]);
      // Si soporte <> 0 la variable aparece sola, en caso contrario
      // el generador es producto de mas de una variable
      if (sub_in <> 0)
        {
          v[sub_in] = 1;
        }
      else
        {
          w[k] = i;
          k++;
        }
    }
  // Ahora solo hay que ver que no aparecen variables distintas de
  // las que tenemos marcadas con 1 en v.
  l = size(w);
  // Si w es cero, quiere decir que todos los generadores del ideal
  // son producto de una sola variable, luego es primario.
  if (l == 1 && w[1] == 0)
    {
      return(1);
    }
  // Estudiamos el exponente de los generadores de I oportunos (los
  // que nos indica w).
  for (i = 1 ; i <= l ; i++)
    {
      exp_gen = leadexp(I[w[i]]);
      for (j = 1 ; j <= nvar ; j++)
        {
          if (v[j] == 0)
            {
              if (exp_gen[j] <> 0)
                {
                  return (0);
                }
            }
        }
    }
  // Si hemos recorrido todo el ideal sin que salte el "for"
  // quiere decir que no se ha contradicho la caracterizacion,
  // luego el ideal es primario.
  return(1);
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^4,x^3,z^2,t^5,x*t,w*x^2*z;
 isprimaryMon (I);
 ideal J = w^4,x^3,z^2,t^5,x*t,w*x^2*z,y^3*t^3;
 isprimaryMon (J);
}
//////////////////////////////////////////////////////////////////////
//
proc isirreducibleMon (ideal I)
"USAGE:    isirreducibleMon(I); I ideal
RETURN:   1, if I is irreducible; 0, otherwise.
          (return -1 if I is not a monomial ideal)
ASSUME:   I is a monomial ideal of the basering.
EXAMPLE:  example isirreducibleMon; shows some examples
"
{
  // Variables
  intvec v;
  int n,i,j,sum,control;
  control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos sistema minimal
      I = minbase(I);
    }
  // Ya tenemos el ideal generado por la BASE MINIMAL
  n = size(I);
  // La funcion soporte devuelve 0 si el monomio es producto de mas
  // de una variable. Chequeamos generador a generador si el ideal
  // esta generado por potencias de variables.
  for (i = 1 ; i <= n ; i ++)
    {
      if (soporte(I[i]) == 0)
        {
          return(0);
        }
    }
  return (1);
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^4,x^3,z^2,t^5;
 isirreducibleMon (I);
 ideal J = w^4*x,x^3,z^2,t^5;
 isirreducibleMon (J);
}
//////////////////////////////////////////////////////////////////////
//
proc isartinianMon (ideal I)
"USAGE:    isartinianMon(I); I ideal.
RETURN:   1, if ideal is artinian; 0, otherwise.
          (return -1 if ideal I is not a monmomial ideal).
ASSUME:   I is a monomial ideal of the basering.
EXAMPLE:  example isartinianMon; shows some examples
"
{
  def R = basering;
  int  nvar  = nvars(R);
  // Declaracion de variables
  int i,j,k,cont,sizI;
  intvec v;
  // COMPROBACIONES
  int control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos sistema minimal
      I = minbase(I);
    }
  // Ya tenemos el ideal generado por la BASE MINIMAL
  sizI = size(I);
  // Comprobamos que entre los generadores minimales aparece una
  // potencia de cada. Cuando encontramos un generador que es potencia
  // de una sola variable aumento contador
  cont = 0;
  for (i = 1 ; i <= sizI ; i++)
    {
      if (soporte(I[i]) <> 0)
        {
          cont ++;
          // Solo volvemos a evaluar en caso de que hayamos aumentado
          if (cont == nvar)
            {
              // Ya hemos encontrado que todas las variables aparrecen en
              // el sistema minimal como potencia pura. No queremos seguir
              // buscando
              return (1);
            }
        }
    }
  // Si ha salido, es que faltan variables
  return(0);
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^4,x^3,y^4,z^2,t^6,w^2*x^2*y,w*z*t^4,x^2*y^3,z^2*t^5;
 isartinianMon (I);
 ideal J = w^4,x^3,y^4,z^2,w^2*x^2*y,w*z*t^4,x^2*y^3,z^2*t^5;
 isartinianMon (J);
}
//////////////////////////////////////////////////////////////////////
//
proc isgenericMon (ideal I)
"USAGE:    isgenericMon(I); I ideal.
RETURN:   1, if ideal is generic; 0, otherwise.
          (return -1 if ideal I is not a monomial ideal)
ASSUME:   I is a monomial ideal of the basering.
EXAMPLE:  example isgenericMon; shows some examples.
"
{
  def R = basering;
  int nvar = nvars(R);
  // Declaracion de variables
  int sizI,i,j,k;
  list exp;
 // COMPROBACIONES
  int control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos sistema minimal
      I = minbase(I);
    }
  // Ya tenemos el ideal generado por la BASE MINIMAL
  sizI = size(I);
  // Creamos una lista que tenga los exponentes de todos los
  // generadores.
  for (i = 1 ; i <= sizI ; i++)
    {
      exp[i] = leadexp(I[i]);
    }
  // Ahora hay que ver si alguno se repite, y si uno de ellos
  // lo hace, ya no es gen?rico.
  for (i = 1 ; i <= nvar ; i++)
    {
      for (j = 1 ; j <= sizI-1 ; j++)
        {
          for (k = j+1 ; k <= sizI ; k++)
            {
              // Notar que no se pueden repetir si la variable realmente
              // aparece en el generador, es decir, exponente >1.
              if (exp[j][i] == exp[k][i] & exp[j][i] <> 0)
                {
                  return (0);
                }
            }
        }
    }
  return (1);
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^4,x^3,y^4,z^2,w^2*x^2*y,w*z*t^4,x*y^3,z*t^5;
 isgenericMon (I);
 ideal J = w^4,x^3,y^4,z^3,w^2*x^2*y,w*z*t^4,x*y^3,z^2*t^5;
 isgenericMon (J);
}
//////////////////////////////////////////////////////////////////////
//
proc dimMon (ideal I)
"USAGE:   dimMon (I); I ideal
RETURN:  an integer, the dimension of the affine variety defined by
         the ideal I.
         (returns -1 if I is not a monomial ideal)
ASSUME:  I is a monomial ideal of the basering.
EXAMPLE: example dimMon; shows some examples.
"
{
  def R = basering;
  int nvar = nvars(R);
  // VARIABLES
  int control,sizSum,sumandos,i,j,k,cont;
  intvec index,indexAux,vaux,w;
  list sum, sumAux;
  // La base del ideal tiene que ser la monomial
 // El ideal trivial no es monomial.
  if ( size(I) == 0 )
    {
      ERROR("the ideal is zero, hence not monomial.");
    }
 // COMPROBACIONES
  control = checkIdeal(I);
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          // Sabemos que I ya viene dado por el sistema minimal de
          // generadores, aunque aqui no sea necesario.
          I = isMon[2];
        }
    }
  // Si el ideal es artiniano, es 0-dimensional
  if (isartinianMon(I) == 1)
    {
      return (0);
    }
  int sizI = size(I);
  // v(i) = vector con sizI entradas y donde cada entrada "j" se corresponde
  //        con el exponente del generador "i" en la variable "j"
  for (i = 1 ; i <= nvar ; i++)
    {
      intvec v(i);
      for (j = 1 ; j <= sizI ;j++ )
        {
          v(i)[j] = leadexp(I[j])[i];
        }
    }
  // Vamos a guardar en el vector "index" la ultima variable que se ha
  // sumado y en cada "sum(i)" el vector suma que se corresponde con la
  // entrada "i" del vector index.
  // Inicializo los valores de index y de cada sum
  w[sizI] = 0;
  sum[1] = w;
  index[1] = 0;
  sizSum = 1;
  while ( 1 )
    {
      cont = 1;
      sumandos ++;
      for (i = 1 ; i <= sizSum ; i ++)
        {
          for (j = index[i] + 1 ; j <= nvar ; j ++)
            {
              w = sum[i];
              vaux = w + v(j);
              // Comprobamos
              for (k = 1 ; k <= sizI ; k ++)
                {
                  if (vaux[k] == 0)
                    {
                      break;
                    }
                }
              if (k == sizI +1)
                {
                  return (nvar - sumandos);
                }
              if (j <> nvar)
                {
                  sumAux[cont] = vaux;
                  indexAux[cont] = j;
                  cont ++;
                }
            }
        
        }
      index = indexAux;
      sum = sumAux;
      sizSum = size(sumAux);
    }
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z,t),lp;
 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
 dimMon (I);
 ideal J = w^4,x^3,y^4,z^2,t^6,w^2*x^2*y,w*z*t^4,x^2*y^3,z*t^5;
 dimMon (J);
}
/////////////////////////////////////////////////////////////////////////////
//-------------------------------------------------------------------------//
//-----------------------  DESCOMPOSICIONES  -----------------------------_//
//-------------------------------------------------------------------------//
/////////////////////////////////////////////////////////////////////////////
//
// METODO 1: Metodo directo para descomp. irreducible (ver Vasconcelos)
//
//////////////////////////////////////////////////////////////////////
// Este procedimiento calcula la descomposicion en irreducibles de  //
// un ideal monomial dado por su base minimal de generadores        //
// haciendo uso de la caracterizacion de ideal monomial irreducible //
// (Vasconcelos)                                                    //
//////////////////////////////////////////////////////////////////////
//
static proc irredDec1 (ideal I)
{
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Variables
  int i,j,n,resp;
  list l1,l2;
  intvec exp,v;
  ideal J,K;
  ideal I = fetch(R,I);
  // ----- DESCOMPOSICION IRREDUCIBLE
  // Inicializamos la lista que va a estar formada por los ideales
  // que tenemos que comprobar son irreducibles.
  I = simplify(I,1);
  l1 = I;
  while (size(l1) > 0)
    {
      for (i = 1 ; i <= size(l1) ; i++)
        {
          J = l1[i];
          n = size(J);
          l1 = delete(l1,i);
          // Llamamos a la funcion que va a detectar si el ideal es o
          // no irreducible, y en caso de no serlo sabemos sobre que
          // generador y con que variables hay que aplicar el
          // el resultado teorico.
          v = irredAux (J);
          // No irreducible
          if (size(v) > 1)
            {
              // En este caso, v[1] nos indica el generador del ideal
              // que debemos eliminar.
              exp = leadexp(J[v[1]]);
              if (v[1] == 1)
                {
                  J = J[2..n];
                }
              if (v[1] > 1 && v[1] < n)
                {
                  J = J[1..v[1]-1],J[v[1]+1..n];
                }
              if (v[1] == n)
                {
                  J = J[1..n-1];
                }
              // Ahora vamos a introducir los nuevos generadores en
              // cada uno de los nuevos ideales que generamos. Los
              // ponemos en la lista en la que comprobamos.
              for (j = 1 ; j <= size(v)-1 ; j++)
                {
                  K = J,x(v[j+1])^exp[v[j+1]];
                  l1 = insert(l1,minbase(K));
                }
            }
          // Si v[1]=0, el ideal es irreducible y lo introducimos en
          // la lista l2, que es la que finalmente devolvera las
          // componentes de la descomposicion.
          else
            {
              l2 = insert(l2,J);
            }
        }
    }
  // ---- IRREDUNDANTE
  l2 = irredundant (l2);
  // La salida es la lista de los ideales irreducibles.
  setring R;
  list l2 = fetch(R1,l2);
  kill R1;
  return (l2);
}
//////////////////////////////////////////////////////////////////////
// La siguiente funcion va a obtener una descomposicion primaria    //
// minimal a partir de la irreducible anterior.                     //
//////////////////////////////////////////////////////////////////////
//
static proc primDec1 (ideal I)  {
  // Variables
  list l1,l2;
// ----- DESCOMPOSICION IRREDUCIBLE
  l1 = irredDec1 (I);
// ----- DESCOMPOSICION PRIMARIA
  l2 = irredPrimary (l1);
  return (l2);
}
//
// METODO 2: Metodo directo para descomp. primaria (ver Vasconcelos)
//
//////////////////////////////////////////////////////////////////////        
// La siguiente funcion va a calcular una descomposicion primaria   //
// minimal de un ideal monomial, pero esta vez usando la            //
// caracterizacion de ideal monomial primario y un resultado que    //
//  hace uso de esta (Vasconcelos).                                 //
//////////////////////////////////////////////////////////////////////
//
static proc primDec2 (ideal I)  {
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Variables en el nuevo anillo
  int i,n,max;
  list l1,l2;
  intvec v;
  ideal J,K;
  ideal I = fetch(R,I);
  // Vamos a tener dos listas: l1 que va a guardar todos los ideales
  // que vayamos generando con el resultado citado antes, que seran
  // los que vamos a comprobar si son primarios; y l2, donde guardamos
  // aquellos de l1 que verifiquemos son primarios.
  I = simplify(I,1);
  l1 = I;
  while (size(l1) > 0)
    {
      for (i = 1 ; i <= size(l1) ; i++)
        {
          J = l1[i];
          n = size(J);
          l1 = delete (l1,i);
          // Usamos la funcion que hemos creado para saber si el ideal
          // es primario. Si no lo es devuelve la variable que crea
          // conflicto y los generadores del ideal que luego hay que
          // usar (los que tienen productos de mas de una vble).
          // Se le llama con el sistema minimal de generadores
          v = primAux (J);
          // En caso de no ser primario, hay que buscar el maximo
          // exponente de la variable que aparece en los generadores
          // del ideal multiplicada siempre por otras variables,
          // nunca por si misma.
          if (v[1] == 0)
            {
              max = maxExp(J,v);
              K = J,x(v[2])^max;
              l1 = insert (l1,minbase(K));
              K = quotientIdealMon(J,x(v[2])^max);
              // quotientidealMon devuelve sistema minimal de generadores
              l1 = insert (l1,minbase(K));
            }
          // En caso de ser primario, lo introducimos en la lista
          // conveniente.
          else
            {
              l2 = insert (l2,J);
            }
        }
    }
// ------   IRREDUNDANTE
  l2 = irredundant (l2);
  setring R;
  list l2 = fetch (R1,l2);
  kill R1;
  return (l2);
}
//
// METODO 3: via dual de Alexander y doble dual (Miller)
//
//////////////////////////////////////////////////////////////////////
// Esta funcion calcula la descomposicion irreducible usando el     //
// dual de Alexander teniendo en cuenta que el dual del dual es el  //
// ideal de partida (Miller)                                        //
//////////////////////////////////////////////////////////////////////
//
static proc irredDec3 (ideal I)  {
  def R = basering;
  int i,n,j;
  poly lcm;
  intvec v,exp_I,exp;
  ideal J;
  list l;
  // Hallamos el m.c.m. de los generadores minimales del ideal.
  n = size (I);
  lcm = I[1];
  for ( i = 2 ; i <= n ; i++ )
    {
      lcm = lcmMon (lcm,I[i]);
    }
  v = leadexp (lcm);
  // Calculamos el dual de Alexander.
  // Hacemos una funcion para este paso porque luego lo volveremos a
  // utilizar.
  l = alexDif (v,I);
  // Tenemos los ideales irreducibles cuya interseccion nos da el dual
  // de Alexander. Notar que tenemos tantos ideales como generadores
  // minimales tiene I.
  // Hallamos la base minimal.
  J = minbase(intersect(l[1..size(l)]));
  // Ya solo queda el ultimo paso: hallar de nuevo el dual de
  // Alexander. Sabemos que este proceso ya devuelve la descomposicion
  // irreducible irredundante.
  l = alexDif (v,J);
  return (l);
}
//////////////////////////////////////////////////////////////////////
// En este caso hallamos una descomposicion primaria minimal usando //
// la irreducible irredundante del procedimiento anterior.          //
//////////////////////////////////////////////////////////////////////
//
static proc primDec3 (ideal I)   {
  // Variables
  list l1,l2;
// ----- DESCOMPOSICION IREDUCIBLE
  l1 = irredDec3 (I);
// ----- DESCOMPOSICION PRIMARIA
  l2 = irredPrimary (l1);
  return (l2);
}
//
// METODO 4: via dual de Alexander y cociente (Miller)
//
//////////////////////////////////////////////////////////////////////
// Vamos a hallar las componentes irreducibles de un ideal monomial //
// dado por sus generadores minimales haciendo uso del dual de      //
// Alexander (con el cociente) (Miller)                             //
//////////////////////////////////////////////////////////////////////
//
static proc irredDec4 (ideal I)  {
  // Cambiamos de anillo.
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
 // Variables
  int n,i,j,m,resp;
  poly lcm;
  intvec v;
  ideal J,K;
  list L;
  ideal I = fetch (R,I);
  // Calculamos el l.c.m. de los generadores minimales del ideal.
  n = size (I);
  lcm = I[1];
  for ( i = 2 ; i <= n ; i++ )
    {
      lcm = lcmMon (lcm,I[i]);
    }
  v = leadexp (lcm);
  // Hallamos el ideal J = (x(1)^(l(1)+1),...,x(n)^(l(n)+1)). Como
  // luego, en el procedimiento quotientMon, vamos a hallar la base
  // minimal de cada ideal, no nos preocupa que tengamos un cero en
  // el ideal J.
  for ( i = 1 ; i <= nvar ; i++ )
    {
      J[i] = (x(i))^(v[i]+1);
    }
  // Ahora hacemos el cociente oportuno.
  K = minbase(quotient (J,I));
  // Buscamos aquellos generadores de K que no son divisibles por los
  // generadores de J. Los generadores que son divisibles los hacemos
  // cero y luego los eliminamos.
  m = size (K);
  for ( i = 1 ; i <= m ; i++ )
    {
      resp = membershipMon(K[i],J);
      if ( resp == 1)
        {
          K[i] = 0;
        }
    }
  K = simplify (K,2);
  // Ahora obtenemos las componentes de la descomposicion irreducible,
  // que estan en correspondencia con los generadores minimales de K.
  L = alexDif (v,K);
  // Volvemos al anillo de partida y devolvemos la lista de las
  // componentes irreducibles irredundantes.
  setring R;
  list L = fetch(R1,L);
  kill R1;
  return (L);
}
//////////////////////////////////////////////////////////////////////
// Ahora hallamos una descomposicion primaria irredundante usando   //
// la ultima funcion para hallar las componentes irreducibles de un //
// ideal monomial dado por sus generadores minimales.               //
//////////////////////////////////////////////////////////////////////
//
static proc primDec4 (ideal I) {
  // Variables
  list l1,l2;
// ----- DESCOMPOSICION IREDUCIBLE
  l1 = irredDec4 (I);
// ----- DESCOMPOSICION PRIMARIA
  l2 = irredPrimary (l1);
  return (l2);
}
//
// METODO 5: un misterio!!
//
//////////////////////////////////////////////////////////////////////
// Este procedimiento halla los elementos maximales de la base      //
// estandar del ideal, que se correspoenden con las componentes     //
// irreducibles del ideal 1-1.                                      //
//////////////////////////////////////////////////////////////////////
//
static proc irredDec5 (ideal I)
{
  // New ring
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  ideal I = fetch (R,I);
  //Variables
  int i,j;
  ideal K;
  // Artinianization
  list artiniano = artinian (I);
  if (artiniano[1] == 0)
    {
      I = artiniano[2];
      intvec elimina = artiniano[3];
    }
  // Quotient
  ideal M = (x(1..nvar));
  ideal J = quotient (I,M);
  // Deleting generators lying in I
  for (i = 1 ; i <= size(J) ; i ++)
    {
      if (membershipMon(J[i],I) == 1)
        {
          if (i == 1)
            {
              J = J[2..size(J)];
              i --;
            }
          else
            {
              if (i == size(J))
                {
                  J = J[1..size(J)-1];
                  i --;
                }
              else
                {
                  J = J[1..i-1],J[i+1..size(J)];
                  i --;
                }
            }
        }
    }
  // Exponents of the ideals are going to form the decomposition
  int sizJ = size(J);
  for (i = 1 ; i <= sizJ ; i ++ )
    {
      intvec exp(i) = leadexp(J[i]) + 1;
    }
  // Deleting artinianization process
  if (artiniano[1] == 0)
   {
     // En elimina estan guardadas las variables que hay que eliminar
     for (i = 1 ; i <= nvar ; i ++)
       {
         if (elimina[i] <> 0)
           {
             for (j = 1 ; j <= sizJ ; j ++)
               {
                 if (exp(j)[i] == elimina[i])
                   {
                     exp(j)[i] = 0;
                   }
               }
           }
       }
   }
  // En exp(i) tengo los exponentes de cada variable de las que aparecen
  // en cada ideal.
  list facets;
  for (i = 1 ; i <= sizJ ; i ++)
    {
      J = 0;
      for (j = 1 ; j <= nvar ; j ++)
        {
          if (exp(i)[j] <> 0)
            {
              J = J,x(j)^exp(i)[j];
            }
        }
      J = simplify(J,2);
      facets[i] = J;
    }
  setring R;
  list facets = fetch (R1,facets);
  kill R1;
  return (facets);
}
//////////////////////////////////////////////////////////////////////
// Ahora hallamos una descomposicion primaria irredundante usando   //
// la ultima funcion para hallar las componentes irreducibles de un //
// ideal monomial dado por sus generadores minimales.               //
//////////////////////////////////////////////////////////////////////
//
static proc primDec5 (ideal I) {
  // Variables
  list l1,l2;
// ----- IRREDUCIBLE DECOMPOSITION
  l1 = irredDec5 (I);
// ----- PRIMARY DECOMPOSITION
  l2 = irredPrimary (l1);
  return (l2);
}
//
// METODO 6: via complejo de Scarf (Milovsky)
//
//////////////////////////////////////////////////////////////////////
// Metodo que usa el complejo de Scarf asociado a un ideal monomial //
// de k[x(1)..x(n)] (Milowski)                                      //
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
// Calcula el maximo exponente de la variable x(i) entre los        //
// generadores del ideal.                                           //
//////////////////////////////////////////////////////////////////////
//
static proc maximoExp(ideal I,int i)
{
  // VARIABLES
  int max,j,k,sizI;
  intvec exp;
  sizI = size(I);
  max = 0;
  for (j = 1 ; j <= sizI ; j ++)
    {
      exp = leadexp(I[j]);
      if ( exp[i] > max)
        {
          max = exp[i];
        }
    }
  return(max);
}
//////////////////////////////////////////////////////////////////////
// Esta funcion estudia si un ideal monomial dado por su sistema    //
// minimal de generadores es o no artiniano. En caso de no serlo,   //
// halla el artiniano mas proximo y ademas devuelve los generadores //
// que han sido introducidos.                                       //
//////////////////////////////////////////////////////////////////////
//
static proc artinian (ideal I)
{
  // Cambiamos de anillo
  def R = basering;
  int  nvar  = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Declaracion de variables
  int i,j,k,cont1,cont2,sizI,marcavar,max;
  intvec v,variablesin,cambio;
  list nuevo;
  ideal J;
  ideal I = fetch(R,I);
  sizI = size(I);
  // Comprobamos que entre los generadores minimales aparece una
  // potencia de cada
  cont2 = 0;
  for (i = 1 ; i <= sizI ; i++)
    {
      v = leadexp(I[i]);
      marcavar  = 0;
      cont1 = 0;
      for (j = 1 ; j <= nvar ; j++)
        {
          if (v[j] <> 0)
            {
              cont1 ++;
              marcavar = j;
            }
        }
      // Si cont1=1 hemos encontrado un generador de los que nos
      // interesa."variablesin" guarda el indice de las que estan.
      if (cont1 == 1)
        {
          cont2 ++;
          variablesin[cont2] = marcavar;
        }
    }
  // Notar que como el sistema de generadores es minimal no se
  // va a repetir la potencia de la misma variable. Por tanto basta
  // comprobar si cont2 es o no nvar.
  if (cont2 == nvar)
    {
      setring R;
      list output;
      output[1] = 1;
      return(output);
    }
  else
    {
      J = I;
      // Buscamos las que no estan
      if (cont2 == 0)
        {
          for (i = 1 ; i <= nvar ; i ++)
            {
              max = maximoExp(I,i);
              cambio[i] = max + 1;
              J = J,x(i)^(max + 1);
            }
        }
      else
        {
          for (i = 1 ; i <= nvar ; i++)
            {
              for (j = 1 ; j <= cont2 ; j ++)
                {
                  if (i == variablesin[j])
                    {
                      cambio[i] = 0;
                      break;
                    }
                }
              if (j == cont2 + 1)
                {
                  max = maximoExp(I,i);
                  cambio[i] = max + 1;
                  J = J,x(i)^(max + 1);
                }
            }
        }
      list output;
      output[1] = 0;
      output[2] = J;
      output[3] = cambio;
      setring R;
      list output = fetch (R1,output);
      kill R1;
      return(output);
    }
}
//////////////////////////////////////////////////////////////////////
// En este caso vamos primero a chequear si el ideal es o no        //
// generico y en caso de no serlo vamos a devolver  los cambios,    //
// pues estos son una aplicacion biyectiva.                         //
//////////////////////////////////////////////////////////////////////
//
static proc generic (ideal I)
{
  // New ring
  def R = basering;
  int nvar = nvars(R);
   string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
   execute(s);
   ideal I = fetch(R,I);
  // VARIABLES
  int i,j,k;
  // Cargamos la matriz con los exponentes
  int sizI = size(I);
  intmat EXP[sizI][nvar];
  for (i = 1 ; i <= sizI ; i ++)
    {
      // Ordenamos el vector de exponentes oportuno
      EXP[i,1..nvar] = leadexp(I[i]);
    }

  // Ahora tenemos que ordenarla segun la variable que este en conficto.
  intvec vari,change;
  intmat NEWEXP = EXP;
  list aux;
  int count = 0;
  for (i = 1 ; i <= nvar ; i ++)
    {
      // Buscamos conflicto en la variable x(i), para ello tengo que ordenar
      // la columna i de EXP
      vari = EXP[1..sizI,i];
      aux = sort(vari);
      vari = aux[1];
      change = aux[2];
      for (j = 1 ; j <= sizI - 1 ; j ++)
        {
          if (vari[j] > 0)
            {
              while (NEWEXP[change[j + count] , i] >= vari[j + 1 + count])
                {
                  NEWEXP[change[j + 1 + count] , i] = NEWEXP[change[j + count] , i] + 1;
                  count ++;
                  if (j + 1 + count > sizI)
                    {
                      break;
                    }
                }
            }
          j = j + count;
          count = 0;
        }
    }
  setring R;
  I = fetch(R1,I);
  // Devolvemos tambien la matriz, pues aqui tengo la biyeccion entre los exponentes
  if (EXP == NEWEXP)
    {
      return (1,I);
    }
  else
    {
      // Hallamos el ideal con los nuevos exponentes
      intvec expI;
      for (i = 1 ; i <= sizI ; i ++)
        {
          expI = NEWEXP[i,1..nvar];
          I[i] = monomial(expI);
        }
      return(0,I,EXP,NEWEXP);
    }
}
//////////////////////////////////////////////////////////////////////
// Esta funci?n obtiene una descomposicion irreducible del ideal    //
// monomial a partir de la descomposicion irreducible del idal      //
// generico que le asociamos.                                       //
//////////////////////////////////////////////////////////////////////
//
static proc nonGeneric (EXP,NEWEXP,Faces,sizI)
{
  def R = basering;
  int nvar = nvars(R);
  int sizFaces = size(Faces);
  // Variables
  int i,j,k;
  // Vamos variable a variable
  intvec exp, newexp;
  list newFaces;
  newFaces = Faces;
  for (i = 1 ; i <= nvar ; i ++)
    {
      // comparamos las matrices de exponentes por columnas
      exp = EXP[1..sizI,i];
      newexp = NEWEXP[1..sizI,i];
      if (exp <> newexp)
        {
          for (j = 1 ; j <= sizI ; j ++)
            {
              if (exp[j] <> newexp[j])
                {
                  for (k = 1 ; k <= sizFaces ; k ++)
                    {
                      if (Faces[k][i] == newexp[j])
                        {
                          newFaces[k][i] = exp[j];
                        }
                    }
                }
            }
        }
    }
  return (newFaces);
}
//////////////////////////////////////////////////////////////////////
// Este procedimiento va a dar una faceta inicial para el complejo  //
// de Scarf asociado a I, donde I es un ideal monomial artiniano    //
// y generico (evidentemente I dado por la bse minimal)             //
//////////////////////////////////////////////////////////////////////
//
static proc initialFacet (ideal I)
{
  // Cambiamos de anillo
  // Queremos usar x(1)..x(n) como variables.
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Declaracion de variables.
  int i,sizJ,j,max,aux,sizK,l, elim;
  intvec v;
  list face;
  // Traemos al nuevo anillo el ideal de entrada.
  ideal I = fetch(R,I);
  // Hacemos una copia de I pues ahora modificaremos el sistema
  // de generadores.
  ideal J = I;
  sizJ = size(J);
  // Para cada variable buscamos el maximo exponente, teniendo en
  // cuenta que un mismo generador no puede dar dos exponentes.
  // Vamos a guardar los exponentes en "expIni"
  intvec expIni;
  for (i = 1 ; i <= nvar ; i++)
    {
      max = 0;
      for (j = 1 ; j <= sizJ ; j++)
        {
          aux = leadexp(J[j])[i];
          if (aux > max)
            {
              max = aux;
              elim = j;
            }
        }
      // Guardamos el exponente
      expIni[i] = max;
      // Ahora tenemos el maximo de la variable x(i), por lo que
      // eliminamos el generador en el que la encontramos.
      // Si queda un generador no hay nada que hacer
      if (sizJ <> 1)
        {
          if (elim <> 1 & elim <> sizJ)
            {
              J = J[1..elim-1],J[elim+1..sizJ];
            }
          else
            {
              if (elim == 1)
                {
                  J = J[2..sizJ];
                }
              else
                {
                  J = J[1..sizJ-1];
                }
            }
          sizJ = size(J);
          // Eliminamos la variable x(i) en todos los generadores.
          for (j = 1 ; j <= sizJ ; j++)
            {
              v = leadexp(J[j]);
              if (v [i] <> 0)
                {
                  v[i] = 0;
                  J[j] = monomial(v);
                }
            }
          // Hallamos el sistema minimal de generadores que
          // corresponde al ideal que nos ha quedado.
          J = minbase(J);
          sizJ = size(J);
        }
    }
  // En expIni tenemos los exponentes de los monomios que dan la cara
  // inicial, ahora hay que buscar en el ideal original a que
  // generador se corresponde (el ideal es generico)
  int sizI = size(I);
  for (i = 1 ; i <= nvar ; i ++)
    {
      for (j = 1 ; j <= sizI ; j ++)
        {
          if (expIni[i] == leadexp(I[j])[i])
            {
              face = insert(face, I[j]);
              // Si lo encontramos buscamos el siguiente
              break;
            }
        }
    }
  setring R;
  list face = fetch (R1,face);
  kill R1;
  return (face);
}
//////////////////////////////////////////////////////////////////////
// La funcion que sigue devuelve las facetas adyacentes a una dada  //
// en el complejo de Scarf asociado a I.                            //
//////////////////////////////////////////////////////////////////////
//
static proc adyacency (list l1, ideal I)
{
  // Cambiamos de anillo
  // Queremos usar x(1)..x(n) como variables.
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Declaracion de variables.
  int i,j,cambio,sizI,k,min,m,generador,max;
  list l2,newl1,w,expI;
  poly LCM,newLCM;
  intvec v,newv,expgen;
  // Traemos al nuevo anillo las dos entradas.
  list l1 = fetch(R,l1);
  ideal I = fetch(R,I);
  sizI = size(I);
  // Hallamos el lcm de los elementos, e. d., la faceta del
  // complejo para luego comparar con los nuevos
  LCM = 1;
  for (i = 1 ; i <= nvar ; i++)
    {
      LCM = lcmMon(LCM,l1[i]);
    }
  v = leadexp(LCM);
  // Calculo los exponentes de los monomios de I
  for (i = 1 ; i <= sizI ; i++)
    {
      expI[i] = leadexp(I[i]);
    }
  // Hay que quitar cada elemento de la lista y comprobar si hay o no
  // cara adyacente al simplice que queda, y en caso de haberla, la
  // calculamos.
  for (i = 1 ; i <= nvar ; i++)
    {
      newl1 = delete(l1,i);
      // Hallamos el lcm de los elementos que hay en la nueva lista.
      newLCM = 1;
      for (j = 1 ; j <= nvar - 1 ; j++)
        {
          newLCM = lcmMon(newLCM,newl1[j]);
        }
      // Ahora hay que detectar si alguna variable ha desaparecido
      // en este LCM.
      newv = leadexp(newLCM);
      for (j = 1 ; j <= nvar ; j++)
        {
          if (newv[j] == 0)
            {
              l2[i] = 0;
              break;
            }
        }
      if (j == nvar+1)
        {
      // Si no ha habido ceros, queremos hallar la cara adyacente,
      // es decir, buscar un generador que introducido en l1 de una
      // faceta del complejo.
      // Comparamos los lcm entre s?, para comprobar en que variable
      // contribu?a el monomio que hemos eliminado.
          for (j = 1 ; j <= nvar ; j++)
            {
              if (v[j] <> newv[j])
                {
                  cambio = j;
                  // Una vez encontrado no hay mas
                  break;
                }
            }
          // Hallamos los exponentes de los monomios que quedan
          // para ver cual da el exponente en "newv"
          for (j = 1 ; j <= nvar - 1 ; j++)
            {
              w[j] = leadexp(newl1[j]);
            }
          for (j = 1 ; j <= nvar ; j++)
            {
              for (k = 1 ; k <= nvar - 1 ; k++)
                {
                  if (newv[cambio] == w[k][cambio])
                  {
                    cambio = k;
                    break;
                  }
                }
              // Si no termino el for con k es que hemos encontrado ya
              // los que son iguales.
              if (k <> nvar)
                {
                  break;
                }
            }
        
          // Donde contribuye antes, e.d., en "v"
          for (j = 1 ; j <= nvar ; j++)
            {
              if (w[cambio][j] == v[j])
                {
                  cambio = j;
                  break;
                }
            }
          // Ahora ya buscamos entre los generadores el nuevo monomio.
          // Ponemos de tope para encontrar el minimo el maximo de
          // las variables en el ideal
          max = 0;
          for (m = 1 ; m <= sizI ; m ++)
            {
              if (expI[m][cambio] > max)
                {
                  max = expI[m][cambio];
                }
            }
          min = max;        
          for (j = 1 ; j <= sizI ; j++)
            {
              for (k = 1 ; k <= nvar ; k ++)
                {
                  if (I[j] == l1[k])
                    {
                      break;
                    }
                }
              // El generador no esta en la lista, es de los que hay que
              // comprobar
              if (k == nvar +1)
                {
                  for (m = 1 ; m <= nvar ; m++)
                    {
                      if (m <> cambio)
                        {
                          if (expI[j][m] > newv[m])
                            {
                              break;
                            }
                        }
                      else
                        {
                          if (expI[j][m] <= newv[m])
                            {
                              break;
                            }
                        }
                    }
                  // Si termina el bucle cumple las condiciones
                  // oportunas, solo hay que comparar con el
                  // otro que tengamos.
                  if (m == nvar + 1)
                    {
                      if (expI[j][cambio] <=  min)
                        {
                          min = expI[j][cambio];
                          generador = j;
                        }
                    }
                }
            }
          // En la lista ponemos en el lugar "i" el generador que
          // hay que introducir cuando eliminamos el generador
          // "i" de la lista de entrada.
          l2[i] = I[generador];
        }
    }
  setring R;
  list l2 = fetch(R1,l2);
  kill R1;
  return(l2);
}
//////////////////////////////////////////////////////////////////////
// Metodo que calcula la descomposicion irreducible de un ideal     //
// monomial usando el complejo de Scarf (Milowski)                  //
//////////////////////////////////////////////////////////////////////
//
static proc ScarfMethod (ideal I)
{
  // Cambiamos de anillo
  // Queremos usar x(1)..x(n) como variables.
  def R = basering;
  int nvar = nvars(R);
  // Sabemos que dp siempre es mejor para trabajar, auqque luego para
  // comparar I y genI vamos a cambiarlo al lexicografico.
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Variables
  int i,j,k,sizl1,sizl,cont1,cont2;
  int sizI;
  list auxl,expl1,l1,l,l2,newLCM,expI,expgenI,Faces;
  poly LCM,mon;
  intvec v,w,betas;
  ideal J,genI,artgenI;
  ideal I = fetch(R,I);
  // Comprobamos si el ideal es generico y artiniano, y, en caso de
  // no serlo, obtenemos una modificacion de este ideal que si
  // verifique estas propiedades
  list genericlist = generic(I);
  if (genericlist[1] == 0)
    {
      genI = genericlist[2];
    }
  else
    {
      genI = I;
    }
  list artinianization = artinian (genI);
  if (artinianization[1] == 0)
    {
      artgenI = artinianization[2];
    }
  else
    {
      artgenI = genI;
    }
  // Una vez tenemos el ideal artiniano y generico, podemos hallar
  // el complejo de Scarf asociado al ideal modificado

  // Hay que obtener una cara inicial del ideal.
  list initial = initialFacet(artgenI);
  // Ahora de cada cara que tengamos en una lista obtenemos sus
  // caras adyacentes. Hay que tener en cuenta que si una cara la
  // obtengo como adyacente de otra, cuando calculemos sus adyacentes
  // sale la anterior, luego hay que evitar repetir.
  // Guardamos la primera faceta, su LCM
  LCM = 1;
  for (i = 1 ; i <= nvar ; i++)
    {
      mon = initial[i];
      LCM = lcmMon(LCM,mon);
    }
  v = leadexp(LCM);
  // Guardamos la primera faceta
  Faces[1] = v;
  int sizfaces = 1;
  // Lista de monomios que dan las facetas para hallar sus caras
  // adyacentes
  l[1] = initial;
  sizl = 1;
  // Ahora hayamos las posibles caras maximales adyacentes
  while (sizl <> 0)
    {
      // Hallamos la lista de monomios que hay que introducir
      // cuando eliminamos cada monomio.
      auxl = adyacency(l[1],artgenI);
      cont1 = 1;
      cont2 = 0;
      l1 = 0;
      for (j = 1 ; j <= nvar ; j++)
      {
        if (auxl[j] <> 0)
          {
            l2 = delete(l[1],j);
            l1[cont1] = insert(l2,auxl[j],cont1 + cont2 - 1);
            cont1 ++;
          }
        else
          {
            cont2 ++;
          }
      }
      // Hallamos los nuevos LCM
      sizl1 = size(l1);
      for (i = 1 ; i <= sizl1 ; i++)
        {
          newLCM[i] = 1;
         for (j = 1 ; j <= nvar ; j++)
           {
             newLCM[i] = lcmMon(newLCM[i],l1[i][j]);
           }
         expl1[i] = leadexp(newLCM[i]);
        }
      // Hallamos los LCM de las nuevas caras y eliminamos las que
      // ya esten en la lista Faces
      cont1 = 0;
      cont2 = 0;
      for (i = 1 ; i <= sizl1 ; i++)
        {
          for (j = 1 ; j <= sizfaces ; j++)
            {
              v = expl1[i];
              w = Faces[j];
              if (v == w)
                {
                  // Si ya esta el LCM en la lista, no queremos
                  // seguir buscando
                  break;
                }
            }
          // Si no ha salido del bucle en "j" es que este LCM
          // no esta en la lista de las caras, la introducimos
          if (j == sizfaces + 1)
            {
              Faces = insert(Faces,expl1[i],sizfaces + cont1);
              l = insert(l,l1[i]);
              cont1 ++;
            }
        }
      l = delete(l,cont1 + 1);
      sizl = size(l);
      sizfaces = size(Faces);
    }
  // En "Faces" ya tengo los exponentes que luego seran los exponentes
  // de los ideales que forman la descomposicion.
  // Deshacemos la artinianizacion
  intvec elimin = artinianization[3];
  if (artinianization[1] == 0)
   {
     // En elimina tenemos las variables que hemos introducido
     // y cual es la potencia
     // Solo miro las que tengan cambio
     for (i = 1 ; i <= nvar ; i ++)
       {
         if (elimin[i] <> 0)
           {
             for (j = 1 ; j <= sizfaces ; j ++)
               {
                 if (Faces[j][i] == elimin[i])
                   {
                     Faces[j][i] = 0;
                   }
               }
           }
       }
   }
  // Generico
  sizI = size(I);
  if (genericlist[1] == 0)
    {
      Faces = nonGeneric(genericlist[3],genericlist[4],Faces,sizI);
    }
  // Ya tenemos en Faces los exponentes de las componentes
  // ahora solo hay que obtener los ideales.
 for (i = 1 ; i <= sizfaces ; i ++)
    {
      J = 0;
      for (j = 1 ; j <= nvar ; j ++)
        {
          if (Faces[i][j] <> 0)
            {
              J = J,x(j)^(Faces[i][j]);
            }
        }
      J = simplify(J,2);
      Faces[i] = J;
    }
 // Esta es la parte LENTA computacionalmente si el ideal de partida
 // no es generico
 if (genericlist[1] == 0)
   {
     Faces = irredundant(Faces);
   }
 setring R;
 list components = fetch(R1,Faces);
 kill R1;
 return(components);
}
//////////////////////////////////////////////////////////////////////
// Devuelve una descomposicion primaria minimal de un ideal         //
// monomial via el complejo de Scarf.                               //
//////////////////////////////////////////////////////////////////////
//
static proc scarfMethodPrim (ideal I)
{
 // VARIABLES
  list l1,l2;
  // Hallamos la despomposicion irreducible del ideal dado usando
  // el complejo de Scarf
  l1 = ScarfMethod (I);
// ----- DESCOMPOSICION PRIMARIA
  l2 = irredPrimary (l1);
  return (l2);
}
//
// METODO 7: algoritmo de etiquetas (Roune)
//
//////////////////////////////////////////////////////////////////////
// Las siguientes funciones calculan la descomposicion en           //
// irreducibles de un ideal monomial. En este caso utilizamos el    //
// algoritmo de etiquetas de B. Roune.                              //
//////////////////////////////////////////////////////////////////////
//
static proc phi (list F)
{
  // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  // Variables
  int sizF,i,j;
  poly f;
  list listphi;
  intvec exp,newexp;
  list F = fetch(R,F);
  // F es una lista de pares, que indica una x(i) etiqueta de una
  // cara del ideal. Suponemos que F tiene ordenados sus elementos
  // segun las x(i)
  sizF = size(F);
  for (i = 1 ; i <= sizF ; i ++)
    {
      f = F[i];
      exp = leadexp(f);
      for (j = 1 ; j <= nvar ; j ++)
        {
          if (j <> i)
            {
              exp[j] = exp[j] + 1;
            }
        }
      listphi[i] = monomial(exp);
    }
  // Ya tenemos la lista de los monomios a los que
  // luego haremos el "lcm"
  setring R;
  list listphi = fetch (R1,listphi);
  kill R1;
  return (listphi);
}
//////////////////////////////////////////////////////////////////////
//
static proc pi(poly f)
{
 // Cambiamos de anillo
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  int i,sizI;
  intvec exp;
  poly f = fetch (R,f);
  exp = leadexp(f);
  for (i = 1 ; i <= nvar  ; i ++)
    {
      if (exp[i] <> 0)
        {
          exp[i] = exp[i] - 1;
        }
    }
  f = monomial(exp);
  setring R;
  f = fetch (R1,f);
  kill R1;
  return (f);
}
//////////////////////////////////////////////////////////////////////
//
static proc conditionComplex (intvec posActual,ideal I,ideal S)
{
  def R = basering;
  int nvar = nvars(R);
  // VARIABLES
  int i,nuevo;
  list F;
  // Vemos cual ha sido la ultima incorporacion al ideal, que es el
  // ultimo dentro de posActual que es distinto de 0.
  for (i = 1 ; i <= nvar ; i ++)
    {
      if (posActual[i] == 0)
        {
          break ;
        }
    }
  nuevo = i - 1;
  // No se pueden repetir generadores, se mira que el ultimo que se ha
 // ha introducido no sea de los que ya tenemos
  for (i = 1 ; i <= nuevo - 1 ; i ++)
    {
      if (posActual[i] == posActual[nuevo])
        {
          return (0);
        }
    }
  // Vemos si la variable oportuna divide al generador
  if (leadexp(I[i]) == 0)
    {
      return (0);
    }
  // Caso de que el LCM sea multiplo de los que ya tenemos
  poly LCM = 1;
  for (i = 1 ; i <= nuevo ; i ++)
    {
      F = insert (F,I[posActual[i]],size(F));
    }
  list phiF = phi(F);
  for (i = 1 ; i <= nuevo ; i ++)
    {
      LCM = lcmMon(phiF[i],LCM);
    }
  // Comprobamos si ya tenemos algun divisor del actual
  if (membershipMon(LCM,S) == 1)
    {
      return (0);
    }
  // Ahora vemos si la lista esta en el complejo simplicial
  if (membershipMon(LCM,I) == 1)
    {
      if (membershipMon(pi(LCM),I) == 0)
        {
          return (1,LCM);
        }
    }
  return (0);
}
//////////////////////////////////////////////////////////////////////
//
static proc findFaces (ideal I)
{
  def R = basering;
  int nvar = nvars(R);
  // Variables
  int i;
  ideal S;
  list condiciones;
  // Inicializamos valores
  list F;
  intvec posActual;
  posActual[nvar] = 0;

  int variable = 1;
  int sizI = size(I);
  while (1)
    {
      while (posActual[variable] == sizI)
        {
          posActual[variable] = 0;
          variable --;
          if (variable == 0)
            {
              break;
            }
        }
      // Comprobamos si hemos recorrido todas las posibilidades. Si
      // es as?, terminamos el while
      if (variable == 0)
        {
          break;
        }
      posActual[variable] = posActual[variable] + 1;
      // Comprobamos las condiciones para saber si los generadores que
      // tenemos est?n o no en el complejo.
      condiciones = conditionComplex (posActual,I,S);

      if (condiciones[1] == 1 )
        {
          if (posActual[nvar] <> 0)
            {
              S = S,condiciones[2];
              F = insert (F,condiciones[2]);
            }
          if (variable < nvar)
            {
              variable ++;
            }
        }
    }
  return (F);
}
//////////////////////////////////////////////////////////////////////
// La siguiente funcion calcula la descomposicion en irreducibles de//
// un ideal monomial artininano usando el algoritmo de etiquetas del//
// metodo de Bjarke Roune.                                          //
//////////////////////////////////////////////////////////////////////
//
static proc labelAlgorithm(ideal I)
{
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  ideal I = fetch(R,I);

  // Variables
  int i,j,sizComponents;
  list components;
  // El ideal tiene que ser artininano, si no lo es hacemos el cambio
  // oportuno para que lo sea (luego se deshace).
  ideal artI;
  list artiniano = artinian (I);
  if (artiniano[1] == 0)
    {
      artI = artiniano[2];
      intvec elimina = artiniano[3];
    }
  else
    {
      artI = I;
    }
  // Llamamos a findFaces para que encuentre las caras maximales del
  // complejo asociado al ideal
  components = findFaces(artI);
  sizComponents = size(components);
  list expComponents;
  poly f;
  for (i = 1 ; i <= sizComponents ; i ++)
    {
      f = components[i];
      expComponents[i] = leadexp(f);
    }
  // Deshacemos la artinianizacion
  if (artiniano[1] == 0)
   {
     // En elimina tenemos las variables que hemos introducido
     // y cual es la potencia
     // Solo miro las que tengan cambio
     for (i = 1 ; i <= nvar ; i ++)
       {
         if (elimina[i] <> 0)
           {
             for (j = 1 ; j <= sizComponents ; j ++)
               {
                 if (expComponents[j][i] == elimina[i])
                   {
                     expComponents[j][i] = 0;
                   }
               }
           }
       }
   }
  // En exp(i) tengo los exponentes de cada variable de las que aparecen
  // en cada ideal.
  ideal J;
  list facets;
  for (i = 1 ; i <= sizComponents ; i ++)
    {
      J = 0;
      for (j = 1 ; j <= nvar ; j ++)
        {
          if (expComponents[i][j] <> 0)
            {
              J = J,x(j)^expComponents[i][j];
            }
        }
      J = simplify(J,2);
      facets[i] = J;
    }
  setring R;
  list facets = fetch (R1,facets);
  kill R1;
  return (facets);
}
//////////////////////////////////////////////////////////////////////
// Devuelve una descomposicion primaria minimal de un ideal monomial//
// dado.                                                            //
//////////////////////////////////////////////////////////////////////
//
static proc  labelAlgPrim (ideal I)
{
  // VARIABLES
  list l1,l2;
  // Hallamos la despomposicion irreducible del ideal dado usando
  // el complejo de Scarf
  l1 = labelAlgorithm (I);
// ----- DESCOMPOSICION PRIMARIA
  l2 = irredPrimary (l1);
  return (l2);
}
//
// METODO 8: Gao-Zhu
//
//////////////////////////////////////////////////////////////////////
//
static proc divide (intvec v, intvec w, int k)
{
  def R = basering;
  int nvar = nvars(R);
  // Variables
  int i;
  for (i = 1 ; i <= nvar ; i ++)
    {
      if (i == k)
        {
          if (v[i] <> w[i])
            {
              return (0);
            }
        }
      else
        {
          if (v[i] >= w[i])
            {
              return (0);
            }
        }
    }
  return (1);
}
//////////////////////////////////////////////////////////////////////
//                                        //
//////////////////////////////////////////////////////////////////////
//
static proc incrementalAlg (ideal I)
{
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  ideal I = fetch(R,I);
  // COMPROBACIONES
  // Variables
  int i,sop,j,k,l,m,cont,cont2;
  intvec beta,dbeta,betaaux,elimina;
  // El ideal tiene que ser artininano, si no lo es hacemos el cambio
  // oportuno para que lo sea (luego se deshace).
  list artiniano = artinian (I);
  ideal artI;
  if (artiniano[1] == 0)
    {
      artI = artiniano[2];
      elimina = artiniano[3];
    }
  else
    {
      artI = I;
      elimina[nvar] = 0;
    }
  // Buscamos la primera componente irreducible o, lo que es lo
  // mismo, aquellos generadores que son potencia de una variable.
  // Si el tama?o de elimina es nvar es que hemos a?adido todos los
  // generadores que son potencia luego estar?n todos al final del
  // ideal.
  list MinI,componentes;
  int sizartI = size(artI);
  int sizelimina = size(elimina);
  for (i = 1 ; i <= nvar ; i ++)
    {
      if (elimina[i] == 0)
        {
          // Buscamos en el ideal los generadores que nos interesan
          for (j = 1 ; j <= sizartI ; j ++)
            {
              sop = soporte(artI[j]);
              if (sop <> 0)
                {
                  beta[sop] = leadexp(artI[j])[sop];
                  MinI = insert(MinI,leadexp(artI[j]));
                  if (j <> 1 and j <> sizartI)
                    {
                      artI = artI[1..j - 1],artI[j + 1..sizartI];
                    }
                  else
                    {
                      if (j == 1)
                        {
                          artI = artI[2..sizartI];
                        }
                      else
                        {
                          artI = artI[1..sizartI - 1];
                        }
                    }
                  sizartI = size(artI);
                  break;
                }
            }
        }
      else
        {
          // Buscamos la que esta al final
          sop = soporte(artI[sizartI]);
          beta[sop] = leadexp(artI[sizartI])[sop];
          MinI = insert(MinI,leadexp(artI[sizartI]));
          if (sizartI <> 1)
            {
              artI = artI[1..sizartI - 1];
            }
          else
            {
              artI = artI[1];
            }
          sizartI = size(artI);
        }
    }
  // En beta tenemos la primera componente
  componentes = insert(componentes,beta);
  int sizcomponents = size(componentes);
  int sizMin = size(MinI);
  // Es mas facil trabajar con los exponentes para nuestro objetivo
  // Se elige un nuevo generador, que en nuestro caso es un nuevo
  // exponente.
  int min,max;
  intvec expartI;
  for(i = 1 ; i <= sizartI ; i ++)
    {
      expartI = leadexp(artI[1]);
      if (size(artI) <> 1)
        {
          artI = artI[2..size(artI)];
        }
      // Hay que distinguir T_1 y T_2. Para ello se comparar vectores
      // de la lista actual de generadores con el que se acaba de
      // introducir.
      cont2 = 0;
      for (j = 1 ; j <= sizcomponents ; j ++)
        {
          beta = componentes[1 + cont2];
          // Si el nuevo generador divide a la componente beta, hay
          // que buscar las nuevas componentes
          for (k = 1 ; k <= nvar ; k ++)
            {
              if (expartI[k] >= beta[k])
                {
                  break;
                }
            }
          // Si el bucle anterior termino, divide y hay que hacer
          // los cambios.
          if (k == nvar + 1)
            {
              componentes = delete (componentes,1 + cont2);
              // Buscamos las nuevas componentes calculando las
              // distancias. Para cada variable busco d(beta,k,l)
              for (k = 1 ; k <= nvar ; k ++)
                {
                  betaaux = beta;
                  max = -1;
                  cont = 0;
                  dbeta = 0;
                  for (l = 1 ; l <= nvar ; l ++)
                    {
                      if (l <> k)
                        {
                          min = 32767;
                          cont ++;
                          for (m = 1 ; m <= sizMin ; m ++)
                            {
                              // Estos son de los buenos
                              if (divide(MinI[m],beta,l) == 1)
                                {
                                  if (MinI[m][k] < min)
                                    {
                                      min = MinI[m][k];
                                    }
                                }
                            }
                          dbeta[cont] = min;
                        }
                    }        
                  // Aqui ya tenemos d(beta,k,l) para cada k
                  // Hallamos el maximo cuando terminemos
                  for (l = 1 ; l <= cont ; l ++)
                    {
                      if (dbeta[l] > max)
                        {
                          max = dbeta[l];
                        }
                    }
                  // Condicion para introducir nueva componente
                  if (max < expartI[k])
                    {
                      betaaux[k] = expartI[k];
                      componentes = insert(componentes,betaaux,size(componentes));
                    }
                }
            }
          else
            {
              cont2 ++;
            }
        }
      MinI = insert(MinI,expartI);
      sizMin = size(MinI);
      sizcomponents = size(componentes);
    }
  // Deahacer los cambios de artiniano si se han hecho
  if (artiniano[1] == 0)
   {
     // En elimina tenemos las variables que hemos introducido
     // y cual es la potencia
     // Solo miro las que tengan cambio
     for (i = 1 ; i <= nvar ; i ++)
       {
         if (elimina[i] <> 0)
           {
             for (j = 1 ; j <= sizcomponents ; j ++)
               {
                 if (componentes[j][i] == elimina[i])
                   {
                     componentes[j][i] = 0;
                   }
               }
           }
       }
   }
  // En exp(i) tengo los exponentes de cada variable de las que aparecen
  // en cada ideal.
  ideal J;
  list facets;
  for (i = 1 ; i <= sizcomponents ; i ++)
    {
      J = 0;
      for (j = 1 ; j <= nvar ; j ++)
        {
          if (componentes[i][j] <> 0)
            {
              J = J,x(j)^componentes[i][j];
            }
        }
      J = simplify(J,2);
      facets[i] = J;
    }
  setring R;
  list facets = fetch (R1,facets);
  kill R1;
  return (facets);
}
//////////////////////////////////////////////////////////////////////
//
static proc  incrementalAlgPrim (ideal I)
{
  // VARIABLES
  list l1,l2;
  // Hallamos la despomposicion irreducible del ideal dado usando
  // el algoritmo de Gao-Zhu
  l1 = incrementalAlg (I);
// ----- DESCOMPOSICION PRIMARIA
  l2 = irredPrimary (l1);
  return (l2);
}
//
// METODO 9: slice algorithm (Roune)
//
//////////////////////////////////////////////////////////////////////
// SLICE ALGORITHM (B.Roune)                                        //
//////////////////////////////////////////////////////////////////////
//
static proc divideMon (poly f , poly g)
{
  def R = basering;
  int nvar = nvars(R);
  intvec expf = leadexp(f);
  intvec expg = leadexp(g);
  for (int i = 1 ; i <= nvar ; i ++)
    {
      if (expf[i] > expg[i])
        {
          return (0);
        }
    }
  return (1);
}
//////////////////////////////////////////////////////////////////////
//
static proc pivot (ideal I , poly lcmMin, ideal S)
{
  // I is monomial ideal
  def R = basering;
  I = fetch (R,I);
  S = fetch(R,S);
  int sizI = size(I);
  int nvar = nvars(R);
  intvec explcmMin = leadexp(lcmMin);
  // Variables
  int i,j;
  // The median estrategy
  poly p;
  int cont, exp, median, sizxi, max;
  intvec xiexp;
  for (i = 1 ; i <= nvar ; i ++)
    {
      if (explcmMin[i] >= 2 )
        {
          // Median exponent of x(i) from intersection(minI,x(i))
          cont = 0;
          for (j = 1 ; j <= sizI ; j ++)
            {
              exp = leadexp(I[j])[i];
              if (exp > 0)
                {
                  cont ++;
                  xiexp[cont] = exp;
                }
            }
          xiexp = sort(xiexp)[1];
          sizxi = size(xiexp);
          if (size(xiexp) == 1)
            {
              median = xiexp[1] - 1;
            }
          else
            {
              if (size(xiexp) == 2)
                {
                  median = xiexp[2] - 1;
                }
              else
                {
                  median = xiexp[(size(xiexp) + 1) / 2];
                }
            }
          p = x(i)^median;
          // valid pivot??
          if ( membershipMon(p,S) == 0)
            {
              return(p);
            }
          else
            {
              max = maximoExp(S,i);
              if ( xiexp[sizxi] == max )
                {
                  return(x(i)^(max-1));
                }
            }
          xiexp = 0;
        }
    }
}
//////////////////////////////////////////////////////////////////////
//
static proc simplification (I)
{
  // VARIABLES
  int i, j, k, cont, numdeleted;
  intvec isMaximal;
  def R = basering;
  I = fetch(R,I);
  int sizI = size(I);
  int nvar = nvars(R);
  poly lcmMinI = 1;
  for (i = 1 ; i <= sizI ; i ++)
    {
      lcmMinI = lcmMon(I[i],lcmMinI);
    }
  intvec explcmMinI = leadexp(lcmMinI);
  // Buscamos los elementos que son x(i) maximales. En caso de que
  // un generador del ideal sea maximal para 2 variables distintas,
  // ese generador se elimina.
  isMaximal[sizI] = 0;
  intvec genexp;
   for (i = 1 ; i <= sizI ; i ++)
    {
      genexp = leadexp(I[i]);
      cont = 0;
      for ( j = 1 ; j <= nvar ; j ++)
        {
          if (genexp[j] <> 0 && genexp[j] == explcmMinI[j])
            {
        
              if (cont == 0)
                {
                  cont ++;
                  isMaximal[i] = j;
                }
              else
                {
                  // Porque cuando encontramos que era maximal para
                  // la primera variable, lo guardamos
                  isMaximal[i] = 0;
                  // Eliminamos del ideal
                  if (i <> 1 && i <> sizI)
                    {
                      I = I[1..i - 1],I[i + 1..sizI];
                    }
                  else
                    {
                      if (i == 1)
                        {
                          I = I[2..sizI];
                        }
                      else
                        {
                          I = I[1..sizI - 1];
                        }
                    }
                  i --;
                  sizI = size(I);
                  // Generador i eliminado, miramos el siguiente
                  break;
                }
            }
        }
    }
   // En isMaximal[i] tenemos 0 si I[i] no es maximal,
   // y j si I[i] es maximal en x(j).                                        
  // Matriz de exponentes de los generadores del ideal
  intmat expI[sizI][nvar];
  for (i = 1 ; i <= sizI ; i++)
    {
      expI[i,1..nvar] = leadexp(I[i]);
    }
  // Buscamos ahora cota inferior
  poly lcmMi = 1;
  poly l,gcdMi;
  intvec Mi, mincol,expgcd;
  for (i = 1 ; i <= nvar ; i ++)
    {
      Mi = 0;
      cont = 0;
      for (j = 1 ; j <= sizI ; j ++)
        {
          // De isMaximal solo se usan las entradas que se corresponden con elementos del ideal
          if (expI[j,i] <> 0)
            {
              if (isMaximal[j] == 0 or isMaximal[j] == i)
                {
                  // Elementos del sistema minimal que estan
                  // en Mi
                  cont ++;
                  Mi[cont] = j;
                }
            }
        }
      // Si solo hay un elemento en Mi, no hay nada que hacer
      if (cont > 1)
        {
          gcdMi = I[Mi[1]];
          // Tenemos los generadores a los que hay que hallar el gcd
          for (j = 2; j <= cont ; j ++)
            {
              gcdMi = gcdMon(gcdMi,I[Mi[j]]);
            }
        }
      else
        {
          if (Mi <> 0)
            {
              gcdMi = I[Mi[1]];
            }
          else
            {
              // Falta alguna variable
              return (0,I);
            }
        }
      l = gcdMi/x(i);
      lcmMi = lcmMon(lcmMi,l);
    }
  // Ahora devolvemos la cota inferior, que luego hay que multiplicar
  // por el monomio que define el corte.
  // Devolvemos tambien el ideal (por si se ha modificado).
  return (lcmMi,I);
}
//////////////////////////////////////////////////////////////////////
//
static proc con (ideal I , ideal S , poly q)
{
  def R = basering;
  int nvar = nvars(R);
  I = fetch(R,I);
  S = fetch(R,S);
  q = fetch(R,q);
  // Variables
  int i;
  poly piI;
  int sizI = size(I);
  // Simplification process
  poly p;
  list sol;
  while (1)
    {
      // (I,S,q) normal slice?
      // Como cada vez que introducimos una cota inferior sabemos
      // que la slice actual es la inner slice (la otra es vacio),
      // hay que volver a verificar si es normal
      if ( S <> 0 )
        {
          // m/rad(m) esta en S, para m generador minimal de I??
          for (i = 1 ; i <= sizI ; i ++)
            {        
              piI = pi(I[i]);
              if (membershipMon(piI,S) == 1)
                {
                  if (i == 1)
                    {
                      I = I[2..sizI];
                    }
                  else
                    {
                      if (i == sizI)
                        {
                          I = I[1..sizI - 1];
                        }
                      else
                        {
                          I = I[1..i - 1],I[i + 1..sizI];
                        }
                    }
                  sizI = size(I);
                  i --;
                }
            }
        }
      // Buscamos cota inferior, y si es distinta de 1, simplificamos
      sol = simplification(I);
      p = sol[1];
      if (p == 1)
        {
          break;
        }
      else
        {
          if (p == 0)
            {
              break;
            }
          else
            {
              if (membershipMon(p,I) == 1 )
                {
                  break;
                }
            }
        }
      // Changing slice by simplification
      I = sol[2];
      I = minbase(quotient(I,p));
      q = p*q;
      S = minbase(quotient(S,p));
      sizI = size(I);
    }
  sizI = size(I);
  // (I,S,q) base case?
  poly lcmMinI;
  lcmMinI = 1;
  for (i = 1 ; i <= sizI ; i ++)
    {
      lcmMinI = lcmMon(lcmMinI,I[i]);
    }
  // a:b generates an intvec of length b with constant entries a
  intvec one = 1:nvar;
  if (divideMon(monomial(one),lcmMinI) == 0)
   {
     return (0);
   }
  if (equal(radicalMon(I),I) == 1)
    {
      if (equal(I, maxideal(1)) == 0)
        {
                return (0);
              }
       else
         {
          for (i = 1 ; i <= nvar ; i ++)
            {
              q = q * x(i);
            }
          return (q);
        }
    }
  // Selecting pivot
  p = pivot(I,lcmMinI,S);
  // New slices
  ideal S1 = minbase(quotient(S,p));
  ideal I1 = minbase(quotient(I,p));
  ideal S2 = S,p;
  S2 = minbase(S2);
  return (con(I1,S1,p*q),con(I,S2,q));
}
//////////////////////////////////////////////////////////////////////
//
static proc irredDecMonSlice (ideal I)
{
  // New ring
  def R = basering;
  int nvar = nvars(R);
  string s = "ring R1 =",charstr(R),",x(1..nvar),dp;";
  execute(s);
  ideal I = fetch(R,I);
  int sizI = size(I);
  int i,j;
  // Artinian ideal
  ideal artI;
  list artinianization = artinian(I);
  if (artinianization[1] == 0)
    {
      artI = artinianization[2];
    }
  else
    {
      artI = I;
    }
  // Easy case: 2 variables
  if (nvar == 2)
    {
      artI = sort(artI)[1];
      int sizartI = size(artI);
      for (i = 1 ; i <= sizartI - 1 ; i ++)
        {
          components[i] = x(1)^(leadexp[artI[i]][1])*x(2)^(leadexp[artI[i + 1]][2]);
        }
      setring R;
      list components = fetch(R1,components);
      kill R1;
      return (components);
    }
  ideal irredDec = con (artI,0,1);
  // Delelting zeros
  irredDec = simplify(irredDec,2);
  // Delting, in case, generators
  intvec elimina;
  if (artinianization[1] == 0)
    {
      elimina = artinianization[3];
    }
  else
    {
      elimina = 0;
    }
 // Each generator (monomial) corresponds to an ideal
  list components;
  poly comp;
  intvec exp;
  int sizIrred = size(irredDec);
  ideal auxIdeal;
  for (i = 1 ; i <= sizIrred ; i ++)
    {
      comp = irredDec[i];
      exp = leadexp(comp);
      for (j = 1 ; j <= nvar ; j ++)
        {
          if (exp[j] <> 0)
            {
              if (elimina <> 0)
                {
                  if (exp[j] == elimina[j])
                    {
                      auxIdeal[j] = 0;
                    }
                  else
                    {
                      auxIdeal[j] = x(j)^exp[j];
                    }
                }
              else
                {
                  auxIdeal[j] = x(j)^exp[j];
                }
            }
        }
      components[i] = simplify(auxIdeal,2);
      auxIdeal = 0;
    }
  setring R;
  list components = fetch(R1,components);
  kill R1;
  return (components);
}
//////////////////////////////////////////////////////////////////////
//
static proc primDecMonSlice (ideal I)
{
 // VARIABLES
  list l1,l2;
  // ---- Irreducible decomposition
  // Slice Method
  l1 = irredDecMonSlice (I);
// ----- Primary decomposition
  l2 = irredPrimary (l1);
  return (l2);
}
//////////////////////////////////////////////////////////////////////
//                                                                  //
//                       DECOMPOSITIONS                             //
//                                                                  //
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//
proc irreddecMon
"USAGE:    irreddecMon (I[,alg]); I ideal, alg string.
RETURN:   list, the irreducible components of the monomial ideal I.
          (returns -1 if I is not a monomial ideal).
ASSUME:   I is a monomial ideal of the basering k[x(1)..x(n)].
NOTE:     This procesure returns the irreducible decomposition of I.
          One may call the procedure with different algorithms using
          the optional argument 'alg':
          - the direct method following Vasconcelos' book (alg=vas)
          - via the Alexander dual and using doble dual (alg=add),
          - via the Alexander dual and quotients following E. Miller
            (alg=ad),
          - the formula of irreducible components (alg=for),
          - via the Scarf complex following Milowski (alg=mil),
          - using the label algorihtm of Roune (alg=lr),
          - using the algorithm of Gao-Zhu (alg=gz).
          - using the slice algorithm of Roune (alg=sr).
EXAMPLE:  example irreddecMon; shows some examples.
"
{
  // COMPROBACIONES
  ideal I = #[1];
  int control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos sistema minimal
      I = minbase(I);

    }
  // Si el ideal es irreducible, devolvemos el mismo
  if (isirreducibleMon(I) == 1)
    {
      return (I);
    }
 // Si no me han dado opcion, elijo una yo.
 if (size(#) == 1)
   {
     return (irredDec3(I));
   }
  // Leo la opcion y llamo al procedimiento oportuno
  else
    {
      if (#[2] == "vas")
        {
          return (irredDec1(I));
        }
      if (#[2] == "add")
        {
          return (irredDec3(I));
        }
      if (#[2] == "ad")
        {
          return (irredDec4(I));
        }
      if (#[2] == "for")
        {
          return (irredDec5(I));
        }
      if (#[2] == "mil")
        {
          return (ScarfMethod(I));
        }
      if (#[2] == "lr")
        {
          return (labelAlgorithm(I));
        }
      if (#[2] == "gz")
        {
          return (incrementalAlg(I));
        }
      if (#[2] == "sr")
        {
          return (irredDecMonSlice(I));
        }
    }
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z),Dp;
ideal I = w^3*x*y,w*x*y*z,x^2*y^2*z^2,x^2*z^4,y^3*z;
// Vasconcelos
 irreddecMon (I,"vas");
// Alexander Dual
 irreddecMon (I,"ad");
 // Scarf Complex
 irreddecMon (I,"mil");
 // slice algorithm
 irreddecMon(I,"sr");
}
//////////////////////////////////////////////////////////////////////
//
proc primdecMon
"USAGE:    primdecMon (I[,alg]); I ideal, alg string
RETURN:   list, the components in a minimal primary decomposition of I.
          (returns -1 if I is not a monomial ideal).
ASSUME:   I is a monomial ideal of the basering k[x(1)..x(n)].
NOTE:     This procesure returns a minimal primary decomposition of I.
          One may call the procedure with different algorithms using
          the optional argument 'alg':
          - the direct method for a primary decomposition following
            Vasconcelos' book (alg=vp),
          - from the irreducible decomposition obtained via the direct
            method following Vasconcelos' book (alg=vi),
          - from the irreducible decomposition obtained via the
            Alexander dual and using doble dual (alg=add),
          - from the irreducible decomposition obtained via the
            Alexander dual and quotients following E. Miller (alg=ad),
          - from the irreducible decomposition obtained
            via ........ (alg=for),
          - from the irreducible decomposition obtained via the Scarf
            complex following Milowski (alg=mil),
          - from the irreducible decomposition obtained using the label
            algorihtm of Roune (alg=lr),
          - from the irreducible decomposition obtained using the
            algorithm of Gao-Zhu (alg=gz),
          - from the irreducible decomposition obtained using the slice
            algorithm of Roune (alg=sr).
EXAMPLE:  example primdecMon; shows some examples.
"
{
  // COMPROBACIONES
  ideal I = #[1];
  int control = checkIdeal(I);
  // Si el sistema de generadores no esta formado por monomios, hay
  // que comprobar si el ideal es monomial.
  if (control == 0)
    {
      list isMon = isMonomialGB (I);
      if (isMon[1] == 0)
        {
          ERROR ("the ideal is not monomial.");
        }
      else
        {
          I = isMon[2];
          // Ya lo tenemos con los generadores minimales
        }
    }
  else
    {
      // Generadores monomiales, hallamos sistema minimal
      I = minbase(I);
    }
  // Estudiamos si el ideal es o no primario
  if (isprimaryMon(I) == 1)
    {
      return (I);
    }
  // Si no me han dado opcion, elijo una yo.
  if (size(#) == 1)
    {
      return(primDec3(I));
    }
  // Leo la opcion y llamo al procedimiento oportuno
  else
    {
      if (#[2] == "vi")
        {
          return (primDec1(I));
        }
      if (#[2] == "vp")
        {
          return (primDec2(I));
        }
      if (#[2] == "add")
        {
          return (primDec3(I));
        }
      if (#[2] == "ad")
        {
          return (primDec4(I));
        }
      if (#[2] == "for")
        {
          return (primDec5(I));
        }
      if (#[2] == "mil")
        {
          return (scarfMethodPrim(I));
        }
      if (#[2] == "lr")
        {
          return (labelAlgPrim(I));
        }
      if (#[2] == "gz")
        {
          return (incrementalAlgPrim(I));
        }
      if (#[2] == "sr")
        {
          return (primDecMonSlice(I));
        }
    }
}
example
{"EXAMPLE:"; echo = 2;
 ring R = 0,(w,x,y,z),Dp;
ideal I = w^3*x*y,w*x*y*z,x^2*y^2*z^2,x^2*z^4,y^3*z;
// Vasconcelos para primaria
 primdecMon(I,"vp");
// Alexander dual
 primdecMon(I,"add");
 // label algorithm
 primdecMon(I,"lr");
 //slice algorithm
 primdecMon(I,"sr");
}