source: git/Singular/LIB/normal.lib @ e8f898

//////////////////////////////////////////////////////////////////////////////
version="version normal.lib 4.0.0.0 Jun_2013 "; // $Id$
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


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 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);     writes a poly in integral closure as element of Quot(R/p)
 iMult(L);           intersection multiplicity of the ideals of the list L

 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 the output of normal, normalP, normalC
 getSmallest(J);     computes the polynomial of smallest degree of J
 getOneVar(J, vari); computes a polynomial of J in the variable vari
 changeDenominator(U1, c1, c2, I); computes ideal U2 such that 1/c1*U1=1/c2*U2

SEE ALSO: locnormal_lib;modnormal_lib
";

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 = list of options. @*
         Optional parameters in list choose (can be entered in any order):@*
         Decomposition:@*
         - \"equidim\" -> computes first an equidimensional decomposition of the
         input ideal, and then the normalization of each component (default).@*
         - \"prim\" -> computes first the minimal associated primes of the input
         ideal, and then the normalization of each prime. (When the input ideal
         is not prime and the minimal associated primes are easy to compute,
         this method is usually faster than \"equidim\".)@*
         - \"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 this assumption
         does not hold, the output might be wrong.@*
         - \"noFac\" -> factorization is avoided in the computation of the
         minimal associated primes;
         Other:@*
         - \"useRing\" -> uses the original ring ordering.@*
         If this option is set and if the ring ordering is not global, normal
         will change to a global ordering only for computing radicals and prime
         or equidimensional decompositions.@*
         If this option is not set, normal changes to dp ordering and performs
         all computations with respect to this ordering.@*
         - \"withDelta\" (or \"wd\") -> returns also the delta invariants.@*
         If the optional parameter choose is not given or empty, only
         \"equidim\" but no other option is used.@*
         - list(\"inputJ\", ideal inputJ) -> takes as initial test ideal the
         ideal inputJ. This option is only for use in other procedures. Using
         this option, the result might not be the normalization.@*
         (Option only valid for global algorithm.)@*
         - list(\"inputC\", ideal inputC) -> takes as initial conductor the
         ideal inputC. This option is only for use in other procedures. Using
         this option, the result might not be the normalization.@*
         (Option only valid for global algorithm.)@*
         Options used for computing integral basis (over rings of two
         variables):@*
         - \"var1\" -> uses a polynomial in the first variable as
         universal denominator.@*
         - \"var2\" -> uses a polynomial in the second variable as universal
         denominator.@*
         If the optional parameter choose is not given or empty, only
         \"equidim\" but no other option is 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 the
         \"prim\" option the minimal associated primes of id are computed first
         and hence normal computes the normalization of the radical of id.@*
NOTE:    \"isPrim\" should only be used if id is known to be prime.
RETURN:  a list, say nor, of size 2 (resp. 3 with option \"withDelta\").
@format  Let R denote the basering and id the input ideal.
         * nor[1] is a list of r rings, where r is the number of associated
         primes P_i with option \"prim\" (resp. >= no of equidimenensional
         components P_i with option \"equidim\").@*
         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 of R/P_i in its field of fractions (as affine ring);
         - @code{normap} gives the normalization map from R/id to
           Ri/norid for each i.@*
         - the direct sum of the rings Ri/norid, i=1,..r, is the normalization
           of R/id as affine algebra; @*
         * nor[2] is a list of size r with information on the normalization of
         the i-th component as module over the basering R:@*
         nor[2][i] is an ideal, say U, in R such that the integral closure
         of basering/P_i is generated as module over R by 1/c * U, with c
         the last element U[size(U)] of U.@*
         * nor[3] (if option \"withDelta\" is set) 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 R/P_i resp. R/id is normal).
@end format
THEORY:  We use here a general algorithm described in [G.-M.Greuel, S.Laplagne,
         F.Seelisch: Normalization of Rings (2009)].@*
         The procedure computes the R-module structure, the algebra structure
         and the delta invariant of the normalization of R/id:@*
         The normalization of R/id is the integral closure of R/id in its total
         ring of fractions. It is a finitely generated R-module and nor[2]
         computes R-module generators of it. More precisely: If U:=nor[2][i]
         and c:=U[size(U)], then c is a non-zero divisor and U/c is an R-module
         in the total ring of fractions, the integral closure of R/P_i. Since
         U[size(U)]/c is equal to 1, R/P_i resp. R/id is contained in the
         integral closure.@*
         The normalization is also an affine algebra over the ground field
         and nor[1] presents it as such. For geometric considerations nor[1] is
         relevant since the variety of the ideal norid in Ri is the
         normalization of the variety of the ideal P_i in R.@*
         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
         For A=K[x1,...,xn]/id we call this number also the delta invariant of
         id. nor[3] returns the delta invariants of the components P_i and of
         id.
NOTE:    To use the i-th ring type e.g.: @code{def R=nor[1][i]; 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 together with the useRing-option (qhweight(id); computes
         weights).
KEYWORDS: normalization; integral closure; delta invariant.
SEE ALSO: normalC, normalP.
EXAMPLE: example normal; shows an example
"
{
  intvec opt = option(get);     // Save current options

  int i,j;
  int decomp;   // Preliminary 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);
  int denomOption;   // Method for choosing the conductor

  ideal inputJ = 0;      // Test ideal given in the input (if any).
  ideal inputC = 0;      // Conductor ideal given in the input (if any).

  list result, resultNew;
  list keepresult;
  list ringStruc;
  ideal U;
  poly c;
  int sp;            // Number of components.

  // Default methods:
  noFac = 0;         // Use facSTD when computing minimal associated primes
  decomp = 2;        // Equidimensional decomposition
  useRing = 0;       // Change first to dp ordering, and perform all
                     // computations there.
  withDelta = 0;     // Do not compute the delta invariant.
  denomOption = 0;   // The default universal denominator is the smallest
                     // degree polynomial.

//--------------------------- define the method ---------------------------
  for ( i=1; i <= size(#); i++ )
  {
    if ( typeof(#[i]) == "string" )
    {
//--------------------------- choosen methods -----------------------
      if ( (#[i]=="isprim") or (#[i]=="isPrim") )
      {decomp = 0;}

      if ( (#[i]=="nodeco") or (#[i]=="noDeco") )
      {decomp = 1;}

      if (#[i]=="prim")
      {decomp = 3;}

      if (#[i]=="equidim")
      {decomp = 2;}

      if ( (#[i]=="nofac") or (#[i]=="noFac") )
      {noFac=1;}

      if ( ((#[i]=="useRing") or (#[i]=="usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
      {useRing = 1;}

      if ( (#[i]=="withDelta") or (#[i]=="wd") or (#[i]=="withdelta"))
      {
        if((decomp == 0) or (decomp == 3))
        {
          withDelta = 1;
        }
        else
        {
          decomp = 3;
          withDelta = 1;
          //Note: the delta invariants cannot be computed with an equidimensional
          //decomposition, hence we compute first the minimal primes
        }
      }
      if (#[i]=="var1")
      {denomOption = 1;}
      if (#[i]=="var2")
      {denomOption = 2;}
    }
    if(typeof(#[i]) == "list"){
      if(size(#[i]) == 2){
        if (#[i][1]=="inputJ"){
          if(typeof(#[i][2]) == "ideal"){
            inputJ = #[i][2];
          }
        }
      }
      if (#[i][1]=="inputC"){
        if(size(#[i]) == 2){
          if(typeof(#[i][2]) == "ideal"){
            inputC = #[i][2];
          }
        }
      }
    }
  }
  kill #;

//------------------------ 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 = attrib(basering,"global");// 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);
  }

//------------------------ trivial checkings ------------------------
  id = groebner(id);
  if((size(id) == 0) or (id[1] == 1))
  {
    // 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 origR;
    list lR = ringlist(origR);
    def ROut = ring(lR);
    setring ROut;
    ideal norid = fetch(origR, id);
    ideal normap = maxideal(1);
    export norid;
    export normap;
    setring origR;
    if(withDelta)
    {
      result = list(list(ROut), list(ideal(1)), list(intvec(0), 0));
    }
    else
    {
      result = list(list(ROut), list(ideal(1)));
    }
    sp = 1;      // number of rings in the output
    option(set, opt);
    normalOutputText(dbg, withDelta, sp);
    return(result);
  }
//------------------------ 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); }
  }
  sp = size(prim);
  if(dbg>=1)
  {
    prim; "";
    "// number of components is", sp;
    "";
  }


//----------------- back to the original ring if required ------------------
// 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;
    ideal inputJ = fetch(origR, inputJ);
    ideal inputC = fetch(origR, inputC);
    if(useRing == 0)
    {
      ideal U;
      poly c;
    }
  }

// ---------------- 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, denomOption, inputJ, inputC);
      printlevel = printlevel - 1;
      for(j = 1; j <= size(norComp); j++)
      {
        newR = norComp[j][3];
        if(!defined(savebasering)) { def savebasering;}
        savebasering=basering;
        setring newR; // must be in a compatible ring to newR
                      // as ringlist may produce ring-dep. stuff
        if(!defined(newRList)) { list newRList;}
        newRList = ringlist(newR);
        setring savebasering;
        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;
          }
        }
        // -- incorporate result for this component to the list of results ---
        if(useRing == 0)
        {
          // We go back to the original ring.
          setring origR;
          U = fetch(R, U);
          c = fetch(R, c);
          newRListO = imap(newR, 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;

  list RL;      // List of rings
  list MG;      // Module generators
  intvec DV;    // Vector of delta's of each component
  for(i = 1; i <= size(result); i++)
  {
    RL[i] = result[i][3];
    MG[i] = lineUpLast(result[i][1], result[i][2]);
    if(withDelta)
    {
      DV[i] = result[i][4];
    }
  }
  if(withDelta)
  {
    resultNew = list(RL, MG, list(DV, deltI));
  }
  else
  {
    resultNew = list(RL, MG);
  }
  sp = size(RL);              //RL = list of rings

  option(set, opt);
  normalOutputText(dbg, withDelta, sp);
  return(resultNew);
}

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", "prim");
  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;
  list nor = normal(i, "withDelta", "prim");  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;
}

///////////////////////////////////////////////////////////////////////////////
// Prints the output text in proc normal.
//
static proc normalOutputText(int dbg, int withDelta, int sp)
// int dbg: printlevel
// int withDelta: output contains information about the delta invariant
// int sp: number of output rings.
{
  if ( dbg >= 0 )
  {
    "";
    if(!withDelta)
    {
      "// 'normal' created a list, say nor, of two elements.";
    }
    else
    {
      "// 'normal' created a list, say nor, of three elements.";
    }
    "// To see the list type";
    "      nor;";
    "";
    "// * nor[1] is a list of", sp, "ring(s).";
    "// To access the i-th ring nor[1][i], give it a name, say Ri, and type";
    "     def R1 = nor[1][1]; setring R1; norid; normap;";
    "// For the other rings type first (if R is the name of your base ring)";
    "     setring R;";
    "// and then continue as for R1.";
    "// Ri/norid is the affine algebra of the normalization of R/P_i where";
    "// P_i is the i-th component of a decomposition of the input ideal id";
    "// and normap the normalization map from R to Ri/norid.";
    "";
    "// * nor[2] is a list of", sp, "ideal(s). Let ci be the last generator";
    "// of the ideal nor[2][i]. Then the integral closure of R/P_i is";
    "// generated as R-submodule of the total ring of fractions by";
    "// 1/ci * nor[2][i].";

    if(withDelta)
    { "";
      "// * nor[3] is a list of an intvec of size", sp, "the delta invariants ";
      "// of the components, and an integer, the total delta invariant ";
      "// of R/id (-1 means infinite, and 0 that R/P_i resp. R/id is normal).";
    }
  }
}


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

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];
   int noRed = 0;
   if(size(Li) > 4)
   {
     if(Li[5] == 1) { noRed = 1; }
   }

   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 zugefuegt
   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];              //### geaendert 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)));


  if(noRed == 0)
  {
    L = substpart(endid,endphi,homo,rw);
    def lastRing=L[1];
    setring lastRing;
    //return(lastRing);
  }
  else
  {
    list RL = ringlist(newRing);
    def lastRing = ring(RL);
    setring lastRing;
    ideal endid = fetch(newRing, endid);
    ideal endphi = fetch(newRing, endphi);
    export(endid);
    export(endphi);
    //def lastRing = newRing;
    //setring R;
    //return(newR);
  }


//   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;
   setring(P);
   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 if and only if 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);
}
///////////////////////////////////////////////////////////////////////////////
static proc deltaP(ideal I)
{
   def R=basering;
   int c,d,i;
   int n=nvars(R);
   list nor;
   if(size(I)>1){ERROR("no hypersurface");}
   ideal J=std(slocus(I));
   if(dim(J)<=0){return(0);}
   poly h;
   d=1;
   while((d)&&(i<n))
   {
      i++;
      h=var(i);
      d=dim(std(J+ideal(h)));
   }
   i=0;
   while(d)
   {
      i++;
      if(i>10){ERROR("delta not found, please inform the authors")};
      h=randomLast(100)[n];
      d=dim(std(J+ideal(h)));
   }
   I=I,h-1;
   if(char(R)<=19)
   {
      nor=normalP(I);
   }
   else
   {
      nor=normal(I);
   }
   return(nor[2][2]);
}

proc genus(ideal I,list #)
"USAGE:   genus(I) or genus(I,<option>); I a 1-dimensional ideal over a perfect field
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:    genus always treats projective curves and takes projective closure if input is affine 1-dim variety.
         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,"nor") uses the normalization to compute delta. Usually genus(I,"nor")
         is slower than genus(I) but sometimes not. @*
         genus(I,"pri") starts with a primary decompsition.
EXAMPLE: example genus; shows an example
"
{
   int w = printlevel-voice+2;  // w=printlevel (default: w=0)

   int ono,rpr,ll;
   if(size(#)>0)
   {
     if(typeof(#[1])=="string")
     {
        if(#[1]=="nor"){ono=1;}
        if(#[1]=="pri"){rpr=1;}
     }
     else { ERROR("invalid option for genus");}
   }
   def R0=basering;
   if((char(basering)>0)||(ono))
   {
     def R1=changeord(list(list("dp",1:nvars(basering))));
     setring R1;
     ideal I=imap(R0,I);
     I=radical(I);
     I=std(I);
     if(dim(I)!=1)
     {
       if(((homog(I))&&(dim(I)!=2))||(!homog(I)))
       {
         ERROR("This is not a curve");
       }
     }
     if(homog(I)&&(dim(I)==2))
     {
       def S=R0;
       setring S;
       ideal J=I;
     }
     else
     {
       def S=changevar(varstr(R0)+",@t");
       setring S;
       ideal J=imap(R1,I);
       J=homog(J,@t);
       J=std(J);
     }

     list LL=normal(J,"prim");
     int pa,i;
     for(i=1;i<=size(LL[1]);i++)
     {
        def T=LL[1][i];
        setring T;
        pa=pa-hilbPoly(std(norid))[1];
        setring S;
        kill T;
     }
     pa=pa+1;
     setring R0;
     return(pa);
   }
   if(!rpr)
   {
      list LZ=minAssGTZ(I);
      if(size(LZ)>1)
      {
         int p_g;
         for(ll=1;ll<=size(LZ);ll++)
         {
            p_g=p_g+genus(LZ[ll],"pri")-1;
         }
         return(p_g+1);
      }
      else
      {
         I=LZ[1];
      }
   }
   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) div 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);

   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,radIDsing);
       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;
          if (inv[2]!=0) { nsing++; }
          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(radIDsing);
           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);
   ring r7=7,(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
"
{
   intvec save_opt=option(get);
   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;
      option(set,save_opt);
      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) div 2;
         option(set,save_opt);
         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) div 2;
         option(set,save_opt);
         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";"";}
        setring R;
        option(set,save_opt);
        return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
      }

      if(w>=1){"Newton-Polygon is degenerated";"";}
/* need to re-consider the degree bound (de):
      // 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];
         int st=w_deg(var(1)^NP[size(NP)][1],v)+1;
         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;
         option(set,save_opt);
         return(list(d*(mu+nb-1) div 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;
         option(set,save_opt);
         return(list(d*(mu+nb-1) div 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);
      }
      setring R;
      option(set,save_opt);
      return(list(d*(mu+nb-1) div 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);
      setring R;
      option(set,save_opt);
      return(list(d*delt,d*tau,d));
   }
   option(set,save_opt);
}
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;
// this code is an exact copy of the proc in paraplanecurves.lib
// (since we do not want to make it non-static)
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 below 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)+");");
      // Rlist may be not defined in this new ring, so we define it again.
      list Rlist2 = ringlist(R);
      def R0 = ring(Rlist2);
      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 singulaeren Orts geaendert (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;
  // In this case we do not eliminate variables, so that the maps
  // are well defined.
  list RR = SM[1],SM[2],J,nzd,1;

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

  list RS = HomJJ(RR);
  //NOTE: HomJJ creates new ring with variables X(i) and T(j)
//-------------------------------------------------------------------------
// 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 polynomial 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,n2,n3;
  intvec V;
  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));
             n2 = size(gnirlist[2]);
             n3 = size(gnirlist[3]);
             for( i=1; i<=ncols(phi); i++)
             {
               gnirlist[2][n2+i] = "Z("+string(i)+")";
             }
             V=0;
             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,n2,n3;
  intvec V;
  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));
         n2 = size(gnirlist[2]);
         n3 = size(gnirlist[3]);
         for( i=1; i<=ncols(phi); i++)
         {
            gnirlist[2][n2+i] = "Z("+string(i)+")";
         }
         V=0;
         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 = radical ideal, choose = optional list of
         strings.
         Optional parameters in list choose (can be entered in any order):@*
         \"withRing\", \"isPrim\", \"noFac\", \"noRed\", where@*
         - \"noFac\" -> factorization is avoided during the computation
         of the minimal associated primes.@*
         - \"isPrim\" -> assumes that the ideal is prime. If the assumption
         does not hold, output might be wrong.@*
         - \"withRing\" -> the ring structure of the normalization is
         computed. The number of variables in the new ring is reduced as much
         as possible.@*
         - \"noRed\" -> when computing the ring structure, no reduction on the
         number of variables is done, it creates one new variable for every
         new module generator of the integral closure in the quotient field.@*
ASSUME:  The characteristic of the ground field must be positive. If the
         option \"isPrim\" is not set, the minimal associated primes of id
         are computed first and hence normalP computes the normalization of
         the radical of id. If option \"isPrim\" is set, the ideal must be
         a prime ideal otherwise the result may be wrong.
RETURN:  a list, say 'nor' of size 2 (resp. 3 if \"withRing\" is set).@*
         ** If option \"withRing\" is not set: @*
         Only the module structure is computed: @*
         * nor[1] is a list of ideals Ii, i=1..r, in the basering R where r
         is the number of minimal associated prime ideals P_i of the input
         ideal id, describing the module structure:@*
         If Ii is given by polynomials g_1,...,g_k in R, then c:=g_k is
         non-zero in the ring R/P_i and g_1/c,...,g_k/c generate the integral
         closure of R/P_i as R-module in the quotient field of R/P_i.@*
         * nor[2] shows the delta invariants: it is a list of an intvec
         of size r, the delta invariants of the r components, and an integer,
         the total delta invariant of R/id
         (-1 means infinite, and 0 that R/P_i resp. R/id is normal). @*
         ** If option \"withRing\" is set: @*
         The ring structure is also computed, and in this case:@*
         * 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 R/P_i, i.e. isomorphic as
           K-algebra (K the ground field) to the integral closure of R/P_i in
           the field of fractions of R/P_i; @*
         - the direct sum of the rings Ri/norid is the normalization
           of R/id; @*
         - @code{normap} gives the normalization map from R to Ri/norid.@*
         * nor[2] gives the module generators of the normalization of R/P_i,
         it is the same as nor[1] if \"withRing\" is not set.@*
         * nor[3] shows the delta invariants, it is the same as nor[2] if
         \"withRing\" is not set.
THEORY:  normalP uses the Leonard-Pellikaan-Singh-Swanson algorithm (using the
         Frobenius) cf. [A. K. Singh, I. Swanson: An algorithm for computing
         the integral closure, arXiv:0901.0871].
         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
         For A=K[x1,...,xn]/id we call this number also the delta invariant of
         id. The procedure returns the delta invariants of the components P_i
         and of id.
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.
         For local or mixed orderings use proc 'normal'.
@*       If the input ideal id is weighted homogeneous a weighted ordering may
         be used (qhweight(id); computes weights).
@*       Works only in characteristic p > 0; use proc normal in char 0.
KEYWORDS: normalization; integral closure; delta invariant.
SEE ALSO: normal, normalC
EXAMPLE: example normalP; shows an example
"
{
   int i,j,y, sr, del, co;
   intvec deli;
   list resu, Resu, prim, Gens, mstdid;
   ideal gens;

   // Default options
   int wring = 0;           // The ring structure is not computed.
   int noRed = 0;           // No reduction is done in the ring structure
   int isPrim = 0;          // Ideal is not assumed to be prime
   int noFac = 0;           // Use facstd when computing the decomposition


   y = printlevel-voice+2;

   if ( attrib(basering,"global") != 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,"noRed") or find(method,"nored") )
   {
     noRed=1;
   }
   if ( find(method,"isPrim") or find(method,"isprim") )
   {
     isPrim=1;
   }
   if ( find(method,"noFac") or 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 -----------
           if(noRed == 0){
             resu = list(computeRing(II,prim[i])) + resu;
           }
           else
           {
             resu = list(computeRing(II,prim[i], "noRed")) + 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:
// To see the result, type
     nor;

// * nor[1] is a list of",sr,"ring(s):
// 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 R1;
// Ri/norid is the affine algebra of the normalization of the i-th
// component R/P_i (where P_i is a min. 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-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:
// To see the result, type
     nor;

// * nor[1] is a list of",sr,"ideal(s), where each ideal nor[1][i] consists
// of elements g1..gk of the basering R such that gj/gk generate the integral
// closure of R/P_i (where P_i is a min. associated prime of the input ideal)
// as R-module in the quotient field of R/P_i;

// * 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 mstd(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(); }
         setring R;
         option(set,op);
         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  --------
      intvec V;
      int n(0);
      list newg2;
      list newg3;
      for ( ii=1; ii<=n3-1; ii++ )
      {
        // If order is a matrix ordering, it is replaced by dp ordering.
        // TODO: replace it only when some of the original
        //       variables are eliminated.
        if(g3[ii][1] == "M"){
          g3[ii][1] = "dp";
          g3[ii][2] = (1..sqroot(size(g3[ii][2])))*0+1;
        }
        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++)
        {
          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
          }
        }
        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;  */

    list gnirlist = ringlist(R);
    def R1 = ring(gnirlist);
    setring R1;
    ideal norid = imap(R, I);
    ideal normap = imap(R, maxid);
    export norid;
    export normap;

    if(noRed == 1){
      setring R;
      return(R1);
    }
    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 geaendert: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(attrib(basering,"global") != 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();"";
    }
    setring R;
    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

// inputJ is used in parametrization of rational curves algorithms, to specify
// a different test ideal.

static proc normalM(ideal I, int decomp, int withDelta, int denomOption, ideal inputJ, ideal inputC){
// 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 universal denominator.
// 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).

// denomOption = 0      -> Uses the smallest degree polynomial
// denomOption = i > 0  -> Uses a polynomial in the i-th variable

  intvec save_opt=option(get);
  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 = attrib(basering,"global");

  poly c;                     // The common denominator.

  def R = basering;

//------------------------ Groebner bases and dimension of I-----------------
  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);

  // ---------------- 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.
  qring Q = I;   // We work in the quotient by the groebner base of the ideal I
  option(redSB);
  option(returnSB);

  // If a conductor ideal was given as input, we use it instead of the
  // singular locus. If a test ideal was given as input, we do not compute the
  // singular locus.
  ideal inputC = fetch(R, inputC);
  ideal inputJ = fetch(R, inputJ);
  if((inputC == 0) && (inputJ == 0))
  {
    // We compute the radical of the ideal of minors modulo the original ideal.
    // This is done only in the first step.
    ideal I = fetch(R, I);
    attrib(I, "isSB", 1);
    ideal IMin = fetch(R, IMin);

    dbprint(dbg, "Computing the jacobian ideal...");

    // If a given conductor ideal is given, we use it.
    // If a given test ideal is given, we don't need to compute the jacobian

    // reduction mod I in 'minor' is not working for local orderings!
    if(attrib(basering,"global"))
    {
      ideal J = minor(jacob(IMin), nvars(basering) - d, I);
    }
    else
    {
      ideal J = minor(jacob(IMin), nvars(basering) - d);
      J = reduce(J, groebner(I));
    }
    J = groebner(J);
  }
  else
  {
    ideal J = fetch(R, inputC);
    J = groebner(J);
  }

  //------------------ We check if the singular locus is empty -------------
  if(J[1] == 1)
  {
    // 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 = fetch(R, I);
    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;
    }
    option(set,save_opt);
    return(list(output));
  }


  // -------------------- election of the universal denominator----------------
  // We first check if a conductor ideal was computed. If not, we don't
  // compute a universal denominator.
  ideal Id1;
  if(J != 0)
  {
    if(denomOption == 0)
    {
      poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
                                     // of J as universal denominator.
    }
    else
    {
      poly condu = getOneVar(J, denomOption);
    }
    if(dbg >= 1)
    {
      "";
      "The universal denominator is ", condu;
    }

    // ----- splitting the ideal by the universal denominator (if possible) -----
    // If the ideal is equidimensional, but not necessarily prime, we check if
    // the universal denominator is a non-zerodivisor of R/I.
    // If not, we split I.
    if((decomp == 1) or (decomp == 2))
    {
      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;
        ideal JDefault = 0; // Now it uses the default J;
        list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault);
        list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault);
        printlevel = printlevel - 1;
        option(set,save_opt);
        list res = nor1 + nor2;
        return(res);
      }
    }
  }
  else
  {
    poly condu = 0;
  }

  // --------------- 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.
  setring R;
  // If a test ideal is given at the input, we use it.
  if(inputJ == 0)
  {
    if(isGlobal == 1)
    {
      ideal J = fetch(Q, J);
      J = J, I;
      if(dbg >= 1)
      {
        "The original singular locus is";
        groebner(J);
        if(dbg >= 2){pause();}
        "";
      }
      // We check if the only singular point is the origin.
      // If so, the radical is the maximal ideal at the origin.
      J = groebner(J);
      if(locAtZero(J))
      {
        J = maxideal(1);
      }
      else
      {
        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);
    }
  }
  else
  {
    ideal J = inputJ;
  }

  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);
  if(denomOption == 0)
  {
    poly D = getSmallest(J);    // Chooses the polynomial of smallest degree as
                                // non-zerodivisor.
  }
  else
  {
    poly D = getOneVar(J, denomOption);
  }
  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 Id2 = quotient(I, Id1);
      // I = Id1 \cap Id2
      printlevel = printlevel + 1;

      ideal JDefault = 0;  // Now it uses the default J;
      list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault);
      list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault);
      printlevel = printlevel - 1;
      option(set,save_opt);
      list res = nor1 + nor2;
      return(res);
    }
  }

  // --------------------- 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, denomOption);
  printlevel = printlevel - 1;
  option(set,save_opt);
  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 or 0 if it was not computed.
// 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 = attrib(basering,"global");
  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);
  intvec save_opt=option(get);
  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 = fetch(R, I);
    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;
    }
    option(set,save_opt);
    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 we have a universal denominator of smaller degree than c,
      // we replace c by it.
      if(condu != 0)
      {
        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;
  }
  option(set,save_opt);
  return(output);
}

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

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

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

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

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

//WARNING - elim is not working here!! Check!!
//It is now replaced by computing an eliminating groebner basis.
proc getOneVar(ideal J, int vari)
"USAGE:   getOneVar(J, vari); J is a 0-dimensional ideal, vari is an integer.
RETURN:  a polynomial of J in the variable indicated by vari of smallest
         degree.@*
NOTE:    Works only over rings of two variables.@*
         It is intended mainly as an auxiliary procedure for computing
         integral bases. @*
EXAMPLE: example getOneVar; shows an example
"
{
  def R = basering;
  list RL = ringlist(R);
  // We keep everything from R but we change the ordering to lp, and we
  // order the variables as needed.
  RL[3] = list(list("lp", 1:2), list("C", 0:1));
  RL[2] = list(var(3-vari), var(vari));
  RL[4]=ideal(0); // does not work with qrings: Ex.7 of paraplanecurves
  def RR = ring(RL);
  setring RR;
  ideal J = imap(R, J);
  J = groebner(J);
  poly g = J[1];
  setring R;
  poly g = imap(RR, g);
  return(g);
}
example
{ "EXAMPLE:";
  printlevel = printlevel+1;
  echo = 2;
  ring s = 0,(x,y),dp;
  ideal J = x3-y, y3;
  getOneVar(J, 1);

  echo = 0;
  printlevel = printlevel-1;
}
///////////////////////////////////////////////////////////////////////////////

proc getSmallest(ideal J)
"USAGE:   getSmallest(J); J is an ideal.
RETURN:  the generator of J of smallest degree. If there are more than one, it
         chooses the one with smallest number of monomials.@*
NOTE:    It looks only at the generator of J, not at all the polynomials in
         the ideal.@*
         It is intended maninly to compute a good universal denominator in the
         normalization algorithms.@*
EXAMPLE: example getSmallest; shows an example
"
{
// Computes the polynomial of smallest degree of J.
//
  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);
}
example
{ "EXAMPLE:";
  printlevel = printlevel+1;
  echo = 2;
  ring s = 0,(x,y),dp;
  ideal J = x3-y, y5, x2-y2+1;
  getSmallest(J);

  echo = 0;
  printlevel = printlevel-1;
}

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

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.
  intvec save_opt=option(get);
  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 = attrib(origEre,"global");      // 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);

  option(set,save_opt);
  return(list(J, ele[1]));
}

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

proc changeDenominator(ideal U1, poly c1, poly c2, ideal I)
"USAGE:   changeDenominator(U1, c1, c2, I); U1 and I ideals, c1 and c2
         polynomials.@*
RETURN:  an ideal U2 such that the A-modules 1/c1 * U1 and 1/c2 * U2 are equal,
         where A = R/I and R is the basering.@*
NOTE:    It assumes that such U2 exists. It is intended maninly as an auxiliary
         procedure in the normalization algorithms.@*
EXAMPLE: example changeDenominator; shows an example
"
{
// Let A = R / I. Given an A-module in the form 1/c1 * U1 (U1 ideal of A), it
// computes a new ideal U2 such that the the A-module 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);
}
example
{
  "EXAMPLE:";
  echo = 2;
  ring s = 0,(x,y),dp;
  ideal I = y5-y4x+4y2x2-x4;
  ideal U1 = normal(I)[2][1];
  poly c1 = U1[4];
  U1;c1;
  // 1/c1 * U1 is the normalization of I.
  ideal U2 = changeDenominator(U1, c1, x3, I);
  U2;
  // 1/x3 * U2 is also the normalization of I, but with a different denominator.
  echo = 0;
}

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

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 w.r.t. 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;
  int size_p=size(p);
  for(i = 1; i <= size_p; i++)
  {
    e = sum(leadexp(p[i]), @v);
    g = g + p[i] * c^(d-e);
  }
  return(g);
}

//                    End procedures for normal
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
//                  Begin 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 list
         of string.
         Optional parameters in list choose (can be entered in any order):@*
         Decomposition:@*
         - \"equidim\" -> computes first an equidimensional decomposition,
         and then the normalization of each component (default).@*
         - \"prim\" -> computes first the minimal associated primes, and then
         the normalization of each prime. @*
         - \"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 the assumption does
         not hold, output might be wrong.@*
         - \"noFac\" -> factorization is avoided in the computation of the
         minimal associated primes;
         Other:@*
         - \"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;@*
         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, if option
         \"prim\" is set the minimal associated primes 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 option \"withGens\" is set).@*
         * nor[1] is always a of r rings, where r is the number of associated
         primes with option \"prim\" (resp. >= no of equidimenensional
         components with option  \"equidim\").@*
         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, i.e. Ri is
          given in the form K[X(1..p),T(1..q)], where K is the ground field;
         - normap gives the normalization map from basering/id to
           Ri/norid for each i (the j-th element of normap is mapped to the
           j-th variable of R).@*
         - the direct sum of the rings Ri/norid is the normalization
           of basering/id; @*
         ** If option \"withGens\" is not set: @*
         * nor[2] shows the delta invariants: nor[2] is a list of an intvec of
         size r, the delta invariants of the r components, and an integer, the
         delta invariant of basering/id. (-1 means infinite, 0 that basering/P_i
         resp. basering/input is normal, -2 means that delta resp. delta of one
         of the components is not computed (which may happen if \"equidim\" is
         given). @*
         ** If option \"withGens\" is set:
         * nor[2] is a list of ideals Ii=nor[2][i], i=1..r, in the basering,
         generating the integral closure of basering/P_i in its quotient field
         as K-algebra (K the ground field):@*
         If Ii is given by polynomials g_1,...,g_k, then c:=g_k is a non-zero
         divisor and the j-th variables of the ring Ri satisfies var(j)=g_j/c,
         j=1..k-1, as element in the quotient field of basering/P_i. The
         g_j/g_k+1 are K-algebra generators  of the integral closure of
         basering/P_i.@*
         * nor[3] shows the delta invariant as above.
THEORY:  We use the Grauert-Remmert-de Jong algorithm [c.f. G.-M. Greuel,
         G. Pfister: A SINGULAR Introduction to Commutative Algebra, 2nd Edition.
         Springer Verlag (2007)].
         The procedure computes the algebra structure and the delta invariant of
         the normalization of R/id:@*
         The normalization is an affine algebra over the ground field K
         and nor[1] presents it as such: Ri = K[X(1..p),T(1..q)] and Ri/norid
         is the integral closure of R/P_i; if option \"withGens\" is set the
         X(j) and T(j) are expressed as quotients in the total ring of
         fractions. Note that the X(j) and T(j) generate the integral closure
         as K-algebra, but not necessarily as R-module (since relations of the
         form X(1)=T(1)*T(2) may have been eliminated). Geometrically the
         algebra structure is relevant since the variety of the ideal norid in
         Ri is the normalization of the variety of the ideal P_i in R.@*
         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
         For A=K[x1,...,xn]/id we call this number also the delta invariant of
         id. nor[3] returns the delta invariants of the components P_i and of
         id.
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.
         For local or mixed orderings use proc 'normal'.
@*       If the input ideal id is weighted homogeneous a weighted ordering may
         be used (qhweight(id); computes weights).
KEYWORDS: normalization; integral closure; delta invariant.
SEE ALSO: normal, normalP.
EXAMPLE: example normalC; shows an example
"
{
   int i,j;
   int withGens, withEqui, withPrim, isPrim, noFac;
   int dbg = printlevel-voice+2;
   int nvar = nvars(basering);
   int chara  = char(basering);
   list result, prim, keepresult;

  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

   // Default methods:
   noFac = 0;         // Use facstd when computing minimal associated primes
   decomp = 2;        // Equidimensional decomposition for nvar > 2
   if (nvar <= 2)
   { decomp = 3; }    // Compute minimal associated primes if nvar <= 2

   if ( attrib(basering,"global") != 1 )
   {
     "";
     "// Not implemented for this ordering,";
     "// please change to global ordering or use proc normal";
     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 -----------------------
   // "withGens": computes algebra generators for each irreducible component
   // ### the extra code for withGens should be incorporated in the general case

   if ( find(method,"withgens") or find(method,"withGens"))
   {
     withGens = 1;
   }

   // the general case: either equidim or minAssGTZ or no decomposition

   if ( find(method,"isprim") or find(method,"isPrim") )
   {decomp = 0; isPrim=1;}

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

   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 the static proc closureGenerators

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

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

         if(dbg >= 2)
         {  prim;""; }
         if(dbg >= 1)
         {
            "// 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(dbg>=1)
         {
            ""; pause(); "";
            "// Computing normalization 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
         // primeClosure is called with list l consisting of the basering
         //### ausprobieren ob primeClosure(l,1) schneller als primeClosure(l)
         // 1 bedeutet: kuerzester nzd
         // l[size(l)] is the delta invariant

         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 zulaessig!
         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 ( dbg >= 0 )
      {"";
"// 'normalC' created a list, say nor, of three lists:
// To see the list type
      nor;

// * nor[1] is a list of",sr,"ring(s)
// 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 (if R is the name of your original basering)
     setring R;
// and then continue as for R1.
// 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 id)
// 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 id
// (-1 means infinite, and 0 that R/P_i resp. R/id is normal).";
      }
      return(Resu);
   }
   //----------------- end method "withGens" --------------------------------

//-------- 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

   if( decomp == 2 )
   {
      withPrim = 0;                 //this is used to check later that prim
                                    //contains equidim but not prime components
      if( dbg >= 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(dbg>=1)
      {  "";
         "// number of equidimensional components:", size(prim);
      }
      if ( !noFac )
      {
        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(dbg>=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( dbg >= 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( dbg >= 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( dbg >= 1 )
      {
         "// We use method 'prim'";
      }

      if( typeof(attrib(id,"isPrim"))=="int" )
      {
         if(attrib(id,"isPrim")==1)
         {
            prim[1]=id;
         }
         else
         {
            if( noFac )
            { prim=minAssGTZ(id,1); }     //does not use factorizing groebner
            else
            { prim=minAssGTZ(id); }       //uses factorizing groebner
         }
      }
      else
      {
            if( noFac )
            { prim=minAssGTZ(id,1); }
            else
            { prim=minAssGTZ(id); }
      }
      if(dbg>=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(dbg>=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 ( dbg>=1 )
      {
           "// delta of component",i; keepresult[skr][1];
      }
   }
   sr = size(result);

   // -------------- Now compute intersection multiplicities -------------
   //intersection multiplicities of list prim, sp=size(prim).
      if ( dbg>=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 ( dbg >= 0)
   {"";
"// 'normalC' created a list, say nor, of two lists:
// To see the result, type
      nor;

// * nor[1] is a list of",sr,"ring(s).
// 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 r/P_i  (where P_i
// is an associated prime or an equidimensional part of the input ideal id)
// and normap the normalization map from the basering to Ri/norid;

// * nor[2] shows the delta-invariant of each component and of id
// (-1 means infinite, 0 that r/P_i resp. r/id is normal, and -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 #)
{
  def r = basering;

  //--------------------------- 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.";
    norTest(I, a1);
    "";
  }

  if(checkP and norP and norM)
  {
    "Comparing with normalP output...";
    if(size(a2) > 0)
    {
      "WARNING: this checking only works if the original ideal was prime.";
      U1 = a1[2][1];
      c1 = U1[size(U1)];
      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;
    }
    else
    {
      "normalP returned no output. Comparison is not possible.";
    }
    "";
  }

  if(checkPC and norM and primCl)
  {
    "Comparing with primeClosure output...";
    if(size(a4) > 0)
    {
      "WARNING: this checking only works if the original ideal was prime.";
      // primeClosure check
      U1 = a1[2][1];
      c1 = U1[size(U1)];
      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;
    }
    else
    {
      "normalC(withGens) returned no output. Comparison is not possible.";
    }
    "";
  }
}

///////////////////////////////////////////////////////////////////////////
static proc sqroot(int n);
{
  int s = 1;
  while(s*s < n) { s++; }
  return(s);
}

///////////////////////////////////////////////////////////////////////////
proc norTest (ideal i, list nor, list #)
"USAGE:   norTest(i,nor,[n]); i=prime ideal, nor=list, n=optional integer
ASSUME:  nor is the output of normal(i) (any options) or
         normalP(i,"withRing") or normalC(i) (any options).
         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 if and only if 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 override 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);
    }
    def R = newnor[1][1];
    qring QAS = std(i);


    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 )
    {
     intvec result = intvec(a);
     setring BAS;
     return (result);
   }
   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 )
    {
       intvec result = intvec(a,b);
       setring BAS;
       return (result);
    }
   b;

    //------------------------ check normality ---------------------------
    list testnor = normal(lnorid[2],"isPrim","noFac", "withDelta");
    //### 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[3][2];             //d = delta
    kill testnor;                              //### sollte ueberfluessig sein
    int d1 = (d==0);                           //d1=1 if delta=0
    dbprint ( prl, "delta: "+string(d) );
    intvec result = intvec(a,b,d1);
    setring BAS;
    return(result);
}
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

   nor = normalC(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

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

//----------------------Examples for normal (new algorithm)------------------
// Timings with Computeserver Dual AMD Opteron 242 1.60GHz.
// Examples from "Normalization of Rings" paper.

// Example 1
// char 0 : normal = 0 secs (7 steps) - normalC = 75 secs
// char 2 : normal = 0 secs (7 steps) - normalP = 0 secs - normalC = 0 secs
// char 5 : normal = 1 secs (7 steps) - normalP = 71 - normalC = 1 secs
// char 11 : normal = 2 secs (7 steps) - normalP = 12 secs - normalC doesn't finish
// char 32003 : normal = 1 secs (7 steps) - normalP doesn't finish - normalC = 1 sec
LIB"normal.lib";
ring r = 2, (x, y), dp;
ideal i = (x-y)*x*(y+x^2)^3-y^3*(x^3+x*y-y^2);
timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");

// Example 2
// char 0  : normal = 1 sec (7 steps) - normalC doesn't finish
// char 3 : normal = 1 secs (8 steps) - normalP = 0 secs - normalC = 4 secs
// char 13 : normal = 1 sec (7 steps) - normalP doesn't finish - normalC = 13 secs
// char 32003 : normal = 1 secs (7 steps) - normalP doesn't finish - normalC = 10 sec
//Example is reducible in char 5 and 7
LIB"normal.lib";
ring r = 3, (x, y), dp;
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;
timeNormal(i, "normal", "normalC", "normalP", "p", "isPrim");

// Example 3
// char 0 : normal = 3 secs (6 steps) - normalC doesn't finish
// char 2 : normal = 1 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
// char 5 : normal = 0 secs (6 steps) - normalP = 8 secs - normalC doesn't finish
LIB"normal.lib";
ring r=5,(x, y),dp;
ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9;
timeNormal(i, "normal", "normalC", "normalP", "isPrim");

// Example 4
// char 0 : normal = 0 secs (1 step) - normalC = 0 secs
// char 5 : normal = 0 secs (1 step) - normalP = 3 secs - normalC = 0 secs
// char 11 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
// char 32003 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
LIB"normal.lib";
ring r=5,(x,y),dp;   // genus 0 4 nodes and 6 cusps im P2
ideal i=(x2+y^2-1)^3 +27x2y2;
timeNormal(i, "normal", "normalC", "normalP", "isPrim");

// Example 5
// char 0 : normal = 0 secs (1 step) - normalC = 0 secs
// char 5 : normal = 1 secs (3 step) - normalP doesn't finish - normalC doesn't finish
// char 11 : normal = 0 secs (1 step) - normalP 0 secs - normalC = 0 secs
// char 32003 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
LIB"normal.lib";
ring r=11,(x,y),dp;    //24 sing, delta 24
ideal i=-x10+x8y2-x6y4-x2y8+2y10-x8+2x6y2+x4y4-x2y6-y8+2x6-x4y2+x2y4+2x4+2x2y2-y4-x2+y2-1;
timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");

// Example 6
// char 2 : normal = 5 secs (2 steps) - normalP = 25 secs - normalC = 166 secs
LIB"normal.lib";
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;
timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");

// Example 7
// char 0 : normal = 11 secs (6 steps) - normalC = 11 secs
// char 2 : normal = 11 secs (6 steps) - normalP = 0 secs - normalC = 11 secs
// char 5 : normal = 11 secs (6 steps) - normalP = 3 secs - normalC = 11 secs
// char 11 : normal = 11 secs (6 steps) - normalP = 43 secs - normalC = 11 secs
// char 32003 : normal = 11 secs (6 steps) - normalP doesn't finish - normalC = 11 secs
LIB"normal.lib";
ring r=11,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3;
timeNormal(i, "normal", "normalC", "normalP", "isPrim");

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

// Other examples with new algorithm

// Example 1
// char 0 : normal = 1 secs (13 steps) - normalC doesn't finish
// char 2 : normal = 1 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
// char 5 : normal = 1 secs (13 steps) - normalP = 29 secs - normalC doesn't finish
ring r=2,(x,y),dp;  //genus 35
ideal i=y30+y13x+x4y5+x3*(x+1)^2;
timeNormal(i, "normal", "normalC", "normalP");

// Example 2
// char 0 : normal = 1 secs (13 steps) - normalC doesn't finish
// char 3 : normal = 2 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
ring r=3,(x,y),dp;  //genus 19, delta 21
ideal i=y20+y13x+x4y5+x3*(x+1)^2;
timeNormal(i, "normal", "normalC", "normalP");

// Example 3
// Very fast with all algorithms
ring r = 3, (x, y), dp;
ideal I = (x-y^2)^2-x*y^3;
timeNormal(I, "normal", "normalC", "normalP", "primCl", "111", "p", "pc");



//----------------------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
//normal(neu)       0     0    1   2    10  155
//### 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=11; 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! haengt 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
//normal(neu) haengt bei return(std(i)) (offensichtlich in eineranderen lib)

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: ### haengt 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

*/