///////////////////////////////////////////////////////////////////////////////
version="$Id: normal.lib,v 1.34.2.17 2002/10/21 14:44:41 lossen Exp $";
category="Commutative Algebra";
info="
LIBRARY:  normal.lib     Normalization of Affine Rings
AUTHORS:  G.-M. Greuel,  greuel@mathematik.uni-kl.de,
@*        G. Pfister,    pfister@mathematik.uni-kl.de

MAIN PROCEDURES:
 normal(I[,wd]);     computes the normalization of basering/I,
                     resp. computes the normalization of basering/I and
                     the delta invariant
 HomJJ(L);           presentation of End_R(J) as affine ring, L a list
 genus(I);           computes genus of the projective curve defined by I

AUXILIARY PROCEDURE: 
 deltaLoc(f,S);      (sum of) delta invariant(s) at conjugated singular points
";

LIB "general.lib";
LIB "poly.lib";
LIB "sing.lib";
LIB "primdec.lib";
LIB "elim.lib";
LIB "presolve.lib";
LIB "inout.lib";
LIB "ring.lib";
LIB "hnoether.lib";
///////////////////////////////////////////////////////////////////////////////

proc HomJJ (list Li)
"USAGE:   HomJJ (Li);  Li = list: ideal SBid, ideal id, ideal J, poly p
ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
@*       SBid = standard basis of id,
@*       J    = ideal of P containing the polynomial p,
@*       p    = nonzero divisor of R
COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
         affine ring, together with the canonical map R --> Hom_R(J,J),
         where R is the quotient ring of P modulo the standard basis SBid.
RETURN:  a list l of two objects
@format
         l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
               such that l[1]/endid = Hom_R(J,J) and
               endphi describes the canonical map R -> Hom_R(J,J)
         l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
         l[3] : an integer, the contribution to delta
@end format
NOTE:    printlevel >=1: display comments (default: printlevel=0)
EXAMPLE: example HomJJ;  shows an example
"
{
//---------- initialisation ---------------------------------------------------

   int isIso,isPr,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
   intvec rw,rw1;
   list L;
   y = printlevel-voice+2;  // y=printlevel (default: y=0)
   def P = basering;
   ideal SBid, id, J = Li[1], Li[2], Li[3];
   poly p = Li[4];
   attrib(SBid,"isSB",1);
   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not

//---- set attributes for special cases where algorithm can be simplified -----
   if( homo==1 )
   {
      rw = ringweights(P);
   }
   if( typeof(attrib(id,"isPrim"))=="int" )
   {
      if(attrib(id,"isPrim")==1)  { isPr=1; }
   }
   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
   {
      if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; }
   }
   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
   {
      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
   }
   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
   {
      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
   }
   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
   {
      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
   }
   if( typeof(attrib(id,"isEquidimensional"))=="int" )
   {
      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
   }
//-------------------------- go to quotient ring ------------------------------
   qring R  = SBid;
   ideal id = fetch(P,id);
   ideal J  = fetch(P,J);
   poly p   = fetch(P,p);
   ideal f,rf,f2;
   module syzf;
//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
   if ( y>=1 )
   {
     "// compute p*Hom(J,J) = p*J:J, p a non-zerodivisor";
     "//   p is equal to:"; "";
     p;
     "";
   }
   f  = quotient(p*J,J);

   if ( y>=1 )
   { "// the module p*Hom(J,J) = p*J:J, p a non-zerodivisor";
      "// p"; p;
      "// f=p*J:J";f;
   }
   f2 = std(p);
  
//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------

   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R

   if ( size(rf) == 0 )
   {
      if ( homog(f) && find(ordstr(basering),"s")==0 )
      {
         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
      }
      else
      {
         ring newR1 = char(P),(X(1..nvars(P))),dp;
      }
      ideal endphi = maxideal(1);
      ideal endid = fetch(P,id);
      L=substpart(endid,endphi,homo,rw);
      def lastRing=L[1];
      setring lastRing;

      attrib(endid,"onlySingularAtZero",oSAZ);
      attrib(endid,"isCohenMacaulay",isCo);
      attrib(endid,"isPrim",isPr);
      attrib(endid,"isIsolatedSingularity",isIso);
      attrib(endid,"isRegInCodim2",isRe);
      attrib(endid,"isEqudimensional",isEq);
      attrib(endid,"isCompleteIntersection",0);
      attrib(endid,"isRad",0);
      L=lastRing;
      L = insert(L,1,1);
      dbprint(y,"// case R = Hom(J,J)");
      if(y>=1)
      {
         "//   R=Hom(J,J)";
         "   ";
         lastRing;
         "   ";
         "//   the new ideal";
         endid;
         "   ";
         "//   the old ring";
         "   ";
         P;
         "   ";
         "//   the old ideal";
         "   ";
         setring P;
         id;
         "   ";
         setring lastRing;
         "//   the map";
         "   ";
         endphi;
         "   ";
         pause();
         newline;
      }
      setring P;
      L[3]=0;
      return(L);
   }
   if(y>=1)
   {
      "// R is not equal to Hom(J,J), we have to try again";
      pause();
      newline;
   }
//---------- Hom(J,J) != R: create new ring and map from old ring -------------
// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module

   f=mstd(f)[2];
   ideal ann=quotient(f2,f);
   int delt;
   if(isIso&&isEq){delt=vdim(std(modulo(f,ideal(p))));}

   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
   q = size(f);

   syzf = syz(f);

   if ( homo==1 )
   {
      rw1 = rw,0;
      for ( ii=2; ii<=q; ii++ )
      {
         rw  = rw, deg(f[ii])-deg(f[1]);
         rw1 = rw1, deg(f[ii])-deg(f[1]);
      }
      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
   }
   else
   {
      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
   }

   map psi1 = P,maxideal(1);
   ideal SBid = psi1(SBid);
   attrib(SBid,"isSB",1);

   qring newR = std(SBid);

   map psi = R,ideal(X(1..nvars(R)));
   ideal id = psi(id);
   ideal f = psi(f);
   module syzf = psi(syzf);
   ideal pf,Lin,Quad,Q;
   matrix T,A;
   list L1;

//---------- computation of Hom(J,J) as affine 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

   pf = f[1]*f;
   T = matrix(ideal(T(1..q)),1,q);
   Lin = ideal(T*syzf);
   if(y>=1)
   {
      "// the ring structure of Hom(J,J) as R-algebra";
      " ";
      "//   the linear relations";
      " ";
      Lin;
      "   ";
   }

   for (ii=2; ii<=q; ii++ )
   {
      for ( jj=2; jj<=ii; jj++ )
      {
         A = lift(pf,f[ii]*f[jj]);
         Quad = Quad, ideal(T(jj)*T(ii) - T*A);          // quadratic relations
      }
   }

   if(y>=1)
   {
      "//   the quadratic relations";
      "   ";
      interred(Quad);
      pause();
      newline;
   }
   Q = Lin,Quad;
   Q = subst(Q,T(1),1);
   Q=mstd(Q)[2];
