proc min_generating_set (matrix P,S)
"USAGE:  min_generating_set(P,S);   P,S matrix
ASSUME: The entries of P,S are homogeneous and ordered by ascending 
        degrees. The first entry of S equals 1. (As satisfied by 
        the first two output matrices of invariant_ring(G).)
RETURN: ideal
NOTE:   The given generators for the output ideal form a minimal 
        generating set for the ring generated by the entries of 
        P,S. The generators are homogeneous and ordered by 
        descending degrees.
"
{
  if (defined(flatten)==0) { LIB "matrix.lib"; }
  ideal I1,I2 = flatten(P),flatten(S); 
  int i1,i2 = size(I1),size(I2);
  // We order the generators by descending degrees
  // (the first generator 1 of I2 is omitted):
  int i,j,s = i1,i2,i1+i2-1;
  ideal I;
  for (int k=1; k<=s; k++)
  { 
    if (i==0) { I[k]=I2[j]; j--; } 
    else
    { 
      if (j==0) { I[k]=I1[i]; i--; } 
      else 
      {
        if (deg(I1[i])>deg(I2[j])) { I[k]=I1[i]; i--; }
        else { I[k]=I2[j]; j--; }
      }
    }
  } 
  intvec deg_I = deg(I[1..s]);
  int n = nvars(basering);
  def BR = basering;

  // Create a new ring with elimination order:
  //---------------------------------------------------------------
  // ****    this part uses the command ringlist which is      ****
  // ****     only available in SINGULAR-3-0-0 or newer        ****
  //---------------------------------------------------------------
  list rData = ringlist(BR);
  intvec wDp;
  for (k=1; k<=n; k++) { 
    rData[2][k] ="x("+string(k)+ ")"; 
    wDp[k]=1;
  }
  for (k=1; k<=s; k++) { rData[2][n+k] ="y("+string(k)+ ")"; } 
  rData[3][1] = list("dp",wDp);
  rData[3][2] = list("wp",deg_I);
  def R_aux = ring(rData);
  setring R_aux;
  //---------------------------------------------------------------

  ideal J;
  map phi = BR, x(1..n);
  ideal I = phi(I);
  for (k=1; k<=s; k++) { J[k] = y(k)-I[k]; } 
  option(redSB);
  J = std(J);

  // Remove all generators that are depending on some x(i) from J: 
  int s_J = size(J);
  for (k=1; k<=s_J; k++) { if (J[k]>=x(n)) {J[k]=0;} }

  // The monomial order on K[y] is chosen such that linear leading 
  // terms in J are in 1-1 correspondence to superfluous generators
  // in I :
  ideal J_1jet = std(jet(lead(J),1));
  intvec to_remove; 
  i=1;
  for (k=1; k<=s; k++)
  { 
    if (reduce(y(k),J_1jet)==0){ to_remove[i]=k; i++; }
  }
  setring BR;
  if (to_remove == 0) { return(ideal(I)); }
  for (i=1; i<=size(to_remove); i++)
  { 
    I[to_remove[i]] = 0;
  } 
  I = simplify(I,2);
  return(I);
}