///////////////////////////////////////////////////////////////////////////////
version="$Id: normal.lib,v 1.37 2001-10-23 10:26:24 pfister 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

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

LIB "general.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 delta;
   if(isIso&&isEq){delta=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]=delta;
   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 choose="wd" the delta-invariant is computed simultaneously
         this may take much more time in the reducible case because the
         factorizing standardbasis 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 choose="wd" an integer
@format
         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)
@end format
NOTE:    to use the i-th ring type: 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).
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 delta, 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]=delta;
         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]=delta;
         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]=delta;
         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]=delta;
      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]=delta;
      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]=delta;
         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]=delta;
      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),delta+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);
         delta=delta+RR[3];
         export norid;
         export normap;
         result=newR7;
         result[size(result)+1]=delta;
         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);
            
          delta=delta+mul+keepresult1[size(keepresult1)]
                         +keepresult1[size(keepresult1)];

          for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
          {
             keepresult1=insert(keepresult1,keepresult2[lauf]);
          }
          keepresult1[size(keepresult1)]=delta;
          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]=delta+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),delta+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]=delta;
         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)]=delta+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 K,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;
   K=std(K);
 
   if(nvars(R)==3)
   {
      if((dim(K)!=2)||(!homog(K))||(size(K)>1)){ERROR("Input is not a curve");}
      execute("ring newR=("+charstr(R)+"),(x,y),dp;");  
      map kappa=R,x,y,1;
      ideal K=kappa(K);  
   }
   if((nvars(R)>3)||(size(K)>1))
   {  // hier geeignet projezieren
      ERROR("This case is not implemented yet");
   }
   if(nvars(R)==2)
   {
      execute("ring newR=("+charstr(R)+"),(x,y),dp;");  
      map kappa=R,x,y;
      ideal K=kappa(K);  
   } 
   
   // assume now that R is a ring with two variables  
   poly p=K[1];   
   ideal I;
   if(homog(p))
   {
      if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");}
      return(-deg(p)+1);
   }
   execute("ring S=("+charstr(R)+"),(x,y,t),dp;");
   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;
   poly F=imap(newR,p);
   F=homog(F,t);
   int d=deg(F);
   int genus=(d-1)*(d-2)/2;
   
