source: git/Singular/LIB/normal.lib @ 1288ef

///////////////////////////////////////////////////////////////////////////////
version="$Id: normal.lib,v 1.50 2009-03-26 16:37:19 Singular Exp $";
category="Commutative Algebra";
info="
LIBRARY:  normal.lib     Normalization of Affine Rings
AUTHORS:  G.-M. Greuel,  greuel@mathematik.uni-kl.de,
@*        S. Laplagne,   slaplagn@dm.uba.ar,
@*        G. Pfister,    pfister@mathematik.uni-kl.de


MAIN PROCEDURES:
 normal(I[...]);     normalization of an affine ring 
 normalP(I,[...];    normalization of an affine ring in positive characteristic
 normalC(I,[...];    normalization of an affine ring through a chain of rings
 HomJJ(L);           presentation of the End_R(J) as affine ring, J an ideal
 genus(I);           computes the geometric genus of a projective curve
 primeClosure(L);    integral closure of R/p, p a prime ideal
 closureFrac(L);     write a poly in integral closure as element of Quot(R/p)
 iMult(L);           intersection multiplicity of the ideals of the list L

AUXILIARY PROCEDURES:
 deltaLoc(f,S);      sum of delta invariants at conjugated singular points
 locAtZero(I);       checks whether the zero set of I is located at 0
 norTest(i,nor;)     checks whether result of normal or mormal_p was correct
";

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";
LIB "reesclos.lib";
LIB "algebra.lib";


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

proc normal(ideal id, list #)
"USAGE:  normal(id [,choose]); id = radical ideal, choose = optional string. @*
         Optional parameters in list choose (can be entered in any order):@*
         Decomposition:@*
         - \"prim\" -> computes first the minimal associated primes, and then 
         the normalization of each prime. (default)@*
         - \"equidim\" -> computes first an equidimensional decomposition,
         and then the normalization of each component. @*
         - \"noDeco\" -> no preliminary decomposition is done. If the ideal is 
         not equidimensional radical, output might be wrong.@*
         - \"isPrim\" -> assumes that the ideal is prime. If assumption does 
         not hold, output might be wrong.@*
         - \"noFac\" -> factorization is avoided in the computation of the 
         minimal associated primes;
         Other:@*
         - \"useRing\" -> uses the original ring ordering.@*
         If the ring ordering is not global, it will change to a global
         ordering only for computing radicals and prime or equidimensional
         decompositions.@*
         If this option is not set, changes to dp ordering and perform all
         computation with respect to this ordering.@*
         - \"withDelta\" (or \"wd\") -> returns also the delta invariants.
         If choose is not given or empty, the default options are used.@*
ASSUME:  The ideal must be radical, for non-radical ideals the output may
         be wrong (id=radical(id); makes id radical). However, when using 
         \"prim\" option the minimal associated primes of id are computed first 
         and hence normal computes the normalization of the radical of id.
         \"isPrim\" should only be used if id is known to be prime.
RETURN:  a list, say nor, of size 1 (resp. 2 with option \"withDelta\"). 
@format 
         The first element nor[1] is a list of r elements, where r is the 
         number of associated primes P_i with option \"prim\" (resp. >= no 
         of equidimenensional components P_i with option \"equidim\").@*
         Each element of nor[1] is a list with the information on the
         normalization of the i-th component, i.e. the integral closure 
         in its field of fractions (as affine ring).@*
         Each list nor[1][i] contains 3 elements (resp. 4 with option
         \"withDelta\"). @*
         The first two element of nor[1][i] give the normalization by
         module generators. The normalization is of the form
         1/c * U, with U = nor[1][i][1] an ideal of R / P_i and
         c = nor[1][i][2] a polynomial in R / P_i.@*
         The third element is a ring Ri=nor[1][i][3], i=1..r, 
         which contains two ideals with given
         names @code{norid} and @code{normap} such that @*
         - Ri/norid is the normalization of the i-th component.
         - the direct sum of the rings Ri/norid is the normalization
           of basering/id; @*
         - @code{normap} gives the normalization map from basering/id to
           Ri/norid for each j.@*
         When option \"withDelta\" is set, the forth element nor[1][i][4] is
         the delta invariant of the i-th component. (-1 means infinite, and 0 
         that basering/P_i is normal)@*
         @*
         Also, When option \"withDelta\" is set the second element of the list
         nor, nor[2], is an integer, the total delta invariant of basering/id
         (-1 means infinite, and 0 that basering/id is normal).
@end format
THEORY:  We use a general algorithm described in [G.-M. Greuel, S. Laplagne, 
         F. Seelisch: Normalization of Rings (2009)].
         The delta invariant of a reduced ring K[x1,...,xn]/id is 
             dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id)
         We call this number also the delta invariant of id.
NOTE:    To use the i-th ring type: @code{def R=nor[1][i][3]; setring R;}.
@*       Increasing/decreasing printlevel displays more/less comments 
         (default: printlevel=0).
@*       Implementation works also for local rings.         
@*       Not implemented for quotient rings.
@*       If the input ideal id is weighted homogeneous a weighted ordering may
         be used (qhweight(id); computes weights).
KEYWORDS: normalization; integral closure; delta invariant.
EXAMPLE: example normal; shows an example
"
{  
   intvec opt = option(get);     // Save current options

   int i,j;
   int decomp;  // Preliminar decomposition:
                // 0 -> no decomposition (id is assumed to be prime)
                // 1 -> no decomposition 
                //      (id is assumed to be equidimensional radical) 
                // 2 -> equidimensional decomposition
                // 3 -> minimal associated primes
   int noFac, useRing, withDelta;
   int dbg = printlevel - voice + 2;
   int nvar = nvars(basering);
   int chara  = char(basering);
   list result;
   list keepresult;
   list ringStruc;
   ideal U;
   poly c;
   
   // Default methods:
   noFac = 0;         // Use facSTD when computing minimal associated primes
   decomp = 3;        // Minimal associated primes
   useRing = 0;       // Change first to dp ordering, and perform all 
                      // computations there.
   withDelta = 0;     // Do not compute the delta invariant.
  

  //--------------------------- define the method --------------------------- 
  string method;                //make all options one string in order to use
                               //all combinations of options simultaneously
  for ( i=1; i <= size(#); i++ )
  { 
    if ( typeof(#[i]) == "string" ) 
    {
      method = method + #[i];
    }
  }

  //--------------------------- choosen methods -----------------------
  if ( find(method,"isprim") or find(method,"isPrim") )
  {decomp = 0;}
  
  if ( find(method,"nodeco") or find(method,"noDeco") )
  {decomp = 1;}
  
  if ( find(method,"equidim") )
  {decomp = 2;}
  
  if ( find(method,"prim") )
  {decomp = 3;}
  
  if ( find(method,"nofac") or find(method,"noFac") )
  {noFac=1;}
  
  if ( (find(method,"useRing") or find(method,"usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
  {useRing = 1;}
  
  if ( find(method,"withDelta") or find(method,"wd") or find(method,"withdelta"))
  {
    if((decomp == 0) or (decomp == 3)){
      withDelta = 1;
    } else {
      "WARNING: the delta invariants cannot be computed with an equidimensional";
      "         decomposition.";
      "         Use option \"prim\" to compute first the minimal associated primes.";
      "";
    }
  }

  kill #;
  list #;

  //------------------------ change ring if required ------------------------
  // If the ordering is not global, we change to dp ordering for computing the 
  // min ass primes.
  // If the ordering is global, but not dp, and useRing = 0, we also change to 
  // dp ordering.
  
  int isGlobal = ord_test(basering);      // Checks if the original ring has
                                          // global ordering.
  
  def origR = basering;   // origR is the original ring
                          // R is the ring where computations will be done
    
  if((useRing  == 1) and (isGlobal == 1)){
    def globR = basering;
  } else {
    // We change to dp ordering.
    list rl = ringlist(origR);
    list origOrd = rl[3];
    list newOrd = list("dp", intvec(1:nvars(origR))), list("C", 0);
    rl[3] = newOrd;
    def globR = ring(rl);
    setring globR;
    ideal id = fetch(origR, id);  
  }

  //------------------------ preliminary decomposition-----------------------
  list prim;
  if(decomp == 2){
    dbprint(dbg, "// Computing the equidimensional decomposition...");
    prim = equidim(id);
  }
  if((decomp == 0) or (decomp == 1)){
    prim = id;
  }
  if(decomp == 3){
    dbprint(dbg, "// Computing the minimal associated primes...");
    if( noFac )
    { prim = minAssGTZ(id,1); }
    else
    { prim = minAssGTZ(id); }
  }

  // if ring was not global and useRing is on, we go back to the original ring
  if((useRing == 1) and (isGlobal != 1)){
    setring origR;
    def R = basering;
    list prim = fetch(globR, prim);
  } else {
    def R = basering;
    if(useRing == 0){
      ideal U;
      poly c;
    }
  }
    
  if(dbg>=1)
  {  
    prim; "";
    "// number of components is", size(prim);
    "";
  }
   
  // ---------------- normalization of the components-------------------------
  // calls normalM to compute the normalization of each component.

  list norComp;       // The normalization of each component.
  int delt;
  int deltI = 0;
  int totalComps = 0;

  setring origR;
  def newROrigOrd;
  list newRListO;
  setring R;
  def newR;
  list newRList;
    
  for(i=1; i<=size(prim); i++)
  {
    if(dbg>=2){pause();}
    if(dbg>=1){
      "// start computation of component",i;
      "   --------------------------------";
    }
    if(groebner(prim[i])[1] != 1){
      if(dbg>=2){
        "We compute the normalization in the ring"; basering;
      }
      printlevel = printlevel + 1;
      norComp = normalM(prim[i], decomp, withDelta);
      printlevel = printlevel - 1;      
      for(j = 1; j <= size(norComp); j++){
        newR = norComp[j][3];
        newRList = ringlist(newR);
        U = norComp[j][1];
        c = norComp[j][2];
        if(withDelta){
          delt = norComp[j][4];
          if((delt >= 0) and (deltI >= 0)){
            deltI = deltI + delt;
          } else {
            deltI = -1;
          }
        }
        if(useRing == 0){
        // We go back to the original ring.
        setring origR;
        U = fetch(R, U);
        c = fetch(R, c);
        newRListO = imap(R, newRList); 
        // We change the ordering in the new ring.
        if(nvars(newR) > nvars(origR)){
          newRListO[3]=insert(origOrd, newRListO[3][1]);
        } else {
          newRListO[3] = origOrd;
        }
        newROrigOrd = ring(newRListO);
        setring newROrigOrd;
        ideal norid = imap(newR, norid);
        ideal normap = imap(newR, normap);
        export norid;
        export normap;
        setring origR;
        totalComps++;
        result[totalComps] = list(U, c, newROrigOrd);
        if(withDelta){
          result[totalComps] = insert(result[totalComps], delt, 3);
        }
        setring R;
        } else {
          setring R;
          totalComps++;
          result[totalComps] = norComp[j];
        }
      }
    }
  }
  
  // -------------------------- delta computation ----------------------------  
  if(withDelta == 1){
    // Intersection multiplicities of list prim, sp=size(prim).
    if ( dbg >= 1 )
    {
      "// Sum of delta for all components: ", deltI;
    }   
    if(size(prim) > 1){
      dbprint(dbg, "// Computing the sum of the intersection multiplicities of the components...");
      int mul = iMult(prim);
      if ( mul < 0 )
      {
        deltI = -1;
      } else {
        deltI = deltI + mul; 
      }
      if ( dbg >= 1 )
      {
        "// Intersection multiplicity is : ", mul;
      }         
    }
  }
  
  // -------------------------- prepare output ------------------------------   
  setring origR;
  if(withDelta){
    result = result, deltI;
  } else {
    result = list(result);
  }
  
  option(set, opt);
  return(result);
}


example
{ "EXAMPLE:";
  printlevel = printlevel+1;
  echo = 2;
  ring s = 0,(x,y),dp;
  ideal i = (x2-y3)*(x2+y2)*x;
  list nor = normal(i, "withDelta");
  nor; 

  // 2 branches have delta = 1, and 1 branch has delta = 0
  // the total delta invariant is 13

  def R2 = nor[1][2][3];  setring R2;
  norid; normap;
   
  echo = 0;
  printlevel = printlevel-1;
  pause("   hit return to continue"); echo=2;

  ring r = 2,(x,y,z),dp;
  ideal i = z3-xy4;
  list nor = normal(i, "withDelta");  nor;  
  // the delta invariant is infinite 
  // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
  // in its quotient field Quot(r/i)

  // the normalization as affine algebra over the ground field:             
  def R = nor[1][1][3]; setring R;
  norid; normap;
}

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

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 map R --> Hom_R(J,J) of affine rings,
         where R is the quotient ring of P modulo the standard basis SBid.
RETURN:  a list l of three 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, = dim_K(Hom_R(J,J)/R) (the contribution to delta)
                if the dimension is finite, -1 otherwise
@end format
NOTE:    printlevel >=1: display comments (default: printlevel=0)
EXAMPLE: example HomJJ;  shows an example
"
{
//---------- initialisation ---------------------------------------------------
   int isIso,isPr,isHy,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";
     "//   the ideal J:";J;
   }
   f  = quotient(p*J,J);

   //### (neu GMG 4.10.08) divide by the greatest common divisor:
   poly gg = gcd( f[1],p );
   for(ii=2; ii <=ncols(f); ii++)
   {
      gg=gcd(gg,f[ii]);
   }
   for(ii=1; ii<=ncols(f); ii++)
   {
      f[ii]=f[ii]/gg;
   }
   p = p/gg;

   if ( y>=1 )
   {  
      "//   the non-zerodivisor p:"; p;
      "//   the module p*Hom(J,J) = p*J:J :"; f;
      "";
   }
   f2 = std(p);

//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------

   //rf = interred(reduce(f,f2)); 
   //### interred hier weggelassen, unten zugefŸgt
   rf = 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);
      endid = simplify(endid,2);
      L = substpart(endid,endphi,homo,rw);   //## hier substpart
      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,"isHypersurface",0);
      attrib(endid,"isCompleteIntersection",0);
      attrib(endid,"isRadical",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 to the new ring";
         endphi;
         "   ";
         pause();
         "";
      }
      setring P;
      L[3]=0;
      return(L);
   }
   if(y>=1)
   {
      "// R is not equal to Hom(J,J), we have to try again";
      pause();
      "";
   }
//---------- 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
// f2=p (i.e. ideal generated by p)

   //f = mstd(f)[2];              //### geŠndert GMG 04.10.08 
   //ideal ann = quotient(f2,f);  //### f durch rf ersetzt  
   rf = mstd(rf)[2];              //rf = NF(f,p), hence <p,rf> = <p,f>
   ideal ann = quotient(f2,rf);   //p:f = p:rf 
 
   //------------- compute the contribution to delta ---------- 
   //delt=dim_K(Hom(JJ)/R (or -1 if infinite)
  
   int 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);          //### psi1 durch fetch ersetzt
   //ideal SBid = psi1(SBid);        
   ideal SBid = fetch(P,SBid);
   attrib(SBid,"isSB",1);

   qring newR = std(SBid);

   //map psi = R,ideal(X(1..nvars(R)));  //### psi durch fetch ersetzt
   //ideal id = psi(id);
   //ideal f = psi(f);
   //module syzf = psi(syzf);
   ideal id = fetch(R,id);
   ideal f = fetch(R,f);
   module syzf = fetch(R,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
// f=p,rf, rf=reduce(f,p), generates pJ:J mod(p), 
// i.e. p*Hom(J,J)/p*R as R-module

   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;
   }

   poly ff;
   for (ii=2; ii<=q; ii++ )
   {
      for ( jj=2; jj<=ii; jj++ )
      {
         ff = NF(f[ii]*f[jj],std(0));       //this makes lift much faster 
         A = lift(pf,ff);                   //ff lin. comb. of elts of pf mod I
         Quad = Quad, ideal(T(jj)*T(ii) - T*A);  //quadratic relations
      }
   }

   if(y>=1)
   {
      "//   the quadratic relations";
      Quad;
      pause();
      newline;
   }
   Q = Lin,Quad;
   Q = subst(Q,T(1),1);
   //Q = mstd(Q)[2];            //### sehr aufwendig, daher weggelassen (GMG)
   //### ev das neue interred
   //mstd dient nur zum verkleinern, die SB-Eigenschaft geht spaeter verloren
   //da in neuen Ring abgebildet und mit id vereinigt
                       
//---------- 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);  
   //hier wird Q weiterverwendet, die SB-Eigenschaft wird nicht verwendet.
   endid = simplify(endid,2);
   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))   //#### ok???
   // {
   //    attrib(endid,"onlySingularAtZero",oSAZ);
   // }
   //kill te;
   attrib(endid,"isCohenMacaulay",isCo);                //#### ok???
   attrib(endid,"isPrim",isPr);
   attrib(endid,"isIsolatedSingularity",isIso);
   attrib(endid,"isRegInCodim2",isRe);
   attrib(endid,"isEquidimensional",isEq);
   attrib(endid,"isHypersurface",0);
   attrib(endid,"isCompleteIntersection",0);
   attrib(endid,"isRadical",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 to the new ring";
      endphi;
      "   ";
      pause();
      "";
   }
   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;
  L[3];              // contribution to delta
}


///////////////////////////////////////////////////////////////////////////////
//compute intersection multiplicities as needed for delta(I) in 
//normalizationPrimes and normalP:

proc iMult (list prim)
"USAGE:   iMult(L);  L a list of ideals
RETURN:  int, the intersection multiplicity of the ideals of L;
         if iMult(L) is infinite, -1 is returned. 
THEORY:  If r=size(L)=2 then iMult(L) = vdim(std(L[1]+L[2])) and in general
         iMult(L) = sum{ iMult(L[j],Lj) | j=1..r-1 } with Lj the intersection
         of L[j+1],...,L[r]. If I is the intersection of all ideals in L then
         we have delta(I) = delta(L[1])+...+delta(L[r]) + iMult(L) where
         delta(I) = vdim (normalisation(R/I)/(R/I)), R the basering. 
EXAMPLE: example iMult; shows an example
"
{    int i,mul,mu;
     int sp = size(prim);
     int y = printlevel-voice+2;
     if ( sp > 1 )
     {
        ideal I(sp-1) = prim[sp];
        mu = vdim(std(I(sp-1)+prim[sp-1]));
        mul = mu;
        if ( y>=1 )
        {
          "// intersection multiplicity of component",sp,"with",sp-1,":"; mu;
        }
        if ( mu >= 0 )
        {
           for (i=sp-2; i>=1 ; i--)
           {
              ideal I(i) = intersect(I(i+1),prim[i+1]);
              mu = vdim(std(I(i)+prim[i]));
              if ( mu < 0 )  
              { 
                break; 
              }
              mul = mul + mu;
              if ( y>=1 )
              {
                "// intersection multiplicity of components",sp,"...",i+1,"with",i; mu;
              }
           }
        }
     }
     return(mul);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s  = 23,(x,y),dp;
   list L = (x-y),(x3+y2);
   iMult(L);
   L = (x-y),(x3+y2),(x3-y4);
   iMult(L);
}  
///////////////////////////////////////////////////////////////////////////////
//check if I has a singularity only at zero, as needed in normalizationPrimes

proc locAtZero (ideal I)
"USAGE:   locAtZero(I);  I = ideal
RETURN:  int, 1 if I has only one point which is located at zero, 0 otherwise
ASSUME:  I is given as a standard bases in the basering
NOTE:    only useful in affine rings, in local rings vdim does the check
EXAMPLE: example locAtZero; shows an example
"
{    
   int ii,jj, caz;                   //caz: conzentrated at zero
   int dbp = printlevel-voice+2;
   int nva = nvars(basering);
   int vdi = vdim(I);     
   if ( vdi < 0 )
   {
      if (dbp >=1) 
      { "// non-isolated singularitiy";""; }         
      return(caz);
   }

   //Now the ideal is 0-dim
   //First an easy test
   //If I is homogenous and not constant it is concentrated at 0
   if( homog(I)==1 && size(jet(I,0))==0)
   {
      caz=1;
      if (dbp >=1) 
      { "// isolated singularity and homogeneous";""; }         
      return(caz);
   }

   //Now the general case with I 0-dim. Choose an appropriate power pot,
   //and check each variable x whether x^pot is in I.
   int mi1 = mindeg1(lead(I));
   int pot = vdi;
   if ( (mi1+(mi1==1))^2 < vdi )
   {
      pot = (mi1+(mi1==1))^2;      //### alternativ: pot = vdi lassen
   }

   while ( 1 )
   {
      caz = 1;
      for ( ii=1; ii<= nva; ii++ )
      {
        if ( NF(var(ii)^pot,I) != 0 )
        {
           caz = 0; break;
        }
      }
      if ( caz == 1 || pot >= vdi )
      {
        if (dbp >=1) 
        {
          "// mindeg, exponent, vdim used in 'locAtZero':", mi1,pot,vdi; "";
        }         
        return(caz);
      }
      else
      { 
        if ( pot^2 < vdi )
        { pot = pot^2; }
        else
        { pot = vdi; }
      }
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),dp;
   poly f = z5+y4+x3+xyz;
   ideal i = jacob(f),f;
   i=std(i);
   locAtZero(i);
   i= std(i*ideal(x-1,y,z));
   locAtZero(i);
}  

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

//The next procedure normalizationPrimes computes the normalization of an 
//irreducible or an equidimensional ideal i. 
//- If i is irreducuble, then the returned list, say nor, has size 2 
//with nor[1] the normalization ring and nor[2] the delta invariant.
//- If i is equidimensional, than the "splitting tools" can create a 
//decomposition of i and nor can have more than 1 ring. 

static proc normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #)
"USAGE:   normalizationPrimes(i,ihp,delt[,si]);  i = equidimensional ideal, 
         ihp = map (partial normalization), delt = partial delta-invariant,
         si = ideal s.t. V(si) contains singular locus (optional)
RETURN:   a list of rings, say nor, and an integer, the delta-invariant
          at the end of the list.
          each ring nor[j], j = 1..size(nor)-1, contains two ideals
          with given names norid and normap such that
           - the direct sum of the rings nor[j]/norid is
             the normalization of basering/i;
           - normap gives the normalization map from basering/id
             to nor[j]/norid (for each j)
          nor[size(nor)] = dim_K(normalisation(P/i) / (P/i)) is the 
          delta-invariant, where P is the basering.
EXAMPLE: example normalizationPrimes; shows an example
"
{
   //Note: this procedure calls itself as long as the test for
   //normality, i.e if R==Hom(J,J), is negative.

   int printlev = printlevel;   //store printlevel in order to reset it later
   int y = printlevel-voice+2;  // y=printlevel (default: y=0)
   if(y>=1)
   {
     "";
     "// START a normalization loop with the ideal";
     i;  "";
     "// in the ring:";
     basering;  "";
     pause();
     "";
   }

   def BAS=basering;
   list result,keepresult1,keepresult2,JM,gnirlist;
   ideal J,SB,MB;
   int depth,lauf,prdim,osaz;
   int ti=timer;
   
   gnirlist = ringlist(BAS);

//----------- the trivial case of a zero ideal as input, RETURN ------------
   if(size(i)==0)
   {
      if(y>=1)
      {
          "// the ideal was the zero-ideal";
      }
     // execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
     //                 +ordstr(basering)+");");
         def newR7 = ring(gnirlist);
         setring newR7;
         ideal norid=ideal(0);
         ideal normap=fetch(BAS,ihp);
         export norid;
         export normap;
         result=newR7;
         result[size(result)+1]=list(delt,delti);
         setring BAS;
         return(result);
   }

//--------------- General NOTATION, compute SB of input ----------------- 
// SM is a list, the result of mstd(i)
// SM[1] = SB of input ideal i, 
// SM[2] = (minimal) generators for i.
// We work with SM and will copy the attributes from i to SM[2]
// JM will be a list, either JM[1]=maxideal(1),JM[2]=maxideal(1)
// in case i has onlySingularAtZero, or JM = mstd(si) where si = #[1],
// or JM = mstd(J) where J is the ideal of the singular locus
// JM[2] must be (made) radical

   if(y>=1)
   {
     "// SB-computation of the ideal";
   }

   list SM = mstd(i);              //Now the work starts
   int dimSM =  dim(SM[1]);        //dimension of variety to normalize
   if(y>=1)
   {
      "// the dimension is:";  dimSM;
   }
//----------------- the general case, set attributes ----------------
   //Note: onlySingularAtZero is NOT preserved under the ring extension
   //basering --> Hom(J,J) (in contrast to isIsolatedSingularity), 
   //therefore we reset it:
  
   attrib(i,"onlySingularAtZero",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,"isHypersurface")==1)
   {
     attrib(SM[2],"isHypersurface",1);
   }
   else
   {
      attrib(SM[2],"isHypersurface",0);
   }

   if(attrib(i,"onlySingularAtZero")==1)
   {
      attrib(SM[2],"onlySingularAtZero",1);
   }
   else
   {
      attrib(SM[2],"onlySingularAtZero",0);
   }

   //------- an easy and cheap test for onlySingularAtZero ---------
   if( (attrib(SM[2],"isIsolatedSingularity")==1) && (homog(SM[2])==1) )
   {
      attrib(SM[2],"onlySingularAtZero",1);
   }

//-------------------- Trivial cases, in each case RETURN ------------------  
// input ideal is the ideal of a partial normalization

   // ------------ Trivial case: input ideal contains 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]=list(delt,delti);
         setring BAS;
         return(result);
   }
   
   // --- Trivial case: input ideal is zero-dimensional and homog ---
   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]=list(delt,delti);
      setring BAS;
      return(result);
   }

   // --- Trivial case: input ideal defines a line ---
   //the one-dimensional, homogeneous case and degree 1 case
   if( (dim(SM[1])==1) && (maxdeg1(SM[2])==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]=list(delt,delti);
      setring BAS;
      return(result);
   }

//---------------------- The non-trivial cases start -------------------  
   //the higher dimensional case
   //we test first hypersurface, CohenMacaulay and complete intersection

   if( ((size(SM[2])+dim(SM[1])) == nvars(basering)) )
   {
      //the 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( size(SM[2]) == 1 )
   {
      attrib(SM[2],"isHypersurface",1);
      if(y>=1)
      {
         "// the ideal is a hypersurface";
      }
   }

   //------------------- compute the singular locus -------------------
   // Computation if singular locus is critical
   // Notation: J ideal of singular locus or (if given) containing it
   // JM = mstd(J) or maxideal(1),maxideal(1)
   // JM[1] SB of singular locus, JM[2] minbasis, dimJ = dim(JM[1]) 
   // SM[1] SB of the input ideal i, SM[2] minbasis
   // Computation if singular locus is critical, because it determines the
   // size of the ring Hom_R(J,J). We only need a test ideal contained in J.

   //----------------------- onlySingularAtZero -------------------------
   if( attrib(SM[2],"onlySingularAtZero") )
   {
       JM = maxideal(1),maxideal(1);
       attrib(JM[1],"isSB",1);
       attrib(JM[2],"isRadical",1);
       if( dim(SM[1]) >=2 )
       {
         attrib(SM[2],"isRegInCodim2",1);
       }
   }

   //-------------------- not onlySingularAtZero -------------------------
   if( attrib(SM[2],"onlySingularAtZero") == 0 )
   {
      //--- the case where an ideal #[1] is given:
      if( size(#)>0 )                 
      {
         J = #[1],SM[2];
         JM = mstd(J);
         if( typeof(attrib(#[1],"isRadical"))!="int" )
         {
            attrib(JM[2],"isRadical",0);
         }
      }

      //--- the case where an ideal #[1] is not given:
      if( (size(#)==0) )
      {
         if(y >=1 )
         {
            "// singular locus will be computed";
         }

         J = SM[1],minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]);
         if( y >=1 )
         {
            "// SB of singular locus will be computed";
         }
         JM = mstd(J);
      }

      int dimJ = dim(JM[1]);
      attrib(JM[1],"isSB",1);
      if( y>=1 )
      {  
         "// the dimension of the singular locus is";  dimJ ; "";
      }
         
      if(dim(JM[1]) <= dim(SM[1])-2)
      {
         attrib(SM[2],"isRegInCodim2",1);
      }

      //------------------ the smooth case, RETURN -------------------
      if( dimJ == -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]=list(delt,delti);
         setring BAS;
         return(result);
      }

      //------- extra check for onlySingularAtZero, relatively cheap ----------
      //it uses the procedure 'locAtZero' from for testing 
      //if an ideal is concentrated at 0
       if(y>=1)
       {
         "// extra test for onlySingularAtZero:";
       }
       if ( locAtZero(JM[1]) )  
       {
           attrib(SM[2],"onlySingularAtZero",1);
           JM = maxideal(1),maxideal(1);
           attrib(JM[1],"isSB",1);
           attrib(JM[2],"isRadical",1);
       }
       else
       {
              attrib(SM[2],"onlySingularAtZero",0);
       }
   }
 
  //displaying the attributes:
   if(y>=2)
   {
      "// the attributes of the ideal are:";  
      "// isCohenMacaulay:", attrib(SM[2],"isCohenMacaulay");
      "// isCompleteIntersection:", attrib(SM[2],"isCompleteIntersection");
      "// isHypersurface:", attrib(SM[2],"isHypersurface");
      "// isEquidimensional:", attrib(SM[2],"isEquidimensional");
      "// isPrim:", attrib(SM[2],"isPrim");
      "// isRegInCodim2:", attrib(SM[2],"isRegInCodim2");
      "// isIsolatedSingularity:", attrib(SM[2],"isIsolatedSingularity");
      "// onlySingularAtZero:", attrib(SM[2],"onlySingularAtZero");
      "// isRad:", attrib(SM[2],"isRad");"";
   }

   //------------- case: CohenMacaulay in codim 2, RETURN ---------------
   if( (attrib(SM[2],"isRegInCodim2")==1) &&   
       (attrib(SM[2],"isCohenMacaulay")==1) )
   {
      if(y>=1)
      {
         "// the ideal was CohenMacaulay and regular in codim 2, hence normal";
      }
      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]=list(delt,delti);
      setring BAS;
      return(result);
   }

//---------- case: isolated singularity only at 0, RETURN ------------
   // In this case things are easier, we can use the maximal ideal as radical
   // of the singular locus; 
   // JM mstd of ideal of singular locus, SM mstd of input ideal

   if( attrib(SM[2],"onlySingularAtZero") )
   {
   //------ check variables for being a non zero-divizor ------
   // SL = ideal of vars not contained in ideal SM[1]:

      attrib(SM[2],"isIsolatedSingularity",1);
      ideal SL = simplify(reduce(maxideal(1),SM[1]),2);
      ideal Ann = quotient(SM[2],SL[1]);
      ideal qAnn = simplify(reduce(Ann,SM[1]),2);
      //NOTE: qAnn=0 iff first var (=SL[1]) not in SM is a nzd of R/SM

   //------------- We found a non-zerodivisor of R/SM -----------------------
   // here the enlarging of the ring via Hom_R(J,J) starts

      if( size(qAnn)==0 )
      {
         if(y>=1)
         {
            "";
            "// the ideal rad(J):"; maxideal(1);
            "";
         }

      // ------------- test for normality, compute Hom_R(J,J) -------------
      // Note:
      // HomJJ (ideal SBid, ideal id, ideal J, poly p) with
      //        SBid = SB of id, J = radical ideal of basering  P with: 
      //        nonNormal(R) is in V(J), J contains the nonzero divisor p
      //        of R = P/id (J = test ideal)
      // returns a list l of three objects
      // l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
      //        s.t. l[1]/endid = Hom_R(J,J) and endphi= 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, = dim_K(Hom_R(J,J)/R) if finite, -1 otherwise

         list RR;
         RR = SM[1],SM[2],maxideal(1),SL[1];
         RR = HomJJ(RR,y);
         // --------------------- non-normal case ------------------
         //RR[2]==0 means that the test for normality is negative
         if( RR[2]==0 )
         {
            def newR=RR[1];
            setring newR;
            map psi=BAS,endphi;
            list JM = psi(JM); //###
            ideal J = JM[2];
            if ( delt>=0 && RR[3]>=0 )
            {
               delt = delt+RR[3];
            }
            else
            { delt = -1; }
            delti[size(delti)]=delt;

            // ---------- recursive call of normalizationPrimes -----------
        //normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #)
        //ihp = (partial) normalisation map from basering
        //#[1] ideal s.t. V(#[1]) contains singular locus of i (test ideal)

            if ( y>=1 )
            {
            "// case: onlySingularAtZero, non-zerodivisor found";
            "// contribution of delta in ringextension R -> Hom_R(J,J):"; delt;
            }

            //intvec atr=getAttrib(endid);
            //"//### case: isolated singularity only at 0, recursive";
            //"size endid:", size(endid), size(string(endid));
            //"interred:";
            //endid = interred(endid);
            //endid = setAttrib(endid,atr);
            //"size endid:", size(endid), size(string(endid));

           printlevel=printlevel+1;
           list tluser = 
                normalizationPrimes(endid,psi(ihp),delt,delti); 
           //list tluser = 
           //     normalizationPrimes(endid,psi(ihp),delt,delti,J); 
           //#### ??? improvement: give also the old ideal of sing locus???

           printlevel = printlev;             //reset printlevel
           setring BAS;   
           return(tluser);          
         }

         // ------------------ the normal case, RETURN -----------------
         // Now RR[2] must be 1, hence the test for normality was positive
         MB=SM[2];
         //execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
         //             +ordstr(basering)+");");
         def newR7 = ring(gnirlist);
         setring newR7;
         ideal norid=fetch(BAS,MB);
         ideal normap=fetch(BAS,ihp);
         if ( delt>=0 && RR[3]>=0 )
         {
               delt = delt+RR[3];
         }
         else
         { delt = -1; }
         delti[size(delti)]=delt;

         intvec atr = getAttrib(norid);

         //"//### case: isolated singularity only at 0, final";
         //"size norid:", size(norid), size(string(norid));
         //"interred:";
         //norid = interred(norid);
         //norid = setAttrib(norid,atr);
         //"size norid:", size(norid), size(string(norid));

         export norid;
         export normap;
         result=newR7;
         result[size(result)+1]=list(delt,delti);
         setring BAS;
         return(result);
      }

   //------ zerodivisor of R/SM was found, gives a splitting ------------
   //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,"isHypersurface",0);
          attrib(id,"isCompleteIntersection",0);
          attrib(id,"isEquidimensional",0);
          attrib(id,"onlySingularAtZero",1);

          ideal id1 = quotient(SM[2],Ann);
          attrib(id1,"isCohenMacaulay",0);
          attrib(id1,"isPrim",0);
          attrib(id1,"isIsolatedSingularity",1);
          attrib(id1,"isRegInCodim2",0);
          attrib(id1,"isHypersurface",0);
          attrib(id1,"isCompleteIntersection",0);
          attrib(id1,"isEquidimensional",0);
          attrib(id1,"onlySingularAtZero",1);

          // ---------- recursive call of normalizationPrimes -----------
          if ( y>=1 )
          {
            "// case: onlySingularAtZero, zerodivisor found, splitting:";
            "// total delta before splitting:", delt;
            "// splitting in two components:";
          }

          printlevel = printlevel+1;  //to see comments in normalizationPrimes
          keepresult1 = normalizationPrimes(id,ihp,0,0);   //1st split factor
          keepresult2 = normalizationPrimes(id1,ihp,0,0);  //2nd split factor
          printlevel = printlev;                           //reset printlevel

          int delt1 = keepresult1[size(keepresult1)][1];
          int delt2 = keepresult2[size(keepresult2)][1];
          intvec delti1 = keepresult1[size(keepresult1)][2];
          intvec delti2 = keepresult2[size(keepresult2)][2];

          if( delt>=0 && delt1>=0 && delt2>=0 )  
          {  ideal idid1=id,id1;
             int mul = vdim(std(idid1));
             if ( mul>=0 )
             {
               delt = delt+mul+delt1+delt2;
             }
             else
             {
               delt = -1;
             }            
          }
         if ( y>=1 )
         {
           "// delta of first component:", delt1;
           "// delta of second componenet:", delt2;
           "// intersection multiplicity of both components:", mul;
           "// total delta after splitting:", delt;
         }

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

          return(keepresult1);
       }
   }
   // Case "onlySingularAtZero" has finished and returned result

//-------------- General case, not onlySingularAtZero, RETURN ---------------
   //test for non-normality, i.e. if Hom(I,I)<>R
   //we can use Hom(I,I) to continue

   //------ check variables for being a non zero-divizor ------
   // SL = ideal of vars not contained in ideal SM[1]:

   ideal SL = simplify(reduce(JM[2],SM[1]),2);
   ideal Ann = quotient(SM[2],SL[1]);
   ideal qAnn = simplify(reduce(Ann,SM[1]),2);
   //NOTE: qAnn=0 <==> first var (=SL[1]) not contained in SM is a nzd of R/SM

   //------------- We found a non-zerodivisor of R/SM -----------------------
   //SM = mstd of ideal of variety, JM = mstd of ideal of singular locus

   if( size(qAnn)==0 )
   {
      list RR;
      list RS;
      // ----------------- Computation of the radical -----------------
      if(y>=1)
      {
         "// radical computation of singular locus";
      }
      J = radical(JM[2]);   //the radical of singular locus
      JM = mstd(J);

      if(y>=1)
      {
        "// radical is equal to:";"";  JM[2]; 
        "";
      }
      // ------------ choose non-zerodivisor of smaller degree ----------
      //### evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen ?
      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];
         }
      }

      // --------------- computation of Hom(rad(J),rad(J)) --------------
      RR=SM[1],SM[2],JM[2],SL[1];

     if(y>=1)
     {
        "// compute Hom(rad(J),rad(J))";
     }

     RS=HomJJ(RR,y);               //most important subprocedure
         
     // ------------------ the normal case, RETURN -----------------
     // RS[2]==1 means that the test for normality was positive
     if(RS[2]==1)
     {
         def lastR=RS[1];
         setring lastR;
         map psi1=BAS,endphi;
         ideal norid=endid;
         ideal normap=psi1(ihp);
         kill endid,endphi;

        intvec atr=getAttrib(norid);

        //"//### general case: not isolated singularity only at 0, final";
        //"size norid:", size(norid), size(string(norid));
        //"interred:";
        //norid = interred(norid);
        //norid = setAttrib(norid,atr);
        //"size norid:", size(norid), size(string(norid));

         export norid;
         export normap;
         result=lastR;
         if ( y>=1 )
         {
            "// case: not onlySingularAtZero, last ring Hom_R(J,J) computed";
            "// delta before last ring:", delt;
         }

         if ( delt>=0 && RS[3]>=0 )
         {
            delt = delt+RS[3];
         }
         else
         { delt = -1; }

        // delti = delti,delt;
         delti[size(delti)]=delt;

         if ( y>=1 )
         {
           "// delta of last ring:", delt;
         }

         result[size(result)+1]=list(delt,delti);
         setring BAS;
         return(result);
     }

    // ----- the non-normal case, recursive call of normalizationPrimes -------
    // RS=HomJJ(RR,y) was computed above, RS[1] contains endid and endphi
    // RS[1] = new ring Hom_R(J,J), RS[2]= 0 or 1, RS[2]=contribution to delta
    // now RS[2]must be 0, i.e. the test for normality was negative

      int n = nvars(basering);
      ideal MJ = JM[2];

      def newR=RS[1];
      setring newR;
      map psi=BAS,endphi;
      if ( y>=1 )
      {
        "// case: not onlySingularAtZero, compute new ring = Hom_R(J,J)";
        "// delta of old ring:", delt;
      }
      if ( delt>=0 && RS[3]>=0 )
      {
         delt = delt+RS[3];
      }
      else
      { delt = -1; }
      if ( y>=1 )
      {
        "// delta of new ring:", delt;
      }

      delti[size(delti)]=delt;
      intvec atr=getAttrib(endid);

      //"//### general case: not isolated singularity only at 0, recursive";
      //"size endid:", size(endid), size(string(endid));
      //"interred:";
      //endid = interred(endid);
      //endid = setAttrib(endid,atr);
      //"size endid:", size(endid), size(string(endid));

      printlevel = printlevel+1;
      list tluser=
          normalizationPrimes(endid,psi(ihp),delt,delti,psi(MJ));
      printlevel = printlev;                //reset printlevel
      setring BAS;
      return(tluser);
   }

   //---- A whole singular component was found, RETURN -----
   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]=lst(delt,delti);
         setring BAS;
         return(result);
   }

   //------ zerodivisor of R/SM was found, gives a splitting ------------
   //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM
   //and we have found a splitting: new1 and new2
   //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
   {
      if(y>=1)
      {
         "// zero-divisor found";
      }
      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];

      //execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
      //                +ordstr(basering)+");");
      def newR1 = ring(gnirlist);
      setring newR1;

      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,"isHypersurface",0);
      attrib(vid,"isCompleteIntersection",0);

      // ---------- recursive call of normalizationPrimes -----------
      if ( y>=1 )
      {
        "// total delta before splitting:", delt;
        "// splitting in two components:";
      }
      printlevel = printlevel+1;
      keepresult1 = 
                  normalizationPrimes(vid,ihp,0,0);  //1st split factor

      list delta1 = keepresult1[size(keepresult1)];

      setring BAS;
      //execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
      //                +ordstr(basering)+");");
      def newR2 = ring(gnirlist);
      setring newR2;

      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,"isHypersurface",0);
      attrib(vid,"isCompleteIntersection",0);
      attrib(vid,"onlySingularAtZero",oSAZ);

      keepresult2 = 
                    normalizationPrimes(vid,ihp,0,0);
      list delta2 = keepresult2[size(keepresult2)];   //2nd split factor
      printlevel = printlev;                          //reset printlevel
 
      setring BAS;

      //compute intersection multiplicity of both components:
      new1 = new1,new2;
      int mul=vdim(std(new1));

     // ----- normalizationPrimes finished, add up results, RETURN --------
      for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
      {
         keepresult1 = insert(keepresult1,keepresult2[lauf]);
      }
      if ( delt >=0 && delta1[1] >=0 && delta2[1] >=0 && mul >=0 )
      {
         delt = delt+mul+delta1[1]+delta2[1];
      }
      else
      {  delt = -1; }
      delti = -2;

      if ( y>=1 )
      {
        "// zero divisor produced a splitting into two components";
        "// delta of first component:", delta1;
        "// delta of second componenet:", delta2;
        "// intersection multiplicity of both components:", mul;
        "// total delta after splitting:", delt;
      }
      keepresult1[size(keepresult1)] = list(delt,delti);
      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";

