proc monomid(ideal I)
"USAGE: monomid(I); I ideal 
RETURN: ideal, the largest monomial subideal in I
NOTE: implements Algorithm 4.4.2
        Saito, Sturmfels, Takayama,
        Groebner Deformations of Hypergeometric Differential Equations
EXAMPLE: example monomid; shows an example
"
{
  def d = basering;
  int n = nvars(basering);
  int i,j,k;
   
  intvec v,w;
 
  // step 0 extend ring

  list rl = ringlist(d);
  for (i=1;i<=n;i++)
  {
    rl[2]= insert(rl[2],"@u("+string(i) +")",n+i-1);
  }
  def rnew = ring(rl);
  setring rnew;
  poly g,h;
  ideal I = fetch(d,I);
 
  // step 1 multi-homogenize
   
  for (i=1;i<=ncols(I);i++)  // each ideal entry
  {
    h = I[i];
    v = 0:2*n;
   
    for (k=1;k<=n;k++)   // determine multidegree
    {
      v[k] = 1;
      w[k] = deg(h,v);   // deg w.r.t. e_k
      v = 0:2*n;
    }
    // w is the multi-degree
   
    for(j=1;j<=size(h);j++)    // each term
    {
      v = leadexp(h[j]);
      v = 0:n,intvec(w[1..n])-intvec(v[1..n]);  // exponent vector
      g = g + h[j]*monomial(v);
    }
    I[i] = g;
  }
  // step 2 saturated quotient

  option(redSB);
  h = 1;
  for (i=1;i<=n;i++) { h = h * var(n+i);}
  I = sat(I,h)[1];
 
  // step 3 select monomials
  ideal J;
  j=0;
 
  for (i=1;i<=ncols(I);i++)
  {
    if (size(I[i])==1 and deg(I[i])>0)
    {
      j++;
      J[j]=I[i];
    }
  }
  setring d;
  ideal J = fetch(rnew,J);
return(J);
}

> LIB "elim.lib";
> ring r=0,(x,y,z),dp;

> // The example in Saito' et al. book on p. 176
>  ideal J = x+y+z,x2+y2+z2,x3+y3+z3;
> momomid(J);
_[1]=z3
_[2]=xyz
_[3]=y3
_[4]=x3
_[5]=y2z2
_[6]=x2z2
_[7]=x2y2

> // The example discussed in this forum:
> ideal I = x^2+y, xy, x^2, z+y;
> monomid(I);
_[1]=z
_[2]=y
_[3]=x2