source: git/Singular/LIB/normal.lib @ 9ebd82

//////////////////////////////////////////////////////////////////////////////
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;
   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,1); I a 1-dimensional ideal
RETURN:  an integer, the geometric genus p_g = p_a - delta of the projective
         curve defined by i, where p_a is the arithmetic genus.
NOTE:    delta is the sum of all local delta-invariants of the singularities,
         i.e. dim(R'/R), R' the normalization of the local ring R of the
         singularity. @*
         genus(i,1) uses the normalization to compute delta. Usually genus(i,1)
         is slower than genus(i) but sometimes not.
EXAMPLE: example genus; shows an example
"
{
   int w = printlevel-voice+2;  // w=printlevel (default: w=0)

   def R0=basering;
   if(char(basering)>0)
   {
     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);
     }
     int pa=1-hilbPoly(J)[1];
     pa=pa-deltaP(J);
     setring R0;
     return(pa);
   }
   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);

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

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

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

   iglob=vdim(std(IDsing));

   if(iglob!=0)//computation of the radical of IDsing
   {
      ideal radIDsing=reduce(IDsing,radID);
      if(size(radIDsing)==0)
      {
         radIDsing=radID;
         attrib(radIDsing,"isSB",1);
      }
      else
      {
         radIDsing=std(radical(IDsing));
      }
      iglob=vdim(radIDsing);
      if((w>=1)&&(iglob))
          {"the non-nodal-cuspidal locus:";radIDsing;pause();"";}
   }
   cusps=vdim(radID)-iglob;
   nsing=nsing+cusps;

   if(iglob==0)
   {
      if(w>=1){"             there are only cusps and nodes";"";}
      tau=vdim(std(I+jacob(I)));
      tauinf=tauinf+tau;
      nodes=tau-2*cusps;
      delt=nodes+cusps;
      nbranch=2*tau-3*cusps;
      nsing=nsing+nodes;
   }
   else
   {
       if(w>=1){"the non-nodal-cuspidal singularities";"";}
       setring C;
       ideal I1=imap(newR,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";"";}
        option(set,save_opt);
        return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
      }

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

      // the following can certainly be made more efficient when replacing
      // 'hnexpansion' (used only for computing number of branches) by
      // successive blowing-up + test if Newton polygon degenerate:
      if(s>2)    //  splitting of f
      {
         if(w>=1){"Newton polygon can be used for splitting";"";}
         intvec v=NP[1][2]-NP[2][2],NP[2][1];
         int de=w_deg(g,v);
         //int st=w_deg(hc,v)+v[1]+v[2];
         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);
      }
      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);
      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(); }
         option(set,op);
         setring R;
         ideal U = imap(Q,U);
         return(U);
      }
      U=K;
   }
}

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

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

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

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

      list gnirlist = ringlist(basering);
      endid = Le[1];
//endphi;"";
      for( ii=1; ii<=n; ii++)
      {
         if( Le[4][ii] == 0 )            //ii=index of substituted var
         {
            endphi = subst(endphi,var(ii),Le[5][ii]);
         }
      }
//endphi;"";
      list g2 = gnirlist[2];             //the varlist
      list g3 = gnirlist[3];             //contains blocks of orderings
      int n3 = size(g3);

   //----------------- first identify module ordering ------------------
      if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
      {
         list gm = g3[n3];              //last blockis module ordering
         g3 = delete(g3,n3);
         int m = 0;
      }
      else
      {
         list gm = g3[1];              //first block is module ordering
         g3 = delete(g3,1);
         int m = 1;
      }
   //---- delete variables which were substituted and weights  --------
      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();"";
    }
    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)[1];
        list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault)[1];
        printlevel = printlevel - 1;
        option(set,save_opt);
        return(list(nor1, nor2));
      }
    }
  }
  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";