Index: ngular/LIB/monomial.lib
===================================================================
--- Singular/LIB/monomial.lib	(revision 06a7af8bc2f204c56e7dc5f5662e372756701089)
+++ 	(revision )
@@ -1,3858 +1,0 @@
-//////////////////////////////////////////////////////////////////////////////
-version = "$Id$";
-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...).
-
-LITERATURE: Miller, Ezra  and Sturmfels, Bernd: Combinatorial Commutative Algebra,
-            Springer 2004
-
-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, i.e.,
-                       no two  minimal generators of it agree in the exponent
-                       of any variable that actually appears in both of them.
- 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 the given generators of I are monomials; 0, otherwise.
-"
-// Aqui NO estoy quitando el caso de que el ideal sea el trivial.
-{
-  int i,n;
-  n = ncols(I);
-  for (i = n ; i >= 1 ; 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 given 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 = ncols(I);
-  for (i = 1 ; i <= sizI ; i++)
-  {
-    g = gcd(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)].
-"
-{
-  // Variables
-  int i,cont,sop;
-  intvec expf;
-  int nvar = nvars(basering);
-  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 = nvar ; i >= 1 ; 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.
-"
-{
-  // Variables
-  int sizI,i,nvar,j,sum;
-  intvec w,exp;
-  sizI = ncols(I);
-  nvar = nvars(basering);
-  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++;
-        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 = ncols(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 (ncols(I) <> ncols(J))
-  {
-    return(0);
-  }
-  // 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
-  return(matrix( sort(I)[1])==matrix(sort(J)[1]));
-  //n = size(I);
-  //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
-  int nvar = nvars(basering);
-  // Variables
-  int i,cont;
-  intvec exp;
-  ideal rad;
-  // 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 = ncols(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] = var(i);
-      cont ++;
-    }
-  }
-  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.
-"
-{
-  // 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(basering);
-  int sizI = ncols(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
-  int nvar = nvars(basering);
-  // Variables
-  int i,j;
-  intvec exp_I,exp;
-  list l;
-  ideal J;
-  int sizI = ncols(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, var(j)^exp[j];
-      }
-    }
-    // Tenemos siempre un cero por la inicializacion de J, entonces
-    // lo quitamos.
-    J = simplify (J,2);
-    l = insert (l,J);
-  }
-  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 given generators of I are the minimal ones;
-          0 & minimal generators of I, otherwise.
-ASSUME:   I is an ideal of the basering given by monomial generators.
-"
-{
-  // VARIABLES
-  int i;
-  ideal J;
- // Quitamos los ceros del sistema de generadores.
-  I = simplify(I,2);
-  int resp = 1;
-  int sizI = ncols(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 = ncols(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 whose given generators are not
-          monomials.
-NOTE:     this procedure is NOT Grobner free and should be used only if the
-          given generators are not monomials. (use first checkIdeal)
-"
-{
-  // Variables
-  int resp;
-  // Si el ideal es cero, no es monomial.
-  if ( size(I) == 0)
-  {
-    return(0);
-  }
-  // Queremos la base de Grobner reducida, para uncidad.
-  intvec save_opt=option(get);
-  option(redSB);
-  // Base de Grobner
-  I = std(I);
-  option(set,save_opt);
-  // Una vez que tenemos la GB, no es mas que comprobar que el ideal
-  // esta generado por monomios.
-  resp = checkIdeal(I);
-  if (resp == 0)
-  {
-    return (0);
-  }
-  else
-  {
-    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)
-{
-  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.
-"
-{
-  // 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);
-  }
-  // Variables
-  int resp,m,k;
-  // Queremos la base de Grobner reducida, para uncidad.
-  intvec save_opt=option(get);
-  option(redSB);
-  // Hallamos GB
-  I = std(I);
-  option(set,save_opt);
-  // 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.
-  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 given 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)
-  {
-    ERROR ("the ideal is not monomial");
-  }
-  // Quitamos los ceros del sistema de generadores.
-  I = simplify(I,2);
-  int sizI = ncols(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 = ncols(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.");
-  }
-  return(f*g/gcd(f,g));
-  //// 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
-"
-{
-  // 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 = ncols(I);
-  int m = size(f);
-  int nvar = nvars(basering);
-  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.
-          USING THE SINGULAR COMMAND 'intersect' IS USUALLY FASTER
-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 = ncols(I);
-  int m = ncols(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.");
-    }
-    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 = ncols(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
-  int nvar = nvars(basering);
-  // Variables
-  int i,m,j,control;
-  poly f;
-  intvec v;
-  ideal rad_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 = ncols(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*var(j);
-      }
-    }
-    rad_I = rad_I,f;
-  }
-  // Hay que devolver el ideal dado por la base minimal, pues se
-  // producen muchos generadores redundantes.
-  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
-"
-{
-  // 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(basering);
-  m = ncols(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 = ncols(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
-"
-{
-  int  nvar  = nvars(basering);
-  // 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 = ncols(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.
-"
-{
-  int nvar = nvars(basering);
-  // 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 = ncols(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.
-"
-{
- // El ideal trivial no es monomial.
-  if ( size(I) == 0 )
-  {
-    ERROR("the ideal is zero, hence not monomial.");
-  }
-  // VARIABLES
-  int control,sizSum,sumandos,i,j,k,cont;
- // 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];
-    }
-  }
-  attrib(I,"isSB",1);
-  return (dim(I));
-//  int nvar = nvars(basering);
-//  intvec index,indexAux,vaux,w;
-//  list sum, sumAux;
-//  // La base del ideal tiene que ser la monomial
-//  // Si el ideal es artiniano, es 0-dimensional
-//  if (isartinianMon(I) == 1)
-//  {
-//    return (0);
-//  }
-//  int sizI = ncols(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)
-{
-  // Variables
-  int i,j,n,resp;
-  list l1,l2;
-  intvec exp,v;
-  ideal J,K;
-  // ----- 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 = ncols(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,var(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.
-  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)
-{
-  // Variables en el nuevo anillo
-  int i,n,max;
-  list l1,l2;
-  intvec v;
-  ideal J,K;
-  // 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 = ncols(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,var(v[2])^max;
-        l1 = insert (l1,minbase(K));
-        K = quotientIdealMon(J,var(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);
-  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)
-{
-  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 = ncols (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.
-  int nvar = nvars(basering);
- // Variables
-  int n,i,j,m,resp;
-  poly lcm;
-  intvec v;
-  ideal J,K;
-  list L;
-  // Calculamos el l.c.m. de los generadores minimales del ideal.
-  n = ncols (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] = (var(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 = ncols (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.
-  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)
-{
-  int nvar = nvars(basering);
-  //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 = maxideal(1);
-  ideal J = quotient (I,M);
-  // Deleting generators lying in I
-  for (i = 1 ; i <= ncols(J) ; i ++)
-  {
-    if (membershipMon(J[i],I) == 1)
-    {
-      if (i == 1)
-      {
-        J = J[2..ncols(J)];
-        i --;
-      }
-      else
-      {
-        if (i == ncols(J))
-        {
-          J = J[1..i-1];
-          i --;
-        }
-        else
-        {
-          J = J[1..i-1],J[i+1..ncols(J)];
-          i --;
-        }
-      }
-    }
-  }
-  // Exponents of the ideals are going to form the decomposition
-  int sizJ = ncols(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,var(j)^exp(i)[j];
-      }
-    }
-    J = simplify(J,2);
-    facets[i] = J;
-  }
-  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 = ncols(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
-  int  nvar  = nvars(basering);
-  // Declaracion de variables
-  int i,j,k,cont1,cont2,sizI,marcavar,max;
-  intvec v,variablesin,cambio;
-  list nuevo;
-  ideal J;
-  sizI = ncols(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)
-  {
-    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,var(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,var(i)^(max + 1);
-        }
-      }
-    }
-    list output;
-    output[1] = 0;
-    output[2] = J;
-    output[3] = cambio;
-    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
-  int nvar = nvars(basering);
-  // VARIABLES
-  int i,j,k;
-  // Cargamos la matriz con los exponentes
-  int sizI = ncols(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;
-    }
-  }
-  // 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)
-{
-  int nvar = nvars(basering);
-  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.
-  int nvar = nvars(basering);
-  // Declaracion de variables.
-  int i,sizJ,j,max,aux,sizK,l, elim;
-  intvec v;
-  list face;
-  // Hacemos una copia de I pues ahora modificaremos el sistema
-  // de generadores.
-  ideal J = I;
-  sizJ = ncols(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 = ncols(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 = ncols(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 = ncols(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;
-      }
-    }
-  }
-  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.
-  int nvar = nvars(basering);
-  // 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;
-  sizI = ncols(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];
-    }
-  }
-  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.
-  int nvar = nvars(basering);
-  // Sabemos que dp siempre es mejor para trabajar, auqque luego para
-  // comparar I y genI vamos a cambiarlo al lexicografico.
-  // 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;
-  // 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 = ncols(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,var(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);
-  }
-  return(Faces);
-}
-//////////////////////////////////////////////////////////////////////
-// 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
-  int nvar = nvars(basering);
-  // Variables
-  int sizF,i,j;
-  poly f;
-  list listphi;
-  intvec exp,newexp;
-  // 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"
-  return (listphi);
-}
-//////////////////////////////////////////////////////////////////////
-//
-static proc pi(poly f)
-{
- // Cambiamos de anillo
-  int nvar = nvars(basering);
-  int i,sizI;
-  intvec exp;
-  exp = leadexp(f);
-  for (i = nvar ; i > 0  ; i --)
-  {
-    if (exp[i] <> 0)
-    {
-      exp[i] = exp[i] - 1;
-    }
-  }
-  f = monomial(exp);
-  return (f);
-}
-//////////////////////////////////////////////////////////////////////
-//
-static proc conditionComplex (intvec posActual,ideal I,ideal S)
-{
-  int nvar = nvars(basering);
-  // 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)
-{
-  int nvar = nvars(basering);
-  // Variables
-  int i;
-  ideal S;
-  list condiciones;
-  // Inicializamos valores
-  list F;
-  intvec posActual;
-  posActual[nvar] = 0;
-
-  int variable = 1;
-  int sizI = ncols(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)
-{
-  int nvar = nvars(basering);
-
-  // 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,var(j)^expComponents[i][j];
-      }
-    }
-    J = simplify(J,2);
-    facets[i] = J;
-  }
-  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)
-{
-  int nvar = nvars(basering);
-  // Variables
-  int i;
-  for (i = nvar ; i > 0 ; 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)
-{
-  int nvar = nvars(basering);
-  // 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 = ncols(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 = ncols(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 = ncols(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 (ncols(artI) <> 1)
-    {
-      artI = artI[2..ncols(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,var(j)^componentes[i][j];
-      }
-    }
-    J = simplify(J,2);
-    facets[i] = J;
-  }
-  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)
-{
-  return (lead(g)/lead(f)!=0);
-  //int nvar = nvars(basering);
-  //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
-  int sizI = ncols(I);
-  int nvar = nvars(basering);
-  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 = var(i)^median;
-      // valid pivot??
-      if ( membershipMon(p,S) == 0)
-      {
-        return(p);
-      }
-      else
-      {
-        max = maximoExp(S,i);
-        if ( xiexp[sizxi] == max )
-        {
-          return(var(i)^(max-1));
-        }
-      }
-      xiexp = 0;
-    }
-  }
-}
-//////////////////////////////////////////////////////////////////////
-//
-static proc simplification (I)
-{
-  // VARIABLES
-  int i, j, k, cont, numdeleted;
-  intvec isMaximal;
-  int sizI = ncols(I);
-  int nvar = nvars(basering);
-  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 = ncols(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 = gcd(gcdMi,I[Mi[j]]);
-      }
-    }
-    else
-    {
-      if (Mi <> 0)
-      {
-        gcdMi = I[Mi[1]];
-      }
-      else
-      {
-        // Falta alguna variable
-        return (0,I);
-      }
-    }
-    l = gcdMi/var(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)
-{
-  int nvar = nvars(basering);
-  // Variables
-  int i;
-  poly piI;
-  int sizI = ncols(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 = ncols(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 = ncols(I);
-  }
-  sizI = ncols(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 * var(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)
-{
-  int nvar = nvars(basering);
-  int sizI = ncols(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 = ncols(artI);
-    for (i = 1 ; i <= sizartI - 1 ; i ++)
-    {
-      components[i] = var(1)^(leadexp[artI[i]][1])*var(2)^(leadexp[artI[i + 1]][2]);
-    }
-    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 = ncols(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] = var(j)^exp[j];
-          }
-        }
-        else
-        {
-          auxIdeal[j] = var(j)^exp[j];
-        }
-      }
-    }
-    components[i] = simplify(auxIdeal,2);
-    auxIdeal = 0;
-  }
-  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 procedure returns the irreducible decomposition of I,
-          i.e., the unique irredundant decomposition of I into irreducible
-          monomial ideals.
-          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 procedure 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");
-}
Index: Singular/LIB/monomialideal.lib
===================================================================
--- Singular/LIB/monomialideal.lib	(revision 2bb08398b9cc77848f151064fed100284feedc24)
+++ Singular/LIB/monomialideal.lib	(revision 2bb08398b9cc77848f151064fed100284feedc24)
@@ -0,0 +1,3858 @@
+//////////////////////////////////////////////////////////////////////////////
+version = "$Id$";
+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...).
+
+LITERATURE: Miller, Ezra  and Sturmfels, Bernd: Combinatorial Commutative Algebra,
+            Springer 2004
+
+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, i.e.,
+                       no two  minimal generators of it agree in the exponent
+                       of any variable that actually appears in both of them.
+ 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 the given generators of I are monomials; 0, otherwise.
+"
+// Aqui NO estoy quitando el caso de que el ideal sea el trivial.
+{
+  int i,n;
+  n = ncols(I);
+  for (i = n ; i >= 1 ; 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 given 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 = ncols(I);
+  for (i = 1 ; i <= sizI ; i++)
+  {
+    g = gcd(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)].
+"
+{
+  // Variables
+  int i,cont,sop;
+  intvec expf;
+  int nvar = nvars(basering);
+  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 = nvar ; i >= 1 ; 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.
+"
+{
+  // Variables
+  int sizI,i,nvar,j,sum;
+  intvec w,exp;
+  sizI = ncols(I);
+  nvar = nvars(basering);
+  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++;
+        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 = ncols(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 (ncols(I) <> ncols(J))
+  {
+    return(0);
+  }
+  // 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
+  return(matrix( sort(I)[1])==matrix(sort(J)[1]));
+  //n = size(I);
+  //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
+  int nvar = nvars(basering);
+  // Variables
+  int i,cont;
+  intvec exp;
+  ideal rad;
+  // 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 = ncols(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] = var(i);
+      cont ++;
+    }
+  }
+  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.
+"
+{
+  // 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(basering);
+  int sizI = ncols(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
+  int nvar = nvars(basering);
+  // Variables
+  int i,j;
+  intvec exp_I,exp;
+  list l;
+  ideal J;
+  int sizI = ncols(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, var(j)^exp[j];
+      }
+    }
+    // Tenemos siempre un cero por la inicializacion de J, entonces
+    // lo quitamos.
+    J = simplify (J,2);
+    l = insert (l,J);
+  }
+  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 given generators of I are the minimal ones;
+          0 & minimal generators of I, otherwise.
+ASSUME:   I is an ideal of the basering given by monomial generators.
+"
+{
+  // VARIABLES
+  int i;
+  ideal J;
+ // Quitamos los ceros del sistema de generadores.
+  I = simplify(I,2);
+  int resp = 1;
+  int sizI = ncols(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 = ncols(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 whose given generators are not
+          monomials.
+NOTE:     this procedure is NOT Grobner free and should be used only if the
+          given generators are not monomials. (use first checkIdeal)
+"
+{
+  // Variables
+  int resp;
+  // Si el ideal es cero, no es monomial.
+  if ( size(I) == 0)
+  {
+    return(0);
+  }
+  // Queremos la base de Grobner reducida, para uncidad.
+  intvec save_opt=option(get);
+  option(redSB);
+  // Base de Grobner
+  I = std(I);
+  option(set,save_opt);
+  // Una vez que tenemos la GB, no es mas que comprobar que el ideal
+  // esta generado por monomios.
+  resp = checkIdeal(I);
+  if (resp == 0)
+  {
+    return (0);
+  }
+  else
+  {
+    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)
+{
+  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.
+"
+{
+  // 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);
+  }
+  // Variables
+  int resp,m,k;
+  // Queremos la base de Grobner reducida, para uncidad.
+  intvec save_opt=option(get);
+  option(redSB);
+  // Hallamos GB
+  I = std(I);
+  option(set,save_opt);
+  // 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.
+  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 given 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)
+  {
+    ERROR ("the ideal is not monomial");
+  }
+  // Quitamos los ceros del sistema de generadores.
+  I = simplify(I,2);
+  int sizI = ncols(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 = ncols(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.");
+  }
+  return(f*g/gcd(f,g));
+  //// 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
+"
+{
+  // 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 = ncols(I);
+  int m = size(f);
+  int nvar = nvars(basering);
+  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.
+          USING THE SINGULAR COMMAND 'intersect' IS USUALLY FASTER
+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 = ncols(I);
+  int m = ncols(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.");
+    }
+    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 = ncols(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
+  int nvar = nvars(basering);
+  // Variables
+  int i,m,j,control;
+  poly f;
+  intvec v;
+  ideal rad_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 = ncols(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*var(j);
+      }
+    }
+    rad_I = rad_I,f;
+  }
+  // Hay que devolver el ideal dado por la base minimal, pues se
+  // producen muchos generadores redundantes.
+  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
+"
+{
+  // 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(basering);
+  m = ncols(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 = ncols(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
+"
+{
+  int  nvar  = nvars(basering);
+  // 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 = ncols(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.
+"
+{
+  int nvar = nvars(basering);
+  // 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 = ncols(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.
+"
+{
+ // El ideal trivial no es monomial.
+  if ( size(I) == 0 )
+  {
+    ERROR("the ideal is zero, hence not monomial.");
+  }
+  // VARIABLES
+  int control,sizSum,sumandos,i,j,k,cont;
+ // 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];
+    }
+  }
+  attrib(I,"isSB",1);
+  return (dim(I));
+//  int nvar = nvars(basering);
+//  intvec index,indexAux,vaux,w;
+//  list sum, sumAux;
+//  // La base del ideal tiene que ser la monomial
+//  // Si el ideal es artiniano, es 0-dimensional
+//  if (isartinianMon(I) == 1)
+//  {
+//    return (0);
+//  }
+//  int sizI = ncols(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)
+{
+  // Variables
+  int i,j,n,resp;
+  list l1,l2;
+  intvec exp,v;
+  ideal J,K;
+  // ----- 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 = ncols(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,var(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.
+  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)
+{
+  // Variables en el nuevo anillo
+  int i,n,max;
+  list l1,l2;
+  intvec v;
+  ideal J,K;
+  // 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 = ncols(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,var(v[2])^max;
+        l1 = insert (l1,minbase(K));
+        K = quotientIdealMon(J,var(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);
+  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)
+{
+  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 = ncols (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.
+  int nvar = nvars(basering);
+ // Variables
+  int n,i,j,m,resp;
+  poly lcm;
+  intvec v;
+  ideal J,K;
+  list L;
+  // Calculamos el l.c.m. de los generadores minimales del ideal.
+  n = ncols (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] = (var(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 = ncols (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.
+  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)
+{
+  int nvar = nvars(basering);
+  //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 = maxideal(1);
+  ideal J = quotient (I,M);
+  // Deleting generators lying in I
+  for (i = 1 ; i <= ncols(J) ; i ++)
+  {
+    if (membershipMon(J[i],I) == 1)
+    {
+      if (i == 1)
+      {
+        J = J[2..ncols(J)];
+        i --;
+      }
+      else
+      {
+        if (i == ncols(J))
+        {
+          J = J[1..i-1];
+          i --;
+        }
+        else
+        {
+          J = J[1..i-1],J[i+1..ncols(J)];
+          i --;
+        }
+      }
+    }
+  }
+  // Exponents of the ideals are going to form the decomposition
+  int sizJ = ncols(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,var(j)^exp(i)[j];
+      }
+    }
+    J = simplify(J,2);
+    facets[i] = J;
+  }
+  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 = ncols(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
+  int  nvar  = nvars(basering);
+  // Declaracion de variables
+  int i,j,k,cont1,cont2,sizI,marcavar,max;
+  intvec v,variablesin,cambio;
+  list nuevo;
+  ideal J;
+  sizI = ncols(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)
+  {
+    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,var(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,var(i)^(max + 1);
+        }
+      }
+    }
+    list output;
+    output[1] = 0;
+    output[2] = J;
+    output[3] = cambio;
+    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
+  int nvar = nvars(basering);
+  // VARIABLES
+  int i,j,k;
+  // Cargamos la matriz con los exponentes
+  int sizI = ncols(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;
+    }
+  }
+  // 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)
+{
+  int nvar = nvars(basering);
+  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.
+  int nvar = nvars(basering);
+  // Declaracion de variables.
+  int i,sizJ,j,max,aux,sizK,l, elim;
+  intvec v;
+  list face;
+  // Hacemos una copia de I pues ahora modificaremos el sistema
+  // de generadores.
+  ideal J = I;
+  sizJ = ncols(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 = ncols(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 = ncols(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 = ncols(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;
+      }
+    }
+  }
+  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.
+  int nvar = nvars(basering);
+  // 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;
+  sizI = ncols(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];
+    }
+  }
+  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.
+  int nvar = nvars(basering);
+  // Sabemos que dp siempre es mejor para trabajar, auqque luego para
+  // comparar I y genI vamos a cambiarlo al lexicografico.
+  // 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;
+  // 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 = ncols(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,var(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);
+  }
+  return(Faces);
+}
+//////////////////////////////////////////////////////////////////////
+// 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
+  int nvar = nvars(basering);
+  // Variables
+  int sizF,i,j;
+  poly f;
+  list listphi;
+  intvec exp,newexp;
+  // 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"
+  return (listphi);
+}
+//////////////////////////////////////////////////////////////////////
+//
+static proc pi(poly f)
+{
+ // Cambiamos de anillo
+  int nvar = nvars(basering);
+  int i,sizI;
+  intvec exp;
+  exp = leadexp(f);
+  for (i = nvar ; i > 0  ; i --)
+  {
+    if (exp[i] <> 0)
+    {
+      exp[i] = exp[i] - 1;
+    }
+  }
+  f = monomial(exp);
+  return (f);
+}
+//////////////////////////////////////////////////////////////////////
+//
+static proc conditionComplex (intvec posActual,ideal I,ideal S)
+{
+  int nvar = nvars(basering);
+  // 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)
+{
+  int nvar = nvars(basering);
+  // Variables
+  int i;
+  ideal S;
+  list condiciones;
+  // Inicializamos valores
+  list F;
+  intvec posActual;
+  posActual[nvar] = 0;
+
+  int variable = 1;
+  int sizI = ncols(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)
+{
+  int nvar = nvars(basering);
+
+  // 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,var(j)^expComponents[i][j];
+      }
+    }
+    J = simplify(J,2);
+    facets[i] = J;
+  }
+  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)
+{
+  int nvar = nvars(basering);
+  // Variables
+  int i;
+  for (i = nvar ; i > 0 ; 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)
+{
+  int nvar = nvars(basering);
+  // 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 = ncols(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 = ncols(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 = ncols(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 (ncols(artI) <> 1)
+    {
+      artI = artI[2..ncols(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,var(j)^componentes[i][j];
+      }
+    }
+    J = simplify(J,2);
+    facets[i] = J;
+  }
+  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)
+{
+  return (lead(g)/lead(f)!=0);
+  //int nvar = nvars(basering);
+  //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
+  int sizI = ncols(I);
+  int nvar = nvars(basering);
+  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 = var(i)^median;
+      // valid pivot??
+      if ( membershipMon(p,S) == 0)
+      {
+        return(p);
+      }
+      else
+      {
+        max = maximoExp(S,i);
+        if ( xiexp[sizxi] == max )
+        {
+          return(var(i)^(max-1));
+        }
+      }
+      xiexp = 0;
+    }
+  }
+}
+//////////////////////////////////////////////////////////////////////
+//
+static proc simplification (I)
+{
+  // VARIABLES
+  int i, j, k, cont, numdeleted;
+  intvec isMaximal;
+  int sizI = ncols(I);
+  int nvar = nvars(basering);
+  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 = ncols(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 = gcd(gcdMi,I[Mi[j]]);
+      }
+    }
+    else
+    {
+      if (Mi <> 0)
+      {
+        gcdMi = I[Mi[1]];
+      }
+      else
+      {
+        // Falta alguna variable
+        return (0,I);
+      }
+    }
+    l = gcdMi/var(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)
+{
+  int nvar = nvars(basering);
+  // Variables
+  int i;
+  poly piI;
+  int sizI = ncols(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 = ncols(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 = ncols(I);
+  }
+  sizI = ncols(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 * var(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)
+{
+  int nvar = nvars(basering);
+  int sizI = ncols(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 = ncols(artI);
+    for (i = 1 ; i <= sizartI - 1 ; i ++)
+    {
+      components[i] = var(1)^(leadexp[artI[i]][1])*var(2)^(leadexp[artI[i + 1]][2]);
+    }
+    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 = ncols(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] = var(j)^exp[j];
+          }
+        }
+        else
+        {
+          auxIdeal[j] = var(j)^exp[j];
+        }
+      }
+    }
+    components[i] = simplify(auxIdeal,2);
+    auxIdeal = 0;
+  }
+  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 procedure returns the irreducible decomposition of I,
+          i.e., the unique irredundant decomposition of I into irreducible
+          monomial ideals.
+          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 procedure 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");
+}
Index: Singular/LIB/normal.lib
===================================================================
--- Singular/LIB/normal.lib	(revision 06a7af8bc2f204c56e7dc5f5662e372756701089)
+++ Singular/LIB/normal.lib	(revision 2bb08398b9cc77848f151064fed100284feedc24)
@@ -2550,4 +2550,5 @@
 "
 {
+   intvec save_opt=option(get);
    option(redSB);
    def R=basering;
@@ -2633,5 +2634,5 @@
       }
    }
-   option(noredSB);
+   option(set,save_opt);
    ideal fstd=std(ideal(f)+jacob(f));
    poly hc=highcorner(fstd);
