LIB "tst.lib";
tst_init();

proc normalisation(ideal i)
"USAGE:  list L=normalisation(i);  i prime ideal
RETURN:  a list L of one ring L[1]=R; R contains the ideal
         norid such that R/norid is the normalisation of 
         basering/i.
NOTE:    to use the ring type def S=L[1];setring S;norid;
"
{
   def BAS=basering;
   list result;
   ideal rf;
   int   ds = -1;
   int isIso;

   if (typeof(attrib(i,"isIsolatedSingularity"))=="int")
   {
     if(attrib(i,"isIsolatedSingularity")==1){isIso=1;}
   }

   if(size(i)!=0)
   {
     list SM=mstd(i);
     i=SM[2];
     ideal SBi=SM[1];
   
     int n=nvars(BAS);
     int d=dim(SBi);
  //------------------- the singular locus ---------------
     if(isIso)
     {
       list singM=maxideal(1), maxideal(1),
       ds=0;
     }
     else
     {
       ideal sing=minor(jacob(i),n-d)+i;
       list singM=mstd(sing);
       ds=dim(singM[1]);
      }
     if(ds!=-1)
     {
  //----------------- computation of the radical ---------
      
     if (isIso)
     {
       ideal J=maxideal(1);
     }
     else
     {
       ideal J=radical(singM[2]);
     }
  //------------------ go to quotient ring ---------------
       qring R=SBi;
       ideal J=fetch(BAS,J);
       ideal i=fetch(BAS,i);
       poly p=J[1];
  //-------- computation of p*Hom(J,J) as R-ideal---------
       ideal f=quotient(p*J,J);
       ideal rf = interred(reduce(f,std(p)));
  // represents p*Hom(J,J)/p*R = Hom(J,J)/R 
     }  
   }  
  //-------- Test: Hom(J,J) == R ?, if yes, go home ------  
   if ( size(rf) == 0 )    
   {
      execute("ring newR="+charstr(basering)+",
         ("+varstr(basering)+"),("+ordstr(basering)+");");
      ideal norid=fetch(BAS,i);
      export norid;
      result=newR;
      setring BAS;
      return(result);
   }
//------------------- Case: Hom(J,J)!= R ------------------
// create new ring and map form old ring, the ring 
// newR/SBi+syzf will be isomorphic to Hom(J,J) as R-module

   f=p,rf; 
   //generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
   int q=size(f);
   module syzf=syz(f);
   
   ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
   map psi1 = BAS,maxideal(1);
   ideal SBi = psi1(SBi);
   attrib(SBi,"isSB",1);

   qring newRq = SBi;
   map psi = R,ideal(X(1..nvars(R)));
   ideal i = psi(i);
   ideal f = psi(f);
   module syzf = psi(syzf);

//------- computation of Hom(J,J) as ring ----------------
// determine kernel of: 
// R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. 
// This is of course the same as the kernel of
//  R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
// It is a fact, that the kernel is generated by the linear 
// and the quadratic relations

   ideal pf=f[1]*f;
   matrix T=matrix(ideal(T(1..q)),1,q);
   ideal Lin = ideal(T*syzf); //the linear relations
   
   int ii,jj;
   matrix A;
   ideal Q;

   for (ii=2; ii<=q; ii++ )
   {
      for ( jj=2; jj<=ii; jj++ )
      {
         A = lift(pf,f[ii]*f[jj]);
         Q = Q, ideal(T(jj)*T(ii) - T*A); 
         //the quadratic relations
      }
   }
   Q = Lin+Q;
   Q = subst(Q,T(1),1);
   Q = interred(reduce(Q,std(0)));

   ring newR = char(R),(X(1..nvars(R)),T(2..q)),dp;
   ideal k=imap(newRq,Q)+imap(newRq,i);
   if(isIso)
   {
     attrib(k,"isIsolatedSingularity",1);
   }   
   result=normalisation(k);
   setring BAS;
   return(result);
}

LIB"primdec.lib";
ring R=0, (x,y,z),dp;
ideal I=zy2-zx3-x6;
list nor=normalisation(I);
def S=nor[1];
setring S;
norid;

tst_status(1);$