//SINGULAR Procedure4.6.8


//=============================== we need ========================
LIB"elim.lib";
LIB"ring.lib";
proc prepareQuotientring(ideal i)
"USAGE:   prepareQuotientring(i); i standard basis
RETURN:   a list l of two strings:
          l[1] to define K[x\u,u ], u a maximal independent 
          set for i
          l[2] to define K(u)[x\u ], u a maximal independent 
          set for i
          both rings with lexicographical ordering
"  
{
  string va,pa;
//v describes the independent set u: var(j) is in 
//u iff v[j]!=0
  intvec v=indepSet(i); 
  int k;

  for(k=1;k<=size(v);k++)
  {
    if(v[k]!=0)
    {
      pa=pa+"var("+string(k)+"),";
    }
    else
    {
      va=va+"var("+string(k)+"),";
    }
  }
  
  pa=pa[1..size(pa)-1];
  va=va[1..size(va)-1];
  
  string newring
  ="ring nring=("+charstr(basering)+"),("+va+","+pa+"),lp;";  
  string quotring
  ="ring quring=("+charstr(basering)+","+pa+"),("+va+"),lp;";
  return(newring,quotring);
}
proc prepareSat(ideal i)
{
   int k;
   poly p=leadcoef(i[1]);
   for(k=2;k<=size(i);k++)
   {
     p=p*leadcoef(i[k]);
   }
   return(p);
}
proc squarefree(poly f, int i)
{
   poly h=gcd(f,diff(f,var(i)));
   poly g=lift(h,f)[1][1];    //  f/h
   return(g);
}

//=================================================================

proc radical(ideal i)
"USAGE:  radical(i); i ideal
RETURN:  ideal = the radical of i
NOTE:    algorithm of Krick/Logar
"  
{
   def  BAS = basering;
   ideal j;
   int n=nvars(BAS);
   int k;

   option(redSB);   
   ideal SBi=std(i);
   option(noredSB);
   int d=dim(SBi);

//-----the trivial cases
   if ((d==-1)||(n==d)||(n==1))
   {
      return(ideal(squarefree(SBi[1],1)));
   }
//-----the zero-dimensional case
   if (d==0)
   {
      j=finduni(SBi);
      for(k=1;k<=size(j);k++)
      {
         i=i,squarefree(cleardenom(j[k]),k);
      }
      return(std(i));
   }
//-----prepare the quotientring with respect to a maximal 
//     independent set
   list quotring=prepareQuotientring(SBi);
   execute (quotring[1]);
//-----we use this ring to compute a standardbasis of 
//     i*quring which is in i
   ideal i=std(imap(BAS,i));
//-----pass to the quotientring with respect to a maximal 
//     independent set
   execute( quotring[2]);
   ideal i=imap(nring,i);
   kill nring;
//-----computation of the zerodimensional radical
   ideal ra=radical(i);
//-----preparation for saturation
   poly p=prepareSat(ra);
   poly q=prepareSat(i);
//-----back to the original ring
   setring BAS;
   poly p=imap(quring,p);
   poly q=imap(quring,q);
   ideal ra=imap(quring,ra);
   kill quring;
//-----compute the intersection of ra with BAS
   ra=sat(ra,p)[1];
//----now we have radical(i)=intersection(ra,radical((i,q)))
   return(intersect(ra,radical(i+q)));
}

   ring  r = 0,(x,y,z),dp;
   ideal i = 
   intersect(ideal(x,y,z)^3,ideal(x-y-z)^2,ideal(x-y,x-z)^2);
   ideal pr= radical(i);
   pr;