//NOTE concerning iteration of maps: Let phi: x->f(y,z), y->g(x,z) then
//phi: x+y+z->f(y,z)+g(x,z)+z, phi(phi):x+y+z->f(g(x,z),z)+g(f(y,z),z)+z
//and so on: none of the x or y will be eliminated
//Now subst: first x and then y: x+y+z->f(g(x,z),z)+g(x,z)+z eliminates y
//further subst replaces x by y, makes no sense (objects more compicated).
//Subst first y and then x eliminates x
//In our situation we have triangular form: x->f(y,z), y->g(z).
//phi: x+y+z->f(y,z)+g(z)+z, phi(phi):x+y+z->f(g(z),z)+g(z)+z eliminates x,y
//subst x,y: x+y+z->f(g(z),z)+g(z)+z, eliminates x,y
//subst y,x: x+y+z->f(y,z)+g(z)+z eliminates only x
//HENCE: substitute vars depending on most other vars first
//However, if the sytem xi-fi is reduced then xi does not appear in any of the 
//fj and hence the order does'nt matter when substitutinp xi by fi

{
   def newRing = basering;
   int ii,jj;
   map phi = newRing,maxideal(1);    //identity map
   list Le = elimpart(endid);
   //this proc and the next loop try to substitute as many variables as 
   //possible indices of substituted variables

   int q = size(Le[2]);    //q vars, stored in Le[2], have been substitutet
   intvec rw1 = 0;         //will become indices of substituted variables
   rw1[nvars(basering)] = 0;
   rw1 = rw1+1;            //rw1=1,..,1 (as many 1 as nvars(basering))

   while( size(Le[2]) != 0 )
   {
      endid = Le[1];
      if ( defined(ps) )
      { kill ps; }
      map ps = newRing,Le[5];
      phi = ps(phi);
      for(ii=1;ii<=size(Le[2]);ii++)
      {
         phi=phi(phi);
      }
      //eingefuegt wegen x2-y2z2+z3

      for( ii=1; ii<=size(rw1); ii++ )
      {
         if( Le[4][ii]==0 )        //ii = index of var which was substituted
         {
            rw1[ii]=0;             //substituted vars have entry 0 in rw1
         }
      }
      Le=elimpart(endid);          //repeated application of elimpart
      q = q + size(Le[2]);
   }
   endphi = phi(endphi);
//---------- return -----------------------------------------------------------
// first the trivial case, where all variable have been eliminated
   if( nvars(newRing) == q )
   {
     ring lastRing = char(basering),T(1),dp;
     ideal endid = T(1);
     ideal endphi = T(1);
     for(ii=2; ii<=q; ii++ )
     {
        endphi[ii] = 0;
     }
     export(endid,endphi);
     list L = lastRing;
     setring newRing;
     return(L);
   }

// 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;
   jj = 1;

   for(ii=1; ii<=size(rw1); ii++ )
   {
      if ( rw1[ii]==1 ) { lastmap[ii] = T(jj); jj=jj+1; }
      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
   }
   map phi1 = newRing,lastmap;
   ideal endid  = phi1(endid);      //### bottelneck
   ideal endphi = phi1(endphi);

/*
Versuch: subst statt phi
   for(ii=1; ii<=size(rw1); ii++ )
   {
      if ( rw1[ii]==1 ) { endid = subst(endid,var(ii),T(jj)); }
      if ( rw1[ii]==0 ) { endid = subst(endid,var(ii),0); }
   }
*/
   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 genus(i,1)
         is slower than genus(i) but sometimes not.
EXAMPLE: example genus; shows an example
"
{
   int w = printlevel-voice+2;  // w=printlevel (default: w=0)

   def R0=basering;
   I=std(I);
   if(dim(I)!=1)
   {
      if(((homog(I))&&(dim(I)!=2))||(!homog(I)))
      {
        // ERROR("This is not a curve");
        if(w==1){"** WARNING: Input does not define a curve **"; "";}
      }
   }
   list L=elimpart(I);
   if(size(L[2])!=0)
   {
      map psi=R0,L[5];
      I=std(psi(I));
   }
   if(size(I)==0)
   {
      return(0);
   }
   list N=findvars(I,0);
   if(size(N[1])==1)
   {

      poly p=I[1];
     // if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");}
      return(-deg(p)+1);
   }
   if(size(N[1]) < nvars(R0))
   {
      string newvar=string(N[1]);
      execute("ring R=("+charstr(R0)+"),("+newvar+"),dp;");
      ideal I =imap(R0,I);
      attrib(I,"isSB",1);
   }
   else
   {
      def R=basering;
   }
   if(dim(I)==2)
   {
      def newR=basering;
   }
   else
   {
      if(dim(I)==0)
      {
         execute("ring Rhelp=("+charstr(R0)+"),(@s,@t),dp;");
      }
      else
      {
         execute("ring Rhelp=("+charstr(R0)+"),(@s),dp;");
      }
      def newR=R+Rhelp;
      setring newR;
      ideal I=imap(R,I);
      I=homog(I,@s);
      attrib(I,"isSB",1);
   }

   if((nvars(basering)<=3)&&(size(I)>1))
   {
       ERROR("This is not equidimensional");
   }

   intvec hp=hilbPoly(I);
   int p_a=1-hp[1];
   int d=hp[2];

   if(w>=1)
   {
      "";"The ideal of the projective curve:";"";I;"";
      "The coefficients of the Hilbert polynomial";hp;
      "arithmetic genus:";p_a;
      "degree:";d;"";
   }

   intvec v = hilb(I,1);
   int i,o;
   if(nvars(basering)>3)
   {
      map phi=newR,maxideal(1);
      int de;
      ideal K,L1;
      matrix M;
      poly m=var(4);
      poly he;
      for(i=5;i<=nvars(basering);i++){m=m*var(i);}
      K=eliminate(I,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(I);
         K=eliminate(K,m,v);
         if(size(K)==1){de=deg(K[1]);}
         if((i==5)&&(d!=de))
         {
            K=reduce(equidimMax(I),I);
            if(size(K)!=0){ERROR("This is not equidimensional");}
         }
         if(i==10)
         {
            J;K;
            ERROR("genus: did not find a good projection for to
                           the plain");
         }
         if(i<5)
         {
            M=sparsetriag(nvars(newR),nvars(newR),80-5*i,i);
         }
         else
         {
            if(i<8)
            {
               M=transpose(sparsetriag(nvars(newR),nvars(newR),80-5*i,i));
            }
            else
            {
               he=0;
               while(he==0)
               {
                  M=randommat(nvars(newR),nvars(newR),ideal(1),20);
                  he=det(M);
               }
            }
         }
         L1=M*transpose(maxideal(1));
         phi=newR,L1;
      }
      I=K;
   }
   poly p=I[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++;}
   if(inv[2]==2){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);
      delt=nor[size(nor)][2];
      genus=genus-delt-deltainf;
      setring R0;
      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 R0;
   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=simplify(std(singL),1);
   attrib(singL,"isSB",1);
   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;
      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
      // 'hnexpansion' (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];
         setring R;
         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";}
      ideal fac=factorize(f,1);
      if(size(fac)>1)
      {
         nb=0;
         for(i=1;i<=size(fac);i++)
         {
            nb=nb+deltaLoc(fac[i],maxideal(1))[3];
         }
         setring R;
         return(list(d*(mu+nb-1)/2,d*tau,d*nb));
      }
      list HNEXP=hnexpansion(f);
      if (typeof(HNEXP[1])=="ring") {
        def altring = basering;
        def HNEring = HNEXP[1]; setring HNEring;
        nb=size(hne);
        setring R;
        kill HNEring;
      }
      else
      {
        nb=size(HNEXP);
      }
      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);
//}


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

proc primeClosure (list L, list #)
"USAGE:    primeClosure(L [,c]); L a list of a ring containing a prime ideal
                                 ker, c an optional integer
RETURN:   a list L (of size n+1) consisting of rings L[1],...,L[n] such that
          - L[1] is a copy of (not a reference to!) the input ring L[1]
          - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi
            such that
                    L[1]/ker --> ... --> L[n]/ker
            are injections given by the corresponding ideals phi, and L[n]/ker
            is the integral closure of L[1]/ker in its quotient field.
          - all rings L[i] contain a polynomial nzd such that elements of
            L[i]/ker are quotients of elements of L[i-1]/ker with denominator
            nzd via the injection phi.
            L[n+1] is the delta invariant 
NOTE:     - L is constructed by recursive calls of primeClosure itself.
          - c determines the choice of nzd:
               - c not given or equal to 0: first generator of the ideal SL,
                 the singular locus of Spec(L[i]/ker)
               - c<>0: the generator of SL with least number of monomials.
EXAMPLE:  example primeClosure; shows an example
"
{
  //---- Start with a consistency check:

  if (!(typeof(L[1])=="ring"))
  {
      "// Parameter must be a ring or a list containing a ring!";
      return(-1);
  }

  int dblvl = printlevel-voice+2;
  list gnirlist = ringlist(basering);

  //---- Some auxiliary variables:
  int delt;                      //finally the delta invariant
  if ( size(L) == 1 )
  {
      L[2] = delt;              //set delta to 0
  }
  int n = size(L)-1;            //L without delta invariant

  //---- How to choose the non-zerodivisor later on?

  int nzdoption=0;
  if (size(#)>0)
  {
      nzdoption=#[1];
  }

// R0 is the ring to work with, if we are in step one, make a copy of the
// input ring, so that all objects are created in the copy, not in the original
// ring (therefore a copy, not a reference is defined).

  if (n==1)
  {
      def R = L[1];
      list Rlist = ringlist(R);
      def BAS = basering;
      setring R;
      if (!(typeof(ker)=="ideal"))
      {
          "// No ideal ker in the input ring!";
          return (-1);
      }
      ker=simplify(interred(ker),15);
      //execute ("ring R0="+charstr(R)+",("+varstr(R)+"),("+ordstr(R)+");");
      def R0 = ring(Rlist);
      setring R0;
      ideal ker=fetch(R,ker);

      // check whether we compute the normalization of the blow up of
      // an isolated singularity at the origin (checked in normalI)

      if (typeof(attrib(L[1],"iso_sing_Rees"))=="int")
      {
        attrib(R0,"iso_sing_Rees",attrib(L[1],"iso_sing_Rees"));
      }
      L[1]=R0;
  }
  else
  {
      def R0 = L[n];
      setring R0;
  }


// In order to apply HomJJ from normal.lib, we need the radical of the singular
// locus of ker, J:=rad(ker):

   list SM=mstd(ker); 

// In the first iteration, we have to compute the singular locus "from
// scratch". 
// In further iterations, we can fetch it from the previous one but
// have to compute its radical
// the next rings R1 contain already the (fetched) ideal

  if (n==1)                              //we are in R0=L[1]
  {
      if (typeof(attrib(R0,"iso_sing_Rees"))=="int")
      {
        ideal J;
        for (int s=1;s<=attrib(R0,"iso_sing_Rees");s++)
        {
          J=J,var(s);
        }
        J = J,SM[2];
        list JM = mstd(J);
      }
      else
      {
        if ( dblvl >= 1 )
        {"";
           "// compute the singular locus";
        }
        //### Berechnung des singulŠren Orts geŠndert (ist so schneller)
        ideal J = minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]);
        J = J,SM[2];
        list JM = mstd(J);
      }

      if ( dblvl >= 1 )
      {"";
         "// dimension of singular locus is", dim(JM[1]);
         if (  dblvl >= 2 )
         {"";
            "// the singular locus is:"; JM[2];
         }
      }

      if ( dblvl >= 1 )
      {"";
         "// compute radical of singular locus";
      }

      J = simplify(radical(JM[2]),2);
      if ( dblvl >= 1 )
      {"";
         "// radical of singular locus is:"; J;
         pause();
      }
  }
  else
  {
      if ( dblvl >= 1 )
      {"";
         "// compute radical of test ideal in ideal of singular locus";
      }
      J = simplify(radical(J),2);
      if ( dblvl >= 1 )
      {"";
         "// radical of test ideal is:"; J;
         pause();
      }
  }

  // having computed the radical J of/in the ideal of the singular locus,
  // we now need to pick an element nzd of J; 
  // NOTE: nzd must be a non-zero divisor mod ker, i.e. not contained in ker

  poly nzd = J[1];
  poly nzd1 = NF(nzd,SM[1]);
  if (nzd1 != 0)
  {
     if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) )
     {
        nzd = nzd1;
     }
  }

  if (nzdoption || nzd1==0)
  {
    for (int ii=2;ii<=ncols(J);ii++)
    {
      nzd1 = NF(J[ii],SM[1]);
      if ( nzd1 != 0 )
      {
        if ( (deg(nzd)>=deg(J[ii])) && (size(nzd)>size(J[ii])) )
        {
          nzd=J[ii];
        }
        if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) )
        {
          nzd = nzd1;
        }
      }
    }
  }

  export nzd;
  list RR = SM[1],SM[2],J,nzd;

  if ( dblvl >= 1 )
  {"";
     "// compute the first ring extension:";
  }

  list RS = HomJJ(RR);

