//SINGULAR Procedure4.6.2

============================= we need ===================================
LIB"ring.lib";
proc primaryTest(ideal i, poly p)
"USAGE:   primaryTest(i,p); i standard basis with respect to
          lp, p irreducible polynomial in K[var(n)], 
          p^a=i[1] for some a;
ASSUME:   i is a zero-dimensional ideal.
RETURN:   an ideal, the radical of i if i is primary and in 
          general position with respect to lp, 
          the zero ideal else.
EXAMPLE:  primaryTest; shows an example
"  
{
   int m,e;
   int n=nvars(basering);
   poly t;
   ideal prm=p;

   for(m=2;m<=size(i);m++)
   {
     if(size(ideal(leadexp(i[m])))==1)
     {
       n--;
//----------------i[m] has a power of var(n) as leading term
       attrib(prm,"isSB",1);
//--- ?? i[m]=(c*var(n)+h)^e modulo prm for h 
//     in K[var(n+1),...], c in K ??
       e=deg(lead(i[m]));
       t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))
       /var(n)^(e-1);
       i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
//---if not (0) is returned, else c*var(n)+h is added to prm
       if (reduce(i[m]-t^e,prm,1) !=0)
       {
         return(ideal(0));
       }
       prm = prm,cleardenom(simplify(t,1));
     }
   }
   return(prm);
}
========================================================================

proc zeroDecomp(ideal i)
"USAGE:  zeroDecomp(i); i zero-dimensional ideal
RETURN:  list l of lists of two ideals such that the
         intersection(l[j][1], j=1..)=i, the l[i][1] are 
         primary and the l[i][2] their radicals
NOTE:    algorithm of Gianni/Trager/Zacharias
EXAMPLE: example zeroDecomp; shows an example
"  
{
   def  BAS = basering;
//----the ordering is changed to the lexicographical one
   changeord("R","lp");
   ideal i=fetch(BAS,i);
   int n=nvars(R);
   int k;
   list result,rest;
   ideal primary,prim;
   option(redSB);
   
//------the random coordinate change and its inverse
   ideal m=maxideal(1);
   m[n]=0;
   poly p=(random(100,1,n)*transpose(m))[1,1]+var(n);
   m[n]=p;
   map phi=R,m;
   m[n]=2*var(n)-p;
   map invphi=R,m;
   ideal j=groebner(phi(i));
//-------------factorization of the first element in i
   list fac=factorize(j[1],2);   
//-------------computation of the primaries and primes
   for(k=1;k<=size(fac[1]);k++)
   {
     p=fac[1][k]^fac[2][k];
     primary=groebner(j+p);
     prim=primaryTest(primary,fac[1][k]);
//---test whether all ideals were primary and in general 
//   position
     if(prim==0)
     {
       rest[size(rest)+1]=i+invphi(p);
     }
     else
     {
       result[size(result)+1]=
         list(std(i+invphi(p)),std(invphi(prim)));
     }
   }
//-------treat the bad cases collected in the rest again
   for(k=1;k<=size(rest);k++)
   {
     result=result+zeroDecomp(rest[k]);
   }
   option(noredSB);
   setring BAS;
   list result=imap(R,result);
   kill R;
   return(result);
}

   ring  r = 32003,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= zeroDecomp(i);
   pr;