//---------- reduce number of variables by substitution, if possible ----------
   if (homo==1)
   {
      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
   }
   else
   {
      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
   }

   ideal endid  = imap(newR,id),imap(newR,Q);
   ideal endphi = ideal(X(1..nvars(R)));

   L=substpart(endid,endphi,homo,rw);

   def lastRing=L[1];
   setring lastRing;

   attrib(endid,"onlySingularAtZero",0);
   map sigma=R,endphi;
   ideal an=sigma(ann);
   export(an);  //noetig?
   ideal te=an,endid;
   if(isIso&&(size(reduce(te,std(maxideal(1))))==0))
   {
      attrib(endid,"onlySingularAtZero",oSAZ);
   }
   kill te;
   attrib(endid,"isCohenMacaulay",isCo);
   attrib(endid,"isPrim",isPr);
   attrib(endid,"isIsolatedSingularity",isIso);
   attrib(endid,"isRegInCodim2",isRe);
   attrib(endid,"isEquidimensional",isEq);
   attrib(endid,"isCompleteIntersection",0);
   attrib(endid,"isRad",0);
   if(y>=1)
   {
      "//   the new ring after reduction of the number of variables";
      "   ";
      lastRing;
      "   ";
      "//   the new ideal";
      "   ";
      endid;
      "   ";
      "//   the old ring";
      "   ";
      P;
      "   ";
      "//   the old ideal";
      "   ";
      setring P;
      id;
      "   ";
      setring lastRing;
      "//   the map";
      "   ";
      endphi;
      "   ";
      pause();
      newline;
   }
   L = lastRing;
   L = insert(L,0,1);
   L[3]=delt;
   return(L);
}
example
{"EXAMPLE:";  echo = 2;
  ring r   = 0,(x,y),wp(2,3);
  ideal id = y^2-x^3;
  ideal J  = x,y;
  poly p   = x;
  list Li = std(id),id,J,p;
  list L   = HomJJ(Li);
  def end = L[1];    // defines ring L[1], containing ideals endid, endphi
  setring end;       // makes end the basering
  end;
  endid;             // end/endid is isomorphic to End(r/id) as ring
  map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
  psi;
}

///////////////////////////////////////////////////////////////////////////////