//-------------------------------------------------------------------------
// If we've reached the integral closure (as determined by the result of
// HomJJ), then we are done, otherwise we have to prepare the next iteration.

  if (RS[2]==1)     // we've reached the integral closure, we are still in R0
    {
      kill J;
      if ( n== 1)
      {
        def R1 = RS[1];
        setring R1;
        ideal norid, normap = endid, endphi;
        kill endid,  endphi;

        //"//### case: primeClosure, final";
        //"size norid:", size(norid), size(string(norid));
        //"interred:";
        //norid = interred(norid);
        //"size norid:", size(norid), size(string(norid));

        export (norid, normap);
        L[1] = R1;
      }
      return(L);
    }
  else                        // prepare the next iteration
    {
      if (n==1)               // In the first iteration: keep only the data
      {                       // needed later on.
         kill RR,SM;
         export(ker);
      }
      if ( dblvl >= 1 )
      {"";
         "// computing the next ring extension, we are in loop"; n+1;
      }

      def R1 = RS[1];         // The data of the next ring R1:
      delt = RS[3];           // the delta invariant of the ring extension
      setring R1;             // keep only what is necessary and kill
      ideal ker=endid;        // everything else.
      export(ker);
      ideal norid=endid;

      //"//### case: primeClosure, loop", n+1;
      //"size norid:", size(norid), size(string(norid));
      //"interred:";
      //norid = interred(norid);        //????
      //"size norid:", size(norid), size(string(norid));

      export(norid);
      kill endid;

      map phi = R0,endphi;                        // fetch the singular locus
      ideal J = mstd(simplify(phi(J)+ker,4))[2];  // ideal J in R1
      export(J);
      if(n>1)
      {
         ideal normap=phi(normap);
      }
      else
      {
         ideal normap=endphi;
      }
      export(normap);
      kill phi;              // we save phi as ideal, not as map, so that
      ideal phi=endphi;      // we have more flexibility in the ring names
      kill endphi;           // later on.
      export(phi);
      L=insert(L,R1,n);       // Add the new ring R1 and go on with the
                              // next iteration
      if ( L[size(L)] >= 0 && delt >= 0 )
      {
         delt = L[size(L)] + delt;
      }
      else
      {
         delt = -1;
      }      
      L[size(L)] = delt;

      if (size(#)>0)
      {
          return (primeClosure(L,#));
      }
      else
      {
          return(primeClosure(L));         // next iteration.
      }
    }
}
example
{
  "EXAMPLE:"; echo=2;
  ring R=0,(x,y),dp;
  ideal I=x4,y4;
  def K=ReesAlgebra(I)[1];        // K contains ker such that K/ker=R[It]
  list L=primeClosure(K);
  def R(1)=L[1];                  // L[4] contains ker, L[4]/ker is the
  def R(4)=L[4];                  // integral closure of L[1]/ker
  setring R(1);
  R(1);
  ker;
  setring R(4);
  R(4);
  ker;
}

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

proc closureFrac(list L)
"USAGE:    closureFrac (L); L a list of size n+1 as in the result of
          primeClosure, L[n] contains an additional poly f
CREATE:   a list fraction of two elements of L[1], such that
          f=fraction[1]/fraction[2] via the injections phi L[i]-->L[i+1].
EXAMPLE:  example closureFrac; shows an example
"
{
// Define some auxiliary variables:

  int n=size(L)-1;
  int i,j,k,l;
  string mapstr;
  for (i=1; i<=n; i++) { def R(i) = L[i]; }

// The quotient representing f is computed as in 'closureGenerators' with
// the differences that
//   - the loop is done twice: for the numerator and for the denominator;
//   - the result is stored in the list fraction and
//   - we have to make sure that no more objects of the rings R(i) survive.

  for (j=1; j<=2; j++)
    {
      setring R(n);
      if (j==1)
      {
         poly p=f;
      }
      else
      {
         p=1;
      }

      for (k=n; k>1; k--)
      {
          if (j==1)
          {
             map phimap=R(k-1),phi;
          }

          p=p*phimap(nzd);

          if (j==2)
          {
            kill phimap;
          }

          if (j==1)
          {
             //### noch abfragen ob Z(i) definiert ist
             list gnirlist = ringlist(R(k));
             int n2 = size(gnirlist[2]);
             int n3 = size(gnirlist[3]);
             for( i=1; i<=ncols(phi); i++)
             {
               gnirlist[2][n2+i] = "Z("+string(i)+")";
             }
             intvec V;
             V[ncols(phi)]=0; V=V+1;
             gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
             def S(k) = ring(gnirlist);
             setring S(k);

             //execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
             //          Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
             //          +"),dp("+string(ncols(phi))+"));");

              ideal phi = imap(R(k),phi);
              ideal J = imap (R(k),ker);
              for (l=1;l<=ncols(phi);l++)
              {
                  J=J+(Z(l)-phi[l]);
              }
              J=groebner(J);
              poly h=NF(imap(R(k),p),J);
          }
          else
          {
              setring S(k);
              h=NF(imap(R(k),p),J);
              setring R(k);
              kill p;
          }

          setring R(k-1);

          if (j==1)
          {
              ideal maxi;
              maxi[nvars(R(k))] = 0;
              maxi = maxi,maxideal(1);
              map backmap = S(k),maxi;
         
              //mapstr=" map backmap = S(k),";
              //for (l=1;l<=nvars(R(k));l++)
              //{
              //  mapstr=mapstr+"0,";
              //}
              //execute (mapstr+"maxideal(1);");
              poly p;
          }
          p=NF(backmap(h),std(ker));
          if (j==2)
          {
            kill backmap;
          }
        }

      if (j==1)
        {
          if (defined(fraction))
            {
              kill fraction;
              list fraction=p;
            }
          else
            {
              list fraction=p;
            }
        }
      else
        {
          fraction=insert(fraction,p,1);
        }
    }
  export(fraction);
  return ();
}
example
{
  "EXAMPLE:"; echo=2;
  ring R=0,(x,y),dp;
  ideal ker=x2+y2;
  export ker;
  list L=primeClosure(R);          // We normalize R/ker
  for (int i=1;i<=size(L);i++) { def R(i)=L[i]; }
  setring R(2);
  kill R;
  phi;                             // The map R(1)-->R(2)
  poly f=T(2);                     // We will get a representation of f
  export f;
  L[2]=R(2);
  closureFrac(L);
  setring R(1);
  kill R(2);
  fraction;                        // f=fraction[1]/fraction[2] via phi
  kill R(1);
}

///////////////////////////////////////////////////////////////////////////////
// closureGenerators is called inside proc normal (option "withGens" )
// 

// INPUT is the output of proc primeClosure (except for the last element, the 
// delta invariant) : hence input is a list L consisting of rings 
// L[1],...,L[n] (denoted R(1)...R(n) below) such that
// - L[1] is a copy of (not a reference to!) the input ring L[1]
// - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi 
// such that
//                L[1]/ker --> ... --> L[n]/ker
// are injections given by the corresponding ideals phi, and L[n]/ker
// is the integral closure of L[1]/ker in its quotient field.
// - all rings L[i] contain a polynomial nzd such that elements of
// L[i]/ker are quotients of elements of L[i-1]/ker with denominator
// nzd via the injection phi.

// COMPUTE: In the list L of rings R(1),...,R(n), compute representations of
// the ring variables of the last ring R(n) as fractions of elements of R(1):
// The proc returns an ideal preim s.t. preim[i]/preim[size(preim)] expresses
// the ith variable of R(n) as fraction of elements of the basering R(1)
// preim[size(preim)] is a non-zero divisor of basering/i.

proc closureGenerators(list L);    
{
  def Rees=basering;         // when called inside normalI (in reesclos.lib) 
                             // the Rees Algebra is the current basering

  // ------- First of all we need some variable declarations -----------
  int n = size(L);                // the number of rings R(1)-->...-->R(n)
  int length = nvars(L[n]);       // the number of variables of the last ring
  int j,k,l;
  string mapstr;   
  list preimages;                 
  //Note: the empty list belongs to no ring, hence preimages can be used 
  //later in R(1)
  //this is not possible for ideals (belong always to some ring)

  for (int i=1; i<=n; i++) 
  { 
     def R(i)=L[i];          //give the rings from L a name
  }

  // For each variable (counter j) and for each intermediate ring (counter k):
  // Find a preimage of var_j*phi(nzd(k-1)) in R(k-1).
  // Finally, do the same for nzd.

  for (j=1; j <= length+1; j++ )
  {
      setring R(n);
      if (j==1)
      {
        poly p;
      }
      if (j <= length )
      {
        p=var(j);
      }
      else
      {
        p=1;
      }
      //i.e. p=j-th var of R(n) for j<=length and p=1 for j=length+1

      for (k=n; k>1; k--)
      {

        if (j==1)
        {
          map phimap=R(k-1),phi;   //phimap:R(k-1)-->R(n), k=2..n, is the map
                                   //belonging to phi in R(n)
        }                             

        p = p*phimap(nzd);

          // Compute the preimage of [p mod ker(k)] under phi in R(k-1):
          // As p is an element of Image(phi), there is a polynomial h such
          // that h is mapped to [p mod ker(k)], and h can be computed as the
          // normal form of p w.r.t. a Groebner basis of
          // J(k) := <ker(k),Z(l)-phi(k)(l)> in R(k)[Z]=:S(k)

        if (j==1)   // In the first iteration: Create S(k), fetch phi and
                    // ker(k) and construct the ideal J(k).
        {
         //### noch abfragen ob Z(i) definiert ist
         list gnirlist = ringlist(R(k));
         int n2 = size(gnirlist[2]);
         int n3 = size(gnirlist[3]);
         for( i=1; i<=ncols(phi); i++)
         {
            gnirlist[2][n2+i] = "Z("+string(i)+")";
         }
         intvec V;
         V[ncols(phi)]=0;
         V=V+1;
         gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
         def S(k) = ring(gnirlist);
         setring S(k);

        // execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
        //           Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
        //           +"),dp("+string(ncols(phi))+"));");

          ideal phi = imap(R(k),phi);
          ideal J = imap (R(k),ker);
          for ( l=1; l<=ncols(phi); l++ )
          {
             J=J+(Z(l)-phi[l]);
          }
          J = groebner(J);
          poly h = NF(imap(R(k),p),J);
        }
        else
        {
           setring S(k);
           h = NF(imap(R(k),p),J);
        }

        setring R(k-1);

        if (j==1)  // In the first iteration: Compute backmap:S(k)-->R(k-1)
        {
           ideal maxi;
           maxi[nvars(R(k))] = 0;
           maxi = maxi,maxideal(1);
           map backmap = S(k),maxi;
         
           //mapstr=" map backmap = S(k),";
           //for (l=1;l<=nvars(R(k));l++)
           //{
           //  mapstr=mapstr+"0,";
           //}
           //execute (mapstr+"maxideal(1);");

           poly p;
        }
        p = NF(backmap(h),std(ker));
     }
     // Whe are down to R(1), store here the result in the list preimages
     preimages = insert(preimages,p,j-1);
  }
  ideal preim;                  //make the list preimages to an ideal preim
  for ( i=1; i<=size(preimages); i++ )
  {
     preim[i] = preimages[i];
  }
  // R(1) was a copy of Rees, so we have to get back to the basering Rees from
  // the beginning and fetch the result (the ideal preim) to this ring.
  setring Rees;
  return (fetch(R(1),preim));
}

///////////////////////////////////////////////////////////////////////////////
//                From here: procedures for char p with Frobenius
///////////////////////////////////////////////////////////////////////////////

proc normalP(ideal id,list #)
"USAGE:   normalP(id [,choose]); id a radical ideal, choose a comma separated 
         list of optional strings: \"withRing\" or \"isPrim\" or \"noFac\".@*
         If choose = \"noFac\", factorization is avoided during the computation
         of the minimal associated primes; choose = \"isPrim\" disables the 
         computation of the minimal associated primes (should only be used
         if the ideal is known to be prime).@*
ASSUME:  The characteristic of the ground field must be positive. The ideal
         id is assumed to be radical. However, if choose does not contain
         \"isPrim\" the minimal associated primes of id are computed first
         and hence normal computes the normalization of the radical of id.
         If choose = \"isPrim\" the ideal must be a prime ideal (this is not
         tested) otherwise the result may be wrong.
RETURN:  a list, say 'nor' of size 2 (default) or, if \"withRing\" is given,
         of size 3. @*         
         nor[1] (resp. nor[2] if \"withRing\" is given) is a list of ideals 
         Ii, i=1..r, in the basering where r is the number of associated prime 
         ideals P_i (irreducible components) of id.@*
         - Ii is an ideal given by polynomials g_1,...,g_k,g_k+1 such that 
           g_1/g_k+1,...,g_k/g_k+1 generate the integral closure of
           basering/P_i as module (over the basering) in its quotient field.@*
         
         nor[2] (resp. nor[3] if choose=\"withRing\") is a list of an intvec 
         of size r, the delta invariants of the r components, and an integer, 
         the total delta invariant of basering/id (-1 means infinite, and 0 
         that basering/P_i resp. basering/input is normal). @* 

         If the optional string \"withRing\" is given, the ring structure of
         the normalization is computed too and nor[1] is a list of r rings.@*
         Each ring Ri = nor[1][i], i=1..r, contains two ideals with given names
         @code{norid} and @code{normap} such that @*
         - Ri/norid is the normalization of the i-th rime component P_i, i.e.
           isomorphic to the integral closure of basering/P_i in its field 
           of fractions; @*
         - the direct sum of the rings Ri/norid is the normalization 
           of basering/id; @*
         - @code{normap} gives the normalization map from basering/P_i to
           Ri/norid (for each i).@*
THEORY:  The delta invariant of a reduced ring K[x1,...,xn]/id is 
             dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id)
         We call this number also the delta invariant of id. normalP uses 
         the qth-power algorithm suggested by Leonard-Pellikaan (using the 
         Frobenius) in part similair to an implementation by Singh-Swanson. 
NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
@*       Increasing/decreasing printlevel displays more/less comments 
         (default: printlevel = 0).
@*       Not implemented for local or mixed orderings or quotient rings.
@*       If the input ideal id is weighted homogeneous a weighted ordering may
         be used (qhweight(id); computes weights).
@*       works only in characteristic p > 0.
KEYWORDS: normalization; integral closure; delta invariant.
SEE ALSO: normal
EXAMPLE: example normalP; shows an example
"
{
   int i,j,y, wring, isPrim, noFac, sr, del, co;;
   intvec deli;
   list resu, Resu, prim, Gens, mstdid;
   ideal gens;
   y = printlevel-voice+2;

   if ( ord_test(basering) != 1)
   {
     "";
     "// Not implemented for this ordering,";
     "// please change to global ordering!";
     return(resu);
   }
   if ( char(basering) <= 0)
   {
     "";
     "// Algorithm works only in positive characteristic,";
     "// use procedure 'normal' if the characteristic is 0";
     return(resu);
   }

//--------------------------- define the method --------------------------- 
   string method;                //make all options one string in order to use
                                 //all combinations of options simultaneously
   for ( i=1; i<= size(#); i++ )
   { 
     if ( typeof(#[i]) == "string" ) 
     {
       method = method + #[i];
     }
   }

   if ( find(method,"withring") or find(method,"withRing") )
   {
     wring=1;
   }
   if ( find(method,"isPrim") or find(method,"isprim") )
   {
     isPrim=1;
   }
   if ( find(method,"noFac") )
   {
     noFac=1;
   }

   kill #;
   list #;
//--------------------------- start computation --------------------------- 
   ideal II,K1,K2;

   //----------- check first (or ignore) if input id is prime -------------

   if ( isPrim )
   {
      prim[1] = id;
      if( y >= 0 )
      { "";
    "// ** WARNING: result is correct if ideal is prime (not checked) **";
    "// disable option \"isPrim\" to decompose ideal into prime components";"";
      }
   }
   else
   {
      if(y>=1)
      {  "// compute minimal associated primes"; }

      if( noFac )
      { prim = minAssGTZ(id,1); }
      else
      { prim = minAssGTZ(id); }

      if(y>=1)
      {  
         prim;"";
         "// number of irreducible components is", size(prim);
      }
   }

   //----------- compute integral closure for every component -------------

      for(i=1; i<=size(prim); i++)
      {
         if(y>=1)
         {
            ""; pause(); "";
            "// start computation of component",i;
            "   --------------------------------";
         }
         if(y>=1)
         {  "// compute SB of ideal"; 
         }
         mstdid = mstd(prim[i]);
         if(y>=1)
         {  "// dimension of component is", dim(mstdid[1]);"";}

      //------- 1-st main subprocedure: compute module generators ----------
         printlevel = printlevel+1;
         II = normalityTest(mstdid); 

      //------ compute also the ringstructure if "withRing" is given -------
         if ( wring )
         {
         //------ 2-nd main subprocedure: compute ring structure -----------
            resu = list(computeRing(II,prim[i])) + resu;
         }
         printlevel = printlevel-1;

      //----- rearrange module generators s.t. denominator comes last ------
         gens=0;
         for( j=2; j<=size(II); j++ )
         {
            gens[j-1]=II[j];
         }
         gens[size(gens)+1]=II[1];
         Gens = list(gens) + Gens;
      //------------------------------ compute delta -----------------------
         K1 = mstdid[1]+II;
         K1 = std(K1);
         K2 = mstdid[1]+II[1];
         K2 = std(K2);       
         // K1 = std(mstdid[1],II);      //### besser 
         // K2 = std(mstdid[1],II[1]);   //### besser: Hannes, fixen!
         co = codim(K1,K2);
         deli = co,deli;
         if ( co >= 0 && del >= 0 )
         {
            del = del + co;
         }
         else
         { del = -1; }
      }

      if ( del >= 0 )
      {
         int mul = iMult(prim);
         del = del + mul;
      }
      else
      { del = -1; }

      deli = deli[1..size(deli)-1];
      if ( wring )
      { Resu = resu,Gens,list(deli,del); }
      else
      { Resu = Gens,list(deli,del); }

   sr = size(prim);

//-------------------- Finally print comments and return --------------------
   if(y >= 0)
   {"";
     if ( wring )
     {
"// 'normalP' created a list, say nor, of three lists: 
// nor[1] resp. nor[2] are lists of",sr,"ring(s) resp. ideal(s)
// and nor[3] is a list of an intvec and an integer.
// To see the result, type
     nor;
// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
     def R1 = nor[1][1]; setring R1;  norid; normap;
// for the other rings type first setring r; (if r is the name of your 
// original basering) and then continue as for the first ring;
// Ri/norid is the affine algebra of the normalization of the i-th 
// component r/P_i (where P_i is an associated prime of the input ideal)
// and normap the normalization map from r to Ri/norid;
//    nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
// elements g1..gk of r such that the gj/gk generate the integral 
// closure of r/P_i as (r/P_i)-module in the quotient field of r/P_i.
//    nor[3] shows the delta-invariant of each component and of the input
// ideal (-1 means infinite, and 0 that r/P_i is normal).";
     }
     else
     {
"// 'normalP' computed a list, say nor, of two lists:
//    nor[1] is a list of",sr,"ideal(s), where each ideal nor[1][i] consists
// of elements g1..gk of the basering such that gj/gk generate the integral
// closure of the i-th component as (basering/P_i)-module in the quotient
// field of basering/P_i (P_i an associated prime of the input ideal);
//    nor[2] shows the delta-invariant of each component and of the input ideal 
// (-1 means infinite, and 0 that r/P_i is normal).";
     }
   }

   return(Resu);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r  = 11,(x,y,z),wp(2,1,2);
   ideal i = x*(z3 - xy4 + x2);
   list nor= normalP(i); nor;
   //the result says that both components of i are normal, but i itself
   //has infinite delta
   pause("hit return to continue");

   ring s = 2,(x,y),dp;
   ideal i = y*((x-y^2)^2 - x^3);
   list nor = normalP(i,"withRing"); nor;

   def R2  = nor[1][2]; setring R2;
   norid; normap;
}

///////////////////////////////////////////////////////////////////////////////
// Assume: mstdid is the result of mstd(prim[i]), prim[i] a prime component of
// the input ideal id of normalP.
// Output is an ideal U s.t. U[i]/U[1] are module generators.
 
static proc normalityTest(list mstdid)
{
   int y = printlevel-voice+2;
   intvec op = option(get);
   option(redSB);
   def R = basering;
   int n, p = nvars(R), char(R);
   int ii;

   ideal J = mstdid[1];         //J is the SB of I
   ideal I = mstdid[2];
   int h = n-dim(J);            //codimension of V(I), I is a prime ideal

   //-------------------------- compute singular locus ----------------------
   qring Q = J;                 //pass to quotient ring
   ideal I = imap(R,I);
   ideal J = imap(R,J);
   attrib(J,"isSB",1);
   if ( y >= 1)
   { "start normality test";  "compute singular locus";}

   ideal M = minor(jacob(I),h,J); //use the command minor modulo J (hence J=0)
   M = std(M);                    //this makes M much smaller
   //keep only minors which are not 0 mod I (!) this is important since we 
   //need a nzd mod I

   //---------------- choose nzd from ideal of singular locus --------------
   ideal D = M[1];
   for( ii=2; ii<=size(M); ii++ )            //look for the shortest one
   {
      if( size(M[ii]) < size(D[1]) )
      { 
          D = M[ii]; 
      }
   }

   //--------------- start p-th power algorithm and return ----------------
   ideal F = var(1)^p;
   for(ii=2; ii<=n; ii++)
   {  
      F=F,var(ii)^p; 
   }

   ideal Dp=D^(p-1);
   ideal U=1;
   ideal K,L;
   map phi=Q,F;
   if ( y >= 1)
   {  "compute module generators of integral closure";
      "denominator D is:";  D;  
      pause();
   }  

   ii=0;
   list LK;
   while(1)
   {
      ii=ii+1;
      if ( y >= 1)
      { "iteration", ii; }
      L = U*Dp + I;             
      //### L=interred(L) oder msdt(L)[2]? 
      //Wird dadurch kleiner aber string(L) wird groesser
      K = preimage(Q,phi,L);    //### Improvement by block ordering?
      option(returnSB);
      K = intersect(U,K);          //K is the new U, it is a SB
      LK = mstd(K);
      K = LK[2];

   //---------------------------- simplify output --------------------------
      if(size(reduce(U,LK[1]))==0)  //previous U coincides with new U
      {                             //i.e. we reached the integral closure 
         U=simplify(reduce(U,groebner(D)),2); 
         U = D,U;
         poly gg = gcd(U[1],U[size(U)]);
         for(ii=2; ii<=size(U)-1 ;ii++)
         {
            gg = gcd(gg,U[ii]);
         }
         for(ii=1; ii<=size(U); ii++)
         {
            U[ii]=U[ii]/gg;
         }
         U = simplify(U,6);
         //if ( y >= 1)
         //{ "module generators are U[i]/U[1], with U:"; U; 
         //  ""; pause(); }
         option(set,op);
         setring R;
         ideal U = imap(Q,U);
         return(U);
      }
      U=K;
   }
}

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

static proc substpartSpecial(ideal endid, ideal endphi)
{
   //Note: newRing is of the form (R the original basering): 
   //char(R),(T(1..N),X(1..nvars(R))),(dp(N),...);

   int ii,jj,kk;
   def BAS = basering;
   int n = nvars(basering);

   list Le = elimpart(endid); 
   int q = size(Le[2]);                   //q variables have been substituted
//Le;"";
   if ( q == 0 ) 
   {
      ideal normap = endphi;
      ideal norid = endid;
      export(norid);
      export(normap);
      list L = BAS;
      return(L);
   }

      list gnirlist = ringlist(basering);
      endid = Le[1];
//endphi;"";
      for( ii=1; ii<=n; ii++)
      {
         if( Le[4][ii] == 0 )            //ii=index of substituted var
         {
            endphi = subst(endphi,var(ii),Le[5][ii]);
         }
      }
//endphi;"";
      list g2 = gnirlist[2];             //the varlist
      list g3 = gnirlist[3];             //contains blocks of orderings
      int n3 = size(g3);  
   //----------------- first identify module ordering ------------------ 
      if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
      {
         list gm = g3[n3];              //last blockis module ordering
         g3 = delete(g3,n3);
         int m = 0;
      }
      else
      {
         list gm = g3[1];              //first block is module ordering
         g3 = delete(g3,1);
         int m = 1;
      }
   //---- delete variables which were substituted and weights  -------- 
//gnirlist;"";
      intvec V;
      int n(0);
      list newg2;
      list newg3;
      for ( ii=1; ii<=n3-1; ii++ )
      {
         V = V,g3[ii][2];           //copy weights for ordering in each block
         if ( ii==1 )               //into one intvector
         { 
            V = V[2..size(V)];
         }
        // int n(ii) = size(g3[ii][2]);
        int n(ii) = size(V);
        intvec V(ii);

        for ( jj = n(ii-1)+1; jj<=n(ii); jj++)
        {
//"ii, jj", ii,jj;
          if( Le[4][jj] !=0 )     //jj=index of var which was not substituted
          {
            kk=kk+1;
            newg2[kk] = g2[jj];   //not substituted var from varlist
            V(ii)=V(ii),V[jj];    //weight of not substituted variable
//"V(",ii,")", V(ii);
          }
        }
        if ( size(V(ii)) >= 2 )
        {
           V(ii) = V(ii)[2..size(V(ii))];
           list g3(ii)=g3[ii][1],V(ii);
           newg3 = insert(newg3,g3(ii),size(newg3));
//"newg3"; newg3;
        }
      }
//"newg3"; newg3;
      //newg3 = delete(newg3,1);    //delete empty list
      
/*
//### neue Ordnung, 1 Block fuer alle vars, aber Gewichte erhalten;
//vorerst nicht realisiert, da bei leonhard1 alte Version (neue Variable T(i)
//ein neuer Block) ein kuerzeres Ergebnis liefert
      kill g3;
      list g3; 
      V=0;
      for ( ii= 1; ii<=n3-1; ii++ )
      {
        V=V,V(ii);
      }
      V = V[2..size(V)];

      if ( V==1 )
      {    
         g3[1] = list("dp",V);
      }
      else
      {
         g3[1] = lis("wp",V);          
      }
      newg3 = g3;

//"newg3";newg3;"";
//### Ende neue Ordnung
*/

      if ( m == 0 )
      {
         newg3 = insert(newg3,gm,size(newg3));
      }
      else
      {
         newg3 = insert(newg3,gm);
      }
      gnirlist[2] = newg2;                    
      gnirlist[3] = newg3;

//gnirlist;
      def newBAS = ring(gnirlist);            //change of ring to less vars
      setring newBAS;
      ideal normap = imap(BAS,endphi);
      //normap = simplify(normap,2);
      ideal norid =  imap(BAS,endid);
      export(norid);
      export(normap);
      list L = newBAS;
      return(L);

   //Hier scheint interred gut zu sein, da es Ergebnis relativ schnell 
   //verkleinert. Hier wird z.B. bei leonard1 size(norid) verkleinert aber
   //size(string(norid)) stark vergroessert, aber es hat keine Auswirkungen
   //da keine map mehr folgt. 
   //### Bei Leonard2 haengt interred (BUG)
   //mstd[2] verkleinert norid nocheinmal um die Haelfte, dauert aber 3.71 sec
   //### Ev. Hinweis auf mstd in der Hilfe?

}

///////////////////////////////////////////////////////////////////////////////
// Computes the ring structure of a ring given by module generators.
// Assume: J[i]/J[1] are the module generators in the quotient field
// with J[1] as universal denominator.
// If option "noRed" is not given, a reduction in the number of variables is 
// attempted.
static proc computeRing(ideal J, ideal I, list #)
{
  int i, ii,jj;
  intvec V;                          // to be used for variable weights 
  int y = printlevel-voice+2;
  def R = basering;
  poly c = J[1];                     // the denominator
  list gnirlist = ringlist(basering);
  string svars = varstr(basering);
  int nva = nvars(basering);
  string svar;
  ideal maxid = maxideal(1);
  
  int noRed = 0;     // By default, we try to reduce the number of generators.
  if(size(#) > 0){
    if ( typeof(#[1]) == "string" ) 
    {
      if (#[1] == "noRed"){noRed = 1;}
    }
  }

  if ( y >= 1){"// computing the ring structure...";}
   
  if(c==1)
  {
    if( defined(norid) )  { kill norid; }
    if( defined(normap) ) { kill normap; }
    ideal norid = I;
    ideal normap =  maxid;
    export norid;
    export normap;
    if(noRed == 1){
      setring R;
      return(R);
    } else {
      list L = substpartSpecial(norid,normap);
      def lastRing = L[1];
      setring R;
      return(lastRing);
    }
  }


  //-------------- Enlarge ring by creating new variables ------------------
  //check first whether variables T(i) and then whether Z(i),...,A(i) exist
  //old variable names are not touched

  if ( find(svars,"T(") == 0 )        
  {
    svar = "T";
  }
  else
  {
    for (ii=90; ii>=65; ii--)  
    {
      if ( find(svars,ASCII(ii)+"(") == 0 )        
      {
        svar = ASCII(ii);  break;
      }
    }
  }

  int q = size(J)-1;
  if ( size(svar) != 0 )
  {
    for ( ii=q; ii>=1; ii-- )
    {
      gnirlist[2] = insert(gnirlist[2],svar+"("+string(ii)+")");
    }
  }
  else
  {       
    for ( ii=q; ii>=1; ii-- )
    {
      gnirlist[2] = insert(gnirlist[2],"T("+string(100*nva+ii)+")");
    }
  }

  V[q]=0;                        //create intvec of variable weights
  V=V+1;
  gnirlist[3] = insert(gnirlist[3],list("dp",V));

  //this is a block ordering with one dp-block (1st block) for new vars
  //the remaining weights and blocks for old vars are kept
  //### perhaps better to make only one block, keeping weights ?
  //this might effect syz below
  //alt: ring newR = char(R),(X(1..nvars(R)),T(1..q)),dp;
  //Reihenfolge geŠndert:neue Variablen kommen zuerst, Namen ev. nicht T(i)

  def newR = ring(gnirlist);   
  setring newR;                //new extended ring
  ideal I = imap(R,I);

  //------------- Compute linear and quadratic relations ---------------
  if(y>=1)
  {
     "// compute linear relations:"; 
  }
  qring newQ = std(I);
  
  ideal f = imap(R,J);
  module syzf = syz(f);
  ideal pf = f[1]*f;               
  //f[1] is the denominator D from normalityTest, a non zero divisor of R/I

  ideal newT = maxideal(1);
  newT = 1,newT[1..q];
  //matrix T = matrix(ideal(1,T(1..q)),1,q+1);   //alt
  matrix T = matrix(newT,1,q+1); 
  ideal Lin = ideal(T*syzf);
  //Lin=interred(Lin); 
  //### interred reduziert ev size aber size(string(LIN)) wird groesser

  if(y>=1)
  {
    if(y>=3)
    {
      "//   the linear relations:";  Lin; pause();"";
    }
      "// the ring structure of the normalization as affine algebra";
      "//   number of linear relations:", size(Lin);
  }

  if(y>=1)
  {
    "// compute quadratic relations:"; 
  }
  matrix A;
  ideal Quad;
  poly ff;
  newT = newT[2..size(newT)];
  matrix u;  // The units for non-global orderings.

  // Quadratic relations
  for (ii=2; ii<=q+1; ii++ )
  {
    for ( jj=2; jj<=ii; jj++ )
    {
      ff = NF(f[ii]*f[jj],std(0));     // this makes lift much faster
      // For non-global orderings, we have to take care of the units.
      if(ord_test(basering) != 1){
        A = lift(pf, ff, u);
        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] * u[1, 1]- T*A);  
      } else {
        A = lift(pf,ff);              // ff lin. comb. of elts of pf mod I
        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] - T*A);
      }
      //A = lift(pf, f[ii]*f[jj]);
      //Quad = Quad, ideal(T(jj-1)*T(ii-1) - T*A);  
    }
  }
  Quad = Quad[2..ncols(Quad)];

  if(y>=1)
  {
    if(y>=3)
    {
      "//   the quadratic relations"; Quad; pause();"";
    }
      "//   number of quadratic relations:", size(Quad);
  }
  ideal Q1 = Lin,Quad;     //elements of Q1 are in NF w.r.t. I

  //Q1 = mstd(Q1)[2];
  //### weglassen, ist sehr zeitaufwendig.
  //Ebenso interred, z.B. bei Leonard1 (1. Komponente von Leonard):
  //"size Q1:", size(Q1), size(string(Q1));   //75 60083
  //Q1 = interred(Q1); 
  //"size Q1:", size(Q1), size(string(Q1));   //73 231956 (!)
  //### Speicherueberlauf bei substpartSpecial bei 'ideal norid  = phi1(endid)'
  //Beispiel fuer Hans um map zu testen!

  setring newR;
  ideal endid  = imap(newQ,Q1),I;
  ideal endphi = imap(R,maxid);

  if(noRed == 0){   
    list L=substpartSpecial(endid,endphi);
    def lastRing=L[1];
    if(y>=1)
    {
      "//   number of substituted variables:", nvars(newR)-nvars(lastRing);
      pause();"";
    }
    return(lastRing);
  } else {
    ideal norid = endid;
    ideal normap = endphi;
    export(norid);
    export(normap);
    setring R;
    return(newR);
  }   
}

//                Up to here: procedures for char p with Frobenius
///////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
//                From here: subprocedures for normal

static proc normalM(ideal I, int decomp, int withDelta){
// Computes the normalization of a ring R / I using the module structure as far
// as possible.
// The ring R is the basering.
// Input: ideal I
// Output: a list of 3 elements (resp 4 if withDelta = 1), say nor.
// - nor[1] = U, an ideal of R.
// - nor[2] = c, an element of R.
// U and c satisfy that 1/c * U is the normalization of R / I in the 
// quotient field Q(R/I).
// - nor[3] = ring say T, containing two ideals norid and normap such that
// normap gives the normalization map from R / I to T / norid.
// - nor[4] = the delta invariant, if withDelta = 1.

// Details:
// --------
// Computes the ideal of the minors in the first step and then reuses it in all 
// steps.
// In step s, the denominator is D^s, where D is a nzd of the original quotient 
// ring, contained in the radical of the singular locus.
// This denominator is used except when the degree of D^i is greater than the 
// degree of a conductor.
// The nzd is taken as the smallest degree polynomial in the radical of the 
// singular locus.

// It assumes that the ideal I is equidimensional radical. This is not checked 
// in the procedure!
// If decomp = 0 or decomp = 3 it assumes further that I is prime. Therefore 
// any non-zero element in the jacobian ideal is assumed to be a 
// non-zerodivisor.

// It works over the given basering.
// If it has a non-global ordering, it changes it to dp global only for 
// computing radical.

// The delta invariant is only computed if withDelta = 1, and decomp = 0 or 
// decomp = 3 (meaning that the ideal is prime).

  option("redSB");
  option("returnSB");
  int step = 0;                       // Number of steps. (for debugging)
  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
  int i;                              // counter
  int isGlobal = ord_test(basering);

  poly c;                     // The common denominator.

  def R = basering;
    
  // ---------------- computation of the singular locus --------------------- 
  // We compute the radical of the ideal of minors modulo the original ideal.
  // This is done only in the first step.

  if(isGlobal == 1){
    list IM = mstd(I);
    I = IM[1];
    ideal IMin = IM[2];   // A minimal set of generators in the groebner basis.
  } else {
    // The minimal set of generators is not computed by mstd for
    // non-global orderings.
    I = groebner(I);
    ideal IMin = I;
  }
  int d = dim(I);
  
  qring Q = I;   // We work in the quotient by the groebner base of the ideal I
  option("redSB");
  option("returnSB");
  ideal I = fetch(R, I);
  attrib(I, "isSB", 1);
  ideal IMin = fetch(R, IMin);

  dbprint(dbg, "Computing the jacobian ideal...");
  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
  J = groebner(J);
  
  // -------------------- election of the conductor -------------------------   
  poly condu = getSmallest(J);   // Choses the polynomial of smallest degree 
                                 // of J as universal denominator.
  if(dbg >= 1){
    "";
    "The universal denominator is ", condu;
  }

  // ---------  splitting the ideal by the conductor (if possible) ----------- 
  // If the ideal is equidimensional, but not necessarily prime, we check if 
  // the conductor is a non-zerodivisor of R/I.
  // If not, we split I.
  if((decomp == 1) or (decomp == 2)){  
    ideal Id1 = quotient(0, condu);
    if(size(Id1) > 0){
      // We have to split.
      if(dbg >= 1){
        "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu;
      }
      setring R;
      ideal Id1 = fetch(Q, Id1), I;
      Id1 = groebner(Id1);
      ideal Id2 = quotient(I, Id1);
      // I = I1 \cap I2
      printlevel = printlevel + 1;
      list nor1 = normalM(Id1, decomp, withDelta)[1];
      list nor2 = normalM(Id2, decomp, withDelta)[1];
      printlevel = printlevel - 1;    
      return(list(nor1, nor2));
    }
  }

  // --------------- computation of the first test ideal ---------------------   
  // To compute the radical we go back to the original ring.
  // If we are using a non-global ordering, we must change to the global 
  // ordering.
  if(isGlobal == 1){
    setring R;
    ideal J = fetch(Q, J);
    J = J, I;
    if(dbg >= 1){
      "The original singular locus is";
      groebner(J);
      if(dbg >= 2){pause();}
       "";
    }
    J = radical(J);
  } else {
    // We change to global dp ordering.
    list rl = ringlist(R);
    list origOrd = rl[3];
    list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0);
    rl[3] = newOrd;
    def globR = ring(rl);
    setring globR;
    ideal J = fetch(Q, J);
    ideal I = fetch(R, I);
    J = J, I;
    if(dbg >= 1){
      "The original singular locus is";
      groebner(J);
      if(dbg>=2){pause();}
      "";
    }
    J = radical(J);
    setring R;
    ideal J = fetch(globR, J);
  }
    
  if(dbg >= 1){
    "The radical of the original singular locus is";
    J;
    if(dbg>=2){pause();}
  }

  // ---------------- election of the non zero divisor --------------------- 
  setring Q;
  J = fetch(R, J);
  J = interred(J);          
  poly D = getSmallest(J);    // Chooses the poly of smallest degree as 
                              // non-zerodivisor.
  if(dbg >= 1){
    "The non zero divisor is ", D;
    "";
  }

  // ------- splitting the ideal by the non-zerodivisor (if possible) --------   
  // If the ideal is equidimensional, but not necessarily prime, we check if D
  // is actually a non-zerodivisor of R/I.
  // If not, we split I.
  if((decomp == 1) or (decomp == 2)){  
    // We check if D is actually a non-zerodivisor of R/I.
    // If not, we split I.
    Id1 = quotient(0, D);
    if(size(Id1) > 0){
      // We have to split.
      if(dbg >= 1){
        "A zerodivisor was found. We split the ideal. The zerodivisor is ", D;
      }    
      setring R;
      ideal Id1 = fetch(Q, Id1), I;
      Id1 = groebner(Id1);
      ideal I2 = quotient(I, Id1);
      // I = Id1 \cap Id2
      printlevel = printlevel + 1;
      list nor1 = normalM(Id1, decomp, withDelta)[1];
      list nor2 = normalM(Id2, decomp, withDelta)[1];
      printlevel = printlevel - 1;    
      return(list(nor1, nor2));    
    }
  }

  // --------------------- normalization ------------------------------------ 
  // We call normalMEqui to compute the normalization.
  setring R;
  poly D = fetch(Q, D);
  poly condu = fetch(Q, condu);
  J = fetch(Q, J);
  printlevel = printlevel + 1;
  list result = normalMEqui(I, J, condu, D, withDelta);
  printlevel = printlevel - 1;
  return(list(result));
}

///////////////////////////////////////////////////////////////////////////////
  
static proc normalMEqui(ideal I, ideal origJ, poly condu, poly D, int withDelta)
// Here is where the normalization is actually computed.

// Computes the normalization of R/I. (basering is R)
// I is assumed to be radical and equidimensional.
// origJ is the first test ideal.
// D is a non-zerodivisor of R/I.
// condu is a non-zerodivisor in the conductor.
// If withDelta = 1, computes the delta invariant.
{
  int step = 0;                       // Number of steps. (for debugging)
  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
  int i;                              // counter
  int isNormal = 0;                   // check for exiting the loop
  int isGlobal = ord_test(basering);
  int delt;
  
  def R = basering;
  poly c = D;
  ideal U;
  ideal cJ;
  list testOut;                 // Output of proc testIdeal
                                // (the test ideal and the ring structure)
    
  qring Q = groebner(I);
  option("redSB");
  option("returnSB");  
  ideal J = imap(R, origJ);
  poly c = imap(R, c);
  poly D = imap(R, D);
  poly condu = imap(R, condu);
  ideal cJ;
  ideal cJMod;

  dbprint(dbg, "Preliminar step begins.");

  // --------------------- computation of A1 ------------------------------- 
  dbprint(dbg, "Computing the quotient (DJ : J)...");
  ideal U = groebner(quotient(D*J, J));  
  ideal oldU = 1;

  if(dbg >= 2){
    "The quotient is"; U;
  }

  // ----------------- Grauer-Remmert criterion check ----------------------- 
  // We check if the equality in Grauert - Remmert criterion is satisfied.
  isNormal = checkInclusions(D*oldU, U);
  if(isNormal == 0){
    if(dbg >= 1){
      "In this step, we have the ring 1/c * U, with c =", c;
      "and U = "; U;
    }
  } else {        
    // The original ring R/I was normal. Nothing to do.
    // We define anyway a new ring, equal to R, to be able to return it.
    setring R;
    list lR = ringlist(R);
    def ROut = ring(lR);
    setring ROut;
    ideal norid = 0;
    ideal normap = maxideal(1);
    export norid;
    export normap;
    setring R;
    if(withDelta){
      list output = ideal(1), poly(1), ROut, 0;
    } else {
      list output = ideal(1), poly(1), ROut;
    }
    return(output);
  }

  // ----- computation of the chain of ideals A1 c A2 c ... c An ------------ 
  while(isNormal == 0){
    step++;
    if(dbg >= 1){
      "";
      "Step ", step, " begins.";
    }
    dbprint(dbg, "Computing the test ideal...");
    
    // --------------- computation of the test ideal ------------------------ 
    // Computes a test ideal for the new ring. 
    // The test ideal will be the radical in the new ring of the original
    // test ideal.
    setring R;
    U = imap(Q, U);
    c = imap(Q, c);
    testOut = testIdeal(I, U, origJ, c, D);
    cJ = testOut[1];
    
    setring Q;
    cJ = imap(R, cJ);
    cJ = groebner(cJ);
    
    // cJ / c is now the ideal mapped back.
    // We have the generators as an ideal in the new ring, 
    // but we want the generators as an ideal in the original ring.
    cJMod = getGenerators(cJ, U, c);

    if(dbg >= 2){
      "The test ideal in this step is ";
      cJMod;
    }
    
    cJ = cJMod;

    // ------------- computation of the quotient (DJ : J)-------------------- 
    oldU = U;
    dbprint(dbg, "Computing the quotient (c*D*cJ : cJ)...");
    U = quotient(c*D*cJ, cJ);
    if(dbg >= 2){"The quotient is "; U;}

    // ------------- Grauert - Remmert criterion check ---------------------- 
    // We check if the equality in Grauert - Remmert criterion is satisfied.
    isNormal = checkInclusions(D*oldU, U);

    if(isNormal == 1){        
      // We go one step back. In the last step we didnt get antyhing new, 
      // we just verified that the ring was already normal.
      dbprint(dbg, "The ring in the previous step was already normal.");
      dbprint(dbg, "");
      U = oldU;
    } else {
      // ------------- preparation for next iteration ---------------------- 
      // We have to go on.
      // The new denominator is chosen.
      c = D * c;
      if(deg(c) > deg(condu)){
        U = changeDenominatorQ(U, c, condu);
        c = condu;
      }
      if(dbg >= 1){
        "In this step, we have the ring 1/c * U, with c =", c;
        "and U = ";
        U;
        if(dbg>=2){pause();}
      }
    }
  }

  // ------------------------- delta computation ---------------------------- 
  if(withDelta){ 
    ideal UD = groebner(U);
    delt = vdim(std(modulo(UD, c)));
  }
  
  // -------------------------- prepare output -----------------------------   
  setring R;
  U = fetch(Q, U);
  c = fetch(Q, c);

  // Ring structure of the new ring
  def ere = testOut[2];
  if(withDelta){
    list output = U, c, ere, delt;
  } else {
    list output = U, c, ere;
  }
  return(output);
 
}

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

static proc lineUp(ideal U, poly c)
// Sets c as the first generator of U.
{
  int i;
  ideal newU = c;
  for (i = 1; i <= ncols(U); i++){
    if(U[i] != c){
      newU = newU, U[i];
    }
  }
  return(newU);
}

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

static proc getSmallest(ideal J){
// Computes the polynomial of smallest degree of J.
// If there are more than one, it chooses the one with smallest number 
// of monomials.
  int i;
  poly p = J[1];
  int d = deg(p);
  int di;
  for(i = 2; i <= ncols(J); i++){
    if(J[i] != 0){
      di = deg(J[i]);
      if(di < d){
        p = J[i];
        d = di;
      } else {
        if(di == d){
          if(size(J[i]) < size(p)){
            p = J[i];
          }
        }
      }
    }
  }
  return(p);
}

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

static proc getGenerators(ideal J, ideal U, poly c){
// Computes the generators of J as an ideal in the original ring,
// where J is given by generators in the new ring.

// The new ring is given by 1/c * U in the total ring of fractions.

  int i, j;                             // counters;
  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)  
  poly p;                               // The lifted polynomial
  ideal JGr = groebner(J);              // Groebner base of J

  if(dbg>1){"Checking for new generators...";}
  for(i = 1; i <= ncols(J); i++){
    for(j = 1; j <= ncols(U); j++){
      p = lift(c, J[i]*U[j])[1,1];
      p = reduce(p, JGr);
      if(p != 0){
        if(dbg>1){
          "New polynoial added:", p;
          if(dbg>4){pause();}
        }
        JGr = JGr, p;
        JGr = groebner(JGr);
        J = J, p;
      }
    }
  }
  return(J);
}

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

static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D){
// Internal procedure, used in normalM.
// Computes the test ideal in the new ring.
// It takes the original test ideal and computes the radical of it in the 
// new ring.

// The new ring is 1/c * U.
// The original test ideal is origJ.
// The original ring is R / I, where R is the basering.
  int i;                                // counter
  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)  
  def R = basering;                      // We dont work in the quo
  ideal J = origJ;
  
  // ---------- computation of the ring structure of 1/c * U ----------------
  U = lineUp(U, c);
  
  if(dbg > 1){"Computing the new ring structure...";}
  list ele = computeRing(U, I, "noRed");

  def origEre = ele[1];
  setring origEre;
  if(dbg > 1){"The relations are"; norid;}

  // ---------------- setting the ring to work in  --------------------------  
  int isGlobal = ord_test(origEre);      // Checks if the original ring has 
                                         // global ordering.
  if(isGlobal != 1){
    list rl = ringlist(origEre);
    list origOrd = rl[3];
    list newOrd = list("dp", intvec(1:nvars(origEre))), list("C", 0);
    rl[3] = newOrd;
    def ere = ring(rl);     // globR is the original ring but 
                            // with a global ordering.
    setring ere;
    ideal norid = imap(origEre, norid);
  } else {
    def ere = origEre;
  }

  ideal I = imap(R, I);
  ideal J = imap(R, J);
  J = J, norid, I;


  // ----- computation of the test ideal using the ring structure of Ai -----  
  option("redSB");
  option("returnSB");

  if(dbg > 1){"Computing the radical of J...";}
  J = radical(J);
  if(dbg > 1){"Computing the interreduction of the radical...";}
  J = groebner(J);  
  //J = interred(J);
  if(dbg > 1){
    "The radical in the generated ring is";
    J;
    if(dbg>4){pause();}
  }
  
  setring ere;
  
  // -------------- map from Ai to the total ring of fractions ---------------  
  // Now we must map back this ideal J to U_i / c in the total ring of
  // fractions.
  // The map sends T_j -> u_j / c.
  // The map is built by the following steps:
  // 1) We compute the degree of the generators of J with respect to the 
  //    new variables T_j.
  // 2) For each generator, we multiply each term by a power of c, as if
  //    taking c^n as a common denominator (considering the new variables as
  //    a polynomial in the old variables divided by c).
  // 3) We replace the new variables T_j by the corresponding numerator u_j.
  // 4) We lift the resulting polynomial to change the denominator 
  //    from c^n to c.
  int nNewVars = nvars(ere) - nvars(R);      // Number of new variables
  poly c = imap(R, c);
  intvec @v = 1..nNewVars;    // Vector of the new variables. 
                              // They must be the first ones.
  if(dbg > 1){"The indices of the new variables are", @v;}

  // ---------------------- step 1 of the mapping ---------------------------  
  intvec degs;
  for(i = 1; i<=ncols(J); i++){
    degs[i] = degSubring(J[i], @v);
  }
  if(dbg > 1){
    "The degrees with respect to the new variables are";
    degs;
  }
  
  // ---------------------- step 2 of the mapping ---------------------------    
  ideal mapJ = mapBackIdeal(J, c, @v);
  
  setring R;

  // ---------------------- step 3 of the mapping ---------------------------  
  ideal z;                    // The variables of the original ring in order.
  for(i = 1; i<=nvars(R); i++){
    z[i] = var(i);
  }
      
  map f = ere, U[2..ncols(U)], z[1..ncols(z)]; // The map to the original ring.
  if(dbg > 1){
    "The map is ";
    f;
    if(dbg>4){pause();}
  }
  
  if(dbg > 1){
    "Computing the map...";
  }
      
  J = f(mapJ);
  if(dbg > 1){
    "The ideal J mapped back (before lifting) is";
    J;
    if(dbg>4){pause();}
  }

  // ---------------------- step 4 of the mapping ---------------------------    
  qring Q = groebner(I);
  ideal J = imap(R, J);
  poly c = imap(R, c);
  for(i = 1; i<=ncols(J); i++){
    if(degs[i]>1){
      J[i] = lift(c^(degs[i]-1), J[i])[1,1];
    } else {
      if(degs[i]==0){
        J[i] = c*J[i];
      }
    }
  }

  if(dbg > 1){
    "The ideal J lifted is";
    J;
    if(dbg>4){pause();}
  }

  // --------------------------- prepare output ----------------------------    
  J = groebner(J);
  
  setring R;
  J = imap(Q, J);
  
  return(list(J, ele[1]));
}

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

static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I){
// Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the 
// ring is 1/c2 * U2.
// The base ring is R, but the computations are to be done in R / I.
  int a;      // counter
  def R = basering;
  qring Q = I;
  ideal U1 = fetch(R, U1);
  poly c1 = fetch(R, c1);
  poly c2 = fetch(R, c2);
  ideal U2 = changeDenominatorQ(U1, c1, c2);
  setring R;
  ideal U2 = fetch(Q, U2);
  return(U2);
}

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

static proc changeDenominatorQ(ideal U1, poly c1, poly c2){
// Given a ring in the form 1/c1 * U, it computes a new U2 st the ring 
// is 1/c2 * U2.
// The base ring is already a quotient ring R / I.
  int a;      // counter
  ideal U2;
  poly p;
  for(a = 1; a <= ncols(U1); a++){
    p = lift(c1, c2*U1[a])[1,1];
    U2[a] = p;
  }
  return(U2);
}

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

static proc checkInclusions(ideal U1, ideal U2){
// Checks if the identity A = Hom(J, J) of Grauert-Remmert criterion is 
// satisfied.
  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
  list reduction1;
  list reduction2;
  
  // ---------------------- inclusion Hom(J, J) c A -------------------------  
  if(dbg > 1){"Checking the inclusion Hom(J, J) c A:";}
  // This interred is used only because a bug in groebner!
  U1 = groebner(U1);
  reduction1 = reduce(U2, U1);
  if(dbg > 1){reduction1[1];}

  // ---------------------- inclusion A c Hom(J, J) -------------------------  
  // The following check should always be satisfied.
  // This is only used for debugging.
  if(dbg > 1){
    "and the inclusion A c Hom(J, J): (this should always be satisfied)";
    // This interred is used only because a bug in groebner!
    U2 = groebner(U2);
    reduction2 = reduce(U1, groebner(U2));
    reduction2[1];
    if(size(reduction2[1]) > 0){
      "Something went wrong... (this inclusion should always be satisfied)";
      ~;
    } else {
      if(dbg>4){pause();}
    }
  }
  
  if(size(reduction1[1]) == 0){
    // We are done! The ring computed in the last step was normal.
    return(1);
  } else {
    return(0);
  }
}

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

static proc degSubring(poly p, intvec @v){
// Computes the degree of a polynomial taking only some variables as variables
// and the others as parameters.

// The degree is taken wrt the variables indicated in v.
  int i;      // Counter
  int d = 0;  // The degree
  int e;      // Degree (auxiliar variable)
  for(i = 1; i <= size(p); i++){
    e = sum(leadexp(p[i]), @v);
    if(e > d){d = e;}
  }
  return(d);
}
    
///////////////////////////////////////////////////////////////////////////////

static proc mapBackIdeal(ideal I, poly c, intvec @v){
// Modifies all polynomials in I so that a map x(i) -> y(i)/c can be 
// carried out.

// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.

  int i;  // counter
  for(i = 1; i <= ncols(I); i++){
    I[i] = mapBackPoly(I[i], c, @v);
  }
  return(I);
}

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

static proc mapBackPoly(poly p, poly c, intvec @v){
// Multiplies each monomial of p by a power of c so that a map x(i) -> y(i)/c 
// can be carried out.

// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.
  int i;  // counter
  int e;  // exponent
  int d = degSubring(p, @v);
  poly g = 0;
  for(i = 1; i <= size(p); i++){
    e = sum(leadexp(p[i]), @v);
    g = g + p[i] * c^(d-e);
  }
  return(g);
}

//                Up to here: procedures for normal
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
//                From here: procedures for normalC

// We first define resp. copy some attributes to be used in proc normal and
// static proc normalizationPrimes, and ..., to speed up computation in
// special cases
//NOTE:  We use the following attributes:
// 1     attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
// 2     attrib(id,"isCompleteIntersection");  //--- complete intersection  
// 3     attrib(id,"isHypersurface");          //--- hypersurface 
// 4     attrib(id,"isEquidimensional");       //--- equidimensional ideal 
// 5     attrib(id,"isPrim");                  //--- prime ideal 
// 6     attrib(id,"isRegInCodim2");           //--- regular in codimension 2 
// 7     attrib(id,"isIsolatedSingularity");   //--- isolated singularities 
// 8     attrib(id,"onlySingularAtZero");      //--- only singular at 0 
// 9     attrib(id,"isRadical");               //--- radical ideal 
//Recall: (attrib(id,"xy"),1) sets attrib xy to TRUE and
//        (attrib(id,"xy"),0) to FALSE

static proc getAttrib (ideal id)
"USAGE:   getAttrib(id);  id=ideal
COMPUTE: check attributes for id. If the attributes above are defined, 
         take its value, otherwise define it and set it to 0 
RETURN:  intvec of size 9, with entries 0 or 1,  values of attributes defined
         above (in this order)
EXAMPLE: no example
"
{
  int isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;

  if( typeof(attrib(id,"isCohenMacaulay"))=="int" )   
  {
    if( attrib(id,"isCohenMacaulay")==1 ) 
    { isCoM=1; isEq=1; }
  }

  if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
  {
    if(attrib(id,"isCompleteIntersection")==1) 
    { isCoI=1; isCoM=1; isEq=1; }
  }
    
  if( typeof(attrib(id,"isHypersurface"))=="int" )
  {
    if(attrib(id,"isHypersurface")==1)
    { isHy=1; isCoI=1; isCoM=1; isEq=1; }
  }
   
  if( typeof(attrib(id,"isEquidimensional"))=="int" )
  {
    if(attrib(id,"isEquidimensional")==1) 
    { isEq=1; }
  }
   
  if( typeof(attrib(id,"isPrim"))=="int" )
  {
    if(attrib(id,"isPrim")==1)  
    { isPr=1; }
  }

  if( typeof(attrib(id,"isRegInCodim2"))=="int" )
  {
    if(attrib(id,"isRegInCodim2")==1)  
    { isReg=1; }
  }
   
  if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
  {
    if(attrib(id,"isIsolatedSingularity")==1)  
    { isIso=1; }
  }
 
  if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
  {
    if(attrib(id,"onlySingularAtZero")==1)
    { oSAZ=1; }
  }

  if( typeof(attrib(id,"isRad"))=="int" )
  {
    if(attrib(id,"isRad")==1)  
    { isRad=1; }
  }

  intvec atr = isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;
  return(atr);
}

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

static proc setAttrib (ideal id, intvec atr)
"USAGE:   setAttrib(id,atr);  id ideal, atr intvec
COMPUTE: set attributes to id specified by atr
RETURN:  id, with assigned attributes from atr
EXAMPLE: no example
"
{
  attrib(id,"isCohenMacaulay",atr[1]);         //--- Cohen Macaulay
  attrib(id,"isCompleteIntersection",atr[2]);  //--- complete intersection  
  attrib(id,"isHypersurface",atr[3]);          //--- hypersurface 
  attrib(id,"isEquidimensional",atr[4]);       //--- equidimensional ideal 
  attrib(id,"isPrim",atr[5]);                  //--- prime ideal 
  attrib(id,"isRegInCodim2",atr[6]);           //--- regular in codimension 2 
  attrib(id,"isIsolatedSingularity",atr[7]);   //--- isolated singularities 
  attrib(id,"onlySingularAtZero",atr[8]);      //--- only singular at 0 
  attrib(id,"isRadical",atr[9]);               //--- radical ideal 

  return(id);
}

///////////////////////////////////////////////////////////////////////////////
// copyAttribs is not used anywhere so far

static proc copyAttribs (ideal id1, ideal id)
"USAGE:   copyAttribs(id1,id);  id1, id ideals
COMPUTE: copy attributes from id1 to id
RETURN:  id, with assigned attributes from id1
EXAMPLE: no example
"
{
  if( typeof(attrib(id1,"isCohenMacaulay"))=="int" )   
  {
    if( attrib(id1,"isCohenMacaulay")==1 )
    {
      attrib(id,"isEquidimensional",1);
    }
  }
  else
  {
    attrib(id,"isCohenMacaulay",0);
  }  

  if( typeof(attrib(id1,"isCompleteIntersection"))=="int" )
  {
    if(attrib(id1,"isCompleteIntersection")==1)
    {
      attrib(id,"isCohenMacaulay",1);
      attrib(id,"isEquidimensional",1);
    }
  }
  else
  {
    attrib(id,"isCompleteIntersection",0);
  }
    
  if( typeof(attrib(id1,"isHypersurface"))=="int" )
  {
    if(attrib(id1,"isHypersurface")==1)
    {
      attrib(id,"isCompleteIntersection",1);
      attrib(id,"isCohenMacaulay",1);
      attrib(id,"isEquidimensional",1);
    }
  }
  else
  {
    attrib(id,"isHypersurface",0);
  }
   
  if( (typeof(attrib(id1,"isEquidimensional"))=="int") )
  {
    if(attrib(id1,"isEquidimensional")==1)
    {
      attrib(id,"isEquidimensional",1);
    }
  }
  else
  {
    attrib(id,"isEquidimensional",0);
  }

  if( typeof(attrib(id1,"isPrim"))=="int" )
  {
    if(attrib(id1,"isPrim")==1)
    {
      attrib(id,"isEquidimensional",1);
    }
  }
  else
  {
    attrib(id,"isPrim",0);
  }
   
  if( (typeof(attrib(id1,"isRegInCodim2"))=="int") )
  {
    if(attrib(id1,"isRegInCodim2")==1)
    {
      attrib(id,"isRegInCodim2",1);
    }
  }
  else
  {
    attrib(id,"isRegInCodim2",0);
  }
    
  if( (typeof(attrib(id1,"isIsolatedSingularity"))=="int") )
  {
    if(attrib(id1,"isIsolatedSingularity")==1)
    {
      attrib(id,"isIsolatedSingularity",1);
    }
  }
  else
  {
    attrib(id,"isIsolatedSingularity",0);
  }
 
  if( typeof(attrib(id1,"onlySingularAtZero"))=="int" )
  {
    if(attrib(id1,"onlySingularAtZero")==1)
    {
      attrib(id,"isIsolatedSingularity",1);
    }
  }
  else
  {
    attrib(id,"onlySingularAtZero",0);
  }

  if( typeof(attrib(id1,"isRad"))=="int" )
  {
    if(attrib(id1,"isRad")==1)  
    {
      attrib(id,"isRad",1);
    }
  }
  else
  {
    attrib(id,"isRad",0);
  }

  return(id);
}

proc normalC(ideal id, list #)
"USAGE:   normalC(id [,choose]);  id = radical ideal, choose = optional string 
         which may be a combination of \"equidim\", \"prim\", \"withGens\",
         \"isPrim\", \"noFac\".@*
         If choose is not given or empty, a default method is choosen@*
         If choose = \"prim\" resp. \"equidim\" normalC computes a prime resp.
         an equidimensional decomposition and then the normalization of each
         component; if choose = \"noFac\" factorization is avoided in the 
         computation of the minimal associated primes; @*
         if choose = \"withGens\" the minimal associated primes P_i of id are 
         computed and for each P_i, algebra generators of the integral closure 
         of basering/P_i are computed as elements of its quotient field;@*
         \"isPrim\" disables \"equidim\" and \"prim\".@*
ASSUME:  The ideal must be radical, for non-radical ideals the output may
         be wrong (id=radical(id); makes id radical). However, if choose = 
         \"prim\" the minimal associated primes of id are computed first 
         and hence normalC computes the normalization of the radical of id.
         \"isPrim\" should only be used if id is known to be irreducible.
RETURN:  a list, say nor, of size 2 (resp. 3 if choose=\"withGens\"). 
         The first list nor[1] consists of r rings, where r is the number 
         of irreducible components if choose = \"prim\" (resp. >= no of 
         equidimenensional components if choose =  \"equidim\"). 
@format 
         Each ring Ri=nor[1][i], i=1..r, contains two ideals with given
         names @code{norid} and @code{normap} such that @*
         - Ri/norid is the normalization of the i-th component, i.e. the 
          integral closure in its field of fractions (as affine ring);
         - the direct sum of the rings Ri/norid is the normalization
           of basering/id; @*
         - @code{normap} gives the normalization map from basering/id to
           Ri/norid for each j.@*

         If choose=\"withGens\", nor[2] is a list of size r of ideals Ii= 
         nor[2][i] in the basering, generating the integral closure of
         basering/P_i in its quotient field, i=1..r, in two different ways:
         - Ii is given by polynomials 
           g_1,...,g_k,g_k+1 such that the variables T(j) of the ring Ri 
           satisfy T(1)=g_1/g_k+1,..,T(k)=g_k/g_k+1. Hence the g_j/g_k+1
           are algebra generators of the integral closure of basering/P_i
           as sub-algebra in the quotient field of basering/P_i. 
 
         nor[2] (resp. nor[3] if choose=\"withGens\") is a list of an intvec 
         of size r, the delta invariants of the r components, and an integer, 
         the total delta invariant of basering/id (-1 means infinite, and 0 
         that basering/P_i resp. basering/input is normal). 
         Return value -2 means that delta resp. delta of one of the components
         is not computed (which may happen if \"equidim\" is given) .
@end format
THEORY:  The delta invariant of a reduced ring K[x1,...,xn]/id is 
             dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id)
         We call this number also the delta invariant of id.
         We use the algorithm described in [G.-M. Greuel, G. Pfister: 
         A SINGULAR Introduction to Commutative Algebra, 2nd Edition. 
         Springer Verlag (2007)].
NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
@*       Increasing/decreasing printlevel displays more/less comments 
         (default: printlevel=0).
@*       Not implemented for local or mixed orderings and for quotient rings.
@*       If the input ideal id is weighted homogeneous a weighted ordering may
         be used (qhweight(id); computes weights).
KEYWORDS: normalization; integral closure; delta invariant.
EXAMPLE: example normalC; shows an example
"
{  
   int i,j, wgens, withEqui, withPrim, isPrim, nfac;
   int y = printlevel-voice+2;
   int nvar = nvars(basering);
   int chara  = char(basering);
   list result,prim,keepresult;

   if ( ord_test(basering) != 1 )
   {
     "";
     "// Not implemented for this ordering,";
     "// please change to global ordering!";
     return(result);
   }

//--------------------------- define the method --------------------------- 
   string method;                //make all options one string in order to use
                                 //all combinations of options simultaneously
   for ( i=1; i <= size(#); i++ )
   { 
     if ( typeof(#[i]) == "string" ) 
     {
       method = method + #[i];
     }
   }

   //--------------------------- choosen methods -----------------------
   //we have the methods 
   // "withGens": computes algebra generators for each irreducible component
   // the general case: either "equidim" or not (then applying minAssGTZ)

   if ( find(method,"withgens") or find(method,"withGens"))
   {
     wgens=1;
   }
   if ( find(method,"equidim") )
   {
     withEqui=1;
   }
   if ( find(method,"prim") )
   {
     withPrim=1;
   }
   if ( find(method,"isprim") or find(method,"isPrim") )
   {
     isPrim=1;
   }
   if ( find(method,"noFac") )
   {
     nfac=1;            
   }

   //------------------ default method, may be changed --------------------
   //Use a prime decomposition only for small no of variables otherwise equidim 
   //Compute delta whenever possible.

   if ( ((chara==0 || chara>30) && nvar>2 && !withPrim && !isPrim) ||
   ((chara==0 || chara>30) && size(#)  ==  0) )
   {
      withEqui=1;
   }
   kill #;
   list #;

//------- Special algorithm with computation of the generators, RETURN -------
   //--------------------- method "withGens" ----------------------------------
   //the integral closure is computed in proc primeClosure. In the general case
   //it is computed in normalizationPrimes. The main difference is that  in
   //primeClosure the singular locus is only computed in the first iteration
   //that no attributes are used and that the generators are computed.
   //In primeClosure the (algebra) generators for each irreducible component
   //are computed in static proc closureGenerators

   if( wgens )
   {
      if( y >= 1 )
      {  "";
         "// We use method 'withGens'";
      }
      if ( isPrim )
      {
         prim[1] = id;
         if( y >= 0 )
         { 
           "";
           "// ** WARNING: result is correct if ideal is prime (not checked) **";
           "// if procedure is called with string \"prim\", primality is checked";
         }
      }
      else
      {
         if(y>=1)
         {  "// compute minimal associated primes"; }

         if( nfac )
         { prim = minAssGTZ(id,1); }
         else
         { prim = minAssGTZ(id); }

         if(y>=1)
         {  
            prim;"";
            "// number of irreducible components is", size(prim);
         }
      }
   //----------- compute integral closure for every component -------------
      int del;
      intvec deli;
      list Gens,l,resu,Resu;
      ideal gens;
      def R = basering;
      poly gg;

      for(i=1; i<=size(prim); i++)
      {
         if(y>=1)
         {
            ""; pause(); "";
            "// start computation of component",i;
            "   --------------------------------";
         }

         if( defined(ker) ) { kill ker; }
         ideal ker = prim[i];
         export(ker);
         l = R;
         l = primeClosure(l,1);               //here the work is done
         //### ausprobieren ob primeClosure(l,1) schneller als primeClosure(l)
         // 1 bedeutet: kŸrzester nzd

         if ( l[size(l)] >= 0 && del >= 0 )
         {
            del = del + l[size(l)];
         }
         else
         { del = -1; }
         deli = l[size(l)],deli;

         l = l[1..size(l)-1];
         resu = list(l[size(l)]) + resu;
         gens = closureGenerators(l);         //computes algebra(!) generators

         //NOTE: gens[i]/gens[size(gens)] expresses the ith variable of resu[1] 
         //(the normalization) as fraction of elements of the basering;
         //the variables of resu[1] are algebra generators.
         //gens[size(gens)] is a non-zero divisor of basering/i

         //divide by the greatest common divisor:
         gg = gcd( gens[1],gens[size(gens)] );
         for(j=2; j<=size(gens)-1; j++)
         {
            gg=gcd(gg,gens[j]);
         }
         for(j=1; j<=size(gens); j++)
         {
            gens[j]=gens[j]/gg;
         }
         Gens = list(gens) + Gens;

/*       ### Da die gens Algebra-Erzeuger sind, ist reduce nach Bestimmung 
         der Algebra-Variablen T(i) nicht zulŠssig!
         for(i=1;i<=size(gens)-1;i++) 
         {
            gens[i]= reduce(gens[i],std(gens[size(gens)]));
         }
         for(i=size(gens)-1; i>=1; i--)
         {
            if(gens[i]==0)
            { gens = delete(gens,i); }
         }
*/
         if( defined(ker) ) { kill ker; }
      }

      if ( del >= 0 )
      {
         int mul = iMult(prim);
         del = del + mul;
      }
      else
      { del = -1; }
      deli = deli[1..size(deli)-1];
      Resu = resu,Gens,list(deli,del); 
      int sr = size(resu);

      if ( y >= 0 )
      {"";
"// 'normalC' created a list, say nor, of three lists: 
// nor[1] resp. nor[2] are lists of",sr,"ring(s) resp. ideal(s)
// and nor[3] is a list of an intvec and an integer.
// To see the result, type
     nor;
// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
     def R1 = nor[1][1]; setring R1;  norid; normap;
// for the other rings type first setring r; (if r is the name of your 
// original basering) and then continue as for the first ring;
// Ri/norid is the affine algebra of the normalization of the i-th 
// component r/P_i (where P_i is an associated prime of the input ideal)
// and normap the normalization map from r to Ri/norid;
//    nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
// elements g1..gk of r such that the gj/gk generate the integral 
// closure of r/P_i as sub-algebra in the quotient field of r/P_i, with
// gj/gk being mapped by normap to the j-th variable of Ri;
//    nor[3] shows the delta-invariant of each component and of the input
// ideal (-1 means infinite, and 0 that r/P_i resp. r/input is normal).";
      }
      return(Resu);
   }

//-------- The general case without computation of the generators -----------
// (attrib(id,"xy"),1) sets attrib xy to TRUE and (attrib(id,"xy"),0) to FALSE
// We use the following attributes:
//   attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
//   attrib(id,"isCompleteIntersection");  //--- complete intersection  
//   attrib(id,"isHypersurface");          //--- hypersurface 
//   attrib(id,"isEquidimensional",-1);    //--- equidimensional ideal 
//   attrib(id,"isPrim");                  //--- prime ideal 
//   attrib(id,"isRegInCodim2");           //--- regular in codimension 2 
//   attrib(id,"isIsolatedSingularity";    //--- isolated singularities 
//   attrib(id,"onlySingularAtZero");      //--- only singular at 0   

 //------------------- first set the attributes ----------------------
   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )   
   {
      if( attrib(id,"isCohenMacaulay")==1 )
      {
         attrib(id,"isEquidimensional",1);
      }
   }
   else
   {
      attrib(id,"isCohenMacaulay",0);
   }  

   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
   {
      if(attrib(id,"isCompleteIntersection")==1)
      {
         attrib(id,"isCohenMacaulay",1);
         attrib(id,"isEquidimensional",1);
      }
   }
   else
   {
      attrib(id,"isCompleteIntersection",0);
   }
    
   if( typeof(attrib(id,"isHypersurface"))=="int" )
   {
      if(attrib(id,"isHypersurface")==1)
      {
         attrib(id,"isCompleteIntersection",1);
         attrib(id,"isCohenMacaulay",1);
         attrib(id,"isEquidimensional",1);
      }
   }
   else
   {
      attrib(id,"isHypersurface",0);
   }
   
   if( ! (typeof(attrib(id,"isEquidimensional"))=="int") )
   {
         attrib(id,"isEquidimensional",0);
   }
   
   if( typeof(attrib(id,"isPrim"))=="int" )
   {
      if(attrib(id,"isPrim")==1)
      {
         attrib(id,"isEquidimensional",1);
      }
   }
   else
   {
      attrib(id,"isPrim",0);
   }
   
   if( ! (typeof(attrib(id,"isRegInCodim2"))=="int") )
   {
         attrib(id,"isRegInCodim2",0);
   }
   
   if( ! (typeof(attrib(id,"isIsolatedSingularity"))=="int") )
   {
         attrib(id,"isIsolatedSingularity",0);
   }
 
   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
   {
      if(attrib(id,"onlySingularAtZero")==1)
      {
         attrib(id,"isIsolatedSingularity",1);
      }
   }
   else
   {
      attrib(id,"onlySingularAtZero",0);
   }

   //-------------- compute equidimensional decomposition --------------------
   //If the method "equidim" is given, compute the equidim decomposition
   //and goto the next step (no normalization 
   //ACHTUNG: equidim berechnet bei nicht reduzierten id die eingebetteten
   //Komponenten als niederdim Komponenten, waehrend diese bei primdecGTZ
   //nicht auftauchen: ideal(x,y)*xy
   //this is default for nvars > 2 and char = 0 or car > 30

   if( withEqui )
   {
      withPrim = 0;                 //this is used to check later that prim
                                    //contains equidim but not prime components
      if( y >= 1 )
      {
         "// We use method 'equidim'";
      }
      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)
      {  "";
         "// number of equidimensional components:", size(prim);
      }
      if ( !nfac )
      {
        intvec opt = option(get);
        option(redSB);
        for(j=1; j<=size(prim); j++)    
        {
           keepresult = keepresult+facstd(prim[j]);
        }
        prim = keepresult;
        if ( size(prim) == 0 )
        {
          prim=ideal(0);     //Bug in facstd, liefert leere Liste bei 0-Ideal
        }
        
        if(y>=1)
        {  "";
         "// number of components after application of facstd:", size(prim);
        }
        option(set,opt);
      }
   }

   //------------------- compute associated primes -------------------------
   //the case where withEqui = 0, here the min. ass. primes are computed
   //start with the computation of the minimal associated primes:

   else                                          
   {
    if( isPrim )
    {
      if( y >= 0 )
      {
         "// ** WARNING: result is correct if ideal is prime";
         "// or equidimensional (not checked) **";
         "// disable option \"isPrim\" to decompose ideal into prime";
         "// or equidimensional components";"";
      }
      if( y >= 1 )
      {
        "// We use method 'isPrim'";"";
      }
      prim[1]=id;
    }
    else
    {
      withPrim = 1;                 //this is used to check later that prim
                                    //contains prime but not equidim components
      if( y >= 1 )
      {
         "// We use method 'prim'";
      }

      if( typeof(attrib(id,"isPrim"))=="int" )
      {
         if(attrib(id,"isPrim")==1)
         {
            prim[1]=id;
         }
         else
         {
            if( nfac )
            { prim=minAssGTZ(id,1); }     //does not use factorizing groebner
            else
            { prim=minAssGTZ(id); }       //uses factorizing groebner
         }
      }
      else
      {
            if( nfac )
            { prim=minAssGTZ(id,1); }
            else
            { prim=minAssGTZ(id); }
      }
      if(y>=1)
      {  "";
         "// number of irreducible components:", size(prim);
      }
    }
   }

   //----- for each component (equidim or irred) compute normalization -----
   int sr, skr, del;
   intvec deli;  
   int sp = size(prim);     //size of list prim (# irred or equidim comp)

   for(i=1; i<=sp; i++)
   {
      if(y>=1)
      {  "";
         "// computing the normalization of component",i;
         "   ----------------------------------------";
      }
      //-------------- first set attributes for components ------------------
      attrib(prim[i],"isEquidimensional",1);
      if( withPrim )
      {
         attrib(prim[i],"isPrim",1);
      }
      else
      { attrib(prim[i],"isPrim",0); }

      if(attrib(id,"onlySingularAtZero")==1)
      { attrib(prim[i],"onlySingularAtZero",1); }
      else
      { attrib(prim[i],"onlySingularAtZero",0); }

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

      if( attrib(id,"isHypersurface")==1 )
      {
         attrib(prim[i],"isHypersurface",1); 
         attrib(prim[i],"isCompleteIntersection",1);
         attrib(prim[i],"isCohenMacaulay",1);
      }
      else
      { attrib(prim[i],"isHypersurface",0); }

      if ( sp == 1)         //the case of one component: copy attribs from id
      {
        if(attrib(id,"isRegInCodim2")==1)
        {attrib(prim[i],"isRegInCodim2",1); }
        else
        {attrib(prim[i],"isRegInCodim2",0); }

        if(attrib(id,"isCohenMacaulay")==1)
        {attrib(prim[i],"isCohenMacaulay",1); }
        else
        {attrib(prim[i],"isCohenMacaulay",0); }

        if(attrib(id,"isCompleteIntersection")==1)
        {attrib(prim[i],"isCompleteIntersection",1); }
        else
        {attrib(prim[i],"isCompleteIntersection",0); }
      }
      else
      {  
        attrib(prim[i],"isRegInCodim2",0);
        attrib(prim[i],"isCohenMacaulay",0);
        attrib(prim[i],"isCompleteIntersection",0);
      }

      //------ Now compute the normalization of each component ---------
      //note: for equidimensional components the "splitting tools" can
      //create further decomposition
      //We now start normalizationPrimes with 
      //ihp = partial normalisation map = identity map = maxideal(1)
      //del = partial delta invariant = 0 
      //deli= intvec of partial delta invariants of components
      //in normalizationPrimes all the work is done:

      keepresult = normalizationPrimes(prim[i],maxideal(1),0,0);

      for(j=1; j<=size(keepresult)-1; j++)
      {
         result=insert(result,keepresult[j]);
      }
      skr = size(keepresult);

      //compute delta:             
      if( del >= 0 && keepresult[skr][1] >=0 )
      {
         del = del + keepresult[skr][1];
      }
      else
      {  
         del = -1; 
      }
      deli = keepresult[skr][2],deli;

      if ( y>=1 )
      {
           "// delta of component",i; keepresult[skr][1];
      } 
   }
   sr = size(result);

   // -------------- Now compute intersection multiplicities -------------
   //intersection multiplicities of list prim, sp=size(prim).
      if ( y>=1 )
      {
        "// Sum of delta for all components"; del;
        if ( sp>1 )
        {
           "// Compute intersection multiplicities of the components";
        }
      }   

      if ( sp > 1 )
      {
        int mul = iMult(prim);
        if ( mul < 0 )
        {
           del = -1;
        }
        else
        {
           del = del + mul; 
        }
      }
   deli = deli[1..size(deli)-1];
   result = result,list(deli,del);

//--------------- Finally print comments and return ------------------
   if ( y >= 0)
   {"";
"// 'normalC' created a list, say nor, of two lists:
// nor[1] is a list of",sr,"ring(s), nor[2] a list of an intvec and an int.
// To see the result, type
      nor;
// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
      def R1 = nor[1][1];  setring R1;  norid;  normap;
// and similair for the other rings nor[1][i];
// Ri/norid is the affine algebra of the normalization of the i-th component 
// r/P_i (where P_i is a prime or an equidimensional ideal of the input ideal) 
// and normap the normalization map from the basering to Ri/norid;
//    nor[2] shows the delta-invariant of each component and of the input
// ideal (-1 means infinite, -2 that delta of a component was not computed).";
   }
   return(result);
}

example
{ "EXAMPLE:";
   printlevel = printlevel+1;
   echo = 2;
   ring s = 0,(x,y),dp;
   ideal i = (x2-y3)*(x2+y2)*x;

   list nor = normalC(i);

   nor; 
   // 2 branches have delta = 1, and 1 branch has delta = 0
   // the total delta invariant is 13

   def R2 = nor[1][2];  setring R2;
   norid; normap;
   
   echo = 0;
   printlevel = printlevel-1;
   pause("   hit return to continue"); echo=2;

   ring r = 2,(x,y,z),dp;
   ideal i = z3-xy4;
   nor = normalC(i);  nor;  
   // the delta invariant is infinite 
   // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
   // in its quotient field Quot(r/i)

   // the normalization as affine algebra over the ground field:             
   def R = nor[1][1]; setring R;
   norid; normap;

   echo = 0;  
   pause("   hit return to continue");echo = 2;

   setring r;
   nor = normalC(i, "withGens", "prim");    // a different algorithm
   nor;
}

//////////////////////////////////////////////////////////////////////////////
//closureRingtower seems not to be used anywhere
static proc closureRingtower(list L)
"USAGE:    closureRingtower(list L); L a list of rings
CREATE:   rings R(1),...,R(n) such that R(i)=L[i] for all i
EXAMPLE:  example closureRingtower; shows an example
"
{
  int n=size(L);
  for (int i=1;i<=n;i++)
    {
      if (defined(R(i))) {
        string s="Fixed name R("+string(i)+") leads to conflict with existing "
              +"object having this name";
        ERROR(s);
      }
      def R(i)=L[i];
      export R(i);
    }

  return();
}
example
{
  "EXAMPLE:"; echo=2;
  ring R=0,(x,y),dp;
  ideal I=x4,y4;
  list L=primeClosure(ReesAlgebra(I)[1]);
  L=delete(L,size(L));
L;
  closureRingtower(L);
  R(1);
  R(4);
  kill R(1),R(2),R(3),R(4);
}

//                Up to here: procedures for normalC
///////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
//                From here: miscellaneous procedures

// Used for timing and comparing the different normalization procedures.
// Option (can be entered in any order)
// "normal"   -> uses the new algortihm (normal)
// "normalP"  -> uses normalP
// "normalC"  -> uses normalC, without "withGens" option
// "primCl"   -> uses normalC, with option "withGens".
// "111"      -> checks the output of normalM using norTest.
// "p"        -> compares the output of norM with the output of normalP 
//               ("normalP" option must also be set).
// "pc"       -> compares the output of norM with the output of normalC with 
//               option "withGens"
//               ("primCl" option must also be set).
proc timeNormal(ideal I, list #){
  //--------------------------- define the method --------------------------- 
  int isPrim, useRing;
  int decomp = -1;
  int norM, norC, norP, primCl;
  int checkP, check111, checkPC;
  int i;
  ideal U1, U2, W;
  poly c1, c2;
  int ch;
  string check;
  string method;                //make all options one string in order to use
                                //all combinations of options simultaneously
  for ( i=1; i <= size(#); i++ )
  { 
    if ( typeof(#[i]) == "string" ) 
    {
      method = method + #[i];
    }
  }
  if ( find(method, "normal"))
  {norM = 1;}  
  if ( find(method, "normalP") and (char(basering) > 0))
  {norP = 1;}  
  if ( find(method, "normalC"))
  {norC = 1;}      
  if ( find(method, "primCl"))
  {primCl = 1;}  
  if ( find(method, "isprim") or find(method,"isPrim") )
  {decomp = 0;}
  if ( find(method, "p") )
  {checkP = 1;}  
  if ( find(method, "pc") )
  {checkPC = 1;}  
  if ( find(method, "111") )
  {check111 = 1;}    
      
  int tt;
  if(norM){
    tt = timer; 
    if(decomp == 0){
      "Running normal(useRing, isPrim)...";
      list a1 = normal(I, "useRing", "isPrim");
      "Time normal(useRing, isPrim): ", timer - tt;
    } else {
      "Running normal(useRing)...";
      list a1 = normal(I, "useRing");
      "Time normal(useRing): ", timer - tt;
    }
    "";
  }
  if(norP){
    tt = timer; 
    if(decomp == 0){
      "Running normalP(isPrim)...";
      list a2 = normalP(I, "isPrim");
      "Time normalP(isPrim): ", timer - tt;
    } else {
      "Running normalP()...";
      list a2 = normalP(I);
      "Time normalP(): ", timer - tt;
    }
    "";
  }
  
  if(norC){
    tt = timer; 
    if(decomp == 0){
      "Running normalC(isPrim)...";
      list a3 = normalC(I, "isPrim");
      "Time normalC(isPrim): ", timer - tt;
    } else {
      "Running normalC()...";
      list a3 = normalC(I);
      "Time normalC(): ", timer - tt;
    }
    "";
  }
  
  if(primCl){
    tt = timer; 
    if(decomp == 0){
      "Running normalC(withGens, isPrim)...";
      list a4 = normalC(I, "isPrim", "withGens");
      "Time normalC(withGens, isPrim): ", timer - tt;
    } else {
      "Running normalC(withGens)...";    
      list a4 = normalC(I, "withGens");
      "Time normalC(withGens): ", timer - tt;
    }
    "";
  }
  
  if(check111 and norM){
    "Checking output with norTest...";
    "WARNING: this checking only works if the original ideal was prime.";
    list norL = list(a1[1][1][3]);
    norTest(I, list(norL));
    "";
  }
  
  if(checkP and norP and norM){
    "Comparing with normalP output...";
    "WARNING: this checking only works if the original ideal was prime.";
    U1 = a1[1][1][1];
    c1 = a1[1][1][2];
    U2 = a2[1][1];
    c2 = a2[1][1][size(a2[1][1])];
    W = changeDenominator(U1, c1, c2, groebner(I));
    qring q = groebner(I);
    ideal U2 = fetch(r, U2);
    ideal W = fetch(r, W);
    ch = 0;
    if(size(reduce(U2, groebner(W))) == 0){
      "U2 c U1";
      ch = 1;
    }
    if(size(reduce(W, groebner(U2))) == 0){
      "U1 c U2";
      ch = ch + 1;
    }
    if(ch == 2){
      "Output of normalP is equal.";
    } else {
      "ERROR: Output of normalP is different.";
    }
    "";
    setring r;
    kill q;
  }  

  if(checkPC and norM and primCl){    
    "Comparing with primeClosure output...";
    "WARNING: this checking only works if the original ideal was prime.";
    // primeClosure check
    U1 = a1[1][1][1];
    c1 = a1[1][1][2];
    U2 = a4[2][1];
    c2 = a4[2][1][size(a4[2][1])];
    W = changeDenominator(U1, c1, c2, groebner(I));
    qring q = groebner(I);
    ideal U2 = fetch(r, U2);
    ideal W = fetch(r, W);
    ch = 0;
    if(size(reduce(U2, groebner(W))) == 0){
      "U2 c U1";
      ch = 1;
    }
    if(size(reduce(W, groebner(U2))) == 0){
      "U1 c U2";
      ch = ch + 1;
    }
    if(ch == 2){
      "Output of normalC(withGens) is equal.";
    } else {
      "ERROR: Output of normalC(withGens) is different.";
    }
    "";
    setring r;
    kill q;    
  }
}

///////////////////////////////////////////////////////////////////////////
proc norTest (ideal i, list nor, list #)
"USAGE:   norTest(i,nor,[n]); i=prime ideal, nor=list, n=optional integer
ASSUME:  nor is the resulting list of normal(i) (any options) or the result
         of normalP(i,"withRing"). In particular, the ring nor[1][1] contains
         the ideal norid and the map normap: basering/i --> nor[1][1]/norid.
RETURN:  an intvec v such that:
@format
         v[1] = 1 if the normap is injective and 0 otherwise
         v[2] = 1 if the normap is finite and 0 otherwise
         v[3] = 1 if nor[1][1]/norid is normal and 0 otherwise 
@end format
         If n=1 (resp n=2) only v[1] (resp. v[2]) is computed and returned
THEORY:  The procedure can be used to test whether the computation of the 
         normalization was correct: basering/i --> nor[1][1]/norid is the 
         normalization of basering/i iff v=1,1,0. 
NOTE:    For big examples it can be hard to fully test correctness; the
         partial test norTest(i,nor,2) is usually fast
EXAMPLE: example norTest; shows an example
"
{
//### Sollte erweitert werden auf den reduziblen Fall: einen neuen affinen
// Ring nor[1][1]+...+nor[1][r] (direkte Summe) erzeugen, map dorthin
// definieren und dann testen.

    int prl = printlevel - voice + 2;
    int a,b,d;
    int n,ii;
    if (size(#) > 0) {  n = #[1];  }

    def BAS = basering;

    //### make a copy of nor to have a cpoy of nor[1][1]  (not a reference to)
    // in order not to overwrite norid and normap. 
    // delete nor[2] (if it contains the module generators, which are not used)
    // s.t. newnor does not belong to a ring.

    list newnor = nor;         
    if ( size(newnor) == 3 )
    {
       newnor = delete(newnor,2);
    }
    qring QAS = std(i);
    
    def R = newnor[1][1]; 
    setring R;
    int nva = nvars(R);
    string svars = varstr(R);
    string svar;

    norid = interred(norid);

    //--------- create new ring with one dp block keeping weights ------------
    list LR = ringlist(R);
    list g3 = LR[3];
    int n3 = size(g3);
    list newg3;
    intvec V;

    //--------- check first whether variables Z(i),...,A(i) exist -----------
    for (ii=90; ii>=65; ii--)  
    {
       if ( find(svars,ASCII(ii)+"(") == 0 )        
       {
          svar = ASCII(ii);  break;
       }
    }
    if ( size(svar) != 0 )
    {
        for ( ii = 1; ii <= nva; ii++ )
        {
            LR[2][ii] = svar+"("+string(ii)+")";
            V[ii] = 1;
        }
    }
    else
    {       
        for ( ii = 1; ii <= nva; ii++ )
        {
           LR[2][ii] = "Z("+string(100*nva+ii)+")";
           V[ii] = 1;
        }
    }
    
    if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
    {
       list gm = g3[n3];       //last blockis module ordering
       newg3[1] = list("dp",V);
       newg3 = insert(newg3,gm,size(newg3));
    }
    else
    {
       list gm = g3[1];              //first block is module ordering
       newg3[1] = list("dp",V);
       newg3 = insert(newg3,gm);
    }
    LR[3] = newg3;
//LR;"";
    def newR = ring(LR);

    setring newR;
    ideal norid = fetch(R,norid);
    ideal normap = fetch(R,normap);
    if( defined(lnorid) )  { kill lnorid; }     //um ** redefinig zu beheben
    if( defined(snorid) )  { kill snorid; }     //sollte nicht noetig sein

    //----------- go to quotient ring for checking injectivity -------------
//"mstd";
    list lnorid = mstd(norid);
    ideal snorid = lnorid[1];
//"size mstdnorid:", size(snorid),size(lnorid[2]);
//"size string mstdnorid:", size(string(snorid)),size(string(lnorid[2]));
    qring QR = snorid;
    ideal qnormap = fetch(newR,normap);
    //ideal qnormap = imap(newR,normap);
    //ideal qnormap = imap(R,normap);
    map Qnormap = QAS,qnormap;    //r/id --> R/norid
    
    //------------------------ check injectivity ---------------------------
//"injective:";
    a = is_injective(Qnormap,QAS);          //a. Test for injectivity of Qnormap
    dbprint ( prl, "injective: "+string(a) );
    if ( n==1 ) { return (a); }
   a;

    //------------------------ check finiteness ---------------------------
    setring newR;
    b = mapIsFinite(normap,BAS,lnorid[2]);  //b. Test for finiteness of normap
    dbprint ( prl, "finite: "+string(b) );
    if ( n==2 ) { return (intvec(a,b)); }
   b;

    //------------------------ check normality ---------------------------
    list testnor = normal(lnorid[2],"isPrim","noFac");
    //### Problem: bei mehrfachem Aufruf von norTest gibt es 
    // ** redefining norid & ** redefining normap
    //Dies produziert Fehler, da alte norid und normap ueberschrieben werden
    //norid und normap werden innnerhalb von proc computeRing ueberschrieben
    //Die Kopie newR scheint das Problem zu loesen

    printlevel=prl;
    d = testnor[size(testnor)][2];             //d = delta
    kill testnor;                              //### sollte ueberfluessig sein
    int d1 = (d==0);                           //d1=1 if delta=0
    dbprint ( prl, "delta: "+string(d) );
    return(intvec(a,b,d1));
}
example
{ "EXAMPLE:"; echo = 2;
   int prl = printlevel;
   printlevel = -1;
   ring r = 0,(x,y),dp;
   ideal i = (x-y^2)^2 - y*x^3;
   list nor = normal(i);
   norTest(i,nor);                //1,1,1 means that normal was correct

   ring s = 2,(x,y),dp;
   ideal i = (x-y^2)^2 - y*x^3;
   nor = normalP(i,"withRing");
   norTest(i,nor);               //1,1,1 means that normalP was correct
   printlevel = prl; 
}

///////////////////////////////////////////////////////////////////////////
//
//                            EXAMPLES
//
///////////////////////////////////////////////////////////////////////////
/*
//commands for computing the normalization:
// options for normal:  "equidim", "prim"
//                      "noDeco", "isPrim", "noFac"
//                       (prim by default)                      
// options for normalP: "withRing", "isPrim" or "noFac"
// options for normalC: "equidim", "prim", "withGens", 
//                      "noDeco", "isPrim", "noFac"

//Commands for testing 'normal' 
 list nor = normal(i); nor;            
 list nor = normal(i,"isPrim");nor; 
 list nor = normal(i,"equidim");nor;
 list nor = normal(i,"prim");nor;
 list nor = normal(i,"equidim","noFac");nor;
 list nor = normal(i,"prim","noFac");nor;

//Commands for testing 'normalP' in positive char
 list nor = normalP(i);nor;              //withGens but no ringstructure
 list nor = normalP(i,"withRing"); nor;  //compute the ringstructure
 list nor = normalP(i,"isPrim"); nor;    //if i is known to be prime

//Commands for testing 'normalC' 
 list nor = normal(i); nor;            
 list nor = normal(i,"withGens");nor;
 list nor = normal(i,"isPrim");nor; 
 list nor = normal(i,"equidim");nor;
 list nor = normal(i,"prim");nor;
 list nor = normal(i,"equidim","noFac");nor;
 list nor = normal(i,"prim","noFac");nor;

//Commands for testing correctness (i must be prime):
list nor = normalP(i,"withRing","isPrim");
list nor = normal(i,"isPrim");
norTest(i,nor);       //full test for not too big examples (1,1,1 => ok)
norTest(i,nor,2);     //partial test for big examples (1,1 => ok)
factorize(i[1]);      //checks for irreducibility

//----------------------Test Example for charp -------------------
//Zu tun:
//### nach minor nur std statt mstd verwenden
//***hat bei keinem Beisp etwas gebracht -> wieder zurueck
//### wenn interred ok, dann wieder einsetzen (am Schluss)
//### bottelnecks bei maps beheben
//### minor verbessern
//### preimage verbessern (Ist imm Kern map oder imap verwendet?)
//### Gleich in Ordnung dp wechseln, ringlist verwenden
//### interred ev nur zum Schluss 
//    (z.B. wenn nacher std; wenn nacher minor: testen )

//Zeiten mit normalV5.lib (mstd aktiv, interred inaktiv)

//SWANSON EXAMPLES: (Macaulay2, icFracP=normalP, icFractions<->normal)
//---------------------------------------------------------------------
//1. Series Fp[x,y,u,v]/(x2v-y2u)
//-------------------------------
//characteristic p   2   3    5    7    11   13   17   37   97 
//icFracP          0.04 0.03 0.04 0.04 0.04 0.05 0.05 0.13 0.59  Mac
//normalP           0   0    0    0     0    0    0    0   1    Sing
//icFractions      0.08 0.09 0.09 0.09 0.14 0.15 0.15 0.15 0.15  Mac
//normal             0   0    0    0     0    0    0    0    0   Sing

2. Series Fp[u, v, w, x, y, z]/u2x4+uvy4+v2z4
//--------------------------------------------
//characteristic p 2    3    5    7   11 
//icFracP         0.07 0.22 9.67 143 12543 
//normalP          0    0    5   42  1566
//icFractions     1.16   *    *   *    *       *: > 6h
//normal            0    0    0   0    0

//3. Series Fp[u, v, w, x, y, z]/(u2xp+uvyp+v2zp) 
//-----------------------------------------------
//characteristic p  2    3    5    7    11   13  17 19 23 
//icFracP          0.06 0.07 0.09 0.27 1.81 4.89 26 56 225 
//normalP          0     0    0    0    1    2  6  10  27
//icFractions      0.16 1.49 75.00 4009  *    *   *  *  * 
//normal            0     0    2    836
//### p=7 normal braucht 807 sec in:
// ideal endid  = phi1(endid);      //### bottelneck'

//1.
int p = 2;  ring r = p,(u,v,x,y,z),dp; ideal i = x2v-y2u;
//2.
int p = 7; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2x4+uvy4+v2z4;
//3.
int p=23; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2*x^p+uv*y^p+v2*z^p;

//IRREDUCIBLE EXAMPLES:
//---------------------  
//timing for MacBookPro 2.2GHz Intel Core 2 Duo, 4GB Ram
//Sing. ix86Mac-darwin version 3-1-0 (3100-2008101314)  Oct 13 2008 14:46:59
//if no time is given: < 1  sec

//Apply:
list nor = normal(i,"isPrim"); nor;
list nor = normalP(i,"withRing","isPrim"); nor; 
def R=nor[1][1]; setring R; norid; normap;
setring r;
norTest(i,nor);

int tt = timer; 
list nor = normalP(i,"withRing","isPrim"); nor; 
timer-tt;
int tt = timer; 
list nor = normal(i,"isPrim");
timer-tt;

ring r=19,(x,y,u,v),dp;    //delta -1
ideal i=x2v-y2u; 
//norTest 2 sec

ring r=2,(y,x2,x1),lp;     //delta -1
ideal i=y^4+y^2*x2*x1+x2^3*x1^2+x2^2*x1^3;
//### norid hat 1 Element nach interred

ring r  = 11,(x,y,z),wp(2,1,2); //alles < 1 sec
ideal i=z3 - xy4 + x2;          //not reduced, delta =0 ok
ideal i=y4+x5+y2x;              //not reduced, delta -1
//interred verkleinert norid

ring r=3,(u,v,x,y,z),dp;   //delta -1
ideal i=u2x3+uvy3+v2z3;

ring r=3,(u,v,x,y,z),dp;   //delta -1
ideal i=u2x4+uvy4+v2z4;
//norTest(i,nor);  0 sec, norTest(i,nor) haengt!

ring r=5,(u,v,x,y,z),dp;   //delta -1
ideal i=u2x6+uvy6+v2z6;  
//normalP 5sec, normalC 1sec
//V5: norTest(i,nor); 45 sec bei normalP, V6 12 sec
//28 sec bei normal

ring r=5,(u,v,x,y,z),dp;   //delta -1
ideal i=u2x5+uvy5+v2z5;
//normalP 1sec, normalC 1 sec, 
//norTest lange: minor(jacob(I),h,J) 193 (308)sec, haengt dann bei M = std(M); 
//norTest(i,nor,2); verwenden!
//Sing 3.0-4 orig  >9h! hŠngt bei Q = mstd(Q)[2];

ring r=2,(y,x),wp(12,5);  //delta 3
ideal i=y5+y2x4+y2x+yx2+x12; 
//normalP 0 sec (Test 0 sec), normalC 2 sec (Test 2 sec)
//normalC withGens (ohne interred) 0sec

ring r=2,(y,x),dp;       //delta= 22
ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9;
//normalP 1sec, interred verkleinert norid betraechtlich
//normalC haengt bei minor, ideal im loop wird zu gross ###
//interred bei normalC vergroeesert string um Faktor 4000!
//withGens haengt bei interred in loop 4 (> 10 h) oder 
//(nach Ausschalten von interred) bei 
//int delt=vdim(std(modulo(f,ideal(p)))); (>?h)

//Leonard1: (1. Komponente von Leonard),  delta -1 
ring r=2,(v,u,z,y,x),dp;
ideal i = z3+zyx+y3x2+y2x3, uyx+z2,uz+z+y2x+yx2, u2+u+zy+zx,
          v3+vux+vz2+vzyx+vzx+uz3+uz2y+z3+z2yx2;
//normalP 5 sec (withRing 9 sec), norTest(i,nor,2); 45 sec
//normalC 102sec, 99sec
//### Zeit wird bei ideal Ann = quotient(SM[2],SL[1]); und bei
// f  = quotient(p*J,J); verbraucht
//withGens (ohne interred) 131sec, norTest(i,nor,2); 2min25sec
//norTest(i,nor,2);  45 sec

 ring r=2,(y,x),wp(25,21); //Leonard2, delta 232
 ring r=2,(y,x),dp;
 ideal i=
 y^21+y^20*x +y^18*(x^3+x+1) +y^17*(x^3+1) +y^16*(x^4+x)
 +y^15*(x^7+x^6+x^3+x+1) +y^14*x^7 +y^13*(x^8+x^7+x^6+x^4+x^3+1)
 +y^12*(x^9+x^8+x^4+1) +y^11*(x^11+x^9+x^8+x^5+x^4+x^3+x^2)
 +y^10*(x^12+x^9+x^8+x^7+x^5+x^3+x+1)
 +y^9*(x^14+x^13+x^10+x^9+x^8+x^7+x^6+x^3+x^2+1)
 +y^8*(x^13+x^9+x^8+x^6+x^4+x^3+x) +y^7*(x^16+x^15+x^13+x^12+x^11+x^7+x^3+x)
 +y^6*(x^17+x^16+x^13+x^9+x^8+x) +y^5*(x^17+x^16+x^12+x^7+x^5+x^2+x+1)
 +y^4*(x^19+x^16+x^15+x^12+x^6+x^5+x^3+1)
 +y^3*(x^18+x^15+x^12+x^10+x^9+x^7+x^4+x)
 +y^2*(x^22+x^21+x^20+x^18+x^13+x^12+x^9+x^8+x^7+x^5+x^4+x^3)
 +y*(x^23+x^22+x^20+x^17+x^15+x^14+x^12+x^9)
 +(x^25+x^23+x^19+x^17+x^15+x^13+x^11+x^5);
//normalP: dp 2sec withRing 8sec, 
//wp 4sec, withRing:51sec Zeit in lin = subst(lin, var(ii), vip); in elimpart ),
//norTest(i,nor,2): haengt bei mstd(norid);
//### normalC: (m. interred): haengt bei endid = interred(endid);
//GEFIXTES INTERRED ABWARTEN. Dann interred aktivieren
//interred(norid) haengt u. mst(norid) zu lange
//(o.interred): haengt bei  haengt bei list SM = mstd(i); 
//ideal in der Mitte zu gross 
//i = Ideal (size 118, 13 var) fuer die neue Normalisierung

REDUCIBLE EXAMPLES:
------------------
//Apply:
int tt = timer; 
list nor=normalP(i,"isPrim","withRing");
timer-tt;

list nor = normal(i); nor;
list nor = normalC(i); nor;
list nor = normalC(i, "withGens"); nor;
list nor = normalP(i,"withRing"); nor;
list nor = normalP(i); nor; 
def R=nor[1][1]; setring R; norid; normap;

//Leonhard 4 Komponenten, dim=2, delta: 0,0,0,-1
ring r=2,(v,u,z,y,x),dp;      //lp zu lange
ideal i=z3+zyx+y3x2+y2x3, uyx+z2, v3+vuyx+vux+vzyx+vzx+uy3x2+uy2x+zy3x+zy2x2;
//normalP: 19 sec, withRing: 22 sec
//normalC ohne (mit) interred: 112 (113)sec, equidim: 99sec
//normalC 1. mal 111 sec, (2.mal) 450sec!! 3.mal 172 sec
//(unterschiedlich lange primdec, mit Auswirkungen)
//char 19: normalC: 15sec , withGens: 14sec (o.interr.)

//----------------------Test Example for special cases -------------------
int tt = timer; 
list nor=normalP(i,"withRing");nor;
//list nor=normalP(i,"withRing", "isPrim");nor;
timer-tt;
def R1 = nor[1][1]; setring R1;  norid; normap; interred(norid);
setring r;

int tt = timer; 
list nor=normal(i,"isPrim");nor;
timer-tt;

ring r = 29,(x,y,z),dp;
ideal i = x2y2,x2z2;       //Nicht equidimensional, equidim reduziert nicht, ok
ideal i  = xyz*(z3-xy4);   //### interred(norid) verkuerzt
//je 0 sec

ideal j = x,y;
ideal i = j*xy;
equidim(i);            
//hat eingebettete Komponente, equidim rechnet wie in Beschreibung (ok)
    
ring r  = 19,(x,y),dp;
   ideal i = x3-y4;                   //delta = 3
   ideal i = y*x*(x3-y4);             //delta = 11; 0,0,3
   ideal i = (x2-y3)*(x3-y4);         //delta = 13; 1,3
   ideal i = (x-y)*(x3+y2)*(x3-y4);   //delta = 23; 0,1,3     
   ideal i = (x-1)*(x3+y2)*(x2-y3);   //delta = 16; 0,1,1
   ideal i = (x-y^2)^2 - y*x^3;       //delta = 3
   //singularities at not only at 0, hier rechnet equidim falsch

// -------------------------- General Examples  ---------------------------//Huneke, irred., delta=2 (Version 3-0-4: < 1sec) 
//Version 3-0-6 default: 1sec, mit gens 2sec, mit delta 5 sec
//(prim,noFac):ca 7 Min, prim:ca 10 min(wg facstd)
//
// "equidim" < 1sec irred. 5sec
// ring r=31991,(a,b,c,d,e),dp;
ring r=2,(a,b,c,d,e),dp;                    //delta=2
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;
//normalC: char 2, 31991: 0 sec (isPrim); char 2, equidim: 7 sec
//norTest(i,nor,2); 1sec
//normalP char 2: 1sec (isPrim) 
//size(norid); size(string(norid));21 1219 interred(norid): 21 1245 (0 sec)

int tt = timer; 
list nor=normalC(i);nor;
timer-tt;

list nor = normalP(i,"isPrim");

//Vasconcelos irred., delta -1 (dauert laenger)
//auf macbook pro = 20 sec mit alter Version, 
//Sing 3-0-6:
// Char 32003: "equidim" 30 sec, "noFac": 30sec
//gens: nach 9 min abgebr (haengt in Lin = ideal(T*syzf);) !!!! Hans zu tun
//Char 2: default (charp) 2 sec, normalC ca 30 sec
//ring r=32003,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
ring r=2,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3;
//normalC: char 2: 22,  sec (mit und ohne isPrim)
//normalP char 2: 0sec (isPrim)      o. interred
//char 32003: ### haengt in ideal endid  = phi1(endid);

//-------------------------------------------------------
//kleine Beispiele:

//Theo1, irred, delta=-1
//normalC: 1sec, normalP: 3 sec
ring r=32003,(x,y,z),wp(2,3,6); //dim 2,dim slocus 1
ideal i=zy2-zx3-x6;
//normalC: char 2,19,32003: 0  sec (isPrim)
//normalP (isPrim) char 2,19: 0sec, char 29: 1sec

//Theo1a, CohenMacaulay regular in codim 2, dim slocus=1, delta=0 
//normalC: 0 sec, normalP: haegt in K=preimage(R,phi,L);
ring r=32003,(x,y,z,u),dp;
ideal i=zy2-zx3-x6+u2;
//normalC: char 2,32003: 0  sec (isPrim)
//normalP (isPrim) char 2: 0sec, char 19: haengt in K = preimage(Q,phi,L);

//Theo2, irreduzibel, reduziert, < 1sec, delta -1
ring r=0,(x,y,z),wp(3,4,12);
ideal i=z*(y3-x4)+x8;
//normalC: char 2,32003,0: 0  sec (isPrim)
//normalP (isPrim) char 2: 0 1sec, char 19: 1sec char 29: 7 sec

//Theo2a, reduiziert, 2-dim, dim_slocus=1, alte Version 3 sec,
//normalP ca 30 sec, normalC ca 4sec, delta -1
//ring r=32003,(T(1..4)),wp(3,4,12,17);
//ring r=11,(T(1..4)),dp;
ring r=11,(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;
//normalC: char 2,32003: 0  sec (isPrim)
//normalP (isPrim) char 2: 0sec, char 11 2se, char 19: 13sec 
//norTest 48sec in char11
//### interred verkuerzt
//char 29: haengt in K = preimage(Q,phi,L);

//Theo3, irred, 2-dim, 1-dim sing, < 1sec
ring r=11,(x,y,z),wp(3,5,15);
ideal i=z*(y3-x5)+x10;
//normalC: char 2,0: 0  sec (withRing)
//normalP (withRing) char 2,11: 0sec, char 19: 13sec norTest 12sec(char 11)

//Theo4 reducible, delta (0,0,0) -1
ring r=29,(x,y,z),dp;
ideal i=(x-y)*(x-z)*(y-z); 
//normalC: char 2,32003: 0  sec 
//normalP char withRing 2, 29: 0sec, 6sec

//Theo6
ring r=32003,(x,y,z),dp;
ideal i=x2y2+x2z2+y2z2;
//normalC: char 2,32003: 0  sec 
//normalP char withRing 2, 29: 0sec, 4sec

//Sturmfels, CM, 15 componenten, alle glatt
ring r=0,(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;
//normalC car 11, 0: 1sec, normalP 0 sec

//riemenschneider, , dim 3, 5 Komp. delta (0,0,0,0,0), -1
ring r=2,(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;
//alles 0 sec in char 2

//4 Komponenten, alle glatt, 0sec
ring r=11,(x,y,z,t),dp;
ideal i=x2z+xzt,xyz,xy2-xyt,x2y+xyt;

//dim 3, 2 Komponenten delta (-1,0), -1
ring r=2,(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;
//alles 0 sec in char 2
//---------------------------------------------------------
int tt = timer; 
list nor=normalP(i,"normalC","withRing");nor;
timer-tt;

//St_S/Y, 3 Komponenten, 2 glatt, 1 normal 
//charp haengt (in char 20) in K=preimage(R,phi,L);
//ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ring r=11,(b,s,t,u,v,w,x,y,z),dp;
ideal i=wy-vz,vx-uy,tv-sw,su-bv,tuy-bvz;
//normalC: char 2,32003: 0  sec 
//normalP char withRing 2: 1sec, char 11: 40sec

//Horrocks: cahr 0: 17 (8 in char 11) Komponenten alle normal, delta 1
//char 11: 8 Komponenten alle normal, delta -1
ring r=0,(a,b,c,d,e,f),dp;
//ring r=11,(a,b,c,d,e,f),dp; //Charp bis p = 200 ca 3sec
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;
//normalC: char 0: 1sec char 11: 0sec
//normalP: char 11: 0sec

//2sec mit normalC, in char 2 ebenfalls (char 20 mit charp >1 min)
//4 Komp. in char 2, delta (0,0,0,0) -1, char 11:delta (-1,0,0,0) -1
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;
//normalC: char 2,32003: 1 sec 
//normalP char withRing 2: 1sec, char 11: 40sec

//---------------------------------------------------------
//genus:
int tt = timer; 
list nor=normal(i, "noFac");nor;
timer-tt;

//Yoshihiko Sakai, irred, 0sec, delta = 8 
ring r=0,(x,y),dp;                    //genus 0 4 nodes and 6 cusps im P2
//ring r=7,(x,y),dp;                  //charp haengt in K = preimage(Q,phi,L)
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=23,(x,y),dp;  //genus 34, delta 2
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;
//normalC, normalP 0 sec

//Rob Koelman
//ring r=0,(x,y,z),dp;      //dim sing = 1 (nach ca 15 min abgebrochen)
ring r=32003,(x,y,z),dp;   
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;
//normalC char 32003: 10 sec, char 0 : 

//ring r=0,(x,y),dp;//genus 10  with 26 cusps (nach ca 4 min abgebrochen)
ring r=32003,(x,y),dp;    //24 sing, delta 24
ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;
//normalC 32003: 4 sec, char 0: abgebrochen bei pr = facstd(i); ###

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

ring r=2,(x,y),dp;  //genus 4, 2 Zweige, delta (13,9) 89
ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
//normalC: char 2 : 1sec, char 0: lange
//normalP char 2 withRing: 0sec

ring r=2,(y,z,w,u),dp; //2 Komp. genus -5
ideal i=y2+z2+w2+u2,w4-u4;
//normalC: char 2 : 0sec, char 0: 1sec
//normalP char 2 withRing: 0sec

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

ring r=0,(x,y,t),dp;  //irred, delta -1
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;
//char 0: 0sec

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

//---------------------------------------------------------
//Probleme mit normalC in char 2 und char 0

int tt = timer; 
list nor=normalC(i,"withRing");nor;
timer-tt;

//Mark van Hoeij
ring r=3,(x,y),dp;  //genus 19, delta 21
ideal i=y20+y13x+x4y5+x3*(x+1)^2;
//normalC: char 2 > 10 min   bei list SM = mstd(i);###
//normalP char 2 withRing: 0sec, char 11: haengt bei K = preimage(Q,phi,L);

ring r=2,(x,y),dp;  //genus 35
ideal i=y30+y13x+x4y5+x3*(x+1)^2;
//char 0 abgebrochen bei list SM = mstd(i); ### 
//char 2 nach ca 30 min
//normalC: char 2: abgebr. bei list SM = mstd(i);  //Now the work starts' 
//normalC, withGens, char 2: abgebrochen bei Q=mstd(Q)[2];
//normalP char 2 withRing: 0sec 

ring r=0,(x,y),dp;   //irred, genus 55, delta 21
ideal i=y40+y13x+x4y5+x3*(x+1)^2;
//normalC: char 2 lange
//normalP char 2 withRing: 0sec

ring r=29,(x,y,t),dp; //char 0: genus -5, 4 Komp, delta (-1,-1,0,0), -1
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;
//normalC: char 29 : 0sec, char 0: 0sec  //char 29 6 Komponenten
//normalP char 29 withRing: 1sec

//-------------------------- problematic examples ------------------------
//ring r=0,(x,y,t),dp;
ring r=32003,(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;
//char 0: lange (es liegt an den grossen Zahlen), char 32003: 0 sec

//dasselbe Beipiel in char 19: irred
ring r=0,(x,y,t),dp;
ideal i=
5x8+6x7y-3x6y2+7x5y3-6x4y4-8x3y5-5x2y6-8y8+4x7t+8x6yt+2x5y2t-6x4y3t+9x3y4t+9x2y5t-xy6t-7x6t2+7x5yt2-x4y2t2+5x3y3t2+7x2y4t2+xy5t2-3y6t2-4x5t3 -3x4yt3+2x3y2t3-7x2y3t3-6xy4t3-3y5t3-5x4t4-3x3yt4-4x2y2t4-8xy3t4 +7y4t4+x3t5+9x2yt5+9xy2t5-8y3t5-2y2t6+4xt7-7yt7-9t8;
//normalP: char 2,3: 0sec, norTest 0,2 sec, char 11 haengt bei peimage
//normalC: char 3: 0 sec, char 0: 1sec

//ring r=0,(x,y),dp;
ring r=32003,(x,y),dp; 
ideal i= 
x30y21+21x29y20+210x28y19+10x27y19+1330x27y18+190x26y18+5985x26y17 
+1710x25y17+20349x25y16+45x24y17+9690x24y16+54264x24y15+765x23y16 
+38760x23y15+116280x23y14+6120x22y15+116280x22y14+120x21y15 
+203490x22y13+30600x21y14+271320x21y13+1799x20y14+293930x21y12+107100x20y13 
+503880x20y12+12586x19y13+352716x20y11+278460x19y12+210x18y13+755820x19y11 
+54509x18y12+352716x19y10+556920x18y11+2723x17y12+923780x18y10+163436x17y11 
+293930x18y9+875160x17y10+16296x16y11+923780x17y9+359359x16y10+252x15y11 
+203490x17y8+1093950x16y9+59598x15y10+755820x16y8+598598x15y9+2751x14y10 
+116280x16y7+1093950x15y8+148610x14y9+503880x15y7+769197x14y8+13650x13y9 
+54264x15y6+875160x14y7+266805x13y8+210x12y9+271320x14y6+768768x13y7 
+40635x12y8+20349x14y5+556920x13y6+354816x12y7+1855x11y8+116280x13y5 
+597597x12y6+80640x11y7+5985x13y4+278460x12y5+353892x11y6+7280x10y7+38760x12y4
+358358x11y5+112014x10y6+120x9y7+1330x12y3+107100x11y4+264726x10y5+16660x9y6 
+9690x11y3+162799x10y4+111132x9y5+805x8y6+210x11y2+30600x10y3+146685x9y4 
+24500x8y5+1710x10y2+54236x9y3+78750x8y4+2310x7y5+21x10y+6120x9y2+58520x8y3 
+24010x7y4+45x6y5+190x9y+12509x8y2+39060x7y3+3675x6y4+x9+765x8y+15918x7y2 
+15680x6y3+204x5y4+10x8+1786x7y+12915x6y2+3500x5y3+45x7+2646x6y+6580x5y2 
+366x4y3+119x6+2562x5y+1995x4y2+10x3y3+203x5+1610x4y+324x3y2+231x4+630x3y 
+23x2y2+175x3+141x2y+85x2+16xy+24x+y+4; 
list nor = normal(i);   
//normalC: char 0: ### hŠngt in SB of singular locus JM = mstd(J);
//normalC: char 32003,"noFac","equidim": 0sec, "noFac": 1sec
// ev neues interred
genus(i);         // haengt bei int mu=vdim(std(jacob(f)));
                  //### ist das noetig?

//Singular rechnet genus richtig, auch im Fall, dass Kurve irreduzibel, 
//aber nicht absolut irreduzibel ist:
ring r = 0,(x,y),dp;
ideal i = x2+y2;      //irreduzibel /Q aber reduzibel /C (x-iy)*(x+iy)
factorize( x2+y2);    //liefert irreduzibel
genus(i);             //sollte 0+0-2+1= -1 sein
genus(i,1);           //beides ist korrekt in Singular

*/