//   if(w>=1){"test for smoothness";}
//   if(dim(std(jacob(F)))==0)  //the smooth case
//   {
//      setring R;
//      return(genus);
//   }
 
   int delta,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,
       tauloc,tausing,k,rat,nbranchinf,nbranch,nodes;
   list inv;

   if(w>=1){"singularities at oo";}
   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);
      for(k=1;k<=size(qr);k++)
      {
         inv=deltaLoc(g,qr[k]);
         deltainf=deltainf+inv[1];
         tauinf=tauinf+inv[2];
         nbranchinf=nbranchinf+inv[3];
      }
   }
   inv=deltaLoc(h,maxideal(1)); 
   deltainf=deltainf+inv[1];
   tauinf=tauinf+inv[2];
   if(inv[2]>0){nbranchinf=nbranchinf+inv[3];}
  
   if(w>=1)
   {
      "branches at oo:";nbranchinf;
      "tau at oo:";tauinf;
      "delta at oo:";deltainf;
      "singularities not at oo";
   }
                          
   setring newR;         //the singularities at the affine part
   map sigma=S,var(1),var(2),1;
   I=sigma(F);

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

   ideal I1=jacob(I);
   matrix Hess[2][2]=jacob(I1);
   ideal ID=I+I1+ideal(det(Hess));

   if(w>=1){"the cusps and nodes";}
   
   ideal radID=std(radical(ID));
   ideal IDsing=minor(jacob(ID),2)+radID;
   iglob=vdim(std(IDsing));

   if(iglob!=0)
   {
      ideal radIDsing=reduce(IDsing,radID);
      if(size(radIDsing)==0)
      {
         radIDsing=radID;
         attrib(radIDsing,"isSB",1);
      }
      else
      {
         radIDsing=std(radical(IDsing));
      }
      iglob=vdim(radIDsing);
   }
   cusps=vdim(radID)-iglob; 
        
   if(w>=1){"the other singularities";}
  
   if(iglob==0)   //only cusps and double points
   {
      tau=vdim(std(I+jacob(I)));
      nodes=tau-2*cusps; 
      delta=nodes+cusps;
      nbranch=2*tau-3*cusps;    
   }
   else
   {
       setring C;                     
       ideal I1=imap(newR,IDsing);
       iloc=vdim(std(I1));
       if(iglob==iloc)  // only cusps and nodes outside 0
       { 
          setring newR;
          tau=vdim(std(I+jacob(I)));
          inv=deltaLoc(I[1],maxideal(1));
          delta=inv[1];
          tauloc=inv[2];
          nodes=tau-tauloc-2*cusps; 
          nbranch=inv[3]+ 2*nodes+cusps;          
          delta=delta+nodes+cusps;
        }
        else
        {
           setring newR;
           list pr=minAssGTZ(IDsing);                   
           for(k=1;k<=size(pr);k++)
           {
              inv=deltaLoc(I[1],pr[k]);
              delta=delta+inv[1];
              tausing=tausing+inv[2];
              nbranch=nbranch+inv[3];
           }
           tau=vdim(std(I+jacob(I)));

           nodes=tau-tausing-2*cusps;
           delta=delta+nodes+cusps;  
           nbranch=nbranch+2*nodes+cusps;            
        }
   }

   if(w>=1)
   {
      "branches :";nbranch;
      "nodes:"; nodes;
      "cusps:";cusps;
      "tau :";tau;
      "delta:";delta;
   }
   genus=genus-delta-deltainf;
   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)
{
   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;

   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);
   }
   else
   {
      poly p;
      poly c;
      map psi;
      while((deg(singL[1])>1)&&(deg(singL[2])>1))
      {
         psi=S,x,y+random(-100,100)*x;
         singL=psi(singL);
         singL=std(singL);
      }
      if(deg(singL[2])==1){p=singL[1];c=singL[2][2];}
      if(deg(singL[1])==1)
      {
         psi=S,y,x;
         f=psi(f);
         singL=psi(singL);
         p=singL[2];
         c=singL[1][2];
      }
      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 alpha=S,x-c,y+a;
     
      poly f=alpha(f);
      f=cleardenom(f);
   }
   ideal fstd=std(ideal(f)+jacob(f));
   poly hc=highcorner(fstd);
   int tau=vdim(fstd);      
   int o=ord(f);
   int delta,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 
      {
         setring R;
         delta=(tau+1)/2;
         return(list(d*delta,d*tau,d*(2*delta-tau+1)));   
      }
      if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2)))
      {//D_k- singularity
         setring R;
         delta=(tau+2)/2;
         return(list(d*delta,d*tau,d*(2*delta-tau+1)));
      }

      int mu=vdim(std(jacob(f)));
      poly g=f+var(1)^mu+var(2)^mu;  //to obtain a convinient Newton-polygon
 
      list NP=newton(g);

      int s=size(NP);
      int nN=-NP[1][2]-NP[s][1]+1;  // computation of the Newton-number
      intmat m[2][2];
      for(i=1;i<=s-1;i++)
      {
         m=NP[i+1],NP[i];
         nN=nN+det(m);
      }

      if(mu==nN)                   // 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]);
        }
        return(list(d*(mu+nb-1)/2,d*tau,d*nb));
      }

   //da reddevelop nur benutzt wird, um die Anzahl der Zweige zu bestimmen
   //kann man das sicher schneller machen:
   //die Aufblasung durchfuehren und stets testen, ob das Newton-polyeder
   //nicht ausgeartet ist.

      if(s>2)    //  splitting of f
      {
         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);
      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);
      list hne=reddevelop(f);
      nb=size(hne);
      if(nb==1)
      {
         delta=invariants(hne[1])[5]/2;
         setring R;
         kill HNEring;
         return(list(d*delta,d*tau,d));
      }
      setring R;
      kill HNEring;
      //delta direkt aus reddevelop zurueckgeben
      ERROR("the case of small characteristic is not fully implemented yet");
   }
}

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 newton (poly f)
{
   def R1=basering;
   execute("ring R2=("+charstr(R1)+"),("+varstr(R1)+"),ls;");
   poly f=imap(R1,f);
   intvec A=(0,ord(subst(f,var(1),0)));
   intvec B=(ord(subst(f,var(2),0)),0);
   intvec C,H; list L;
   int abbruch,i;
   poly hilf;
   L[1]=A;
   f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
   map xytausch=R2,var(2),var(1);
   for (i=2; f!=0; i++) 
   {
      abbruch=0;
      while (abbruch==0) 
      {
         C=leadexp(f);          
         if(jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
         { 
            abbruch=1; 
         }       
         else 
         { 
            f=jet(f,-C[1]-1,intvec(-1,0)); 
         }
     }
     hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));
     H=leadexp(xytausch(hilf));
     A=H[2],H[1]; 
     L[i]=A;
     f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]-A[1]));
   }
   L[i]=B;
   setring R1;
   return(L);
}

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

/*
                           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;   // genuss 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);

*/