proc normal(ideal id, list #)
"USAGE:   normal(i [,choose]);  i a radical ideal, choose empty, 1 or \"wd\"
         if choose=1 the normalization of the associated primes is computed
         (which is sometimes more efficient);
         if @code{choose=\"wd\"} the delta invariant is computed
         simultaneously; this may take much more time in the reducible case,
         since the factorizing standard basis algorithm cannot be used.
ASSUME:  The ideal must be radical, for non-radical ideals the output may
         be wrong (i=radical(i); makes i radical)
RETURN:  a list of rings, say nor and in case of @code{choose=\"wd\"} an 
         integer at the end of the list. 
         Each ring @code{nor[i]} contains two ideals with given names 
         @code{norid} and @code{normap} such that@*
         - the direct sum of the rings @code{nor[i]/norid} is the 
         normalization of basering/id;@*
         - @code{normap} gives the normalization map from basering/id to 
         @code{nor[i]/norid} (for each i).
NOTE:    to use the i-th ring type: @code{def R=nor[i]; setring R;}.
@*       Increasing printlevel displays more comments (default: printlevel=0).
@*       Not implemented for local or mixed orderings.
@*       If the input ideal i is weighted homogeneous a weighted ordering may
         be used (qhweight(i); computes weights).
KEYWORDS: normalization; delta invariant.
EXAMPLE: example normal; shows an example
"
{
   int i,j,y,withdelta;
   string sr;
   list result,prim,keepresult;
   y = printlevel-voice+2;
   if(size(#)>0)
   {
     if(typeof(#[1])=="string")
     {
       kill #;
       list #;
       withdelta=1;
     }
   }
   attrib(id,"isRadical",1);
   if ( ord_test(basering) != 1)
   {
     "";
     "// Not implemented for this ordering,";
     "// please change to global ordering!";
     return(result);
   }
   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
   {
      if(attrib(id,"isCompleteIntersection")==1)
      {
         attrib(id,"isCohenMacaulay",1);
         attrib(id,"isEquidimensional",1);
      }
   }
   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
   {
      if(attrib(id,"isCohenMacaulay")==1)
      {
         attrib(id,"isEquidimensional",1);
      }
   }
   if( typeof(attrib(id,"isPrim"))=="int" )
   {
      if(attrib(id,"isPrim")==1)
      {
         attrib(id,"isEquidimensional",1);
      }
   }
   if(size(#)==0)
   {
      if( typeof(attrib(id,"isEquidimensional"))=="int" )
      {
         if(attrib(id,"isEquidimensional")==1)
         {
            prim[1]=id;
         }
         else
         {
            prim=equidim(id);
         }
      }
      else
      {
         prim=equidim(id);
      }
      if(y>=1)
      {
         "// we have ",size(prim),"equidimensional components";
      }
      if(withdelta &&(size(prim)>1))
      {
         withdelta=0;
         "WARNING!  cannot compute delta,";
         "the ideal is not equidimensional";
      }
      if(!withdelta)
      {
         option(redSB);
         for(j=1;j<=size(prim);j++)
         {
            keepresult=keepresult+facstd(prim[j]);
         }
         prim=keepresult; 
         option(noredSB); 
      }
   }
   else
   {
      if( typeof(attrib(id,"isPrim"))=="int" )
      {
         if(attrib(id,"isPrim")==1)
         {
            prim[1]=id;
         }
         else
         {
            prim=minAssPrimes(id);
         }
      }
      else
      {
         prim=minAssPrimes(id);
      }
      if(y>=1)
      {
         "// we have ",size(prim),"irreducible components";
      }
   }
   for(i=1; i<=size(prim); i++)
   {
      if(y>=1)
      {
         "// we are in loop ",i;
      }
      attrib(prim[i],"isCohenMacaulay",0);
      if(size(#)!=0)
      {
         attrib(prim[i],"isPrim",1);
      }
      else
      {
         attrib(prim[i],"isPrim",0);
      }
      attrib(prim[i],"isRegInCodim2",0);
      attrib(prim[i],"isIsolatedSingularity",0);
      attrib(prim[i],"isEquidimensional",1);
      attrib(prim[i],"isCompleteIntersection",0);
      attrib(prim[i],"onlySingularAtZero",0);
      if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
      {
            if(attrib(id,"onlySingularAtZero")==1)
             {attrib(prim[i],"onlySingularAtZero",1); }
      }

      if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
      {
            if(attrib(id,"isIsolatedSingularity")==1)
             {attrib(prim[i],"isIsolatedSingularity",1); }
      }

      if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
      {
            if((attrib(id,"isIsolatedSingularity")==1)&&(size(#)==0))
             {attrib(prim[i],"isIsolatedSingularity",1); }
      }
      keepresult=normalizationPrimes(prim[i],maxideal(1),0);

      for(j=1;j<=size(keepresult)-1;j++)
      {
         result=insert(result,keepresult[j]);
      }
      sr = string(size(result));      
   }
   if(withdelta)
   {
      sr = string(size(keepresult)-1);
      result=keepresult;
   }
      dbprint(y+1,"
// 'normal' created a list of "+sr+" ring(s).
// nor["+sr+"+1] is the delta-invariant in case of choose=wd.
// To see the rings, type (if the name of your list is nor):
     show( nor);
// To access the 1-st ring and map (similar for the others), type:
     def R = nor[1]; setring R;  norid; normap;
// R/norid is the 1-st ring of the normalization and
// normap the map from the original basering to R/norid");

      //kill endphi,endid;
      return(result);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),wp(2,1,2);
   ideal i=z3-xy4;
   list nor=normal(i);
   show(nor);
   def r1=nor[1];
   setring r1;
   norid;
   normap;

   ring s=0,(x,y),dp;
   ideal i=(x-y^2)^2 - y*x^3;
   nor=normal(i,"wd");
   //the delta-invariant
   nor[size(nor)];
}

///////////////////////////////////////////////////////////////////////////////
static proc normalizationPrimes(ideal i,ideal ihp,int delt, list #)
"USAGE:   normalizationPrimes(i,ihp[,si]);  i equidimensional ideal, ihp map
         (partial normalization),delta partial delta-invariant,
          # ideal of singular locus
RETURN:   a list of rings, say nor, and an integer, the delta-invariant
          at the end of the list.
          each ring nor[i] contains two ideals
          with given names norid and normap such that
           - the direct sum of the rings nor[i]/norid is
             the normalization of basering/id;
           - normap gives the normalization map from basering/id
             to nor[i]/norid (for each i)
EXAMPLE: example normalizationPrimes; shows an example
"
{

   int y = printlevel-voice+2;  // y=printlevel (default: y=0)

   if(y>=1)
   {
     "";
     "// START a normalization loop with the ideal";  "";
     i;  "";
     basering;  "";
     pause();
     newline;
   }

   def BAS=basering;
   list result,keepresult1,keepresult2;
   ideal J,SB,MB;
   int depth,lauf,prdim;
   int ti=timer;

   if(size(i)==0)
   {
      if(y>=1)
      {
          "// the ideal was the zero-ideal";
      }
         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");
         ideal norid=ideal(0);
         ideal normap=fetch(BAS,ihp);
         export norid;
         export normap;
         result=newR7;
         result[size(result)+1]=delt;
         setring BAS;
         return(result);
   }

   if(y>=1)
   {
     "// SB-computation of the input ideal";
   }
   list SM=mstd(i);                //here the work starts
   int dimSM =  dim(SM[1]);        //dimension of variety to normalize
  // Case: Get an ideal containing a unit
   if( dimSM == -1)
   {  "";
      "      // A unit ideal was found.";
      "      // Stop with partial result computed so far";"";

         MB=SM[2];
         intvec rw;
         list LL=substpart(MB,ihp,0,rw);
         def newR6=LL[1];
         setring newR6;
         ideal norid=endid;
         ideal normap=endphi;
         kill endid,endphi;
         export norid;
         export normap;
         result=newR6;
         result[size(result)+1]=delt;
         setring BAS;
         return(result);
   }

   if(y>=1)
   {
      "//   the dimension is:"; "";
      dimSM;"";
   }
   if(size(#)>0)
   {
      if(attrib(i,"onlySingularAtZero"))
      {
         list JM=maxideal(1),maxideal(1);
         attrib(JM[1],"isSB",1);
         attrib(JM[2],"isRad",1);
      }
      else
      {
         ideal te=#[1],SM[2];
         list JM=mstd(te);
         kill te;
         if( typeof(attrib(#[1],"isRad"))!="int" )
         {
            attrib(JM[2],"isRad",0);
         }
      }
   }

   if(attrib(i,"isPrim")==1)
   {
      attrib(SM[2],"isPrim",1);
   }
   else
   {
      attrib(SM[2],"isPrim",0);
   }
   if(attrib(i,"isIsolatedSingularity")==1)
   {
      attrib(SM[2],"isIsolatedSingularity",1);
   }
   else
   {
      attrib(SM[2],"isIsolatedSingularity",0);
   }
   if(attrib(i,"isCohenMacaulay")==1)
   {
      attrib(SM[2],"isCohenMacaulay",1);
   }
   else
   {
      attrib(SM[2],"isCohenMacaulay",0);
   }
   if(attrib(i,"isRegInCodim2")==1)
   {
      attrib(SM[2],"isRegInCodim2",1);
   }
   else
   {
      attrib(SM[2],"isRegInCodim2",0);
   }
   if(attrib(i,"isEquidimensional")==1)
   {
      attrib(SM[2],"isEquidimensional",1);
   }
   else
   {
      attrib(SM[2],"isEquidimensional",0);
   }
   if(attrib(i,"isCompleteIntersection")==1)
   {
      attrib(SM[2],"isCompleteIntersection",1);
   }
   else
   {
      attrib(SM[2],"isCompleteIntersection",0);
   }
   if(attrib(i,"onlySingularAtZero")==1)
   {
      attrib(SM[2],"onlySingularAtZero",1);
   }
   else
   {
      attrib(SM[2],"onlySingularAtZero",0);
   }
   if((attrib(SM[2],"isIsolatedSingularity")==1)&&(homog(SM[2])==1))
   {
      attrib(SM[2],"onlySingularAtZero",1);
   }

   //the smooth case
   if(size(#)>0)
   {
      if(dim(JM[1])==-1)
      {
         if(y>=1)
         {
            "// the ideal was smooth";
         }
         MB=SM[2];
         intvec rw;;
         list LL=substpart(MB,ihp,0,rw);
         def newR6=LL[1];
         setring newR6;
         ideal norid=endid;
         ideal normap=endphi;
         kill endid,endphi;
         export norid;
         export normap;
         result=newR6;
         result[size(result)+1]=delt;
         setring BAS;
         return(result);
     }
   }

   //the zero-dimensional case
   if((dim(SM[1])==0)&&(homog(SM[2])==1))
   {
      if(y>=1)
      {
         "// the ideal was zero-dimensional and homogeneous";
      }
      MB=maxideal(1);
      intvec rw;
      list LL=substpart(MB,ihp,0,rw);
      def newR5=LL[1];
      setring newR5;
      ideal norid=endid;
      ideal normap=endphi;
      kill endid,endphi;
      export norid;
      export normap;
      result=newR5;
      result[size(result)+1]=delt;
      setring BAS;
      return(result);
   }

   //the one-dimensional case
   //in this case it is a line because
   //it is irreducible and homogeneous
   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
        &&(homog(SM[2])==1))
   {
      if(y>=1)
      {
         "// the ideal defines a line";
      }
      MB=SM[2];
      intvec rw;
      list LL=substpart(MB,ihp,0,rw);
      def newR4=LL[1];
      setring newR4;
      ideal norid=endid;
      ideal normap=endphi;
      kill endid,endphi;
      export norid;
      export normap;
      result=newR4;
      result[size(result)+1]=delt;
      setring BAS;
      return(result);
   }

   //the higher dimensional case
   //we test first of all CohenMacaulay and
   //complete intersection
   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
   {
      //test for complete intersection
      attrib(SM[2],"isCohenMacaulay",1);
      attrib(SM[2],"isCompleteIntersection",1);
      attrib(SM[2],"isEquidimensional",1);
      if(y>=1)
      {
         "// the ideal is a complete intersection";
      }
   }
   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(i)==1))
   {
      int tau=vdim(std(i+jacob(i)));
      if(tau>=0)
      {
        execute("ring BASL="+charstr(BAS)+",("+varstr(BAS)+"),ds;");
        ideal i=imap(BAS,i);
        int tauloc=vdim(std(i+jacob(i)));
        setring BAS;
        attrib(SM[2],"onlySingularAtZero",(tau==tauloc));
        kill BASL;
      }
   }

   //compute the singular locus
   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(#)==0))
   {
      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
      if(y >=1 )
      {
         "// SB of singular locus will be computed";
      }
      ideal sin=J,SM[1];
      list JM=mstd(sin);
      attrib(JM[1],"isSB",1);

      //JM[1] SB of singular locus, JM[2]=minbasis of singular locus
      //SM[1] SB of the input ideal, SM[2] minbasis
      if(y>=1)
      {
         "//   the dimension of the singular locus is:";"";
         dim(JM[1]); "";
      }

      if(dim(JM[1])==-1)
      {
         if(y>=1)
         {
            "// the ideal is smooth";
         }
         MB=SM[2];
         intvec rw;
         list LL=substpart(MB,ihp,0,rw);
         def newR3=LL[1];
         setring newR3;
         ideal norid=endid;
         ideal normap=endphi;
         kill endid,endphi;
         export norid;
         export normap;
         result=newR3;
         result[size(result)+1]=delt;
         setring BAS;
         return(result);
      }
      if(dim(JM[1])==0)
      {
         attrib(SM[2],"isIsolatedSingularity",1);
         if(homog(SM[2])){attrib(SM[2],"onlySingularAtZero",1);}
      }
      if(dim(JM[1])<=dim(SM[1])-2)
      {
         attrib(SM[2],"isRegInCodim2",1);
      }
   }
   else
   {
     if(size(#)==0)
     {
        list JM=maxideal(1),maxideal(1);

        attrib(JM[1],"isSB",1);
        if(dim(SM[1])>=2){attrib(SM[2],"isRegInCodim2",1);}
     }
   }
   if((attrib(SM[2],"isRegInCodim2")==1)&&(attrib(SM[2],"isCohenMacaulay")==1))
   {
      if(y>=1)
      {
            "// the ideal was CohenMacaulay and regular in codimension 2";
      }
      MB=SM[2];
      intvec rw;
      list LL=substpart(MB,ihp,0,rw);
      def newR6=LL[1];
      setring newR6;
      ideal norid=endid;
      ideal normap=endphi;
      kill endid,endphi;
      export norid;
      export normap;
      result=newR6;
      result[size(result)+1]=delt;
      setring BAS;
      return(result);
   }

   //if it is an isolated singularity only at 0 things are easier
   //JM ideal of singular locus, SM ideal of variety
   if(attrib(SM[2],"onlySingularAtZero"))          
   {
      attrib(SM[2],"isIsolatedSingularity",1);
      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
      //vars not contained in ideal
      ideal Ann=quotient(SM[2],SL[1]);
      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
      //qAnn=0 ==> the first var(=SL[1]) not contained in SM is a nzd of R/SM
      if(size(qAnn)==0)
      {
         if(y>=1)
         {
            "";
            "//   the ideal rad(J):";
            "";
            maxideal(1);
            newline;
         }
         //again test for normality
         //Hom(I,R)=R
         list RR;
         RR=SM[1],SM[2],maxideal(1),SL[1];

         RR=HomJJ(RR,y);
         if(RR[2]==0)
         {
            def newR=RR[1];
            setring newR;
            map psi=BAS,endphi;
            list tluser=normalizationPrimes(endid,psi(ihp),delt+RR[3],an);
            setring BAS;
            return(tluser);
         }
         MB=SM[2];
         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");
         ideal norid=fetch(BAS,MB);
         ideal normap=fetch(BAS,ihp);
         delt=delt+RR[3];
         export norid;
         export normap;
         result=newR7;
         result[size(result)+1]=delt;
         setring BAS;
         return(result);

       }
      //Now the case where qAnn!=0, i.e.SL[1] is a zero divisor of R/SM
      //and we have found a splitting: id and id1
      //id=Ann defines components of R/SM in the complement of V(SL[1])
      //id1 defines components of R/SM in the complement of V(id)
       else
       {
          ideal id=Ann;
          attrib(id,"isCohenMacaulay",0);
          attrib(id,"isPrim",0);
          attrib(id,"isIsolatedSingularity",1);
          attrib(id,"isRegInCodim2",0);
          attrib(id,"isCompleteIntersection",0);
          attrib(id,"isEquidimensional",0);
          attrib(id,"onlySingularAtZero",1);

          ideal id1=quotient(SM[2],Ann);
          //ideal id=SL[1],SM[2];
          attrib(id1,"isCohenMacaulay",0);
          attrib(id1,"isPrim",0);
          attrib(id1,"isIsolatedSingularity",1);
          attrib(id1,"isRegInCodim2",0);
          attrib(id1,"isCompleteIntersection",0);
          attrib(id1,"isEquidimensional",0);
          attrib(id1,"onlySingularAtZero",1);
          
          ideal t1=id,id1;
          int mul=vdim(std(t1));
          kill t1;

          keepresult1=normalizationPrimes(id,ihp,0);

          keepresult2=normalizationPrimes(id1,ihp,0);
            
          delt=delt+mul+keepresult1[size(keepresult1)]
                         +keepresult1[size(keepresult1)];

          for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
          {
             keepresult1=insert(keepresult1,keepresult2[lauf]);
          }
          keepresult1[size(keepresult1)]=delt;
          return(keepresult1);
       }
   }

   //test for non-normality
   //Hom(I,I)<>R
   //we can use Hom(I,I) to continue
   ideal SL=simplify(reduce(JM[2],SM[1]),2);
   ideal Ann=quotient(SM[2],SL[1]);
   ideal qAnn=simplify(reduce(Ann,SM[1]),2);

   if(size(qAnn)==0)
   {
      list RR;
      list RS;
      //now we have to compute the radical
      if(y>=1)
      {
         "// radical computation of singular locus";
      }

      J=radical(JM[2]);   //the singular locus
      JM=mstd(J);
      
      if((vdim(JM[1])==1)&&(size(reduce(J,std(maxideal(1))))==0))
      {
         attrib(SM[2],"onlySingularAtZero",1);
      }
      if(y>=1)
      {
        "//   radical is equal to:";"";
        J;
        "";
      }

      if(deg(SL[1])>deg(J[1]))
      {
         Ann=quotient(SM[2],J[1]);
         qAnn=simplify(reduce(Ann,SM[1]),2);
         if(size(qAnn)==0){SL[1]=J[1];}
      }

      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
      RR=SM[1],SM[2],JM[2],SL[1];

      //   evtl eine geeignete Potenz von JM?
     if(y>=1)
     {
        "// compute Hom(rad(J),rad(J))";
     }

     RS=HomJJ(RR,y);

      if(RS[2]==1)
      {
         def lastR=RS[1];
         setring lastR;
         map psi1=BAS,endphi;
         ideal norid=endid;
         ideal normap=psi1(ihp);
         kill endid,endphi;
         export norid;
         export normap;
         result=lastR;
         result[size(result)+1]=delt+RS[3];
         setring BAS;
         return(result);
      }
      int n=nvars(basering);
      ideal MJ=JM[2];

      def newR=RS[1];
      setring newR;

      map psi=BAS,endphi;
      list tluser=
             normalizationPrimes(endid,psi(ihp),delt+RS[3],psi(MJ));
      setring BAS;
      return(tluser);
   }
    // A component with singular locus the whole component found
   if( Ann == 1)
   {
      "// Input appeared not to be a radical ideal!";
      "// A (everywhere singular) component with ideal";
      "// equal to its Jacobian ideal was found";
      "// Procedure will stop with partial result computed so far";"";

         MB=SM[2];
         intvec rw;
         list LL=substpart(MB,ihp,0,rw);
         def newR6=LL[1];
         setring newR6;
         ideal norid=endid;
         ideal normap=endphi;
         kill endid,endphi;
         export norid;
         export normap;
         result=newR6;
         result[size(result)+1]=delt;
         setring BAS;
         return(result);
   }
   else
   {
      int equi=attrib(SM[2],"isEquidimensional");
      int oSAZ=attrib(SM[2],"onlySingularAtZero");
      int isIs=attrib(SM[2],"isIsolatedSingularity");

      ideal new1=Ann;
      ideal new2=quotient(SM[2],Ann);
      //ideal new2=SL[1],SM[2];

      int mul;
      if(equi&&isIs)
      {
          ideal t2=new1,new2;
          mul=vdim(std(t2));
      }
      execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");
      if(y>=1)
      {
         "// zero-divisor found";
      }
      ideal vid=fetch(BAS,new1);
      ideal ihp=fetch(BAS,ihp);
      attrib(vid,"isCohenMacaulay",0);
      attrib(vid,"isPrim",0);
      attrib(vid,"isIsolatedSingularity",isIs);
      attrib(vid,"isRegInCodim2",0);
      attrib(vid,"onlySingularAtZero",oSAZ);
      attrib(vid,"isEquidimensional",equi);      
      attrib(vid,"isCompleteIntersection",0);

      keepresult1=normalizationPrimes(vid,ihp,0);
      int delta1=keepresult1[size(keepresult1)];
      setring BAS;

      execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");

      ideal vid=fetch(BAS,new2);
      ideal ihp=fetch(BAS,ihp);
      attrib(vid,"isCohenMacaulay",0);
      attrib(vid,"isPrim",0);
      attrib(vid,"isIsolatedSingularity",isIs);
      attrib(vid,"isRegInCodim2",0);
      attrib(vid,"isEquidimensional",equi);     
      attrib(vid,"isCompleteIntersection",0);
      attrib(vid,"onlySingularAtZero",oSAZ);

      keepresult2=normalizationPrimes(vid,ihp,0);
      int delta2=keepresult2[size(keepresult2)];

      setring BAS;
      
      for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
      {
         keepresult1=insert(keepresult1,keepresult2[lauf]);
      }
      keepresult1[size(keepresult1)]=delt+mul+delta1+delta2;
      return(keepresult1);
   }
}
example
{ "EXAMPLE:";echo = 2;
   // Huneke
   ring qr=31991,(a,b,c,d,e),dp;
   ideal i=
   5abcde-a5-b5-c5-d5-e5,
   ab3c+bc3d+a3be+cd3e+ade3,
   a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
   abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
   ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
   a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
   a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;

   list pr=normalizationPrimes(i);
   def r1=pr[1];
   setring r1;
   norid;
   normap;
}
///////////////////////////////////////////////////////////////////////////////
static proc substpart(ideal endid, ideal endphi, int homo, intvec rw)

"//Repeated application of elimpart to endid, until no variables can be
//directy substituded. homo=1 if input is homogeneous, rw contains
//original weights, endphi (partial) normalization map";

{
   def newRing=basering;
   int ii,jj;
   map phi = basering,maxideal(1);

   list Le = elimpart(endid);
                                     //this proc and the next loop try to
   int q = size(Le[2]);              //substitute as many variables as possible
   intvec rw1 = 0;                   //indices of substituted variables
   rw1[nvars(basering)] = 0;
   rw1 = rw1+1;

   while( size(Le[2]) != 0 )
   {
      endid = Le[1];
      map ps = newRing,Le[5];

      phi = ps(phi);
      for(ii=1;ii<=size(Le[2])-1;ii++)
      {
         phi=phi(phi);
      }
      //eingefuegt wegen x2-y2z2+z3
      kill ps;

      for( ii=1; ii<=size(rw1); ii++ )
      {
         if( Le[4][ii]==0 )
         {
            rw1[ii]=0;                             //look for substituted vars
         }
      }
      Le=elimpart(endid);
      q = q + size(Le[2]);
   }
   endphi = phi(endphi);

//---------- return -----------------------------------------------------------
// in the homogeneous case put weights for the remaining vars correctly, i.e.
// delete from rw those weights for which the corresponding entry of rw1 is 0

   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
   {
      jj=1;
      for( ii=2; ii<size(rw1); ii++)
      {
         jj++;
         if( rw1[ii]==0 )
         {
            rw=rw[1..jj-1],rw[jj+1..size(rw)];
            jj=jj-1;
         }
      }
      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }

      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),(a(rw),dp);
   }
   else
   {
      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
   }

   ideal lastmap;
   q = 1;
   for(ii=1; ii<=size(rw1); ii++ )
   {
      if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; }
      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
   }
   map phi1 = newRing,lastmap;
   ideal endid  = phi1(endid);
   ideal endphi = phi1(endphi);
   export(endid);
   export(endphi);
   list L = lastRing;
   setring newRing;
   return(L);
}
/////////////////////////////////////////////////////////////////////////////

proc genus(ideal I,list #)
"USAGE:   genus(I) or genus(i,1); I a 1-dimensional ideal
RETURN:  an integer, the geometric genus p_g = p_a - delta of the projective 
         curve defined by I, where p_a is the arithmetic genus.
NOTE:    delta is the sum of all local delta-invariants of the singularities,
         i.e. dim(R'/R), R' the normalization of the local ring R of the
         singularity.
         genus(i,1) uses the normalization to compute delta. Usually this
         is slow but sometimes not.
EXAMPLE: example genus; shows an example
"
{
   int w = printlevel-voice+2;  // w=printlevel (default: w=0)

   def R=basering;

   if(homog(I))
   {
      execute("ring newR=("+charstr(R)+"),("+varstr(R)+"),dp;"); 
      ideal I=imap(R,I);  
   }
   else
   {
      execute("ring newR=("+charstr(R)+"),("+varstr(R)+",@t),dp;"); 
      ideal I=imap(R,I);  
      I=homog(std(I),@t);
   }
   ideal J=std(I);
   if((dim(J)!=2)||((nvars(basering)==2)&&(dim(J)!=1)))
   {
      ERROR("This is not a curve");
   }
   if((nvars(basering)<=3)&&(size(J)>1))
   {
       ERROR("This is not equidimensional");
   } 
   
   intvec hp=hilbPoly(J);
   int p_a=1-hp[1];
   int d=hp[2];

   if(w>=1)
   {
      "";"The ideal of the projective curve:";"";J;"";
      "The coefficients of the Hilbert polynomial";hp;
      "arithmetic genus:";p_a;
      "degree:";d;"";
   }
   
   intvec v = hilb(J,1);
   int o,i;
   
   if(nvars(basering)>3)
   {
      map phi=newR,maxideal(1);
      int de;
      ideal K,L;
      poly m=var(4);
      for(i=5;i<=nvars(basering);i++){m=m*var(i);}
      K=eliminate(J,m,v);
      if(size(K)==1){de=deg(K[1]);}
      m=var(1);
      for(i=2;i<=nvars(basering)-3;i++){m=m*var(i);}
      i=0;
      while(d!=de)
      {
         o=1;
         i++;
         K=phi(J);
         K=eliminate(K,m,v);
         if(size(K)==1){de=deg(K[1]);}
         if(i==5)
         {
            K=reduce(equidimMax(J),J);
            if(size(K)!=0){ERROR("This is not equidimensional");}
          }
          if(i==10){ERROR("did not find a good projection");}
          L=sparsetriag(nvars(newR),nvars(newR),80-5*i,i)*transpose(maxideal(1));
          phi=newR,L;
      }
      J=K;
   }
   poly p=J[1];
   
   if(nvars(basering)==2)
   {
      if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");}
      return(-deg(p)+1);
   }
   
   execute("ring S=("+charstr(R)+"),(x,y,t),dp;");
   ideal L=maxideal(1);
   execute("ring C=("+charstr(R)+"),(x,y),ds;");
   ideal I;
   execute("ring A=("+charstr(R)+"),(x,t),dp;");
   map phi=S,1,x,t;
   map psi=S,x,1,t;
   poly g,h;
   ideal I,I1;
   execute("ring B=("+charstr(R)+"),(x,t),ds;");
   
   setring S;
   if(o)
   {
     for(i=1;i<=nvars(newR)-3;i++){L[i]=0;}
     L=L,maxideal(1);
   }   
   map sigma=newR,L;
   poly F=sigma(p);
   if(w>=1){"the projected curve:";"";F;"";}
   
   kill newR;
   
   int genus=(d-1)*(d-2)/2;
   if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();}
 
   int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing,
       tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf;
   list inv;

   if(w>=1)
     {"";"analyse the singularities at oo";"";"singular locus at (1,x,0):";"";}
   setring A;        
   g=phi(F);
   h=psi(F);
   I=g,jacob(g),var(2);
   I=std(I);
   if(deg(I[1])>0)
   {
      list qr=minAssGTZ(I);
      if(w>=1){qr;"";}
      
      for(k=1;k<=size(qr);k++)
      {
         if(w>=1){ nsing=nsing+vdim(std(qr[k]));}
         inv=deltaLoc(g,qr[k]);
         deltainf=deltainf+inv[1];
         tauinf=tauinf+inv[2];
         l=vdim(std(qr[k]));
         if(inv[2]==l){nodesinf=nodesinf+l;}
         if(inv[2]==2*l){cuspsinf=cuspsinf+l;}         
         nbranchinf=nbranchinf+inv[3];
      }
   }
   else
   {
     if(w>=1){"            the curve is smooth at (1,x,0)";"";}
   }
   if(w>=1){"singular locus at (0,1,0):";"";}
   inv=deltaLoc(h,maxideal(1)); 
   if((w>=1)&&(inv[2]!=0)){ nsing++;}
   deltainf=deltainf+inv[1];
   tauinf=tauinf+inv[2];
   if(inv[2]==1){nodesinf=nodeainf++;}
   if(inv[2]==2){cuspsinf=cuspsinf++;}         
   
   if((w>=1)&&(inv[2]==0)){"            the curve is smooth at (0,1,0)";"";}
   if(inv[2]>0){nbranchinf=nbranchinf+inv[3];}
  
   if(w>=1)
   {
      if(tauinf==0)
      {
        "            the curve is smooth at oo";"";
      }
      else
      {
         "number of singularities at oo:";nsing;
         "nodes at oo:";nodesinf;
         "cusps at oo:";cuspsinf;
         "branches at oo:";nbranchinf;
         "Tjurina number at oo:";tauinf;
         "delta at oo:";deltainf;
         "Milnor number at oo:";2*deltainf-nbranchinf+nsing;
         pause();
      }
      "singularities at (x,y,1):";"";
   }
   execute("ring newR=("+charstr(R)+"),(x,y),dp;");
   //the singularities at the affine part
   map sigma=S,var(1),var(2),1;
   ideal I=sigma(F);

   if(size(#)!=0)     
   {                              //uses the normalization to compute delta
      list nor=normal(I,"wd");
      delt=nor[size(nor)];
      genus=genus-delt-deltainf;
      setring R;
      return(genus);
   }

   ideal I1=jacob(I);
   matrix Hess[2][2]=jacob(I1);
   ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1
   
   ideal radID=std(radical(ID));//the non-nodal locus
   if(w>=1){"the non-nodal locus:";"";radID;pause();"";}
   if(deg(radID[1])==0)
   {
     ideal IDsing=1;
   }
   else
   {
     ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID 
   }
   
   iglob=vdim(std(IDsing));

   if(iglob!=0)//computation of the radical of IDsing
   {
      ideal radIDsing=reduce(IDsing,radID);
      if(size(radIDsing)==0)
      {
         radIDsing=radID;
         attrib(radIDsing,"isSB",1);
      }
      else
      {
         radIDsing=std(radical(IDsing));
      }
      iglob=vdim(radIDsing);
      if((w>=1)&&(iglob))
          {"the non-nodal-cuspidal locus:";radIDsing;pause();"";}
   }
   cusps=vdim(radID)-iglob;
   nsing=nsing+cusps;
   
   if(iglob==0)
   {
      if(w>=1){"             there are only cusps and nodes";"";}
      tau=vdim(std(I+jacob(I)));
      tauinf=tauinf+tau;
      nodes=tau-2*cusps; 
      delt=nodes+cusps;
      nbranch=2*tau-3*cusps; 
      nsing=nsing+nodes;
   }
   else
   {
       if(w>=1){"the non-nodal-cuspidal singularities";"";}
       setring C;                     
       ideal I1=imap(newR,IDsing);
       iloc=vdim(std(I1));
       if(iglob==iloc)
       { 
          if(w>=1){"only cusps and nodes outside (0,0,1)";}
          setring newR;
          tau=vdim(std(I+jacob(I)));
          tauinf=tauinf+tau;
          inv=deltaLoc(I[1],maxideal(1));
          delt=inv[1];
          tauloc=inv[2];
          nodes=tau-tauloc-2*cusps; 
          nsing=nsing+nodes;
          nbranch=inv[3]+ 2*nodes+cusps;         
          delt=delt+nodes+cusps;
          if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";}
        }
        else
        {
           setring newR;
           list pr=minAssGTZ(IDsing);
           if(w>=1){pr;}
           
           for(k=1;k<=size(pr);k++)
           {
              if(w>=1){nsing=nsing+vdim(std(pr[k]));}
              inv=deltaLoc(I[1],pr[k]);
              delt=delt+inv[1];
              tausing=tausing+inv[2];
              nbranch=nbranch+inv[3];
           }
           tau=vdim(std(I+jacob(I)));
           tauinf=tauinf+tau;
           nodes=tau-tausing-2*cusps;
           nsing=nsing+nodes;
           delt=delt+nodes+cusps;  
           nbranch=nbranch+2*nodes+cusps;            
        }
   }
   genus=genus-delt-deltainf;
   if(w>=1)
   {
      "The projected plane curve has locally:";"";
      "singularities:";nsing;
      "branches:";nbranch+nbranchinf;
      "nodes:"; nodes+nodesinf;
      "cusps:";cusps+cuspsinf;
      "Tjurina number:";tauinf;
      "Milnor number:";2*(delt+deltainf)-nbranch-nbranchinf+nsing;
      "delta of the projected curve:";delt+deltainf;
      "delta of the curve:";p_a-genus;
      "genus:";genus;
      "====================================================";
      "";
   }
   setring R;
   return(genus);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y),dp;
   ideal i=y^9 - x^2*(x - 1)^9;
   genus(i);
}

///////////////////////////////////////////////////////////////////////////////
proc deltaLoc(poly f,ideal singL)
"USAGE:  deltaLoc(f,J);  f poly, J ideal 
ASSUME: f is reduced bivariate polynomial; basering has exactly two variables;
        J is irreducible prime component of the singular locus of f (e.g., one
        entry of the output of @code{minAssGTZ(I);}, I = <f,jacob(f)>). 
RETURN:  list L: 
@texinfo
@table @asis
@item @code{L[1]}; int:
         the sum of (local) delta invariants of f at the (conjugated) singular
         points given by J.
@item @code{L[2]}; int:
         the sum of (local) Tjurina numbers of f at the (conjugated) singular
         points given by J.
@item @code{L[3]}; int:
         the sum of (local) number of branches of f at the (conjugated) 
         singular points given by J.
@end table
@end texinfo
NOTE:    procedure makes use of @code{execute}; increasing printlevel displays 
         more comments (default: printlevel=0).
SEE ALSO: delta, tjurina
KEYWORDS: delta invariant; Tjurina number
EXAMPLE: example deltaLoc;  shows an example
"
{
   option(redSB);
   def R=basering;
   execute("ring S=("+charstr(R)+"),(x,y),lp;");
   map phi=R,x,y;
   ideal singL=phi(singL);
   singL=std(singL);
   int d=vdim(singL);
   poly f=phi(f);
   int i;
   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
   if(d==1)
   {
      map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2];
      f=alpha(f);

      execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;");
      poly f=imap(S,f);
      ideal singL=imap(S,singL);
      if((w>=1)&&(ord(f)>=2))
      {
        "local analysis of the singularities";"";
        basering;
        singL;
        f;
        pause();
      }
   }
   else
   {
      poly p;
      poly c;
      map psi;
      number co;
      
      while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1))
      {
         psi=S,x,y+random(-100,100)*x;
         singL=psi(singL);
         singL=std(singL);
          f=psi(f);
      }
      
      if(deg(lead(singL[2]))==1)
      {
         p=singL[1];
         c=singL[2]-lead(singL[2]);
         co=leadcoef(singL[2]);
      }
      if(deg(lead(singL[1]))==1)
      {
         psi=S,y,x;
         f=psi(f);
         singL=psi(singL);
         p=singL[2];
         c=singL[1]-lead(singL[1]);;
         co=leadcoef(singL[1]);
      }
      
      execute("ring B=("+charstr(S)+"),a,dp;");
      map beta=S,a,a;
      poly p=beta(p);
      
      execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;");
      number p=number(imap(B,p));
      
      minpoly=p;
      //number c=number(imap(S,c));
      map iota=S,a,a;
      number c=number(iota(c));
      number co=iota(co);
      
      map alpha=S,x-c/co,y+a;
      poly f=alpha(f);
      f=cleardenom(f);
      if((w>=1)&&(ord(f)>=2))
      {
        "local analysis of the singularities";"";
        basering;
        alpha;
        f;
        pause();
        "";
      } 
   }
   option(noredSB);
   ideal fstd=std(ideal(f)+jacob(f));
   poly hc=highcorner(fstd);
   int tau=vdim(fstd);      
   int o=ord(f);
   int delt,nb;
   
   if(tau==0)                 //smooth case
   {
      setring R;
      return(list(0,0,1));
   }
   if((char(basering)>=181)||(char(basering)==0))
   {
      if(o==2)                //A_k-singularity 
      {
        if(w>=1){"A_k-singularity";"";}        
         setring R;
         delt=(tau+1)/2;
         return(list(d*delt,d*tau,d*(2*delt-tau+1)));   
      }
      if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2)))
      { 
        if(w>=1){"D_k- singularity";"";}
        
         setring R;
         delt=(tau+2)/2;
         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
      }

      int mu=vdim(std(jacob(f)));
      poly g=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon
 
      list NP=newtonpoly(g);
      if(w>=1){"Newton-Polygon:";NP;"";}
      int s=size(NP);
       
      if(is_NND(f,mu,NP))     
      { // the Newton-polygon is non-degenerate                            
        // compute nb, the number of branches
        for(i=1;i<=s-1;i++)
        {
          nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]);
        }
        if(w>=1){"Newton-Polygon is non-degenerated";"";}
        return(list(d*(mu+nb-1)/2,d*tau,d*nb));
      }

      if(w>=1){"Newton-Polygon is degenerated";"";}

      // the following can certainly be made more efficient when replacing
      // 'reddevelop' (used only for computing number of branches) by 
      // successive blowing-up + test if Newton polygon degenerate:
      if(s>2)    //  splitting of f
      {
         if(w>=1){"Newton polygon can be used for splitting";"";}
         intvec v=NP[1][2]-NP[2][2],NP[2][1];
         int de=w_deg(g,v);
         int st=w_deg(hc,v)+v[1]+v[2];
         poly f1=var(2)^NP[2][2];
         poly f2=jet(g,de,v)/var(2)^NP[2][2];
         poly h=g-f1*f2;
         de=w_deg(h,v);
         poly k;
         ideal wi=var(2)^NP[2][2],f2;
         matrix li;
         while(de<st)
         { 
           k=jet(h,de,v);
           li=lift(wi,k);
           f1=f1+li[2,1];
           f2=f2+li[1,1];
           h=g-f1*f2;
           de=w_deg(h,v);
         }
         nb=deltaLoc(f1,maxideal(1))[3]+deltaLoc(f2,maxideal(1))[3];
         return(list(d*(mu+nb-1)/2,d*tau,d*nb));
      }

      f=jet(f,deg(hc)+2);
      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
      list hne=reddevelop(f);
      nb=size(hne);
      setring R;
      kill HNEring;       
      return(list(d*(mu+nb-1)/2,d*tau,d*nb));
   }
   else             //the case of small characteristic
   {
      f=jet(f,deg(hc)+2);
      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
      delt=delta(f);
      return(list(d*delt,d*tau,d));
   }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x,y),dp;
  poly f=(x2+y^2-1)^3 +27x2y2;
  ideal I=f,jacob(f);
  I=std(I);
  list qr=minAssGTZ(I);
  size(qr);
  // each component of the singular locus either describes a cusp or a pair
  // of conjugated nodes:
  deltaLoc(f,qr[1]); 
  deltaLoc(f,qr[2]); 
  deltaLoc(f,qr[3]); 
  deltaLoc(f,qr[4]); 
  deltaLoc(f,qr[5]); 
  deltaLoc(f,qr[6]);
}
///////////////////////////////////////////////////////////////////////////////
// compute the weighted degree of p 
static proc w_deg(poly p, intvec v)
{
   if(p==0){return(-1);}
   int d=0;
   while(jet(p,d,v)==0){d++;}
   d=(transpose(leadexp(jet(p,d,v)))*v)[1];
   return(d);  
}

//proc hilbPoly(ideal J)
//{  
//   poly hp;
//   int i;
//   if(!attrib(J,"isSB")){J=std(J);}
//   intvec v = hilb(J,2);
//   for(i=1; i<=size(v); i++){ hp=hp+v[i]*(var(1)-i+2);}
//   return(hp);
//}

///////////////////////////////////////////////////////////////////////////

/*
                           Examples:
LIB"normal.lib";
//Huneke
ring qr=31991,(a,b,c,d,e),dp;
ideal i=
5abcde-a5-b5-c5-d5-e5,
ab3c+bc3d+a3be+cd3e+ade3,
a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;


//Vasconcelos (dauert laenger: 70 sec)
ring r=32003,(x,y,z,w,t),dp;
ideal i=
x2+zw,
y3+xwt,
xw3+z3t+ywt2,
y2w4-xy2z2t-w3t3;

//Theo1
ring r=32003,(x,y,z),wp(2,3,6);
ideal i=zy2-zx3-x6;

//Theo1a (CohenMacaulay and regular in codimension 2)
ring r=32003,(x,y,z,u),wp(2,3,6,6);
ideal i=zy2-zx3-x6+u2;


//Theo2
ring r=32003,(x,y,z),wp(3,4,12);
ideal i=z*(y3-x4)+x8;

//Theo2a
ring r=32003,(T(1..4)),wp(3,4,12,17);
ideal i=
T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
T(1)^6*T(2)*T(3)+T(1)^2*T(2)^4*T(3)+T(1)^3*T(2)^2*T(4)-T(1)^2*T(2)*T(3)^2+T(4)^2;

//Theo3
ring r=32003,(x,y,z),wp(3,5,15);
ideal i=z*(y3-x5)+x10;


//Theo4
ring r=32003,(x,y,z),dp;
ideal i=(x-y)*(x-z)*(y-z);

//Theo5
ring r=32003,(x,y,z),wp(2,1,2);
ideal i=z3-xy4;

//Theo6
ring r=32003,(x,y,z),dp;
ideal i=x2y2+x2z2+y2z2;

ring r=32003,(a,b,c,d,e,f),dp;
ideal i=
bf,
af,
bd,
ad;

//ist CM
//Sturmfels
ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ideal i=
bv+su,
bw+tu,
sw+tv,
by+sx,
bz+tx,
sz+ty,
uy+vx,
uz+wx,
vz+wy,
bvz;

//J S/Y
ring r=32003,(x,y,z,t),dp;
ideal i=
x2z+xzt,
xyz,
xy2-xyt,
x2y+xyt;

//St_S/Y
ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ideal i=
wy-vz,
vx-uy,
tv-sw,
su-bv,
tuy-bvz;

//Horrocks:
ring r=32003,(a,b,c,d,e,f),dp;
ideal i=
adef-16000be2f+16001cef2,
ad2f+8002bdef+8001cdf2,
abdf-16000b2ef+16001bcf2,
a2df+8002abef+8001acf2,
ad2e-8000bde2-7999cdef,
acde-16000bce2+16001c2ef,
a2de-8000abe2-7999acef,
acd2+8002bcde+8001c2df,
abd2-8000b2de-7999bcdf,
a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
abde-8000b2e2-acdf-16001bcef-8001c2f2,
abcd-16000b2ce+16001bc2f,
a2cd+8002abce+8001ac2f,
a2bd-8000ab2e-7999abcf,
ab3f-3bdf3,
a2b2f-2adf3-16000bef3+16001cf4,
a3bf+4aef3,
ac3e-10668cde3,
a2c2e+10667ade3+16001be4+5334ce3f,
a3ce+10669ae3f,
bc3d+8001cd3e,
ac3d+8000bc3e+16001cd2e2+8001c4f,
b2c2d+16001ad4+4000bd3e+12001cd3f,
b2c2e-10668bc3f-10667cd2ef,
abc2e-cde2f,
b3cd-8000bd3f,
b3ce-10668b2c2f-10667bd2ef,
abc2f-cdef2,
a2bce-16000be3f+16001ce2f2,
ab3d-8000b4e-8001b3cf+16000bd2f2,
ab2cf-bdef2,
a2bcf-16000be2f2+16001cef3,
a4d-8000a3be+8001a3cf-2ae2f2;


ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ideal k=
wy-vz,
vx-uy,
tv-sw,
su-bv,
tuy-bvz;
ideal j=x2y2+x2z2+y2z2;
ideal i=mstd(intersect(j,k))[2];


ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ideal i=
wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;


//riemenschneider  
ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
ideal i=
xz,
vx,
ux,
su,
qu,
txy,
stx,
qtx,
uv2z-uwz,
uv3-uvw,
puv2-puw;

ring r=0,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;

//Yoshihiko Sakai
ring r=0,(x,y),dp;   //genus 0 4 nodes and 6 cusps   
ideal i=(x2+y^2-1)^3 +27x2y2;

ring r=0,(x,y),dp;   //genus 0
ideal i=(x-y^2)^2 - y*x^3;

ring r=0,(x,y),dp;  //genus 4
ideal i=y3-x6+1;

int m=9;           // q=9: genus 0
int p=2;
int q=9;//2,...,9
ring r=0,(x,y),dp;  
ideal i=y^m - x^p*(x - 1)^q;

ring r=0,(x,y),dp;  //genus 19
ideal i=55*x^8+66*y^2*x^9+837*x^2*y^6-75*y^4*x^2-70*y^6-97*y^7*x^2;


ring r=0,(x,y),dp;  //genus 34
ideal i=y10+(-2494x2+474)*y8+(84366+2042158x4-660492)*y6
        +(128361096x4-47970216x2+6697080-761328152x6)*y4
        +(-12024807786x4-506101284x2+15052058268x6+202172841-3212x8)*y2
        +34263110700x4-228715574724x6+5431439286x2+201803238
        -9127158539954x10-3212722859346x8;

//Rob Koelman
ring r=0,(x,y,z),dp;//genus 10  with 26 cusps
ideal i=
761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+
 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2-
 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4-
 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+
 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-
 z^10-474*z^8*y^2-84366*z^6*y^4;

ring r=0,(x,y),dp;//genus 10  with 26 cusps
ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;


ring r=0,(x,y),dp;   // genus 1  with 5 cusps
ideal i=57y5+516x4y-320x4+66y4-340x2y3+73y3+128x2-84x2y2-96x2y;

//Mark van Hoeij
ring r=0,(x,y),dp;  //genus 19
ideal i=y20+y13x+x4y5+x3*(x+1)^2;

ring r=0,(x,y),dp;  //genus 35
ideal i=y30+y13x+x4y5+x3*(x+1)^2;

ring r=0,(x,y),dp;   //genus 55
ideal i=y40+y13x+x4y5+x3*(x+1)^2;

ring r=0,(x,y),dp;  //genus 4
ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);

ring r=0,(y,z,w,u),dp; //genus -5
ideal i=y2+z2+w2+u2,w4-u4;

ring r=0,(x,y,t),dp; //genus -5
ideal i=x8+8x7y+32x6y2+80x5y3+136x4y4+160x3y5+128x2y6+64xy7+16y8+4x6t2+24x5yt2+72x4y2t2+128x3y3t2+144x2y4t2+96xy5t2+32y6t2+14x4t4+56x3yt4+112x2y2t4+112xy3t4+40y4t4+20x2t6+40xyt6+8y2t6+9t8;

ring r=0,(y,z,w,u),dp; //genus 9
ideal i=y2+z2+w2+u2,z4+w4+u4;

ring r=0,(x,y,t),dp;
ideal i=
25x8+200x7y+720x6y2+1520x5y3+2064x4y4+1856x3y5+1088x2y6+384xy7+64y8-12x6t2-72x5yt2-184x4y2t2-256x3y3t2-192x2y4t2-64xy5t2-2x4t4-8x3yt4+16xy3t4+16y4t4+4x2t6+8xyt6+8y2t6+t8;

ring r=0,(x,y,t),dp;
ideal i=
32761x8+786264x7y+8314416x6y2+50590224x5y3+193727376x4y4+478146240x3y5+742996800x2y6+664848000xy7+262440000y8+524176x7t+11007696x6yt+99772992x5y2t+505902240x4y3t+1549819008x3y4t+2868877440x2y5t+2971987200xy6t+1329696000y7t+3674308x6t2+66137544x5yt2+499561128x4y2t2+2026480896x3y3t2+4656222144x2y4t2+5746386240xy5t2+2976652800y6t2+14737840x5t3+221067600x4yt3+1335875904x3y2t3+4064449536x2y3t3+6226336512xy4t3+3842432640y5t3+36997422x4t4+443969064x3yt4+2012198112x2y2t4+4081745520xy3t4+3126751632y4t4+59524208x3t5+535717872x2yt5+1618766208xy2t5+1641991392y3t5+59938996x2t6+359633976xyt6+543382632y2t6+34539344xt7+103618032yt7+8720497t8;

ring r=32003,(x,y,z,w,u),dp;
ideal i=x2+y2+z2+w2+u2,x3+y3+z3,z4+w4+u4;

*/