Changeset a9431c in git

Timestamp:

Mar 26, 2009, 5:37:19 PM (14 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')

Children:

53e33d93f45a825808d84cade536120482dc46c3

Parents:

54248e343a5d77bfd7f19c65855707d3950619e8

Message:

*hannes: new normal


git-svn-id: file:///usr/local/Singular/svn/trunk@11594 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/normal.lib (modified) (63 diffs)

Singular/LIB/normal.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: normal.lib,v 1.49 2006-07-25 17:54:28 Singular Exp $";
+version="$Id: normal.lib,v 1.50 2009-03-26 16:37:19 Singular Exp $";
 category="Commutative Algebra";
 info="
 LIBRARY:  normal.lib     Normalization of Affine Rings
 AUTHORS:  G.-M. Greuel,  greuel@mathematik.uni-kl.de,
+@*        S. Laplagne,   slaplagn@dm.uba.ar,
 @*        G. Pfister,    pfister@mathematik.uni-kl.de
 
+
 MAIN PROCEDURES:
- normal(I[,wd]);     computes the normalization of basering/I,
-                     resp. computes the normalization of basering/I and
-                     the delta invariant
- HomJJ(L);           presentation of End_R(J) as affine ring, L a list
- genus(I);           computes genus of the projective curve defined by I
- primeClosure(L);    integral closure of R/p (p prime), L a list
- closureFrac(L)      write poly in integral closure as element of Quot(R/p)
-
-AUXILIARY PROCEDURE:
- deltaLoc(f,S);      (sum of) delta invariant(s) at conjugated singular points
+ normal(I[...]);     normalization of an affine ring 
+ normalP(I,[...];    normalization of an affine ring in positive characteristic
+ normalC(I,[...];    normalization of an affine ring through a chain of rings
+ HomJJ(L);           presentation of the End_R(J) as affine ring, J an ideal
+ genus(I);           computes the geometric genus of a projective curve
+ primeClosure(L);    integral closure of R/p, p a prime ideal
+ closureFrac(L);     write a poly in integral closure as element of Quot(R/p)
+ iMult(L);           intersection multiplicity of the ideals of the list L
+
+AUXILIARY PROCEDURES:
+ deltaLoc(f,S);      sum of delta invariants at conjugated singular points
+ locAtZero(I);       checks whether the zero set of I is located at 0
+ norTest(i,nor;)     checks whether result of normal or mormal_p was correct
 ";
 
@@ -30 +35 @@
 LIB "hnoether.lib";
 LIB "reesclos.lib";
+LIB "algebra.lib";
+
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc normal(ideal id, list #)
+"USAGE:  normal(id [,choose]); id = radical ideal, choose = optional string. @*
+         Optional parameters in list choose (can be entered in any order):@*
+         Decomposition:@*
+         - \"prim\" -> computes first the minimal associated primes, and then 
+         the normalization of each prime. (default)@*
+         - \"equidim\" -> computes first an equidimensional decomposition,
+         and then the normalization of each component. @*
+         - \"noDeco\" -> no preliminary decomposition is done. If the ideal is 
+         not equidimensional radical, output might be wrong.@*
+         - \"isPrim\" -> assumes that the ideal is prime. If assumption does 
+         not hold, output might be wrong.@*
+         - \"noFac\" -> factorization is avoided in the computation of the 
+         minimal associated primes;
+         Other:@*
+         - \"useRing\" -> uses the original ring ordering.@*
+         If the ring ordering is not global, it will change to a global
+         ordering only for computing radicals and prime or equidimensional
+         decompositions.@*
+         If this option is not set, changes to dp ordering and perform all
+         computation with respect to this ordering.@*
+         - \"withDelta\" (or \"wd\") -> returns also the delta invariants.
+         If choose is not given or empty, the default options are used.@*
+ASSUME:  The ideal must be radical, for non-radical ideals the output may
+         be wrong (id=radical(id); makes id radical). However, when using 
+         \"prim\" option the minimal associated primes of id are computed first 
+         and hence normal computes the normalization of the radical of id.
+         \"isPrim\" should only be used if id is known to be prime.
+RETURN:  a list, say nor, of size 1 (resp. 2 with option \"withDelta\"). 
+@format 
+         The first element nor[1] is a list of r elements, where r is the 
+         number of associated primes P_i with option \"prim\" (resp. >= no 
+         of equidimenensional components P_i with option \"equidim\").@*
+         Each element of nor[1] is a list with the information on the
+         normalization of the i-th component, i.e. the integral closure 
+         in its field of fractions (as affine ring).@*
+         Each list nor[1][i] contains 3 elements (resp. 4 with option
+         \"withDelta\"). @*
+         The first two element of nor[1][i] give the normalization by
+         module generators. The normalization is of the form
+         1/c * U, with U = nor[1][i][1] an ideal of R / P_i and
+         c = nor[1][i][2] a polynomial in R / P_i.@*
+         The third element is a ring Ri=nor[1][i][3], i=1..r, 
+         which contains two ideals with given
+         names @code{norid} and @code{normap} such that @*
+         - Ri/norid is the normalization of the i-th component.
+         - the direct sum of the rings Ri/norid is the normalization
+           of basering/id; @*
+         - @code{normap} gives the normalization map from basering/id to
+           Ri/norid for each j.@*
+         When option \"withDelta\" is set, the forth element nor[1][i][4] is
+         the delta invariant of the i-th component. (-1 means infinite, and 0 
+         that basering/P_i is normal)@*
+         @*
+         Also, When option \"withDelta\" is set the second element of the list
+         nor, nor[2], is an integer, the total delta invariant of basering/id
+         (-1 means infinite, and 0 that basering/id is normal).
+@end format
+THEORY:  We use a general algorithm described in [G.-M. Greuel, S. Laplagne, 
+         F. Seelisch: Normalization of Rings (2009)].
+         The delta invariant of a reduced ring K[x1,...,xn]/id is 
+             dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id)
+         We call this number also the delta invariant of id.
+NOTE:    To use the i-th ring type: @code{def R=nor[1][i][3]; setring R;}.
+@*       Increasing/decreasing printlevel displays more/less comments 
+         (default: printlevel=0).
+@*       Implementation works also for local rings.         
+@*       Not implemented for quotient rings.
+@*       If the input ideal id is weighted homogeneous a weighted ordering may
+         be used (qhweight(id); computes weights).
+KEYWORDS: normalization; integral closure; delta invariant.
+EXAMPLE: example normal; shows an example
+"
+{  
+   intvec opt = option(get);     // Save current options
+
+   int i,j;
+   int decomp;  // Preliminar decomposition:
+                // 0 -> no decomposition (id is assumed to be prime)
+                // 1 -> no decomposition 
+                //      (id is assumed to be equidimensional radical) 
+                // 2 -> equidimensional decomposition
+                // 3 -> minimal associated primes
+   int noFac, useRing, withDelta;
+   int dbg = printlevel - voice + 2;
+   int nvar = nvars(basering);
+   int chara  = char(basering);
+   list result;
+   list keepresult;
+   list ringStruc;
+   ideal U;
+   poly c;
+   
+   // Default methods:
+   noFac = 0;         // Use facSTD when computing minimal associated primes
+   decomp = 3;        // Minimal associated primes
+   useRing = 0;       // Change first to dp ordering, and perform all 
+                      // computations there.
+   withDelta = 0;     // Do not compute the delta invariant.
+  
+
+  //--------------------------- define the method --------------------------- 
+  string method;                //make all options one string in order to use
+                               //all combinations of options simultaneously
+  for ( i=1; i <= size(#); i++ )
+  { 
+    if ( typeof(#[i]) == "string" ) 
+    {
+      method = method + #[i];
+    }
+  }
+
+  //--------------------------- choosen methods -----------------------
+  if ( find(method,"isprim") or find(method,"isPrim") )
+  {decomp = 0;}
+  
+  if ( find(method,"nodeco") or find(method,"noDeco") )
+  {decomp = 1;}
+  
+  if ( find(method,"equidim") )
+  {decomp = 2;}
+  
+  if ( find(method,"prim") )
+  {decomp = 3;}
+  
+  if ( find(method,"nofac") or find(method,"noFac") )
+  {noFac=1;}
+  
+  if ( (find(method,"useRing") or find(method,"usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
+  {useRing = 1;}
+  
+  if ( find(method,"withDelta") or find(method,"wd") or find(method,"withdelta"))
+  {
+    if((decomp == 0) or (decomp == 3)){
+      withDelta = 1;
+    } else {
+      "WARNING: the delta invariants cannot be computed with an equidimensional";
+      "         decomposition.";
+      "         Use option \"prim\" to compute first the minimal associated primes.";
+      "";
+    }
+  }
+
+  kill #;
+  list #;
+
+  //------------------------ change ring if required ------------------------
+  // If the ordering is not global, we change to dp ordering for computing the 
+  // min ass primes.
+  // If the ordering is global, but not dp, and useRing = 0, we also change to 
+  // dp ordering.
+  
+  int isGlobal = ord_test(basering);      // Checks if the original ring has
+                                          // global ordering.
+  
+  def origR = basering;   // origR is the original ring
+                          // R is the ring where computations will be done
+    
+  if((useRing  == 1) and (isGlobal == 1)){
+    def globR = basering;
+  } else {
+    // We change to dp ordering.
+    list rl = ringlist(origR);
+    list origOrd = rl[3];
+    list newOrd = list("dp", intvec(1:nvars(origR))), list("C", 0);
+    rl[3] = newOrd;
+    def globR = ring(rl);
+    setring globR;
+    ideal id = fetch(origR, id);  
+  }
+
+  //------------------------ preliminary decomposition-----------------------
+  list prim;
+  if(decomp == 2){
+    dbprint(dbg, "// Computing the equidimensional decomposition...");
+    prim = equidim(id);
+  }
+  if((decomp == 0) or (decomp == 1)){
+    prim = id;
+  }
+  if(decomp == 3){
+    dbprint(dbg, "// Computing the minimal associated primes...");
+    if( noFac )
+    { prim = minAssGTZ(id,1); }
+    else
+    { prim = minAssGTZ(id); }
+  }
+
+  // if ring was not global and useRing is on, we go back to the original ring
+  if((useRing == 1) and (isGlobal != 1)){
+    setring origR;
+    def R = basering;
+    list prim = fetch(globR, prim);
+  } else {
+    def R = basering;
+    if(useRing == 0){
+      ideal U;
+      poly c;
+    }
+  }
+    
+  if(dbg>=1)
+  {  
+    prim; "";
+    "// number of components is", size(prim);
+    "";
+  }
+   
+  // ---------------- normalization of the components-------------------------
+  // calls normalM to compute the normalization of each component.
+
+  list norComp;       // The normalization of each component.
+  int delt;
+  int deltI = 0;
+  int totalComps = 0;
+
+  setring origR;
+  def newROrigOrd;
+  list newRListO;
+  setring R;
+  def newR;
+  list newRList;
+    
+  for(i=1; i<=size(prim); i++)
+  {
+    if(dbg>=2){pause();}
+    if(dbg>=1){
+      "// start computation of component",i;
+      "   --------------------------------";
+    }
+    if(groebner(prim[i])[1] != 1){
+      if(dbg>=2){
+        "We compute the normalization in the ring"; basering;
+      }
+      printlevel = printlevel + 1;
+      norComp = normalM(prim[i], decomp, withDelta);
+      printlevel = printlevel - 1;      
+      for(j = 1; j <= size(norComp); j++){
+        newR = norComp[j][3];
+        newRList = ringlist(newR);
+        U = norComp[j][1];
+        c = norComp[j][2];
+        if(withDelta){
+          delt = norComp[j][4];
+          if((delt >= 0) and (deltI >= 0)){
+            deltI = deltI + delt;
+          } else {
+            deltI = -1;
+          }
+        }
+        if(useRing == 0){
+        // We go back to the original ring.
+        setring origR;
+        U = fetch(R, U);
+        c = fetch(R, c);
+        newRListO = imap(R, newRList); 
+        // We change the ordering in the new ring.
+        if(nvars(newR) > nvars(origR)){
+          newRListO[3]=insert(origOrd, newRListO[3][1]);
+        } else {
+          newRListO[3] = origOrd;
+        }
+        newROrigOrd = ring(newRListO);
+        setring newROrigOrd;
+        ideal norid = imap(newR, norid);
+        ideal normap = imap(newR, normap);
+        export norid;
+        export normap;
+        setring origR;
+        totalComps++;
+        result[totalComps] = list(U, c, newROrigOrd);
+        if(withDelta){
+          result[totalComps] = insert(result[totalComps], delt, 3);
+        }
+        setring R;
+        } else {
+          setring R;
+          totalComps++;
+          result[totalComps] = norComp[j];
+        }
+      }
+    }
+  }
+  
+  // -------------------------- delta computation ----------------------------  
+  if(withDelta == 1){
+    // Intersection multiplicities of list prim, sp=size(prim).
+    if ( dbg >= 1 )
+    {
+      "// Sum of delta for all components: ", deltI;
+    }   
+    if(size(prim) > 1){
+      dbprint(dbg, "// Computing the sum of the intersection multiplicities of the components...");
+      int mul = iMult(prim);
+      if ( mul < 0 )
+      {
+        deltI = -1;
+      } else {
+        deltI = deltI + mul; 
+      }
+      if ( dbg >= 1 )
+      {
+        "// Intersection multiplicity is : ", mul;
+      }         
+    }
+  }
+  
+  // -------------------------- prepare output ------------------------------   
+  setring origR;
+  if(withDelta){
+    result = result, deltI;
+  } else {
+    result = list(result);
+  }
+  
+  option(set, opt);
+  return(result);
+}
+
+
+example
+{ "EXAMPLE:";
+  printlevel = printlevel+1;
+  echo = 2;
+  ring s = 0,(x,y),dp;
+  ideal i = (x2-y3)*(x2+y2)*x;
+  list nor = normal(i, "withDelta");
+  nor; 
+
+  // 2 branches have delta = 1, and 1 branch has delta = 0
+  // the total delta invariant is 13
+
+  def R2 = nor[1][2][3];  setring R2;
+  norid; normap;
+   
+  echo = 0;
+  printlevel = printlevel-1;
+  pause("   hit return to continue"); echo=2;
+
+  ring r = 2,(x,y,z),dp;
+  ideal i = z3-xy4;
+  list nor = normal(i, "withDelta");  nor;  
+  // the delta invariant is infinite 
+  // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
+  // in its quotient field Quot(r/i)
+
+  // the normalization as affine algebra over the ground field:             
+  def R = nor[1][1][3]; setring R;
+  norid; normap;
+}
+
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -39 +400 @@
 @*       p    = nonzero divisor of R
 COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
-         affine ring, together with the canonical map R --> Hom_R(J,J),
+         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 two objects
+RETURN:  a list l of three objects
 @format
          l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
@@ -47 +408 @@
                endphi describes the canonical map R -> Hom_R(J,J)
          l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
-         l[3] : an integer, the contribution to delta
+         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)
@@ -54 +416 @@
 {
 //---------- initialisation ---------------------------------------------------
-   int isIso,isPr,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
+   int isIso,isPr,isHy,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
    intvec rw,rw1;
    list L;
@@ -103 +465 @@
    if ( y>=1 )
    {
-     "// compute p*Hom(J,J) = p*J:J, p a non-zerodivisor";
-     "//   p is equal to:"; "";
-     p;
-     "";
+     "// 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 module p*Hom(J,J) = p*J:J, p a non-zerodivisor";
-      "// p"; p;
-      "// f=p*J:J";f;
+   {  
+      "//   the non-zerodivisor p:"; p;
+      "//   the module p*Hom(J,J) = p*J:J :"; f;
+      "";
    }
    f2 = std(p);
@@ -119 +492 @@
 //---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
 
-   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R
+   //rf = interred(reduce(f,f2)); 
+   //### interred hier weggelassen, unten zugefŸgt
+   rf = reduce(f,f2);       //represents p*Hom(J,J)/p*R = Hom(J,J)/R
    if ( size(rf) == 0 )
    {
@@ -132 +507 @@
       ideal endphi = maxideal(1);
       ideal endid = fetch(P,id);
-      L=substpart(endid,endphi,homo,rw);
-      def lastRing=L[1];
+      endid = simplify(endid,2);
+      L = substpart(endid,endphi,homo,rw);   //## hier substpart
+      def lastRing = L[1];
       setring lastRing;
 
@@ -142 +518 @@
       attrib(endid,"isRegInCodim2",isRe);
       attrib(endid,"isEqudimensional",isEq);
+      attrib(endid,"isHypersurface",0);
       attrib(endid,"isCompleteIntersection",0);
-      attrib(endid,"isRad",0);
+      attrib(endid,"isRadical",0);
       L=lastRing;
       L = insert(L,1,1);
@@ -150 +527 @@
       {
          "//   R=Hom(J,J)";
-         "   ";
          lastRing;
-         "   ";
          "//   the new ideal";
          endid;
          "   ";
          "//   the old ring";
-         "   ";
          P;
-         "   ";
          "//   the old ideal";
-         "   ";
          setring P;
          id;
          "   ";
          setring lastRing;
-         "//   the map";
-         "   ";
+         "//   the map to the new ring";
          endphi;
          "   ";
          pause();
-         newline;
+         "";
       }
       setring P;
@@ -181 +552 @@
       "// R is not equal to Hom(J,J), we have to try again";
       pause();
-      newline;
+      "";
    }
 //---------- Hom(J,J) != R: create new ring and map from old ring -------------
 // the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
-
-   f=mstd(f)[2];
-   ideal ann=quotient(f2,f);
-   int delt;
-   if(isIso&&isEq){delt=vdim(std(modulo(f,ideal(p))));}
+// f2=p (i.e. ideal generated by p)
+
+   //f = mstd(f)[2];              //### geŠndert GMG 04.10.08 
+   //ideal ann = quotient(f2,f);  //### f durch rf ersetzt  
+   rf = mstd(rf)[2];              //rf = NF(f,p), hence <p,rf> = <p,f>
+   ideal ann = quotient(f2,rf);   //p:f = p:rf 
+ 
+   //------------- compute the contribution to delta ---------- 
+   //delt=dim_K(Hom(JJ)/R (or -1 if infinite)
+  
+   int delt=vdim(std(modulo(f,ideal(p))));    
 
    f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
    q = size(f);
-
    syzf = syz(f);
 
@@ -211 +587 @@
    }
 
-   map psi1 = P,maxideal(1);
-   ideal SBid = psi1(SBid);
+   //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)));
-   ideal id = psi(id);
-   ideal f = psi(f);
-   module syzf = psi(syzf);
+   //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;
@@ -231 +611 @@
 // 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;
@@ -238 +620 @@
    {
       "// the ring structure of Hom(J,J) as R-algebra";
-      " ";
-      "//   the linear relations";
-      " ";
+      "//   the linear relations:";
       Lin;
-      "   ";
-   }
-
+   }
+
+   poly ff;
    for (ii=2; ii<=q; ii++ )
    {
       for ( jj=2; jj<=ii; jj++ )
       {
-         A = lift(pf,f[ii]*f[jj]);
-         Quad = Quad, ideal(T(jj)*T(ii) - T*A);          // quadratic relations
+         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
       }
    }
@@ -257 +638 @@
    {
       "//   the quadratic relations";
-      "   ";
-      interred(Quad);
+      Quad;
       pause();
       newline;
@@ -264 +644 @@
    Q = Lin,Quad;
    Q = subst(Q,T(1),1);
-   Q=mstd(Q)[2];
+   //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)
@@ -275 +659 @@
    }
 
-   ideal endid  = imap(newR,id),imap(newR,Q);
+   ideal endid  = imap(newR,id),imap(newR,Q);  
+   //hier wird Q weiterverwendet, die SB-Eigenschaft wird nicht verwendet.
+   endid = simplify(endid,2);
    ideal endphi = ideal(X(1..nvars(R)));
 
-   L=substpart(endid,endphi,homo,rw);
+   L = substpart(endid,endphi,homo,rw);
 
    def lastRing=L[1];
@@ -287 +673 @@
    ideal an=sigma(ann);
    export(an);  //noetig?
-   ideal te=an,endid;
-   if(isIso&&(size(reduce(te,std(maxideal(1))))==0))
-   {
-      attrib(endid,"onlySingularAtZero",oSAZ);
-   }
-   kill te;
-   attrib(endid,"isCohenMacaulay",isCo);
+   //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,"isRad",0);
+   attrib(endid,"isRadical",0);
    if(y>=1)
    {
-      "//   the new ring after reduction of the number of variables";
-      "   ";
+      "// the new ring after reduction of the number of variables";
       lastRing;
-      "   ";
       "//   the new ideal";
-      "   ";
-      endid;
-      "   ";
-      "//   the old ring";
-      "   ";
+      endid;  "";
+      "// the old ring";
       P;
-      "   ";
       "//   the old ideal";
-      "   ";
       setring P;
       id;
       "   ";
       setring lastRing;
-      "//   the map";
-      "   ";
+      "//   the map to the new ring";
       endphi;
       "   ";
       pause();
-      newline;
+      "";
    }
    L = lastRing;
    L = insert(L,0,1);
-   L[3]=delt;
+   L[3] = delt;
    return(L);
 }
@@ -338 +717 @@
   ideal J  = x,y;
   poly p   = x;
-  list Li = std(id),id,J,p;
+  list Li  = std(id),id,J,p;
   list L   = HomJJ(Li);
   def end = L[1];    // defines ring L[1], containing ideals endid, endphi
@@ -346 +725 @@
   map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
   psi;
+  L[3];              // contribution to delta
 }
 
+
 ///////////////////////////////////////////////////////////////////////////////
-
-proc normal(ideal id, list #)
-"USAGE:   normal(i [,choose]);  i a radical ideal, choose empty, 1 ,\"g\" \"wd\"
-         if choose=1 the normalization of the associated primes is computed
-         (which is sometimes more efficient);
-         if @code{choose=\"wd\"} the delta invariant is computed
-         simultaneously; this may take much more time in the reducible case,
-         since the factorizing standard basis algorithm cannot be used.@*
-         if @code{choose=\"g\"} the generators of the integral closure
-         are computed.
-ASSUME:  The ideal must be radical, for non-radical ideals the output may
-         be wrong (i=radical(i); makes i radical)
-RETURN:  a list of rings, say nor and in case of @code{choose=\"wd\"} an
-         integer at the end of the list.
-         In case of @code{choose=\"g\"} it returns a list of polynomials
-         g_1,...,g_k+1 at the end of the list such that g_1/g_k+1,...
-         g_k/g_k+1 generate the integral closure.
-         Each ring @code{nor[i]} contains two ideals with given names
-         @code{norid} and @code{normap} such that@*
-         - the direct sum of the rings @code{nor[i]/norid} is the
-         normalization of basering/id;@*
-         - @code{normap} gives the normalization map from basering/id to
-         @code{nor[i]/norid} (for each i).
-NOTE:    to use the i-th ring type: @code{def R=nor[i]; setring R;}.
-@*       Increasing printlevel displays more comments (default: printlevel=0).
-@*       Not implemented for local or mixed orderings.
-@*       If the input ideal i is weighted homogeneous a weighted ordering may
-         be used (qhweight(i); computes weights).
-KEYWORDS: normalization; delta invariant.
-EXAMPLE: example normal; shows an example
+//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,j,y,withdelta,withgens;
-   string sr;
-   list result,prim,keepresult;
-   y = printlevel-voice+2;
-   if(size(#)>0)
-   {
-     if(typeof(#[1])=="string")
+{    int i,mul,mu;
+     int sp = size(prim);
+     int y = printlevel-voice+2;
+     if ( sp > 1 )
      {
-       if(#[1]=="g")
-       {
-          withgens=1;
-       }
-       else
-       {
-          withdelta=1;
-       }
-       kill #;
-       list #;
+        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;
+              }
+           }
+        }
      }
-   }
-   attrib(id,"isRadical",1);
-   if ( ord_test(basering) != 1)
-   {
-     "";
-     "// Not implemented for this ordering,";
-     "// please change to global ordering!";
-     return(result);
-   }
-   if(withgens)
-   {
-      list pr=minAssGTZ(id);
-      if(size(pr)>1){ERROR("ideal is not prime");}
-      def R=basering;
-      if(defined(ker)){kill ker;}
-      ideal ker=id;
-      export(ker);
-      list l=R;
-      l=primeClosure(l);
-      list gens=closureGenerators(l);
-      list resu=l[size(l)],gens;
-      dbprint(y+1,"
-// 'normal' created a list of one ring.
-// nor[2] is a list of elements of the basering such that
-// nor[2][i]/nor[2][size(nor[2])] generate the integral closure.
-// To see the rings, type (if the name of your list is nor):
-     show( nor);
-// To access the 1-st ring and map (similar for the others), type:
-     def R = nor[1]; setring R;  norid; normap;
-// R/norid is the 1-st ring of the normalization and
-// normap the map from the original basering to R/norid");
-
-      return(resu);
-   }
-   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
-   {
-      if(attrib(id,"isCompleteIntersection")==1)
-      {
-         attrib(id,"isCohenMacaulay",1);
-         attrib(id,"isEquidimensional",1);
-      }
-   }
-   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
-   {
-      if(attrib(id,"isCohenMacaulay")==1)
-      {
-         attrib(id,"isEquidimensional",1);
-      }
-   }
-   if( typeof(attrib(id,"isPrim"))=="int" )
-   {
-      if(attrib(id,"isPrim")==1)
-      {
-         attrib(id,"isEquidimensional",1);
-      }
-   }
-   if(size(#)==0)
-   {
-      if( typeof(attrib(id,"isEquidimensional"))=="int" )
-      {
-         if(attrib(id,"isEquidimensional")==1)
-         {
-            prim[1]=id;
-         }
-         else
-         {
-            prim=equidim(id);
-         }
-      }
-      else
-      {
-         prim=equidim(id);
-      }
-      if(y>=1)
-      {
-         "// we have ",size(prim),"equidimensional components";
-      }
-      if(withdelta &&(size(prim)>1))
-      {
-         withdelta=0;
-         "WARNING!  cannot compute delta,";
-         "the ideal is not equidimensional";
-      }
-      if(!withdelta)
-      {
-         option(redSB);
-         for(j=1;j<=size(prim);j++)
-         {
-            keepresult=keepresult+facstd(prim[j]);
-         }
-         prim=keepresult;
-         option(noredSB);
-      }
-   }
-   else
-   {
-      if( typeof(attrib(id,"isPrim"))=="int" )
-      {
-         if(attrib(id,"isPrim")==1)
-         {
-            prim[1]=id;
-         }
-         else
-         {
-            prim=minAssGTZ(id);
-         }
-      }
-      else
-      {
-         prim=minAssGTZ(id);
-      }
-      if(y>=1)
-      {
-         "// we have ",size(prim),"irreducible components";
-      }
-   }
-   for(i=1; i<=size(prim); i++)
-   {
-      if(y>=1)
-      {
-         "// we are in loop ",i;
-      }
-      attrib(prim[i],"isCohenMacaulay",0);
-      if(size(#)!=0)
-      {
-         attrib(prim[i],"isPrim",1);
-      }
-      else
-      {
-         attrib(prim[i],"isPrim",0);
-      }
-      attrib(prim[i],"isRegInCodim2",0);
-      attrib(prim[i],"isIsolatedSingularity",0);
-      attrib(prim[i],"isEquidimensional",1);
-      attrib(prim[i],"isCompleteIntersection",0);
-      attrib(prim[i],"onlySingularAtZero",0);
-      if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
-      {
-            if(attrib(id,"onlySingularAtZero")==1)
-             {attrib(prim[i],"onlySingularAtZero",1); }
-      }
-
-      if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
-      {
-            if(attrib(id,"isIsolatedSingularity")==1)
-             {attrib(prim[i],"isIsolatedSingularity",1); }
-      }
-
-      if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
-      {
-            if((attrib(id,"isIsolatedSingularity")==1)&&(size(#)==0))
-             {attrib(prim[i],"isIsolatedSingularity",1); }
-      }
-      keepresult=normalizationPrimes(prim[i],maxideal(1),0);
-
-      for(j=1;j<=size(keepresult)-1;j++)
-      {
-         result=insert(result,keepresult[j]);
-      }
-      sr = string(size(result));
-   }
-   if(withdelta)
-   {
-      sr = string(size(keepresult)-1);
-      result=keepresult;
-   }
-      dbprint(y+1,"
-// 'normal' created a list of "+sr+" ring(s).");
-      if(withdelta)
-      {
-        dbprint(y+1,"// nor["+sr+"+1] is the delta-invariant.");
-      }
-      dbprint(y+1,"// To see the rings, type (if the name of your list is nor):
-     show( nor);
-// To access the 1-st ring and map (similar for the others), type:
-     def R = nor[1]; setring R;  norid; normap;
-// R/norid is the 1-st ring of the normalization and
-// normap the map from the original basering to R/norid");
-
-      //kill endphi,endid;
-      return(result);
+     return(mul);
 }
 example
 { "EXAMPLE:"; echo = 2;
-   ring r=32003,(x,y,z),wp(2,1,2);
-   ideal i=z3-xy4;
-   list nor=normal(i);
-   show(nor);
-   def r1=nor[1];
-   setring r1;
-   norid;
-   normap;
-
-   ring s=0,(x,y),dp;
-   ideal i=(x-y^2)^2 - y*x^3;
-   nor=normal(i,"wd");
-   //the delta-invariant
-   nor[size(nor)];
-   nor=normal(i,"g");
-   nor[size(nor)];
+   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);
+}  
 
 ///////////////////////////////////////////////////////////////////////////////
-static proc normalizationPrimes(ideal i,ideal ihp,int delt, list #)
-"USAGE:   normalizationPrimes(i,ihp[,si]);  i equidimensional ideal, ihp map
-         (partial normalization),delta partial delta-invariant,
-          # ideal of singular locus
+
+//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[i] contains two ideals
+          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[i]/norid is
-             the normalization of basering/id;
+           - 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[i]/norid (for each i)
+             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";  "";
+     "// START a normalization loop with the ideal";
      i;  "";
+     "// in the ring:";
      basering;  "";
      pause();
-     newline;
+     "";
    }
 
    def BAS=basering;
-   list result,keepresult1,keepresult2;
+   list result,keepresult1,keepresult2,JM,gnirlist;
    ideal J,SB,MB;
-   int depth,lauf,prdim;
+   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)
    {
@@ -643 +921 @@
           "// the ideal was the zero-ideal";
       }
-         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
-                      +ordstr(basering)+");");
+     // execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
+     //                 +ordstr(basering)+");");
+         def newR7 = ring(gnirlist);
+         setring newR7;
          ideal norid=ideal(0);
          ideal normap=fetch(BAS,ihp);
@@ -650 +930 @@
          export normap;
          result=newR7;
-         result[size(result)+1]=delt;
+         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 input ideal";
-   }
-   list SM=mstd(i);                //here the work starts
+     "// SB-computation of the ideal";
+   }
+
+   list SM = mstd(i);              //Now the work starts
    int dimSM =  dim(SM[1]);        //dimension of variety to normalize
-  // Case: Get an ideal containing a unit
+   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)
    {  "";
@@ -678 +1055 @@
          export normap;
          result=newR6;
-         result[size(result)+1]=delt;
+         result[size(result)+1]=list(delt,delti);
          setring BAS;
          return(result);
    }
-
-   if(y>=1)
-   {
-      "//   the dimension is:"; "";
-      dimSM;"";
-   }
-   if(size(#)>0)
-   {
-      if(attrib(i,"onlySingularAtZero"))
-      {
-         list JM=maxideal(1),maxideal(1);
-         attrib(JM[1],"isSB",1);
-         attrib(JM[2],"isRad",1);
-      }
-      else
-      {
-         ideal te=#[1],SM[2];
-         list JM=mstd(te);
-         kill te;
-         if( typeof(attrib(#[1],"isRad"))!="int" )
-         {
-            attrib(JM[2],"isRad",0);
-         }
-      }
-   }
-
-   if(attrib(i,"isPrim")==1)
-   {
-      attrib(SM[2],"isPrim",1);
-   }
-   else
-   {
-      attrib(SM[2],"isPrim",0);
-   }
-   if(attrib(i,"isIsolatedSingularity")==1)
-   {
-      attrib(SM[2],"isIsolatedSingularity",1);
-   }
-   else
-   {
-      attrib(SM[2],"isIsolatedSingularity",0);
-   }
-   if(attrib(i,"isCohenMacaulay")==1)
-   {
-      attrib(SM[2],"isCohenMacaulay",1);
-   }
-   else
-   {
-      attrib(SM[2],"isCohenMacaulay",0);
-   }
-   if(attrib(i,"isRegInCodim2")==1)
-   {
-      attrib(SM[2],"isRegInCodim2",1);
-   }
-   else
-   {
-      attrib(SM[2],"isRegInCodim2",0);
-   }
-   if(attrib(i,"isEquidimensional")==1)
-   {
-      attrib(SM[2],"isEquidimensional",1);
-   }
-   else
-   {
-      attrib(SM[2],"isEquidimensional",0);
-   }
-   if(attrib(i,"isCompleteIntersection")==1)
-   {
-      attrib(SM[2],"isCompleteIntersection",1);
-   }
-   else
-   {
-      attrib(SM[2],"isCompleteIntersection",0);
-   }
-   if(attrib(i,"onlySingularAtZero")==1)
-   {
-      attrib(SM[2],"onlySingularAtZero",1);
-   }
-   else
-   {
-      attrib(SM[2],"onlySingularAtZero",0);
-   }
-   if((attrib(SM[2],"isIsolatedSingularity")==1)&&(homog(SM[2])==1))
-   {
-      attrib(SM[2],"onlySingularAtZero",1);
-   }
-
-   //the smooth case
-   if(size(#)>0)
-   {
-      if(dim(JM[1])==-1)
-      {
-         if(y>=1)
-         {
-            "// the ideal was smooth";
-         }
-         MB=SM[2];
-         intvec rw;;
-         list LL=substpart(MB,ihp,0,rw);
-         def newR6=LL[1];
-         setring newR6;
-         ideal norid=endid;
-         ideal normap=endphi;
-         kill endid,endphi;
-         export norid;
-         export normap;
-         result=newR6;
-         result[size(result)+1]=delt;
-         setring BAS;
-         return(result);
-     }
-   }
-
-   //the zero-dimensional case
-   if((dim(SM[1])==0)&&(homog(SM[2])==1))
+   
+   // --- Trivial case: input ideal is zero-dimensional and homog ---
+   if( (dim(SM[1])==0) && (homog(SM[2])==1) )
    {
       if(y>=1)
@@ -803 +1068 @@
       }
       MB=maxideal(1);
-      intvec rw;
-      list LL=substpart(MB,ihp,0,rw);
+      intvec rw; 
+      list LL=substpart(MB,ihp,0,rw);     
       def newR5=LL[1];
       setring newR5;
@@ -813 +1078 @@
       export normap;
       result=newR5;
-      result[size(result)+1]=delt;
+      result[size(result)+1]=list(delt,delti);
       setring BAS;
       return(result);
    }
 
-   //the one-dimensional case
-   //in this case it is a line because
-   //it is irreducible and homogeneous
-   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
-        &&(homog(SM[2])==1))
+   // --- 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)
@@ -839 +1102 @@
       export normap;
       result=newR4;
-      result[size(result)+1]=delt;
+      result[size(result)+1]=list(delt,delti);
       setring BAS;
       return(result);
    }
 
+//---------------------- The non-trivial cases start -------------------  
    //the higher dimensional case
-   //we test first of all CohenMacaulay and
-   //complete intersection
-   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
-   {
-      //test for complete intersection
+   //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);
@@ -858 +1122 @@
       }
    }
-   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(i)==1))
-   {
-      int tau=vdim(std(i+jacob(i)));
-      if(tau>=0)
-      {
-        execute("ring BASL="+charstr(BAS)+",("+varstr(BAS)+"),ds;");
-        ideal i=imap(BAS,i);
-        int tauloc=vdim(std(i+jacob(i)));
-        setring BAS;
-        attrib(SM[2],"onlySingularAtZero",(tau==tauloc));
-        kill BASL;
-      }
-   }
-   //compute the singular locus
-   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(#)==0))
-   {
-      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
-      if(y >=1 )
-      {
-         "// SB of singular locus will be computed";
-      }
-      ideal sin=J,SM[1];
-      list JM=mstd(sin);
+   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);
-
-      //JM[1] SB of singular locus, JM[2]=minbasis of singular locus
-      //SM[1] SB of the input ideal, SM[2] minbasis
-      if(y>=1)
-      {
-         "//   the dimension of the singular locus is:";"";
-         dim(JM[1]); "";
-      }
-
-      if(dim(JM[1])==-1)
+      if( y>=1 )
+      {  
+         "// the 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)
@@ -908 +1212 @@
          export normap;
          result=newR3;
-         result[size(result)+1]=delt;
+         result[size(result)+1]=list(delt,delti);
          setring BAS;
          return(result);
       }
-      if(dim(JM[1])==0)
-      {
-         attrib(SM[2],"isIsolatedSingularity",1);
-         if(homog(SM[2])){attrib(SM[2],"onlySingularAtZero",1);}
-      }
-      if(dim(JM[1])<=dim(SM[1])-2)
-      {
-         attrib(SM[2],"isRegInCodim2",1);
-      }
-   }
-   else
-   {
-     if(size(#)==0)
-     {
-        list JM=maxideal(1),maxideal(1);
-
-        attrib(JM[1],"isSB",1);
-        if(dim(SM[1])>=2){attrib(SM[2],"isRegInCodim2",1);}
-     }
-   }
-   if((attrib(SM[2],"isRegInCodim2")==1)&&(attrib(SM[2],"isCohenMacaulay")==1))
+
+      //------- 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 codimension 2";
+         "// the ideal was CohenMacaulay and regular in codim 2, hence normal";
       }
       MB=SM[2];
@@ -949 +1271 @@
       export normap;
       result=newR6;
-      result[size(result)+1]=delt;
+      result[size(result)+1]=list(delt,delti);
       setring BAS;
       return(result);
    }
 
-   //if it is an isolated singularity only at 0 things are easier
-   //JM ideal of singular locus, SM ideal of variety
-   if(attrib(SM[2],"onlySingularAtZero"))
-   {
+//---------- 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);
-      //vars not contained in ideal
-      ideal Ann=quotient(SM[2],SL[1]);
-      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
-      //qAnn=0 ==> the first var(=SL[1]) not contained in SM is a nzd of R/SM
-      if(size(qAnn)==0)
+      ideal SL = simplify(reduce(maxideal(1),SM[1]),2);
+      ideal Ann = quotient(SM[2],SL[1]);
+      ideal qAnn = simplify(reduce(Ann,SM[1]),2);
+      //NOTE: qAnn=0 iff first var (=SL[1]) not in SM is a nzd of R/SM
+
+   //------------- We found a non-zerodivisor of R/SM -----------------------
+   // here the enlarging of the ring via Hom_R(J,J) starts
+
+      if( size(qAnn)==0 )
       {
          if(y>=1)
          {
             "";
-            "//   the ideal rad(J):";
+            "// the ideal rad(J):"; maxideal(1);
             "";
-            maxideal(1);
-            newline;
          }
-         //again test for normality
-         //Hom(I,R)=R
+
+      // ------------- 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);
-         if(RR[2]==0)
+         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 tluser=normalizationPrimes(endid,psi(ihp),delt+RR[3],an);
-            setring BAS;
-            return(tluser);
+            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)+");");
+         //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);
-         delt=delt+RR[3];
+         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]=delt;
+         result[size(result)+1]=list(delt,delti);
          setring BAS;
          return(result);
-
-       }
-      //Now the case where qAnn!=0, i.e.SL[1] is a zero divisor of R/SM
-      //and we have found a splitting: id and id1
-      //id=Ann defines components of R/SM in the complement of V(SL[1])
-      //id1 defines components of R/SM in the complement of V(id)
-       else
+      }
+
+   //------ 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;
+          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);
-          //ideal id=SL[1],SM[2];
+          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);
 
-          ideal t1=id,id1;
-          int mul=vdim(std(t1));
-          kill t1;
-
-          keepresult1=normalizationPrimes(id,ihp,0);
-
-          keepresult2=normalizationPrimes(id1,ihp,0);
-
-          delt=delt+mul+keepresult1[size(keepresult1)]
-                         +keepresult1[size(keepresult1)];
-
+          // ---------- 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)]=delt;
+          keepresult1[size(keepresult1)]=list(delt,delti);
+
           return(keepresult1);
        }
    }
-
-   //test for non-normality
-   //Hom(I,I)<>R
+   // 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
-   ideal SL=simplify(reduce(JM[2],SM[1]),2);
-   ideal Ann=quotient(SM[2],SL[1]);
-   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
-
-   if(size(qAnn)==0)
+
+   //------ 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;
-      //now we have to compute the radical
+      // ----------------- Computation of the radical -----------------
       if(y>=1)
       {
          "// radical computation of singular locus";
       }
-      J=radical(JM[2]);   //the singular locus
-      JM=mstd(J);
-
-      if((vdim(JM[1])==1)&&(size(reduce(J,std(maxideal(1))))==0))
-      {
-         attrib(SM[2],"onlySingularAtZero",1);
-      }
+      J = radical(JM[2]);   //the radical of singular locus
+      JM = mstd(J);
+
       if(y>=1)
       {
-        "//   radical is equal to:";"";
-        J;
+        "// radical is equal to:";"";  JM[2]; 
         "";
       }
-      if(deg(SL[1])>deg(J[1]))
+      // ------------ 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];}
-      }
-
-      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
+         if(size(qAnn)==0)
+         {
+           SL[1]=J[1];
+         }
+      }
+
+      // --------------- computation of Hom(rad(J),rad(J)) --------------
       RR=SM[1],SM[2],JM[2],SL[1];
 
-      //   evtl eine geeignete Potenz von JM?
      if(y>=1)
      {
@@ -1093 +1534 @@
      }
 
-     RS=HomJJ(RR,y);
-
-      if(RS[2]==1)
-      {
+     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;
@@ -1103 +1546 @@
          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;
-         result[size(result)+1]=delt+RS[3];
+         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);
-      }
-      int n=nvars(basering);
-      ideal MJ=JM[2];
+     }
+
+    // ----- 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+RS[3],psi(MJ));
+          normalizationPrimes(endid,psi(ihp),delt,delti,psi(MJ));
+      printlevel = printlev;                //reset printlevel
       setring BAS;
       return(tluser);
    }
-    // A component with singular locus the whole component found
+
+   //---- A whole singular component was found, RETURN -----
    if( Ann == 1)
    {
@@ -1141 +1649 @@
          export normap;
          result=newR6;
-         result[size(result)+1]=delt;
+         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
    {
-      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);
+      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];
 
-      int mul;
-      if(equi&&isIs)
-      {
-          ideal t2=new1,new2;
-          mul=vdim(std(t2));
-      }
-      execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
-                      +ordstr(basering)+");");
-      if(y>=1)
-      {
-         "// zero-divisor found";
-      }
-      ideal vid=fetch(BAS,new1);
-      ideal ihp=fetch(BAS,ihp);
+      //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);
@@ -1175 +1687 @@
       attrib(vid,"onlySingularAtZero",oSAZ);
       attrib(vid,"isEquidimensional",equi);
+      attrib(vid,"isHypersurface",0);
       attrib(vid,"isCompleteIntersection",0);
 
-      keepresult1=normalizationPrimes(vid,ihp,0);
-      int delta1=keepresult1[size(keepresult1)];
+      // ---------- 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)+");");
-
-      ideal vid=fetch(BAS,new2);
-      ideal ihp=fetch(BAS,ihp);
+      //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);
@@ -1191 +1715 @@
       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);
-      int delta2=keepresult2[size(keepresult2)];
-
+      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]);
-      }
-      keepresult1[size(keepresult1)]=delt+mul+delta1+delta2;
+         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);
    }
@@ -1222 +1770 @@
 
    list pr=normalizationPrimes(i);
-   def r1=pr[1];
+   def r1 = pr[1];
    setring r1;
    norid;
    normap;
 }
+
 ///////////////////////////////////////////////////////////////////////////////
-proc substpart(ideal endid, ideal endphi, int homo, intvec rw)
+static proc substpart(ideal endid, ideal endphi, int homo, intvec rw)
 
 "//Repeated application of elimpart to endid, until no variables can be
@@ -1234 +1783 @@
 //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;
+   def newRing = basering;
    int ii,jj;
-   map phi = basering,maxideal(1);
-
+   map phi = newRing,maxideal(1);    //identity map
    list Le = elimpart(endid);
-                                     //this proc and the next loop try to
-   int q = size(Le[2]);              //substitute as many variables as possible
-   intvec rw1 = 0;                   //indices of substituted variables
+   //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 = 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])-1;ii++)
+      for(ii=1;ii<=size(Le[2]);ii++)
       {
          phi=phi(phi);
       }
       //eingefuegt wegen x2-y2z2+z3
-      kill ps;
 
       for( ii=1; ii<=size(rw1); ii++ )
       {
-         if( Le[4][ii]==0 )
+         if( Le[4][ii]==0 )        //ii = index of var which was substituted
          {
-            rw1[ii]=0;                             //look for substituted vars
+            rw1[ii]=0;             //substituted vars have entry 0 in rw1
          }
       }
-      Le=elimpart(endid);
+      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
@@ -1296 +1875 @@
       ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
    }
-
    ideal lastmap;
-   q = 1;
+   jj = 1;
+
    for(ii=1; ii<=size(rw1); ii++ )
    {
-      if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; }
+      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);
+   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);
@@ -1313 +1901 @@
    return(L);
 }
-/////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
 
 proc genus(ideal I,list #)
-"USAGE:   genus(I) or genus(i,1); I a 1-dimensional ideal
+"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.
+         curve defined by i, where p_a is the arithmetic genus.
 NOTE:    delta is the sum of all local delta-invariants of the singularities,
          i.e. dim(R'/R), R' the normalization of the local ring R of the
-         singularity.
-         genus(i,1) uses the normalization to compute delta. Usually this
-         is slow but sometimes not.
+         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
 "
@@ -1336 +1924 @@
       {
         // ERROR("This is not a curve");
-        if(w==1){"WARNING:This is not a curve";}
+        if(w==1){"** WARNING: Input does not define a curve **"; "";}
       }
    }
@@ -1357 +1945 @@
       return(-deg(p)+1);
    }
-   if(size(N[1])<nvars(R0))
+   if(size(N[1]) < nvars(R0))
    {
       string newvar=string(N[1]);
@@ -1562 +2150 @@
    if(size(#)!=0)
    {                              //uses the normalization to compute delta
-      list nor=normal(I,"wd");
-      delt=nor[size(nor)];
+      list nor=normal(I);
+      delt=nor[size(nor)][2];
       genus=genus-delt-deltainf;
       setring R0;
@@ -1953 +2541 @@
 "USAGE:    primeClosure(L [,c]); L a list of a ring containing a prime ideal
                                  ker, c an optional integer
-RETURN:   a list L consisting of rings L[1],...,L[n] such that
+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
@@ -1963 +2551 @@
             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:
@@ -1971 +2560 @@
 "
 {
-// Start with a consistency check:
+  //---- 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;
-
-
-// Some auxiliary variables:
-
-  int n=size(L);
-
-// How to choose the non-zerodivisor later on?
+  }
+
+  int dblvl = printlevel-voice+2;
+  list gnirlist = ringlist(basering);
+
+  //---- Some auxiliary variables:
+  int delt;                      //finally the delta invariant
+  if ( size(L) == 1 )
+  {
+      L[2] = delt;              //set delta to 0
+  }
+  int n = size(L)-1;            //L without delta invariant
+
+  //---- How to choose the non-zerodivisor later on?
 
   int nzdoption=0;
   if (size(#)>0)
-    {
+  {
       nzdoption=#[1];
-    }
+  }
 
 // R0 is the ring to work with, if we are in step one, make a copy of the
@@ -1999 +2592 @@
 
   if (n==1)
-    {
-      def R=L[1];
-      def BAS=basering;
+  {
+      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)+");");
+      //execute ("ring R0="+charstr(R)+",("+varstr(R)+"),("+ordstr(R)+");");
+      def R0 = ring(Rlist);
+      setring R0;
       ideal ker=fetch(R,ker);
 
@@ -2019 +2615 @@
         attrib(R0,"iso_sing_Rees",attrib(L[1],"iso_sing_Rees"));
       }
-      L=R0;
-    }
+      L[1]=R0;
+  }
   else
-    {
-      def R0=L[n];
+  {
+      def R0 = L[n];
       setring R0;
-    }
+  }
 
 
@@ -2031 +2627 @@
 // locus of ker, J:=rad(ker):
 
-//  if (homog(ker)==1)
-//    {
-      list SM=mstd(ker);
-//    }
-/*  else
-    {
-      list SM=groebner(ker),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.
-
-  if (n==1)
-    {
+// 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")
       {
@@ -2053 +2644 @@
           J=J,var(s);
         }
+        J = J,SM[2];
+        list JM = mstd(J);
       }
       else
       {
-        ideal J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
-      }
-      J=J+SM[2];
-      if (homog(J)==1)
-      {
-        J=mstd(J)[2];
-      }
-      J=radical(interred(J));
-    }
+        if ( dblvl >= 1 )
+        {"";
+           "// compute the singular locus";
+        }
+        //### Berechnung des singulŠren Orts geŠndert (ist so schneller)
+        ideal J = minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]);
+        J = J,SM[2];
+        list JM = mstd(J);
+      }
+
+      if ( dblvl >= 1 )
+      {"";
+         "// dimension of singular locus is", dim(JM[1]);
+         if (  dblvl >= 2 )
+         {"";
+            "// the singular locus is:"; JM[2];
+         }
+      }
+
+      if ( dblvl >= 1 )
+      {"";
+         "// compute radical of singular locus";
+      }
+
+      J = simplify(radical(JM[2]),2);
+      if ( dblvl >= 1 )
+      {"";
+         "// radical of singular locus is:"; J;
+         pause();
+      }
+  }
   else
-    {
-      J=radical(interred(J));
-    }
-
-  // having computed the radical J of the ideal of the singular locus,
-  // we now need to pick an element nzd of J
-
-  poly nzd=J[1];
-  if (nzdoption)
+  {
+      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++)
     {
-      if ( (deg(nzd)>=deg(J[ii])) && (size(nzd)>size(J[ii])) )
-      {
-        nzd=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;
+        }
       }
     }
@@ -2086 +2728 @@
 
   export nzd;
-  list RR=SM[1],SM[2],J,nzd;
-  list RS=HomJJ(RR);
-
-
-//////////////////////////////////////////////////////////////////////////////
-
-
+  list RR = SM[1],SM[2],J,nzd;
+
+  if ( dblvl >= 1 )
+  {"";
+     "// compute the first ring extension:";
+  }
+
+  list RS = HomJJ(RR);
+
+//-------------------------------------------------------------------------
 // If we've reached the integral closure (as determined by the result of
 // HomJJ), then we are done, otherwise we have to prepare the next iteration.
 
-  if (RS[2]==1)               // we've reached the integral closure
+  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);
     }
@@ -2104 +2765 @@
     {
       if (n==1)               // In the first iteration: keep only the data
-        {                     // needed later on.
-          kill RR,SM;
-
-          export(ker);
-        }
-
-      def R1=RS[1];           // The data of the next ring R1:
+      {                       // 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];
+      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)
@@ -2131 +2803 @@
       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
+      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.
-        }
+      }
     }
 }
@@ -2162 +2845 @@
 }
 
-
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-
-proc closureRingtower(list L)
-"USAGE:    closureRingtower(list L); L a list of rings
-CREATE:   rings R(1),...,R(n) such that R(i)=L[i] for all i
-EXAMPLE:  example closureRingtower; shows an example
-"
-{
-  int n=size(L);
-
-  for (int i=1;i<=n;i++)
-    {
-      if (defined(R(i))) {
-        string s="Fixed name R("+string(i)+") leads to conflict with existing "
-              +"object having this name";
-        ERROR(s);
-      }
-      def R(i)=L[i];
-      export R(i);
-    }
-
-  return();
-}
-example
-{
-  "EXAMPLE:"; echo=2;
-  ring R=0,(x,y),dp;
-  ideal I=x4,y4;
-  list L=primeClosure(ReesAlgebra(I)[1]);
-  closureRingtower(L);
-  R(1);
-  R(4);
-  kill R(1),R(2),R(3),R(4);
-}
-
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
 
 proc closureFrac(list L)
-"USAGE:    closureFrac (L); L a list as in the result of primeClosure, L[n]
-                           containing an additional poly f
+"USAGE:    closureFrac (L); L a list of size n+1 as in the result of
+          primeClosure, L[n] contains an additional poly f
 CREATE:   a list fraction of two elements of L[1], such that
           f=fraction[1]/fraction[2] via the injections phi L[i]-->L[i+1].
@@ -2212 +2857 @@
 // Define some auxiliary variables:
 
-  int n=size(L);
-  int j,k,l;
+  int n=size(L)-1;
+  int i,j,k,l;
   string mapstr;
-
-  for (int i=1;i<=n;i++) { def R(i)=L[i]; }
+  for (i=1; i<=n; i++) { def R(i) = L[i]; }
 
 // The quotient representing f is computed as in 'closureGenerators' with
@@ -2224 +2868 @@
 //   - we have to make sure that no more objects of the rings R(i) survive.
 
-  for (j=1;j<=2;j++)
+  for (j=1; j<=2; j++)
     {
       setring R(n);
       if (j==1)
-        {
-          poly p=f;
-        }
+      {
+         poly p=f;
+      }
       else
-        {
-          p=1;
-        }
-
-      for (k=n;k>1;k--)
-        {
-
+      {
+         p=1;
+      }
+
+      for (k=n; k>1; k--)
+      {
           if (j==1)
-            {
-              map phimap=R(k-1),phi;
-            }
+          {
+             map phimap=R(k-1),phi;
+          }
 
           p=p*phimap(nzd);
@@ -2252 +2895 @@
 
           if (j==1)
-            {
-              execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
-                       Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
-                       +"),dp("+string(ncols(phi))+"));");
+          {
+             //### noch abfragen ob Z(i) definiert ist
+             list gnirlist = ringlist(R(k));
+             int n2 = size(gnirlist[2]);
+             int n3 = size(gnirlist[3]);
+             for( i=1; i<=ncols(phi); i++)
+             {
+               gnirlist[2][n2+i] = "Z("+string(i)+")";
+             }
+             intvec V;
+             V[ncols(phi)]=0; V=V+1;
+             gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
+             def S(k) = ring(gnirlist);
+             setring S(k);
+
+             //execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
+             //          Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
+             //          +"),dp("+string(ncols(phi))+"));");
+
               ideal phi = imap(R(k),phi);
-              ideal J= imap (R(k),ker);
+              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)
-            {
-              mapstr=" map backmap = S(k),";
-              for (l=1;l<=nvars(R(k));l++)
-                {
-                  mapstr=mapstr+"0,";
-                }
-              execute (mapstr+"maxideal(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)
@@ -2334 +2996 @@
 }
 
-//////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////
-
-static
-proc closureGenerators(list L);    // called inside normalI
+///////////////////////////////////////////////////////////////////////////////
+// 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);    
 {
-  // In the list L of rings R(1),...,R(n), compute representations of the ring
-  // ring variables of the last ring R(n) as fractions of elements of R(1).
-
-  def Rees=basering;              // when called inside normalI, the Rees
-                                  // Algebra is the current basering.
-
-  // First of all we need some variable declarations...
-
-  list preimages;
-  int n=size(L);                  // the number of rings R(1)-->...-->R(n)
+  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;
-  int length=nvars(L[n]);         // the number of variables of the last ring
-  string mapstr;
-
-  for (int i=1;i<=n;i++) { def R(i)=L[i]; }
+  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):
@@ -2360 +3043 @@
   // Finally, do the same for nzd.
 
-  for (j=1;j<=length+1;j++)
-    {
+  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)
         {
-          poly p;
-        }
-
-      if (j<=nvars(R(n)))
-        {
-          p=var(j);
-        }
-      else
-        {
-          p=1;
-        }
-
-      for (k=n;k>1;k--)
-        {
-
-          if (j==1)
-          {
-            map phimap=R(k-1),phi;      // phimap is the map corresponding
-          }                             // to phi
-
-          p=p*phimap(nzd);
+          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 Gröbner 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).
-            {
-              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)
-            {
-              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));
+          // normal form of p w.r.t. a Groebner basis of
+          // J(k) := <ker(k),Z(l)-phi(k)(l)> in R(k)[Z]=:S(k)
+
+        if (j==1)   // In the first iteration: Create S(k), fetch phi and
+                    // ker(k) and construct the ideal J(k).
+        {
+         //### noch abfragen ob Z(i) definiert ist
+         list gnirlist = ringlist(R(k));
+         int n2 = size(gnirlist[2]);
+         int n3 = size(gnirlist[3]);
+         for( i=1; i<=ncols(phi); i++)
+         {
+            gnirlist[2][n2+i] = "Z("+string(i)+")";
+         }
+         intvec V;
+         V[ncols(phi)]=0;
+         V=V+1;
+         gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
+         def S(k) = ring(gnirlist);
+         setring S(k);
+
+        // execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
+        //           Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
+        //           +"),dp("+string(ncols(phi))+"));");
+
+          ideal phi = imap(R(k),phi);
+          ideal J = imap (R(k),ker);
+          for ( l=1; l<=ncols(phi); l++ )
+          {
+             J=J+(Z(l)-phi[l]);
+          }
+          J = groebner(J);
+          poly h = NF(imap(R(k),p),J);
         }
-
-      // When down to R(1), store the result in the list preimages
-
-      preimages = insert(preimages,p,j-1);
+        else
+        {
+           setring S(k);
+           h = NF(imap(R(k),p),J);
+        }
+
+        setring R(k-1);
+
+        if (j==1)  // In the first iteration: Compute backmap:S(k)-->R(k-1)
+        {
+           ideal maxi;
+           maxi[nvars(R(k))] = 0;
+           maxi = maxi,maxideal(1);
+           map backmap = S(k),maxi;
+         
+           //mapstr=" map backmap = S(k),";
+           //for (l=1;l<=nvars(R(k));l++)
+           //{
+           //  mapstr=mapstr+"0,";
+           //}
+           //execute (mapstr+"maxideal(1);");
+
+           poly p;
+        }
+        p = NF(backmap(h),std(ker));
+     }
+     // Whe are down to R(1), store here the result in the list preimages
+     preimages = insert(preimages,p,j-1);
+  }
+  ideal preim;                  //make the list preimages to an ideal preim
+  for ( i=1; i<=size(preimages); i++ )
+  {
+     preim[i] = preimages[i];
+  }
+  // R(1) was a copy of Rees, so we have to get back to the basering Rees from
+  // the beginning and fetch the result (the ideal preim) to this ring.
+  setring Rees;
+  return (fetch(R(1),preim));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+//                From here: procedures for char p with Frobenius
+///////////////////////////////////////////////////////////////////////////////
+
+proc normalP(ideal id,list #)
+"USAGE:   normalP(id [,choose]); id a radical ideal, choose a comma separated 
+         list of optional strings: \"withRing\" or \"isPrim\" or \"noFac\".@*
+         If choose = \"noFac\", factorization is avoided during the computation
+         of the minimal associated primes; choose = \"isPrim\" disables the 
+         computation of the minimal associated primes (should only be used
+         if the ideal is known to be prime).@*
+ASSUME:  The characteristic of the ground field must be positive. The ideal
+         id is assumed to be radical. However, if choose does not contain
+         \"isPrim\" the minimal associated primes of id are computed first
+         and hence normal computes the normalization of the radical of id.
+         If choose = \"isPrim\" the ideal must be a prime ideal (this is not
+         tested) otherwise the result may be wrong.
+RETURN:  a list, say 'nor' of size 2 (default) or, if \"withRing\" is given,
+         of size 3. @*         
+         nor[1] (resp. nor[2] if \"withRing\" is given) is a list of ideals 
+         Ii, i=1..r, in the basering where r is the number of associated prime 
+         ideals P_i (irreducible components) of id.@*
+         - Ii is an ideal given by polynomials g_1,...,g_k,g_k+1 such that 
+           g_1/g_k+1,...,g_k/g_k+1 generate the integral closure of
+           basering/P_i as module (over the basering) in its quotient field.@*
+         
+         nor[2] (resp. nor[3] if choose=\"withRing\") is a list of an intvec 
+         of size r, the delta invariants of the r components, and an integer, 
+         the total delta invariant of basering/id (-1 means infinite, and 0 
+         that basering/P_i resp. basering/input is normal). @* 
+
+         If the optional string \"withRing\" is given, the ring structure of
+         the normalization is computed too and nor[1] is a list of r rings.@*
+         Each ring Ri = nor[1][i], i=1..r, contains two ideals with given names
+         @code{norid} and @code{normap} such that @*
+         - Ri/norid is the normalization of the i-th rime component P_i, i.e.
+           isomorphic to the integral closure of basering/P_i in its field 
+           of fractions; @*
+         - the direct sum of the rings Ri/norid is the normalization 
+           of basering/id; @*
+         - @code{normap} gives the normalization map from basering/P_i to
+           Ri/norid (for each i).@*
+THEORY:  The delta invariant of a reduced ring K[x1,...,xn]/id is 
+             dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id)
+         We call this number also the delta invariant of id. normalP uses 
+         the qth-power algorithm suggested by Leonard-Pellikaan (using the 
+         Frobenius) in part similair to an implementation by Singh-Swanson. 
+NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
+@*       Increasing/decreasing printlevel displays more/less comments 
+         (default: printlevel = 0).
+@*       Not implemented for local or mixed orderings or quotient rings.
+@*       If the input ideal id is weighted homogeneous a weighted ordering may
+         be used (qhweight(id); computes weights).
+@*       works only in characteristic p > 0.
+KEYWORDS: normalization; integral closure; delta invariant.
+SEE ALSO: normal
+EXAMPLE: example normalP; shows an example
+"
+{
+   int i,j,y, wring, isPrim, noFac, sr, del, co;;
+   intvec deli;
+   list resu, Resu, prim, Gens, mstdid;
+   ideal gens;
+   y = printlevel-voice+2;
+
+   if ( ord_test(basering) != 1)
+   {
+     "";
+     "// Not implemented for this ordering,";
+     "// please change to global ordering!";
+     return(resu);
+   }
+   if ( char(basering) <= 0)
+   {
+     "";
+     "// Algorithm works only in positive characteristic,";
+     "// use procedure 'normal' if the characteristic is 0";
+     return(resu);
+   }
+
+//--------------------------- define the method --------------------------- 
+   string method;                //make all options one string in order to use
+                                 //all combinations of options simultaneously
+   for ( i=1; i<= size(#); i++ )
+   { 
+     if ( typeof(#[i]) == "string" ) 
+     {
+       method = method + #[i];
+     }
+   }
+
+   if ( find(method,"withring") or find(method,"withRing") )
+   {
+     wring=1;
+   }
+   if ( find(method,"isPrim") or find(method,"isprim") )
+   {
+     isPrim=1;
+   }
+   if ( find(method,"noFac") )
+   {
+     noFac=1;
+   }
+
+   kill #;
+   list #;
+//--------------------------- start computation --------------------------- 
+   ideal II,K1,K2;
+
+   //----------- check first (or ignore) if input id is prime -------------
+
+   if ( isPrim )
+   {
+      prim[1] = id;
+      if( y >= 0 )
+      { "";
+    "// ** WARNING: result is correct if ideal is prime (not checked) **";
+    "// disable option \"isPrim\" to decompose ideal into prime components";"";
+      }
+   }
+   else
+   {
+      if(y>=1)
+      {  "// compute minimal associated primes"; }
+
+      if( noFac )
+      { prim = minAssGTZ(id,1); }
+      else
+      { prim = minAssGTZ(id); }
+
+      if(y>=1)
+      {  
+         prim;"";
+         "// number of irreducible components is", size(prim);
+      }
+   }
+
+   //----------- compute integral closure for every component -------------
+
+      for(i=1; i<=size(prim); i++)
+      {
+         if(y>=1)
+         {
+            ""; pause(); "";
+            "// start computation of component",i;
+            "   --------------------------------";
+         }
+         if(y>=1)
+         {  "// compute SB of ideal"; 
+         }
+         mstdid = mstd(prim[i]);
+         if(y>=1)
+         {  "// dimension of component is", dim(mstdid[1]);"";}
+
+      //------- 1-st main subprocedure: compute module generators ----------
+         printlevel = printlevel+1;
+         II = normalityTest(mstdid); 
+
+      //------ compute also the ringstructure if "withRing" is given -------
+         if ( wring )
+         {
+         //------ 2-nd main subprocedure: compute ring structure -----------
+            resu = list(computeRing(II,prim[i])) + resu;
+         }
+         printlevel = printlevel-1;
+
+      //----- rearrange module generators s.t. denominator comes last ------
+         gens=0;
+         for( j=2; j<=size(II); j++ )
+         {
+            gens[j-1]=II[j];
+         }
+         gens[size(gens)+1]=II[1];
+         Gens = list(gens) + Gens;
+      //------------------------------ compute delta -----------------------
+         K1 = mstdid[1]+II;
+         K1 = std(K1);
+         K2 = mstdid[1]+II[1];
+         K2 = std(K2);       
+         // K1 = std(mstdid[1],II);      //### besser 
+         // K2 = std(mstdid[1],II[1]);   //### besser: Hannes, fixen!
+         co = codim(K1,K2);
+         deli = co,deli;
+         if ( co >= 0 && del >= 0 )
+         {
+            del = del + co;
+         }
+         else
+         { del = -1; }
+      }
+
+      if ( del >= 0 )
+      {
+         int mul = iMult(prim);
+         del = del + mul;
+      }
+      else
+      { del = -1; }
+
+      deli = deli[1..size(deli)-1];
+      if ( wring )
+      { Resu = resu,Gens,list(deli,del); }
+      else
+      { Resu = Gens,list(deli,del); }
+
+   sr = size(prim);
+
+//-------------------- Finally print comments and return --------------------
+   if(y >= 0)
+   {"";
+     if ( wring )
+     {
+"// 'normalP' created a list, say nor, of three lists: 
+// nor[1] resp. nor[2] are lists of",sr,"ring(s) resp. ideal(s)
+// and nor[3] is a list of an intvec and an integer.
+// To see the result, type
+     nor;
+// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
+     def R1 = nor[1][1]; setring R1;  norid; normap;
+// for the other rings type first setring r; (if r is the name of your 
+// original basering) and then continue as for the first ring;
+// Ri/norid is the affine algebra of the normalization of the i-th 
+// component r/P_i (where P_i is an associated prime of the input ideal)
+// and normap the normalization map from r to Ri/norid;
+//    nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
+// elements g1..gk of r such that the gj/gk generate the integral 
+// closure of r/P_i as (r/P_i)-module in the quotient field of r/P_i.
+//    nor[3] shows the delta-invariant of each component and of the input
+// ideal (-1 means infinite, and 0 that r/P_i is normal).";
+     }
+     else
+     {
+"// 'normalP' computed a list, say nor, of two lists:
+//    nor[1] is a list of",sr,"ideal(s), where each ideal nor[1][i] consists
+// of elements g1..gk of the basering such that gj/gk generate the integral
+// closure of the i-th component as (basering/P_i)-module in the quotient
+// field of basering/P_i (P_i an associated prime of the input ideal);
+//    nor[2] shows the delta-invariant of each component and of the input ideal 
+// (-1 means infinite, and 0 that r/P_i is normal).";
+     }
+   }
+
+   return(Resu);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r  = 11,(x,y,z),wp(2,1,2);
+   ideal i = x*(z3 - xy4 + x2);
+   list nor= normalP(i); nor;
+   //the result says that both components of i are normal, but i itself
+   //has infinite delta
+   pause("hit return to continue");
+
+   ring s = 2,(x,y),dp;
+   ideal i = y*((x-y^2)^2 - x^3);
+   list nor = normalP(i,"withRing"); nor;
+
+   def R2  = nor[1][2]; setring R2;
+   norid; normap;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Assume: mstdid is the result of mstd(prim[i]), prim[i] a prime component of
+// the input ideal id of normalP.
+// Output is an ideal U s.t. U[i]/U[1] are module generators.
+ 
+static proc normalityTest(list mstdid)
+{
+   int y = printlevel-voice+2;
+   intvec op = option(get);
+   option(redSB);
+   def R = basering;
+   int n, p = nvars(R), char(R);
+   int ii;
+
+   ideal J = mstdid[1];         //J is the SB of I
+   ideal I = mstdid[2];
+   int h = n-dim(J);            //codimension of V(I), I is a prime ideal
+
+   //-------------------------- compute singular locus ----------------------
+   qring Q = J;                 //pass to quotient ring
+   ideal I = imap(R,I);
+   ideal J = imap(R,J);
+   attrib(J,"isSB",1);
+   if ( y >= 1)
+   { "start normality test";  "compute singular locus";}
+
+   ideal M = minor(jacob(I),h,J); //use the command minor modulo J (hence J=0)
+   M = std(M);                    //this makes M much smaller
+   //keep only minors which are not 0 mod I (!) this is important since we 
+   //need a nzd mod I
+
+   //---------------- choose nzd from ideal of singular locus --------------
+   ideal D = M[1];
+   for( ii=2; ii<=size(M); ii++ )            //look for the shortest one
+   {
+      if( size(M[ii]) < size(D[1]) )
+      { 
+          D = M[ii]; 
+      }
+   }
+
+   //--------------- start p-th power algorithm and return ----------------
+   ideal F = var(1)^p;
+   for(ii=2; ii<=n; ii++)
+   {  
+      F=F,var(ii)^p; 
+   }
+
+   ideal Dp=D^(p-1);
+   ideal U=1;
+   ideal K,L;
+   map phi=Q,F;
+   if ( y >= 1)
+   {  "compute module generators of integral closure";
+      "denominator D is:";  D;  
+      pause();
+   }  
+
+   ii=0;
+   list LK;
+   while(1)
+   {
+      ii=ii+1;
+      if ( y >= 1)
+      { "iteration", ii; }
+      L = U*Dp + I;             
+      //### L=interred(L) oder msdt(L)[2]? 
+      //Wird dadurch kleiner aber string(L) wird groesser
+      K = preimage(Q,phi,L);    //### Improvement by block ordering?
+      option(returnSB);
+      K = intersect(U,K);          //K is the new U, it is a SB
+      LK = mstd(K);
+      K = LK[2];
+
+   //---------------------------- simplify output --------------------------
+      if(size(reduce(U,LK[1]))==0)  //previous U coincides with new U
+      {                             //i.e. we reached the integral closure 
+         U=simplify(reduce(U,groebner(D)),2); 
+         U = D,U;
+         poly gg = gcd(U[1],U[size(U)]);
+         for(ii=2; ii<=size(U)-1 ;ii++)
+         {
+            gg = gcd(gg,U[ii]);
+         }
+         for(ii=1; ii<=size(U); ii++)
+         {
+            U[ii]=U[ii]/gg;
+         }
+         U = simplify(U,6);
+         //if ( y >= 1)
+         //{ "module generators are U[i]/U[1], with U:"; U; 
+         //  ""; pause(); }
+         option(set,op);
+         setring R;
+         ideal U = imap(Q,U);
+         return(U);
+      }
+      U=K;
+   }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc substpartSpecial(ideal endid, ideal endphi)
+{
+   //Note: newRing is of the form (R the original basering): 
+   //char(R),(T(1..N),X(1..nvars(R))),(dp(N),...);
+
+   int ii,jj,kk;
+   def BAS = basering;
+   int n = nvars(basering);
+
+   list Le = elimpart(endid); 
+   int q = size(Le[2]);                   //q variables have been substituted
+//Le;"";
+   if ( q == 0 ) 
+   {
+      ideal normap = endphi;
+      ideal norid = endid;
+      export(norid);
+      export(normap);
+      list L = BAS;
+      return(L);
+   }
+
+      list gnirlist = ringlist(basering);
+      endid = Le[1];
+//endphi;"";
+      for( ii=1; ii<=n; ii++)
+      {
+         if( Le[4][ii] == 0 )            //ii=index of substituted var
+         {
+            endphi = subst(endphi,var(ii),Le[5][ii]);
+         }
+      }
+//endphi;"";
+      list g2 = gnirlist[2];             //the varlist
+      list g3 = gnirlist[3];             //contains blocks of orderings
+      int n3 = size(g3);  
+   //----------------- first identify module ordering ------------------ 
+      if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
+      {
+         list gm = g3[n3];              //last blockis module ordering
+         g3 = delete(g3,n3);
+         int m = 0;
+      }
+      else
+      {
+         list gm = g3[1];              //first block is module ordering
+         g3 = delete(g3,1);
+         int m = 1;
+      }
+   //---- delete variables which were substituted and weights  -------- 
+//gnirlist;"";
+      intvec V;
+      int n(0);
+      list newg2;
+      list newg3;
+      for ( ii=1; ii<=n3-1; ii++ )
+      {
+         V = V,g3[ii][2];           //copy weights for ordering in each block
+         if ( ii==1 )               //into one intvector
+         { 
+            V = V[2..size(V)];
+         }
+        // int n(ii) = size(g3[ii][2]);
+        int n(ii) = size(V);
+        intvec V(ii);
+
+        for ( jj = n(ii-1)+1; jj<=n(ii); jj++)
+        {
+//"ii, jj", ii,jj;
+          if( Le[4][jj] !=0 )     //jj=index of var which was not substituted
+          {
+            kk=kk+1;
+            newg2[kk] = g2[jj];   //not substituted var from varlist
+            V(ii)=V(ii),V[jj];    //weight of not substituted variable
+//"V(",ii,")", V(ii);
+          }
+        }
+        if ( size(V(ii)) >= 2 )
+        {
+           V(ii) = V(ii)[2..size(V(ii))];
+           list g3(ii)=g3[ii][1],V(ii);
+           newg3 = insert(newg3,g3(ii),size(newg3));
+//"newg3"; newg3;
+        }
+      }
+//"newg3"; newg3;
+      //newg3 = delete(newg3,1);    //delete empty list
+      
+/*
+//### neue Ordnung, 1 Block fuer alle vars, aber Gewichte erhalten;
+//vorerst nicht realisiert, da bei leonhard1 alte Version (neue Variable T(i)
+//ein neuer Block) ein kuerzeres Ergebnis liefert
+      kill g3;
+      list g3; 
+      V=0;
+      for ( ii= 1; ii<=n3-1; ii++ )
+      {
+        V=V,V(ii);
+      }
+      V = V[2..size(V)];
+
+      if ( V==1 )
+      {    
+         g3[1] = list("dp",V);
+      }
+      else
+      {
+         g3[1] = lis("wp",V);          
+      }
+      newg3 = g3;
+
+//"newg3";newg3;"";
+//### Ende neue Ordnung
+*/
+
+      if ( m == 0 )
+      {
+         newg3 = insert(newg3,gm,size(newg3));
+      }
+      else
+      {
+         newg3 = insert(newg3,gm);
+      }
+      gnirlist[2] = newg2;                    
+      gnirlist[3] = newg3;
+
+//gnirlist;
+      def newBAS = ring(gnirlist);            //change of ring to less vars
+      setring newBAS;
+      ideal normap = imap(BAS,endphi);
+      //normap = simplify(normap,2);
+      ideal norid =  imap(BAS,endid);
+      export(norid);
+      export(normap);
+      list L = newBAS;
+      return(L);
+
+   //Hier scheint interred gut zu sein, da es Ergebnis relativ schnell 
+   //verkleinert. Hier wird z.B. bei leonard1 size(norid) verkleinert aber
+   //size(string(norid)) stark vergroessert, aber es hat keine Auswirkungen
+   //da keine map mehr folgt. 
+   //### Bei Leonard2 haengt interred (BUG)
+   //mstd[2] verkleinert norid nocheinmal um die Haelfte, dauert aber 3.71 sec
+   //### Ev. Hinweis auf mstd in der Hilfe?
+
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Computes the ring structure of a ring given by module generators.
+// Assume: J[i]/J[1] are the module generators in the quotient field
+// with J[1] as universal denominator.
+// If option "noRed" is not given, a reduction in the number of variables is 
+// attempted.
+static proc computeRing(ideal J, ideal I, list #)
+{
+  int i, ii,jj;
+  intvec V;                          // to be used for variable weights 
+  int y = printlevel-voice+2;
+  def R = basering;
+  poly c = J[1];                     // the denominator
+  list gnirlist = ringlist(basering);
+  string svars = varstr(basering);
+  int nva = nvars(basering);
+  string svar;
+  ideal maxid = maxideal(1);
+  
+  int noRed = 0;     // By default, we try to reduce the number of generators.
+  if(size(#) > 0){
+    if ( typeof(#[1]) == "string" ) 
+    {
+      if (#[1] == "noRed"){noRed = 1;}
     }
-
-  // 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 list preimages) to this ring.
-
-  setring Rees;
-  return (fetch(R(1),preimages));
+  }
+
+  if ( y >= 1){"// computing the ring structure...";}
+   
+  if(c==1)
+  {
+    if( defined(norid) )  { kill norid; }
+    if( defined(normap) ) { kill normap; }
+    ideal norid = I;
+    ideal normap =  maxid;
+    export norid;
+    export normap;
+    if(noRed == 1){
+      setring R;
+      return(R);
+    } else {
+      list L = substpartSpecial(norid,normap);
+      def lastRing = L[1];
+      setring R;
+      return(lastRing);
+    }
+  }
+
+
+  //-------------- Enlarge ring by creating new variables ------------------
+  //check first whether variables T(i) and then whether Z(i),...,A(i) exist
+  //old variable names are not touched
+
+  if ( find(svars,"T(") == 0 )        
+  {
+    svar = "T";
+  }
+  else
+  {
+    for (ii=90; ii>=65; ii--)  
+    {
+      if ( find(svars,ASCII(ii)+"(") == 0 )        
+      {
+        svar = ASCII(ii);  break;
+      }
+    }
+  }
+
+  int q = size(J)-1;
+  if ( size(svar) != 0 )
+  {
+    for ( ii=q; ii>=1; ii-- )
+    {
+      gnirlist[2] = insert(gnirlist[2],svar+"("+string(ii)+")");
+    }
+  }
+  else
+  {       
+    for ( ii=q; ii>=1; ii-- )
+    {
+      gnirlist[2] = insert(gnirlist[2],"T("+string(100*nva+ii)+")");
+    }
+  }
+
+  V[q]=0;                        //create intvec of variable weights
+  V=V+1;
+  gnirlist[3] = insert(gnirlist[3],list("dp",V));
+
+  //this is a block ordering with one dp-block (1st block) for new vars
+  //the remaining weights and blocks for old vars are kept
+  //### perhaps better to make only one block, keeping weights ?
+  //this might effect syz below
+  //alt: ring newR = char(R),(X(1..nvars(R)),T(1..q)),dp;
+  //Reihenfolge geŠndert:neue Variablen kommen zuerst, Namen ev. nicht T(i)
+
+  def newR = ring(gnirlist);   
+  setring newR;                //new extended ring
+  ideal I = imap(R,I);
+
+  //------------- Compute linear and quadratic relations ---------------
+  if(y>=1)
+  {
+     "// compute linear relations:"; 
+  }
+  qring newQ = std(I);
+  
+  ideal f = imap(R,J);
+  module syzf = syz(f);
+  ideal pf = f[1]*f;               
+  //f[1] is the denominator D from normalityTest, a non zero divisor of R/I
+
+  ideal newT = maxideal(1);
+  newT = 1,newT[1..q];
+  //matrix T = matrix(ideal(1,T(1..q)),1,q+1);   //alt
+  matrix T = matrix(newT,1,q+1); 
+  ideal Lin = ideal(T*syzf);
+  //Lin=interred(Lin); 
+  //### interred reduziert ev size aber size(string(LIN)) wird groesser
+
+  if(y>=1)
+  {
+    if(y>=3)
+    {
+      "//   the linear relations:";  Lin; pause();"";
+    }
+      "// the ring structure of the normalization as affine algebra";
+      "//   number of linear relations:", size(Lin);
+  }
+
+  if(y>=1)
+  {
+    "// compute quadratic relations:"; 
+  }
+  matrix A;
+  ideal Quad;
+  poly ff;
+  newT = newT[2..size(newT)];
+  matrix u;  // The units for non-global orderings.
+
+  // Quadratic relations
+  for (ii=2; ii<=q+1; ii++ )
+  {
+    for ( jj=2; jj<=ii; jj++ )
+    {
+      ff = NF(f[ii]*f[jj],std(0));     // this makes lift much faster
+      // For non-global orderings, we have to take care of the units.
+      if(ord_test(basering) != 1){
+        A = lift(pf, ff, u);
+        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] * u[1, 1]- T*A);  
+      } else {
+        A = lift(pf,ff);              // ff lin. comb. of elts of pf mod I
+        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] - T*A);
+      }
+      //A = lift(pf, f[ii]*f[jj]);
+      //Quad = Quad, ideal(T(jj-1)*T(ii-1) - T*A);  
+    }
+  }
+  Quad = Quad[2..ncols(Quad)];
+
+  if(y>=1)
+  {
+    if(y>=3)
+    {
+      "//   the quadratic relations"; Quad; pause();"";
+    }
+      "//   number of quadratic relations:", size(Quad);
+  }
+  ideal Q1 = Lin,Quad;     //elements of Q1 are in NF w.r.t. I
+
+  //Q1 = mstd(Q1)[2];
+  //### weglassen, ist sehr zeitaufwendig.
+  //Ebenso interred, z.B. bei Leonard1 (1. Komponente von Leonard):
+  //"size Q1:", size(Q1), size(string(Q1));   //75 60083
+  //Q1 = interred(Q1); 
+  //"size Q1:", size(Q1), size(string(Q1));   //73 231956 (!)
+  //### Speicherueberlauf bei substpartSpecial bei 'ideal norid  = phi1(endid)'
+  //Beispiel fuer Hans um map zu testen!
+
+  setring newR;
+  ideal endid  = imap(newQ,Q1),I;
+  ideal endphi = imap(R,maxid);
+
+  if(noRed == 0){   
+    list L=substpartSpecial(endid,endphi);
+    def lastRing=L[1];
+    if(y>=1)
+    {
+      "//   number of substituted variables:", nvars(newR)-nvars(lastRing);
+      pause();"";
+    }
+    return(lastRing);
+  } else {
+    ideal norid = endid;
+    ideal normap = endphi;
+    export(norid);
+    export(normap);
+    setring R;
+    return(newR);
+  }   
 }
+
+//                Up to here: procedures for char p with Frobenius
+///////////////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////////////
+//                From here: subprocedures for normal
+
+static proc normalM(ideal I, int decomp, int withDelta){
+// Computes the normalization of a ring R / I using the module structure as far
+// as possible.
+// The ring R is the basering.
+// Input: ideal I
+// Output: a list of 3 elements (resp 4 if withDelta = 1), say nor.
+// - nor[1] = U, an ideal of R.
+// - nor[2] = c, an element of R.
+// U and c satisfy that 1/c * U is the normalization of R / I in the 
+// quotient field Q(R/I).
+// - nor[3] = ring say T, containing two ideals norid and normap such that
+// normap gives the normalization map from R / I to T / norid.
+// - nor[4] = the delta invariant, if withDelta = 1.
+
+// Details:
+// --------
+// Computes the ideal of the minors in the first step and then reuses it in all 
+// steps.
+// In step s, the denominator is D^s, where D is a nzd of the original quotient 
+// ring, contained in the radical of the singular locus.
+// This denominator is used except when the degree of D^i is greater than the 
+// degree of a conductor.
+// The nzd is taken as the smallest degree polynomial in the radical of the 
+// singular locus.
+
+// It assumes that the ideal I is equidimensional radical. This is not checked 
+// in the procedure!
+// If decomp = 0 or decomp = 3 it assumes further that I is prime. Therefore 
+// any non-zero element in the jacobian ideal is assumed to be a 
+// non-zerodivisor.
+
+// It works over the given basering.
+// If it has a non-global ordering, it changes it to dp global only for 
+// computing radical.
+
+// The delta invariant is only computed if withDelta = 1, and decomp = 0 or 
+// decomp = 3 (meaning that the ideal is prime).
+
+  option("redSB");
+  option("returnSB");
+  int step = 0;                       // Number of steps. (for debugging)
+  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
+  int i;                              // counter
+  int isGlobal = ord_test(basering);
+
+  poly c;                     // The common denominator.
+
+  def R = basering;
+    
+  // ---------------- computation of the singular locus --------------------- 
+  // We compute the radical of the ideal of minors modulo the original ideal.
+  // This is done only in the first step.
+
+  if(isGlobal == 1){
+    list IM = mstd(I);
+    I = IM[1];
+    ideal IMin = IM[2];   // A minimal set of generators in the groebner basis.
+  } else {
+    // The minimal set of generators is not computed by mstd for
+    // non-global orderings.
+    I = groebner(I);
+    ideal IMin = I;
+  }
+  int d = dim(I);
+  
+  qring Q = I;   // We work in the quotient by the groebner base of the ideal I
+  option("redSB");
+  option("returnSB");
+  ideal I = fetch(R, I);
+  attrib(I, "isSB", 1);
+  ideal IMin = fetch(R, IMin);
+
+  dbprint(dbg, "Computing the jacobian ideal...");
+  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
+  J = groebner(J);
+  
+  // -------------------- election of the conductor -------------------------   
+  poly condu = getSmallest(J);   // Choses the polynomial of smallest degree 
+                                 // of J as universal denominator.
+  if(dbg >= 1){
+    "";
+    "The universal denominator is ", condu;
+  }
+
+  // ---------  splitting the ideal by the conductor (if possible) ----------- 
+  // If the ideal is equidimensional, but not necessarily prime, we check if 
+  // the conductor is a non-zerodivisor of R/I.
+  // If not, we split I.
+  if((decomp == 1) or (decomp == 2)){  
+    ideal Id1 = quotient(0, condu);
+    if(size(Id1) > 0){
+      // We have to split.
+      if(dbg >= 1){
+        "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu;
+      }
+      setring R;
+      ideal Id1 = fetch(Q, Id1), I;
+      Id1 = groebner(Id1);
+      ideal Id2 = quotient(I, Id1);
+      // I = I1 \cap I2
+      printlevel = printlevel + 1;
+      list nor1 = normalM(Id1, decomp, withDelta)[1];
+      list nor2 = normalM(Id2, decomp, withDelta)[1];
+      printlevel = printlevel - 1;    
+      return(list(nor1, nor2));
+    }
+  }
+
+  // --------------- computation of the first test ideal ---------------------   
+  // To compute the radical we go back to the original ring.
+  // If we are using a non-global ordering, we must change to the global 
+  // ordering.
+  if(isGlobal == 1){
+    setring R;
+    ideal J = fetch(Q, J);
+    J = J, I;
+    if(dbg >= 1){
+      "The original singular locus is";
+      groebner(J);
+      if(dbg >= 2){pause();}
+       "";
+    }
+    J = radical(J);
+  } else {
+    // We change to global dp ordering.
+    list rl = ringlist(R);
+    list origOrd = rl[3];
+    list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0);
+    rl[3] = newOrd;
+    def globR = ring(rl);
+    setring globR;
+    ideal J = fetch(Q, J);
+    ideal I = fetch(R, I);
+    J = J, I;
+    if(dbg >= 1){
+      "The original singular locus is";
+      groebner(J);
+      if(dbg>=2){pause();}
+      "";
+    }
+    J = radical(J);
+    setring R;
+    ideal J = fetch(globR, J);
+  }
+    
+  if(dbg >= 1){
+    "The radical of the original singular locus is";
+    J;
+    if(dbg>=2){pause();}
+  }
+
+  // ---------------- election of the non zero divisor --------------------- 
+  setring Q;
+  J = fetch(R, J);
+  J = interred(J);          
+  poly D = getSmallest(J);    // Chooses the poly of smallest degree as 
+                              // non-zerodivisor.
+  if(dbg >= 1){
+    "The non zero divisor is ", D;
+    "";
+  }
+
+  // ------- splitting the ideal by the non-zerodivisor (if possible) --------   
+  // If the ideal is equidimensional, but not necessarily prime, we check if D
+  // is actually a non-zerodivisor of R/I.
+  // If not, we split I.
+  if((decomp == 1) or (decomp == 2)){  
+    // We check if D is actually a non-zerodivisor of R/I.
+    // If not, we split I.
+    Id1 = quotient(0, D);
+    if(size(Id1) > 0){
+      // We have to split.
+      if(dbg >= 1){
+        "A zerodivisor was found. We split the ideal. The zerodivisor is ", D;
+      }    
+      setring R;
+      ideal Id1 = fetch(Q, Id1), I;
+      Id1 = groebner(Id1);
+      ideal I2 = quotient(I, Id1);
+      // I = Id1 \cap Id2
+      printlevel = printlevel + 1;
+      list nor1 = normalM(Id1, decomp, withDelta)[1];
+      list nor2 = normalM(Id2, decomp, withDelta)[1];
+      printlevel = printlevel - 1;    
+      return(list(nor1, nor2));    
+    }
+  }
+
+  // --------------------- normalization ------------------------------------ 
+  // We call normalMEqui to compute the normalization.
+  setring R;
+  poly D = fetch(Q, D);
+  poly condu = fetch(Q, condu);
+  J = fetch(Q, J);
+  printlevel = printlevel + 1;
+  list result = normalMEqui(I, J, condu, D, withDelta);
+  printlevel = printlevel - 1;
+  return(list(result));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+  
+static proc normalMEqui(ideal I, ideal origJ, poly condu, poly D, int withDelta)
+// Here is where the normalization is actually computed.
+
+// Computes the normalization of R/I. (basering is R)
+// I is assumed to be radical and equidimensional.
+// origJ is the first test ideal.
+// D is a non-zerodivisor of R/I.
+// condu is a non-zerodivisor in the conductor.
+// If withDelta = 1, computes the delta invariant.
+{
+  int step = 0;                       // Number of steps. (for debugging)
+  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
+  int i;                              // counter
+  int isNormal = 0;                   // check for exiting the loop
+  int isGlobal = ord_test(basering);
+  int delt;
+  
+  def R = basering;
+  poly c = D;
+  ideal U;
+  ideal cJ;
+  list testOut;                 // Output of proc testIdeal
+                                // (the test ideal and the ring structure)
+    
+  qring Q = groebner(I);
+  option("redSB");
+  option("returnSB");  
+  ideal J = imap(R, origJ);
+  poly c = imap(R, c);
+  poly D = imap(R, D);
+  poly condu = imap(R, condu);
+  ideal cJ;
+  ideal cJMod;
+
+  dbprint(dbg, "Preliminar step begins.");
+
+  // --------------------- computation of A1 ------------------------------- 
+  dbprint(dbg, "Computing the quotient (DJ : J)...");
+  ideal U = groebner(quotient(D*J, J));  
+  ideal oldU = 1;
+
+  if(dbg >= 2){
+    "The quotient is"; U;
+  }
+
+  // ----------------- Grauer-Remmert criterion check ----------------------- 
+  // We check if the equality in Grauert - Remmert criterion is satisfied.
+  isNormal = checkInclusions(D*oldU, U);
+  if(isNormal == 0){
+    if(dbg >= 1){
+      "In this step, we have the ring 1/c * U, with c =", c;
+      "and U = "; U;
+    }
+  } else {        
+    // The original ring R/I was normal. Nothing to do.
+    // We define anyway a new ring, equal to R, to be able to return it.
+    setring R;
+    list lR = ringlist(R);
+    def ROut = ring(lR);
+    setring ROut;
+    ideal norid = 0;
+    ideal normap = maxideal(1);
+    export norid;
+    export normap;
+    setring R;
+    if(withDelta){
+      list output = ideal(1), poly(1), ROut, 0;
+    } else {
+      list output = ideal(1), poly(1), ROut;
+    }
+    return(output);
+  }
+
+  // ----- computation of the chain of ideals A1 c A2 c ... c An ------------ 
+  while(isNormal == 0){
+    step++;
+    if(dbg >= 1){
+      "";
+      "Step ", step, " begins.";
+    }
+    dbprint(dbg, "Computing the test ideal...");
+    
+    // --------------- computation of the test ideal ------------------------ 
+    // Computes a test ideal for the new ring. 
+    // The test ideal will be the radical in the new ring of the original
+    // test ideal.
+    setring R;
+    U = imap(Q, U);
+    c = imap(Q, c);
+    testOut = testIdeal(I, U, origJ, c, D);
+    cJ = testOut[1];
+    
+    setring Q;
+    cJ = imap(R, cJ);
+    cJ = groebner(cJ);
+    
+    // cJ / c is now the ideal mapped back.
+    // We have the generators as an ideal in the new ring, 
+    // but we want the generators as an ideal in the original ring.
+    cJMod = getGenerators(cJ, U, c);
+
+    if(dbg >= 2){
+      "The test ideal in this step is ";
+      cJMod;
+    }
+    
+    cJ = cJMod;
+
+    // ------------- computation of the quotient (DJ : J)-------------------- 
+    oldU = U;
+    dbprint(dbg, "Computing the quotient (c*D*cJ : cJ)...");
+    U = quotient(c*D*cJ, cJ);
+    if(dbg >= 2){"The quotient is "; U;}
+
+    // ------------- Grauert - Remmert criterion check ---------------------- 
+    // We check if the equality in Grauert - Remmert criterion is satisfied.
+    isNormal = checkInclusions(D*oldU, U);
+
+    if(isNormal == 1){        
+      // We go one step back. In the last step we didnt get antyhing new, 
+      // we just verified that the ring was already normal.
+      dbprint(dbg, "The ring in the previous step was already normal.");
+      dbprint(dbg, "");
+      U = oldU;
+    } else {
+      // ------------- preparation for next iteration ---------------------- 
+      // We have to go on.
+      // The new denominator is chosen.
+      c = D * c;
+      if(deg(c) > deg(condu)){
+        U = changeDenominatorQ(U, c, condu);
+        c = condu;
+      }
+      if(dbg >= 1){
+        "In this step, we have the ring 1/c * U, with c =", c;
+        "and U = ";
+        U;
+        if(dbg>=2){pause();}
+      }
+    }
+  }
+
+  // ------------------------- delta computation ---------------------------- 
+  if(withDelta){ 
+    ideal UD = groebner(U);
+    delt = vdim(std(modulo(UD, c)));
+  }
+  
+  // -------------------------- prepare output -----------------------------   
+  setring R;
+  U = fetch(Q, U);
+  c = fetch(Q, c);
+
+  // Ring structure of the new ring
+  def ere = testOut[2];
+  if(withDelta){
+    list output = U, c, ere, delt;
+  } else {
+    list output = U, c, ere;
+  }
+  return(output);
+ 
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc lineUp(ideal U, poly c)
+// Sets c as the first generator of U.
+{
+  int i;
+  ideal newU = c;
+  for (i = 1; i <= ncols(U); i++){
+    if(U[i] != c){
+      newU = newU, U[i];
+    }
+  }
+  return(newU);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc getSmallest(ideal J){
+// Computes the polynomial of smallest degree of J.
+// If there are more than one, it chooses the one with smallest number 
+// of monomials.
+  int i;
+  poly p = J[1];
+  int d = deg(p);
+  int di;
+  for(i = 2; i <= ncols(J); i++){
+    if(J[i] != 0){
+      di = deg(J[i]);
+      if(di < d){
+        p = J[i];
+        d = di;
+      } else {
+        if(di == d){
+          if(size(J[i]) < size(p)){
+            p = J[i];
+          }
+        }
+      }
+    }
+  }
+  return(p);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc getGenerators(ideal J, ideal U, poly c){
+// Computes the generators of J as an ideal in the original ring,
+// where J is given by generators in the new ring.
+
+// The new ring is given by 1/c * U in the total ring of fractions.
+
+  int i, j;                             // counters;
+  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)  
+  poly p;                               // The lifted polynomial
+  ideal JGr = groebner(J);              // Groebner base of J
+
+  if(dbg>1){"Checking for new generators...";}
+  for(i = 1; i <= ncols(J); i++){
+    for(j = 1; j <= ncols(U); j++){
+      p = lift(c, J[i]*U[j])[1,1];
+      p = reduce(p, JGr);
+      if(p != 0){
+        if(dbg>1){
+          "New polynoial added:", p;
+          if(dbg>4){pause();}
+        }
+        JGr = JGr, p;
+        JGr = groebner(JGr);
+        J = J, p;
+      }
+    }
+  }
+  return(J);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D){
+// Internal procedure, used in normalM.
+// Computes the test ideal in the new ring.
+// It takes the original test ideal and computes the radical of it in the 
+// new ring.
+
+// The new ring is 1/c * U.
+// The original test ideal is origJ.
+// The original ring is R / I, where R is the basering.
+  int i;                                // counter
+  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)  
+  def R = basering;                      // We dont work in the quo
+  ideal J = origJ;
+  
+  // ---------- computation of the ring structure of 1/c * U ----------------
+  U = lineUp(U, c);
+  
+  if(dbg > 1){"Computing the new ring structure...";}
+  list ele = computeRing(U, I, "noRed");
+
+  def origEre = ele[1];
+  setring origEre;
+  if(dbg > 1){"The relations are"; norid;}
+
+  // ---------------- setting the ring to work in  --------------------------  
+  int isGlobal = ord_test(origEre);      // Checks if the original ring has 
+                                         // global ordering.
+  if(isGlobal != 1){
+    list rl = ringlist(origEre);
+    list origOrd = rl[3];
+    list newOrd = list("dp", intvec(1:nvars(origEre))), list("C", 0);
+    rl[3] = newOrd;
+    def ere = ring(rl);     // globR is the original ring but 
+                            // with a global ordering.
+    setring ere;
+    ideal norid = imap(origEre, norid);
+  } else {
+    def ere = origEre;
+  }
+
+  ideal I = imap(R, I);
+  ideal J = imap(R, J);
+  J = J, norid, I;
+
+
+  // ----- computation of the test ideal using the ring structure of Ai -----  
+  option("redSB");
+  option("returnSB");
+
+  if(dbg > 1){"Computing the radical of J...";}
+  J = radical(J);
+  if(dbg > 1){"Computing the interreduction of the radical...";}
+  J = groebner(J);  
+  //J = interred(J);
+  if(dbg > 1){
+    "The radical in the generated ring is";
+    J;
+    if(dbg>4){pause();}
+  }
+  
+  setring ere;
+  
+  // -------------- map from Ai to the total ring of fractions ---------------  
+  // Now we must map back this ideal J to U_i / c in the total ring of
+  // fractions.
+  // The map sends T_j -> u_j / c.
+  // The map is built by the following steps:
+  // 1) We compute the degree of the generators of J with respect to the 
+  //    new variables T_j.
+  // 2) For each generator, we multiply each term by a power of c, as if
+  //    taking c^n as a common denominator (considering the new variables as
+  //    a polynomial in the old variables divided by c).
+  // 3) We replace the new variables T_j by the corresponding numerator u_j.
+  // 4) We lift the resulting polynomial to change the denominator 
+  //    from c^n to c.
+  int nNewVars = nvars(ere) - nvars(R);      // Number of new variables
+  poly c = imap(R, c);
+  intvec @v = 1..nNewVars;    // Vector of the new variables. 
+                              // They must be the first ones.
+  if(dbg > 1){"The indices of the new variables are", @v;}
+
+  // ---------------------- step 1 of the mapping ---------------------------  
+  intvec degs;
+  for(i = 1; i<=ncols(J); i++){
+    degs[i] = degSubring(J[i], @v);
+  }
+  if(dbg > 1){
+    "The degrees with respect to the new variables are";
+    degs;
+  }
+  
+  // ---------------------- step 2 of the mapping ---------------------------    
+  ideal mapJ = mapBackIdeal(J, c, @v);
+  
+  setring R;
+
+  // ---------------------- step 3 of the mapping ---------------------------  
+  ideal z;                    // The variables of the original ring in order.
+  for(i = 1; i<=nvars(R); i++){
+    z[i] = var(i);
+  }
+      
+  map f = ere, U[2..ncols(U)], z[1..ncols(z)]; // The map to the original ring.
+  if(dbg > 1){
+    "The map is ";
+    f;
+    if(dbg>4){pause();}
+  }
+  
+  if(dbg > 1){
+    "Computing the map...";
+  }
+      
+  J = f(mapJ);
+  if(dbg > 1){
+    "The ideal J mapped back (before lifting) is";
+    J;
+    if(dbg>4){pause();}
+  }
+
+  // ---------------------- step 4 of the mapping ---------------------------    
+  qring Q = groebner(I);
+  ideal J = imap(R, J);
+  poly c = imap(R, c);
+  for(i = 1; i<=ncols(J); i++){
+    if(degs[i]>1){
+      J[i] = lift(c^(degs[i]-1), J[i])[1,1];
+    } else {
+      if(degs[i]==0){
+        J[i] = c*J[i];
+      }
+    }
+  }
+
+  if(dbg > 1){
+    "The ideal J lifted is";
+    J;
+    if(dbg>4){pause();}
+  }
+
+  // --------------------------- prepare output ----------------------------    
+  J = groebner(J);
+  
+  setring R;
+  J = imap(Q, J);
+  
+  return(list(J, ele[1]));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I){
+// Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the 
+// ring is 1/c2 * U2.
+// The base ring is R, but the computations are to be done in R / I.
+  int a;      // counter
+  def R = basering;
+  qring Q = I;
+  ideal U1 = fetch(R, U1);
+  poly c1 = fetch(R, c1);
+  poly c2 = fetch(R, c2);
+  ideal U2 = changeDenominatorQ(U1, c1, c2);
+  setring R;
+  ideal U2 = fetch(Q, U2);
+  return(U2);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc changeDenominatorQ(ideal U1, poly c1, poly c2){
+// Given a ring in the form 1/c1 * U, it computes a new U2 st the ring 
+// is 1/c2 * U2.
+// The base ring is already a quotient ring R / I.
+  int a;      // counter
+  ideal U2;
+  poly p;
+  for(a = 1; a <= ncols(U1); a++){
+    p = lift(c1, c2*U1[a])[1,1];
+    U2[a] = p;
+  }
+  return(U2);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc checkInclusions(ideal U1, ideal U2){
+// Checks if the identity A = Hom(J, J) of Grauert-Remmert criterion is 
+// satisfied.
+  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
+  list reduction1;
+  list reduction2;
+  
+  // ---------------------- inclusion Hom(J, J) c A -------------------------  
+  if(dbg > 1){"Checking the inclusion Hom(J, J) c A:";}
+  // This interred is used only because a bug in groebner!
+  U1 = groebner(U1);
+  reduction1 = reduce(U2, U1);
+  if(dbg > 1){reduction1[1];}
+
+  // ---------------------- inclusion A c Hom(J, J) -------------------------  
+  // The following check should always be satisfied.
+  // This is only used for debugging.
+  if(dbg > 1){
+    "and the inclusion A c Hom(J, J): (this should always be satisfied)";
+    // This interred is used only because a bug in groebner!
+    U2 = groebner(U2);
+    reduction2 = reduce(U1, groebner(U2));
+    reduction2[1];
+    if(size(reduction2[1]) > 0){
+      "Something went wrong... (this inclusion should always be satisfied)";
+      ~;
+    } else {
+      if(dbg>4){pause();}
+    }
+  }
+  
+  if(size(reduction1[1]) == 0){
+    // We are done! The ring computed in the last step was normal.
+    return(1);
+  } else {
+    return(0);
+  }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc degSubring(poly p, intvec @v){
+// Computes the degree of a polynomial taking only some variables as variables
+// and the others as parameters.
+
+// The degree is taken wrt the variables indicated in v.
+  int i;      // Counter
+  int d = 0;  // The degree
+  int e;      // Degree (auxiliar variable)
+  for(i = 1; i <= size(p); i++){
+    e = sum(leadexp(p[i]), @v);
+    if(e > d){d = e;}
+  }
+  return(d);
+}
+    
+///////////////////////////////////////////////////////////////////////////////
+
+static proc mapBackIdeal(ideal I, poly c, intvec @v){
+// Modifies all polynomials in I so that a map x(i) -> y(i)/c can be 
+// carried out.
+
+// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.
+
+  int i;  // counter
+  for(i = 1; i <= ncols(I); i++){
+    I[i] = mapBackPoly(I[i], c, @v);
+  }
+  return(I);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc mapBackPoly(poly p, poly c, intvec @v){
+// Multiplies each monomial of p by a power of c so that a map x(i) -> y(i)/c 
+// can be carried out.
+
+// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.
+  int i;  // counter
+  int e;  // exponent
+  int d = degSubring(p, @v);
+  poly g = 0;
+  for(i = 1; i <= size(p); i++){
+    e = sum(leadexp(p[i]), @v);
+    g = g + p[i] * c^(d-e);
+  }
+  return(g);
+}
+
+//                Up to here: procedures for normal
+///////////////////////////////////////////////////////////////////////////////
+
+
+///////////////////////////////////////////////////////////////////////////////
+//                From here: procedures for normalC
+
+// We first define resp. copy some attributes to be used in proc normal and
+// static proc normalizationPrimes, and ..., to speed up computation in
+// special cases
+//NOTE:  We use the following attributes:
+// 1     attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
+// 2     attrib(id,"isCompleteIntersection");  //--- complete intersection  
+// 3     attrib(id,"isHypersurface");          //--- hypersurface 
+// 4     attrib(id,"isEquidimensional");       //--- equidimensional ideal 
+// 5     attrib(id,"isPrim");                  //--- prime ideal 
+// 6     attrib(id,"isRegInCodim2");           //--- regular in codimension 2 
+// 7     attrib(id,"isIsolatedSingularity");   //--- isolated singularities 
+// 8     attrib(id,"onlySingularAtZero");      //--- only singular at 0 
+// 9     attrib(id,"isRadical");               //--- radical ideal 
+//Recall: (attrib(id,"xy"),1) sets attrib xy to TRUE and
+//        (attrib(id,"xy"),0) to FALSE
+
+static proc getAttrib (ideal id)
+"USAGE:   getAttrib(id);  id=ideal
+COMPUTE: check attributes for id. If the attributes above are defined, 
+         take its value, otherwise define it and set it to 0 
+RETURN:  intvec of size 9, with entries 0 or 1,  values of attributes defined
+         above (in this order)
+EXAMPLE: no example
+"
+{
+  int isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;
+
+  if( typeof(attrib(id,"isCohenMacaulay"))=="int" )   
+  {
+    if( attrib(id,"isCohenMacaulay")==1 ) 
+    { isCoM=1; isEq=1; }
+  }
+
+  if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
+  {
+    if(attrib(id,"isCompleteIntersection")==1) 
+    { isCoI=1; isCoM=1; isEq=1; }
+  }
+    
+  if( typeof(attrib(id,"isHypersurface"))=="int" )
+  {
+    if(attrib(id,"isHypersurface")==1)
+    { isHy=1; isCoI=1; isCoM=1; isEq=1; }
+  }
+   
+  if( typeof(attrib(id,"isEquidimensional"))=="int" )
+  {
+    if(attrib(id,"isEquidimensional")==1) 
+    { isEq=1; }
+  }
+   
+  if( typeof(attrib(id,"isPrim"))=="int" )
+  {
+    if(attrib(id,"isPrim")==1)  
+    { isPr=1; }
+  }
+
+  if( typeof(attrib(id,"isRegInCodim2"))=="int" )
+  {
+    if(attrib(id,"isRegInCodim2")==1)  
+    { isReg=1; }
+  }
+   
+  if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
+  {
+    if(attrib(id,"isIsolatedSingularity")==1)  
+    { isIso=1; }
+  }
+ 
+  if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
+  {
+    if(attrib(id,"onlySingularAtZero")==1)
+    { oSAZ=1; }
+  }
+
+  if( typeof(attrib(id,"isRad"))=="int" )
+  {
+    if(attrib(id,"isRad")==1)  
+    { isRad=1; }
+  }
+
+  intvec atr = isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;
+  return(atr);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc setAttrib (ideal id, intvec atr)
+"USAGE:   setAttrib(id,atr);  id ideal, atr intvec
+COMPUTE: set attributes to id specified by atr
+RETURN:  id, with assigned attributes from atr
+EXAMPLE: no example
+"
+{
+  attrib(id,"isCohenMacaulay",atr[1]);         //--- Cohen Macaulay
+  attrib(id,"isCompleteIntersection",atr[2]);  //--- complete intersection  
+  attrib(id,"isHypersurface",atr[3]);          //--- hypersurface 
+  attrib(id,"isEquidimensional",atr[4]);       //--- equidimensional ideal 
+  attrib(id,"isPrim",atr[5]);                  //--- prime ideal 
+  attrib(id,"isRegInCodim2",atr[6]);           //--- regular in codimension 2 
+  attrib(id,"isIsolatedSingularity",atr[7]);   //--- isolated singularities 
+  attrib(id,"onlySingularAtZero",atr[8]);      //--- only singular at 0 
+  attrib(id,"isRadical",atr[9]);               //--- radical ideal 
+
+  return(id);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// copyAttribs is not used anywhere so far
+
+static proc copyAttribs (ideal id1, ideal id)
+"USAGE:   copyAttribs(id1,id);  id1, id ideals
+COMPUTE: copy attributes from id1 to id
+RETURN:  id, with assigned attributes from id1
+EXAMPLE: no example
+"
+{
+  if( typeof(attrib(id1,"isCohenMacaulay"))=="int" )   
+  {
+    if( attrib(id1,"isCohenMacaulay")==1 )
+    {
+      attrib(id,"isEquidimensional",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isCohenMacaulay",0);
+  }  
+
+  if( typeof(attrib(id1,"isCompleteIntersection"))=="int" )
+  {
+    if(attrib(id1,"isCompleteIntersection")==1)
+    {
+      attrib(id,"isCohenMacaulay",1);
+      attrib(id,"isEquidimensional",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isCompleteIntersection",0);
+  }
+    
+  if( typeof(attrib(id1,"isHypersurface"))=="int" )
+  {
+    if(attrib(id1,"isHypersurface")==1)
+    {
+      attrib(id,"isCompleteIntersection",1);
+      attrib(id,"isCohenMacaulay",1);
+      attrib(id,"isEquidimensional",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isHypersurface",0);
+  }
+   
+  if( (typeof(attrib(id1,"isEquidimensional"))=="int") )
+  {
+    if(attrib(id1,"isEquidimensional")==1)
+    {
+      attrib(id,"isEquidimensional",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isEquidimensional",0);
+  }
+
+  if( typeof(attrib(id1,"isPrim"))=="int" )
+  {
+    if(attrib(id1,"isPrim")==1)
+    {
+      attrib(id,"isEquidimensional",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isPrim",0);
+  }
+   
+  if( (typeof(attrib(id1,"isRegInCodim2"))=="int") )
+  {
+    if(attrib(id1,"isRegInCodim2")==1)
+    {
+      attrib(id,"isRegInCodim2",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isRegInCodim2",0);
+  }
+    
+  if( (typeof(attrib(id1,"isIsolatedSingularity"))=="int") )
+  {
+    if(attrib(id1,"isIsolatedSingularity")==1)
+    {
+      attrib(id,"isIsolatedSingularity",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isIsolatedSingularity",0);
+  }
+ 
+  if( typeof(attrib(id1,"onlySingularAtZero"))=="int" )
+  {
+    if(attrib(id1,"onlySingularAtZero")==1)
+    {
+      attrib(id,"isIsolatedSingularity",1);
+    }
+  }
+  else
+  {
+    attrib(id,"onlySingularAtZero",0);
+  }
+
+  if( typeof(attrib(id1,"isRad"))=="int" )
+  {
+    if(attrib(id1,"isRad")==1)  
+    {
+      attrib(id,"isRad",1);
+    }
+  }
+  else
+  {
+    attrib(id,"isRad",0);
+  }
+
+  return(id);
+}
+
+proc normalC(ideal id, list #)
+"USAGE:   normalC(id [,choose]);  id = radical ideal, choose = optional string 
+         which may be a combination of \"equidim\", \"prim\", \"withGens\",
+         \"isPrim\", \"noFac\".@*
+         If choose is not given or empty, a default method is choosen@*
+         If choose = \"prim\" resp. \"equidim\" normalC computes a prime resp.
+         an equidimensional decomposition and then the normalization of each
+         component; if choose = \"noFac\" factorization is avoided in the 
+         computation of the minimal associated primes; @*
+         if choose = \"withGens\" the minimal associated primes P_i of id are 
+         computed and for each P_i, algebra generators of the integral closure 
+         of basering/P_i are computed as elements of its quotient field;@*
+         \"isPrim\" disables \"equidim\" and \"prim\".@*
+ASSUME:  The ideal must be radical, for non-radical ideals the output may
+         be wrong (id=radical(id); makes id radical). However, if choose = 
+         \"prim\" the minimal associated primes of id are computed first 
+         and hence normalC computes the normalization of the radical of id.
+         \"isPrim\" should only be used if id is known to be irreducible.
+RETURN:  a list, say nor, of size 2 (resp. 3 if choose=\"withGens\"). 
+         The first list nor[1] consists of r rings, where r is the number 
+         of irreducible components if choose = \"prim\" (resp. >= no of 
+         equidimenensional components if choose =  \"equidim\"). 
+@format 
+         Each ring Ri=nor[1][i], i=1..r, contains two ideals with given
+         names @code{norid} and @code{normap} such that @*
+         - Ri/norid is the normalization of the i-th component, i.e. the 
+          integral closure in its field of fractions (as affine ring);
+         - the direct sum of the rings Ri/norid is the normalization
+           of basering/id; @*
+         - @code{normap} gives the normalization map from basering/id to
+           Ri/norid for each j.@*
+
+         If choose=\"withGens\", nor[2] is a list of size r of ideals Ii= 
+         nor[2][i] in the basering, generating the integral closure of
+         basering/P_i in its quotient field, i=1..r, in two different ways:
+         - Ii is given by polynomials 
+           g_1,...,g_k,g_k+1 such that the variables T(j) of the ring Ri 
+           satisfy T(1)=g_1/g_k+1,..,T(k)=g_k/g_k+1. Hence the g_j/g_k+1
+           are algebra generators of the integral closure of basering/P_i
+           as sub-algebra in the quotient field of basering/P_i. 
+ 
+         nor[2] (resp. nor[3] if choose=\"withGens\") is a list of an intvec 
+         of size r, the delta invariants of the r components, and an integer, 
+         the total delta invariant of basering/id (-1 means infinite, and 0 
+         that basering/P_i resp. basering/input is normal). 
+         Return value -2 means that delta resp. delta of one of the components
+         is not computed (which may happen if \"equidim\" is given) .
+@end format
+THEORY:  The delta invariant of a reduced ring K[x1,...,xn]/id is 
+             dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id)
+         We call this number also the delta invariant of id.
+         We use the algorithm described in [G.-M. Greuel, G. Pfister: 
+         A SINGULAR Introduction to Commutative Algebra, 2nd Edition. 
+         Springer Verlag (2007)].
+NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
+@*       Increasing/decreasing printlevel displays more/less comments 
+         (default: printlevel=0).
+@*       Not implemented for local or mixed orderings and for quotient rings.
+@*       If the input ideal id is weighted homogeneous a weighted ordering may
+         be used (qhweight(id); computes weights).
+KEYWORDS: normalization; integral closure; delta invariant.
+EXAMPLE: example normalC; shows an example
+"
+{  
+   int i,j, wgens, withEqui, withPrim, isPrim, nfac;
+   int y = printlevel-voice+2;
+   int nvar = nvars(basering);
+   int chara  = char(basering);
+   list result,prim,keepresult;
+
+   if ( ord_test(basering) != 1 )
+   {
+     "";
+     "// Not implemented for this ordering,";
+     "// please change to global ordering!";
+     return(result);
+   }
+
+//--------------------------- define the method --------------------------- 
+   string method;                //make all options one string in order to use
+                                 //all combinations of options simultaneously
+   for ( i=1; i <= size(#); i++ )
+   { 
+     if ( typeof(#[i]) == "string" ) 
+     {
+       method = method + #[i];
+     }
+   }
+
+   //--------------------------- choosen methods -----------------------
+   //we have the methods 
+   // "withGens": computes algebra generators for each irreducible component
+   // the general case: either "equidim" or not (then applying minAssGTZ)
+
+   if ( find(method,"withgens") or find(method,"withGens"))
+   {
+     wgens=1;
+   }
+   if ( find(method,"equidim") )
+   {
+     withEqui=1;
+   }
+   if ( find(method,"prim") )
+   {
+     withPrim=1;
+   }
+   if ( find(method,"isprim") or find(method,"isPrim") )
+   {
+     isPrim=1;
+   }
+   if ( find(method,"noFac") )
+   {
+     nfac=1;            
+   }
+
+   //------------------ default method, may be changed --------------------
+   //Use a prime decomposition only for small no of variables otherwise equidim 
+   //Compute delta whenever possible.
+
+   if ( ((chara==0 || chara>30) && nvar>2 && !withPrim && !isPrim) ||
+   ((chara==0 || chara>30) && size(#)  ==  0) )
+   {
+      withEqui=1;
+   }
+   kill #;
+   list #;
+
+//------- Special algorithm with computation of the generators, RETURN -------
+   //--------------------- method "withGens" ----------------------------------
+   //the integral closure is computed in proc primeClosure. In the general case
+   //it is computed in normalizationPrimes. The main difference is that  in
+   //primeClosure the singular locus is only computed in the first iteration
+   //that no attributes are used and that the generators are computed.
+   //In primeClosure the (algebra) generators for each irreducible component
+   //are computed in static proc closureGenerators
+
+   if( wgens )
+   {
+      if( y >= 1 )
+      {  "";
+         "// We use method 'withGens'";
+      }
+      if ( isPrim )
+      {
+         prim[1] = id;
+         if( y >= 0 )
+         { 
+           "";
+           "// ** WARNING: result is correct if ideal is prime (not checked) **";
+           "// if procedure is called with string \"prim\", primality is checked";
+         }
+      }
+      else
+      {
+         if(y>=1)
+         {  "// compute minimal associated primes"; }
+
+         if( nfac )
+         { prim = minAssGTZ(id,1); }
+         else
+         { prim = minAssGTZ(id); }
+
+         if(y>=1)
+         {  
+            prim;"";
+            "// number of irreducible components is", size(prim);
+         }
+      }
+   //----------- compute integral closure for every component -------------
+      int del;
+      intvec deli;
+      list Gens,l,resu,Resu;
+      ideal gens;
+      def R = basering;
+      poly gg;
+
+      for(i=1; i<=size(prim); i++)
+      {
+         if(y>=1)
+         {
+            ""; pause(); "";
+            "// start computation of component",i;
+            "   --------------------------------";
+         }
+
+         if( defined(ker) ) { kill ker; }
+         ideal ker = prim[i];
+         export(ker);
+         l = R;
+         l = primeClosure(l,1);               //here the work is done
+         //### ausprobieren ob primeClosure(l,1) schneller als primeClosure(l)
+         // 1 bedeutet: kŸrzester nzd
+
+         if ( l[size(l)] >= 0 && del >= 0 )
+         {
+            del = del + l[size(l)];
+         }
+         else
+         { del = -1; }
+         deli = l[size(l)],deli;
+
+         l = l[1..size(l)-1];
+         resu = list(l[size(l)]) + resu;
+         gens = closureGenerators(l);         //computes algebra(!) generators
+
+         //NOTE: gens[i]/gens[size(gens)] expresses the ith variable of resu[1] 
+         //(the normalization) as fraction of elements of the basering;
+         //the variables of resu[1] are algebra generators.
+         //gens[size(gens)] is a non-zero divisor of basering/i
+
+         //divide by the greatest common divisor:
+         gg = gcd( gens[1],gens[size(gens)] );
+         for(j=2; j<=size(gens)-1; j++)
+         {
+            gg=gcd(gg,gens[j]);
+         }
+         for(j=1; j<=size(gens); j++)
+         {
+            gens[j]=gens[j]/gg;
+         }
+         Gens = list(gens) + Gens;
+
+/*       ### Da die gens Algebra-Erzeuger sind, ist reduce nach Bestimmung 
+         der Algebra-Variablen T(i) nicht zulŠssig!
+         for(i=1;i<=size(gens)-1;i++) 
+         {
+            gens[i]= reduce(gens[i],std(gens[size(gens)]));
+         }
+         for(i=size(gens)-1; i>=1; i--)
+         {
+            if(gens[i]==0)
+            { gens = delete(gens,i); }
+         }
+*/
+         if( defined(ker) ) { kill ker; }
+      }
+
+      if ( del >= 0 )
+      {
+         int mul = iMult(prim);
+         del = del + mul;
+      }
+      else
+      { del = -1; }
+      deli = deli[1..size(deli)-1];
+      Resu = resu,Gens,list(deli,del); 
+      int sr = size(resu);
+
+      if ( y >= 0 )
+      {"";
+"// 'normalC' created a list, say nor, of three lists: 
+// nor[1] resp. nor[2] are lists of",sr,"ring(s) resp. ideal(s)
+// and nor[3] is a list of an intvec and an integer.
+// To see the result, type
+     nor;
+// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
+     def R1 = nor[1][1]; setring R1;  norid; normap;
+// for the other rings type first setring r; (if r is the name of your 
+// original basering) and then continue as for the first ring;
+// Ri/norid is the affine algebra of the normalization of the i-th 
+// component r/P_i (where P_i is an associated prime of the input ideal)
+// and normap the normalization map from r to Ri/norid;
+//    nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
+// elements g1..gk of r such that the gj/gk generate the integral 
+// closure of r/P_i as sub-algebra in the quotient field of r/P_i, with
+// gj/gk being mapped by normap to the j-th variable of Ri;
+//    nor[3] shows the delta-invariant of each component and of the input
+// ideal (-1 means infinite, and 0 that r/P_i resp. r/input is normal).";
+      }
+      return(Resu);
+   }
+
+//-------- The general case without computation of the generators -----------
+// (attrib(id,"xy"),1) sets attrib xy to TRUE and (attrib(id,"xy"),0) to FALSE
+// We use the following attributes:
+//   attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
+//   attrib(id,"isCompleteIntersection");  //--- complete intersection  
+//   attrib(id,"isHypersurface");          //--- hypersurface 
+//   attrib(id,"isEquidimensional",-1);    //--- equidimensional ideal 
+//   attrib(id,"isPrim");                  //--- prime ideal 
+//   attrib(id,"isRegInCodim2");           //--- regular in codimension 2 
+//   attrib(id,"isIsolatedSingularity";    //--- isolated singularities 
+//   attrib(id,"onlySingularAtZero");      //--- only singular at 0   
+
+ //------------------- first set the attributes ----------------------
+   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )   
+   {
+      if( attrib(id,"isCohenMacaulay")==1 )
+      {
+         attrib(id,"isEquidimensional",1);
+      }
+   }
+   else
+   {
+      attrib(id,"isCohenMacaulay",0);
+   }  
+
+   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
+   {
+      if(attrib(id,"isCompleteIntersection")==1)
+      {
+         attrib(id,"isCohenMacaulay",1);
+         attrib(id,"isEquidimensional",1);
+      }
+   }
+   else
+   {
+      attrib(id,"isCompleteIntersection",0);
+   }
+    
+   if( typeof(attrib(id,"isHypersurface"))=="int" )
+   {
+      if(attrib(id,"isHypersurface")==1)
+      {
+         attrib(id,"isCompleteIntersection",1);
+         attrib(id,"isCohenMacaulay",1);
+         attrib(id,"isEquidimensional",1);
+      }
+   }
+   else
+   {
+      attrib(id,"isHypersurface",0);
+   }
+   
+   if( ! (typeof(attrib(id,"isEquidimensional"))=="int") )
+   {
+         attrib(id,"isEquidimensional",0);
+   }
+   
+   if( typeof(attrib(id,"isPrim"))=="int" )
+   {
+      if(attrib(id,"isPrim")==1)
+      {
+         attrib(id,"isEquidimensional",1);
+      }
+   }
+   else
+   {
+      attrib(id,"isPrim",0);
+   }
+   
+   if( ! (typeof(attrib(id,"isRegInCodim2"))=="int") )
+   {
+         attrib(id,"isRegInCodim2",0);
+   }
+   
+   if( ! (typeof(attrib(id,"isIsolatedSingularity"))=="int") )
+   {
+         attrib(id,"isIsolatedSingularity",0);
+   }
+ 
+   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
+   {
+      if(attrib(id,"onlySingularAtZero")==1)
+      {
+         attrib(id,"isIsolatedSingularity",1);
+      }
+   }
+   else
+   {
+      attrib(id,"onlySingularAtZero",0);
+   }
+
+   //-------------- compute equidimensional decomposition --------------------
+   //If the method "equidim" is given, compute the equidim decomposition
+   //and goto the next step (no normalization 
+   //ACHTUNG: equidim berechnet bei nicht reduzierten id die eingebetteten
+   //Komponenten als niederdim Komponenten, waehrend diese bei primdecGTZ
+   //nicht auftauchen: ideal(x,y)*xy
+   //this is default for nvars > 2 and char = 0 or car > 30
+
+   if( withEqui )
+   {
+      withPrim = 0;                 //this is used to check later that prim
+                                    //contains equidim but not prime components
+      if( y >= 1 )
+      {
+         "// We use method 'equidim'";
+      }
+      if( typeof(attrib(id,"isEquidimensional"))=="int" )
+      {
+         if(attrib(id,"isEquidimensional")==1)
+         {
+            prim[1] = id;
+         }
+         else
+         {
+            prim = equidim(id);
+         }
+      }
+      else
+      {
+         prim = equidim(id);
+      }
+      if(y>=1)
+      {  "";
+         "// number of equidimensional components:", size(prim);
+      }
+      if ( !nfac )
+      {
+        intvec opt = option(get);
+        option(redSB);
+        for(j=1; j<=size(prim); j++)    
+        {
+           keepresult = keepresult+facstd(prim[j]);
+        }
+        prim = keepresult;
+        if ( size(prim) == 0 )
+        {
+          prim=ideal(0);     //Bug in facstd, liefert leere Liste bei 0-Ideal
+        }
+        
+        if(y>=1)
+        {  "";
+         "// number of components after application of facstd:", size(prim);
+        }
+        option(set,opt);
+      }
+   }
+
+   //------------------- compute associated primes -------------------------
+   //the case where withEqui = 0, here the min. ass. primes are computed
+   //start with the computation of the minimal associated primes:
+
+   else                                          
+   {
+    if( isPrim )
+    {
+      if( y >= 0 )
+      {
+         "// ** WARNING: result is correct if ideal is prime";
+         "// or equidimensional (not checked) **";
+         "// disable option \"isPrim\" to decompose ideal into prime";
+         "// or equidimensional components";"";
+      }
+      if( y >= 1 )
+      {
+        "// We use method 'isPrim'";"";
+      }
+      prim[1]=id;
+    }
+    else
+    {
+      withPrim = 1;                 //this is used to check later that prim
+                                    //contains prime but not equidim components
+      if( y >= 1 )
+      {
+         "// We use method 'prim'";
+      }
+
+      if( typeof(attrib(id,"isPrim"))=="int" )
+      {
+         if(attrib(id,"isPrim")==1)
+         {
+            prim[1]=id;
+         }
+         else
+         {
+            if( nfac )
+            { prim=minAssGTZ(id,1); }     //does not use factorizing groebner
+            else
+            { prim=minAssGTZ(id); }       //uses factorizing groebner
+         }
+      }
+      else
+      {
+            if( nfac )
+            { prim=minAssGTZ(id,1); }
+            else
+            { prim=minAssGTZ(id); }
+      }
+      if(y>=1)
+      {  "";
+         "// number of irreducible components:", size(prim);
+      }
+    }
+   }
+
+   //----- for each component (equidim or irred) compute normalization -----
+   int sr, skr, del;
+   intvec deli;  
+   int sp = size(prim);     //size of list prim (# irred or equidim comp)
+
+   for(i=1; i<=sp; i++)
+   {
+      if(y>=1)
+      {  "";
+         "// computing the normalization of component",i;
+         "   ----------------------------------------";
+      }
+      //-------------- first set attributes for components ------------------
+      attrib(prim[i],"isEquidimensional",1);
+      if( withPrim )
+      {
+         attrib(prim[i],"isPrim",1);
+      }
+      else
+      { attrib(prim[i],"isPrim",0); }
+
+      if(attrib(id,"onlySingularAtZero")==1)
+      { attrib(prim[i],"onlySingularAtZero",1); }
+      else
+      { attrib(prim[i],"onlySingularAtZero",0); }
+
+      if(attrib(id,"isIsolatedSingularity")==1)
+      { attrib(prim[i],"isIsolatedSingularity",1); }
+      else
+      { attrib(prim[i],"isIsolatedSingularity",0); }
+
+      if( attrib(id,"isHypersurface")==1 )
+      {
+         attrib(prim[i],"isHypersurface",1); 
+         attrib(prim[i],"isCompleteIntersection",1);
+         attrib(prim[i],"isCohenMacaulay",1);
+      }
+      else
+      { attrib(prim[i],"isHypersurface",0); }
+
+      if ( sp == 1)         //the case of one component: copy attribs from id
+      {
+        if(attrib(id,"isRegInCodim2")==1)
+        {attrib(prim[i],"isRegInCodim2",1); }
+        else
+        {attrib(prim[i],"isRegInCodim2",0); }
+
+        if(attrib(id,"isCohenMacaulay")==1)
+        {attrib(prim[i],"isCohenMacaulay",1); }
+        else
+        {attrib(prim[i],"isCohenMacaulay",0); }
+
+        if(attrib(id,"isCompleteIntersection")==1)
+        {attrib(prim[i],"isCompleteIntersection",1); }
+        else
+        {attrib(prim[i],"isCompleteIntersection",0); }
+      }
+      else
+      {  
+        attrib(prim[i],"isRegInCodim2",0);
+        attrib(prim[i],"isCohenMacaulay",0);
+        attrib(prim[i],"isCompleteIntersection",0);
+      }
+
+      //------ Now compute the normalization of each component ---------
+      //note: for equidimensional components the "splitting tools" can
+      //create further decomposition
+      //We now start normalizationPrimes with 
+      //ihp = partial normalisation map = identity map = maxideal(1)
+      //del = partial delta invariant = 0 
+      //deli= intvec of partial delta invariants of components
+      //in normalizationPrimes all the work is done:
+
+      keepresult = normalizationPrimes(prim[i],maxideal(1),0,0);
+
+      for(j=1; j<=size(keepresult)-1; j++)
+      {
+         result=insert(result,keepresult[j]);
+      }
+      skr = size(keepresult);
+
+      //compute delta:             
+      if( del >= 0 && keepresult[skr][1] >=0 )
+      {
+         del = del + keepresult[skr][1];
+      }
+      else
+      {  
+         del = -1; 
+      }
+      deli = keepresult[skr][2],deli;
+
+      if ( y>=1 )
+      {
+           "// delta of component",i; keepresult[skr][1];
+      } 
+   }
+   sr = size(result);
+
+   // -------------- Now compute intersection multiplicities -------------
+   //intersection multiplicities of list prim, sp=size(prim).
+      if ( y>=1 )
+      {
+        "// Sum of delta for all components"; del;
+        if ( sp>1 )
+        {
+           "// Compute intersection multiplicities of the components";
+        }
+      }   
+
+      if ( sp > 1 )
+      {
+        int mul = iMult(prim);
+        if ( mul < 0 )
+        {
+           del = -1;
+        }
+        else
+        {
+           del = del + mul; 
+        }
+      }
+   deli = deli[1..size(deli)-1];
+   result = result,list(deli,del);
+
+//--------------- Finally print comments and return ------------------
+   if ( y >= 0)
+   {"";
+"// 'normalC' created a list, say nor, of two lists:
+// nor[1] is a list of",sr,"ring(s), nor[2] a list of an intvec and an int.
+// To see the result, type
+      nor;
+// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
+      def R1 = nor[1][1];  setring R1;  norid;  normap;
+// and similair for the other rings nor[1][i];
+// Ri/norid is the affine algebra of the normalization of the i-th component 
+// r/P_i (where P_i is a prime or an equidimensional ideal of the input ideal) 
+// and normap the normalization map from the basering to Ri/norid;
+//    nor[2] shows the delta-invariant of each component and of the input
+// ideal (-1 means infinite, -2 that delta of a component was not computed).";
+   }
+   return(result);
+}
+
+example
+{ "EXAMPLE:";
+   printlevel = printlevel+1;
+   echo = 2;
+   ring s = 0,(x,y),dp;
+   ideal i = (x2-y3)*(x2+y2)*x;
+
+   list nor = normalC(i);
+
+   nor; 
+   // 2 branches have delta = 1, and 1 branch has delta = 0
+   // the total delta invariant is 13
+
+   def R2 = nor[1][2];  setring R2;
+   norid; normap;
+   
+   echo = 0;
+   printlevel = printlevel-1;
+   pause("   hit return to continue"); echo=2;
+
+   ring r = 2,(x,y,z),dp;
+   ideal i = z3-xy4;
+   nor = normalC(i);  nor;  
+   // the delta invariant is infinite 
+   // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
+   // in its quotient field Quot(r/i)
+
+   // the normalization as affine algebra over the ground field:             
+   def R = nor[1][1]; setring R;
+   norid; normap;
+
+   echo = 0;  
+   pause("   hit return to continue");echo = 2;
+
+   setring r;
+   nor = normalC(i, "withGens", "prim");    // a different algorithm
+   nor;
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//closureRingtower seems not to be used anywhere
+static proc closureRingtower(list L)
+"USAGE:    closureRingtower(list L); L a list of rings
+CREATE:   rings R(1),...,R(n) such that R(i)=L[i] for all i
+EXAMPLE:  example closureRingtower; shows an example
+"
+{
+  int n=size(L);
+  for (int i=1;i<=n;i++)
+    {
+      if (defined(R(i))) {
+        string s="Fixed name R("+string(i)+") leads to conflict with existing "
+              +"object having this name";
+        ERROR(s);
+      }
+      def R(i)=L[i];
+      export R(i);
+    }
+
+  return();
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+  ring R=0,(x,y),dp;
+  ideal I=x4,y4;
+  list L=primeClosure(ReesAlgebra(I)[1]);
+  L=delete(L,size(L));
+L;
+  closureRingtower(L);
+  R(1);
+  R(4);
+  kill R(1),R(2),R(3),R(4);
+}
+
+//                Up to here: procedures for normalC
+///////////////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////////////
+//                From here: miscellaneous procedures
+
+// Used for timing and comparing the different normalization procedures.
+// Option (can be entered in any order)
+// "normal"   -> uses the new algortihm (normal)
+// "normalP"  -> uses normalP
+// "normalC"  -> uses normalC, without "withGens" option
+// "primCl"   -> uses normalC, with option "withGens".
+// "111"      -> checks the output of normalM using norTest.
+// "p"        -> compares the output of norM with the output of normalP 
+//               ("normalP" option must also be set).
+// "pc"       -> compares the output of norM with the output of normalC with 
+//               option "withGens"
+//               ("primCl" option must also be set).
+proc timeNormal(ideal I, list #){
+  //--------------------------- define the method --------------------------- 
+  int isPrim, useRing;
+  int decomp = -1;
+  int norM, norC, norP, primCl;
+  int checkP, check111, checkPC;
+  int i;
+  ideal U1, U2, W;
+  poly c1, c2;
+  int ch;
+  string check;
+  string method;                //make all options one string in order to use
+                                //all combinations of options simultaneously
+  for ( i=1; i <= size(#); i++ )
+  { 
+    if ( typeof(#[i]) == "string" ) 
+    {
+      method = method + #[i];
+    }
+  }
+  if ( find(method, "normal"))
+  {norM = 1;}  
+  if ( find(method, "normalP") and (char(basering) > 0))
+  {norP = 1;}  
+  if ( find(method, "normalC"))
+  {norC = 1;}      
+  if ( find(method, "primCl"))
+  {primCl = 1;}  
+  if ( find(method, "isprim") or find(method,"isPrim") )
+  {decomp = 0;}
+  if ( find(method, "p") )
+  {checkP = 1;}  
+  if ( find(method, "pc") )
+  {checkPC = 1;}  
+  if ( find(method, "111") )
+  {check111 = 1;}    
+      
+  int tt;
+  if(norM){
+    tt = timer; 
+    if(decomp == 0){
+      "Running normal(useRing, isPrim)...";
+      list a1 = normal(I, "useRing", "isPrim");
+      "Time normal(useRing, isPrim): ", timer - tt;
+    } else {
+      "Running normal(useRing)...";
+      list a1 = normal(I, "useRing");
+      "Time normal(useRing): ", timer - tt;
+    }
+    "";
+  }
+  if(norP){
+    tt = timer; 
+    if(decomp == 0){
+      "Running normalP(isPrim)...";
+      list a2 = normalP(I, "isPrim");
+      "Time normalP(isPrim): ", timer - tt;
+    } else {
+      "Running normalP()...";
+      list a2 = normalP(I);
+      "Time normalP(): ", timer - tt;
+    }
+    "";
+  }
+  
+  if(norC){
+    tt = timer; 
+    if(decomp == 0){
+      "Running normalC(isPrim)...";
+      list a3 = normalC(I, "isPrim");
+      "Time normalC(isPrim): ", timer - tt;
+    } else {
+      "Running normalC()...";
+      list a3 = normalC(I);
+      "Time normalC(): ", timer - tt;
+    }
+    "";
+  }
+  
+  if(primCl){
+    tt = timer; 
+    if(decomp == 0){
+      "Running normalC(withGens, isPrim)...";
+      list a4 = normalC(I, "isPrim", "withGens");
+      "Time normalC(withGens, isPrim): ", timer - tt;
+    } else {
+      "Running normalC(withGens)...";    
+      list a4 = normalC(I, "withGens");
+      "Time normalC(withGens): ", timer - tt;
+    }
+    "";
+  }
+  
+  if(check111 and norM){
+    "Checking output with norTest...";
+    "WARNING: this checking only works if the original ideal was prime.";
+    list norL = list(a1[1][1][3]);
+    norTest(I, list(norL));
+    "";
+  }
+  
+  if(checkP and norP and norM){
+    "Comparing with normalP output...";
+    "WARNING: this checking only works if the original ideal was prime.";
+    U1 = a1[1][1][1];
+    c1 = a1[1][1][2];
+    U2 = a2[1][1];
+    c2 = a2[1][1][size(a2[1][1])];
+    W = changeDenominator(U1, c1, c2, groebner(I));
+    qring q = groebner(I);
+    ideal U2 = fetch(r, U2);
+    ideal W = fetch(r, W);
+    ch = 0;
+    if(size(reduce(U2, groebner(W))) == 0){
+      "U2 c U1";
+      ch = 1;
+    }
+    if(size(reduce(W, groebner(U2))) == 0){
+      "U1 c U2";
+      ch = ch + 1;
+    }
+    if(ch == 2){
+      "Output of normalP is equal.";
+    } else {
+      "ERROR: Output of normalP is different.";
+    }
+    "";
+    setring r;
+    kill q;
+  }  
+
+  if(checkPC and norM and primCl){    
+    "Comparing with primeClosure output...";
+    "WARNING: this checking only works if the original ideal was prime.";
+    // primeClosure check
+    U1 = a1[1][1][1];
+    c1 = a1[1][1][2];
+    U2 = a4[2][1];
+    c2 = a4[2][1][size(a4[2][1])];
+    W = changeDenominator(U1, c1, c2, groebner(I));
+    qring q = groebner(I);
+    ideal U2 = fetch(r, U2);
+    ideal W = fetch(r, W);
+    ch = 0;
+    if(size(reduce(U2, groebner(W))) == 0){
+      "U2 c U1";
+      ch = 1;
+    }
+    if(size(reduce(W, groebner(U2))) == 0){
+      "U1 c U2";
+      ch = ch + 1;
+    }
+    if(ch == 2){
+      "Output of normalC(withGens) is equal.";
+    } else {
+      "ERROR: Output of normalC(withGens) is different.";
+    }
+    "";
+    setring r;
+    kill q;    
+  }
+}
+
 ///////////////////////////////////////////////////////////////////////////
-
+proc norTest (ideal i, list nor, list #)
+"USAGE:   norTest(i,nor,[n]); i=prime ideal, nor=list, n=optional integer
+ASSUME:  nor is the resulting list of normal(i) (any options) or the result
+         of normalP(i,"withRing"). In particular, the ring nor[1][1] contains
+         the ideal norid and the map normap: basering/i --> nor[1][1]/norid.
+RETURN:  an intvec v such that:
+@format
+         v[1] = 1 if the normap is injective and 0 otherwise
+         v[2] = 1 if the normap is finite and 0 otherwise
+         v[3] = 1 if nor[1][1]/norid is normal and 0 otherwise 
+@end format
+         If n=1 (resp n=2) only v[1] (resp. v[2]) is computed and returned
+THEORY:  The procedure can be used to test whether the computation of the 
+         normalization was correct: basering/i --> nor[1][1]/norid is the 
+         normalization of basering/i iff v=1,1,0. 
+NOTE:    For big examples it can be hard to fully test correctness; the
+         partial test norTest(i,nor,2) is usually fast
+EXAMPLE: example norTest; shows an example
+"
+{
+//### Sollte erweitert werden auf den reduziblen Fall: einen neuen affinen
+// Ring nor[1][1]+...+nor[1][r] (direkte Summe) erzeugen, map dorthin
+// definieren und dann testen.
+
+    int prl = printlevel - voice + 2;
+    int a,b,d;
+    int n,ii;
+    if (size(#) > 0) {  n = #[1];  }
+
+    def BAS = basering;
+
+    //### make a copy of nor to have a cpoy of nor[1][1]  (not a reference to)
+    // in order not to overwrite norid and normap. 
+    // delete nor[2] (if it contains the module generators, which are not used)
+    // s.t. newnor does not belong to a ring.
+
+    list newnor = nor;         
+    if ( size(newnor) == 3 )
+    {
+       newnor = delete(newnor,2);
+    }
+    qring QAS = std(i);
+    
+    def R = newnor[1][1]; 
+    setring R;
+    int nva = nvars(R);
+    string svars = varstr(R);
+    string svar;
+
+    norid = interred(norid);
+
+    //--------- create new ring with one dp block keeping weights ------------
+    list LR = ringlist(R);
+    list g3 = LR[3];
+    int n3 = size(g3);
+    list newg3;
+    intvec V;
+
+    //--------- check first whether variables Z(i),...,A(i) exist -----------
+    for (ii=90; ii>=65; ii--)  
+    {
+       if ( find(svars,ASCII(ii)+"(") == 0 )        
+       {
+          svar = ASCII(ii);  break;
+       }
+    }
+    if ( size(svar) != 0 )
+    {
+        for ( ii = 1; ii <= nva; ii++ )
+        {
+            LR[2][ii] = svar+"("+string(ii)+")";
+            V[ii] = 1;
+        }
+    }
+    else
+    {       
+        for ( ii = 1; ii <= nva; ii++ )
+        {
+           LR[2][ii] = "Z("+string(100*nva+ii)+")";
+           V[ii] = 1;
+        }
+    }
+    
+    if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
+    {
+       list gm = g3[n3];       //last blockis module ordering
+       newg3[1] = list("dp",V);
+       newg3 = insert(newg3,gm,size(newg3));
+    }
+    else
+    {
+       list gm = g3[1];              //first block is module ordering
+       newg3[1] = list("dp",V);
+       newg3 = insert(newg3,gm);
+    }
+    LR[3] = newg3;
+//LR;"";
+    def newR = ring(LR);
+
+    setring newR;
+    ideal norid = fetch(R,norid);
+    ideal normap = fetch(R,normap);
+    if( defined(lnorid) )  { kill lnorid; }     //um ** redefinig zu beheben
+    if( defined(snorid) )  { kill snorid; }     //sollte nicht noetig sein
+
+    //----------- go to quotient ring for checking injectivity -------------
+//"mstd";
+    list lnorid = mstd(norid);
+    ideal snorid = lnorid[1];
+//"size mstdnorid:", size(snorid),size(lnorid[2]);
+//"size string mstdnorid:", size(string(snorid)),size(string(lnorid[2]));
+    qring QR = snorid;
+    ideal qnormap = fetch(newR,normap);
+    //ideal qnormap = imap(newR,normap);
+    //ideal qnormap = imap(R,normap);
+    map Qnormap = QAS,qnormap;    //r/id --> R/norid
+    
+    //------------------------ check injectivity ---------------------------
+//"injective:";
+    a = is_injective(Qnormap,QAS);          //a. Test for injectivity of Qnormap
+    dbprint ( prl, "injective: "+string(a) );
+    if ( n==1 ) { return (a); }
+   a;
+
+    //------------------------ check finiteness ---------------------------
+    setring newR;
+    b = mapIsFinite(normap,BAS,lnorid[2]);  //b. Test for finiteness of normap
+    dbprint ( prl, "finite: "+string(b) );
+    if ( n==2 ) { return (intvec(a,b)); }
+   b;
+
+    //------------------------ check normality ---------------------------
+    list testnor = normal(lnorid[2],"isPrim","noFac");
+    //### Problem: bei mehrfachem Aufruf von norTest gibt es 
+    // ** redefining norid & ** redefining normap
+    //Dies produziert Fehler, da alte norid und normap ueberschrieben werden
+    //norid und normap werden innnerhalb von proc computeRing ueberschrieben
+    //Die Kopie newR scheint das Problem zu loesen
+
+    printlevel=prl;
+    d = testnor[size(testnor)][2];             //d = delta
+    kill testnor;                              //### sollte ueberfluessig sein
+    int d1 = (d==0);                           //d1=1 if delta=0
+    dbprint ( prl, "delta: "+string(d) );
+    return(intvec(a,b,d1));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   int prl = printlevel;
+   printlevel = -1;
+   ring r = 0,(x,y),dp;
+   ideal i = (x-y^2)^2 - y*x^3;
+   list nor = normal(i);
+   norTest(i,nor);                //1,1,1 means that normal was correct
+
+   ring s = 2,(x,y),dp;
+   ideal i = (x-y^2)^2 - y*x^3;
+   nor = normalP(i,"withRing");
+   norTest(i,nor);               //1,1,1 means that normalP was correct
+   printlevel = prl; 
+}
+
+///////////////////////////////////////////////////////////////////////////
+//
+//                            EXAMPLES
+//
+///////////////////////////////////////////////////////////////////////////
 /*
-                           Examples:
-LIB"normal.lib";
-//Huneke
-ring qr=31991,(a,b,c,d,e),dp;
+//commands for computing the normalization:
+// options for normal:  "equidim", "prim"
+//                      "noDeco", "isPrim", "noFac"
+//                       (prim by default)                      
+// options for normalP: "withRing", "isPrim" or "noFac"
+// options for normalC: "equidim", "prim", "withGens", 
+//                      "noDeco", "isPrim", "noFac"
+
+//Commands for testing 'normal' 
+ list nor = normal(i); nor;            
+ list nor = normal(i,"isPrim");nor; 
+ list nor = normal(i,"equidim");nor;
+ list nor = normal(i,"prim");nor;
+ list nor = normal(i,"equidim","noFac");nor;
+ list nor = normal(i,"prim","noFac");nor;
+
+//Commands for testing 'normalP' in positive char
+ list nor = normalP(i);nor;              //withGens but no ringstructure
+ list nor = normalP(i,"withRing"); nor;  //compute the ringstructure
+ list nor = normalP(i,"isPrim"); nor;    //if i is known to be prime
+
+//Commands for testing 'normalC' 
+ list nor = normal(i); nor;            
+ list nor = normal(i,"withGens");nor;
+ list nor = normal(i,"isPrim");nor; 
+ list nor = normal(i,"equidim");nor;
+ list nor = normal(i,"prim");nor;
+ list nor = normal(i,"equidim","noFac");nor;
+ list nor = normal(i,"prim","noFac");nor;
+
+//Commands for testing correctness (i must be prime):
+list nor = normalP(i,"withRing","isPrim");
+list nor = normal(i,"isPrim");
+norTest(i,nor);       //full test for not too big examples (1,1,1 => ok)
+norTest(i,nor,2);     //partial test for big examples (1,1 => ok)
+factorize(i[1]);      //checks for irreducibility
+
+//----------------------Test Example for charp -------------------
+//Zu tun:
+//### nach minor nur std statt mstd verwenden
+//***hat bei keinem Beisp etwas gebracht -> wieder zurueck
+//### wenn interred ok, dann wieder einsetzen (am Schluss)
+//### bottelnecks bei maps beheben
+//### minor verbessern
+//### preimage verbessern (Ist imm Kern map oder imap verwendet?)
+//### Gleich in Ordnung dp wechseln, ringlist verwenden
+//### interred ev nur zum Schluss 
+//    (z.B. wenn nacher std; wenn nacher minor: testen )
+
+//Zeiten mit normalV5.lib (mstd aktiv, interred inaktiv)
+
+//SWANSON EXAMPLES: (Macaulay2, icFracP=normalP, icFractions<->normal)
+//---------------------------------------------------------------------
+//1. Series Fp[x,y,u,v]/(x2v-y2u)
+//-------------------------------
+//characteristic p   2   3    5    7    11   13   17   37   97 
+//icFracP          0.04 0.03 0.04 0.04 0.04 0.05 0.05 0.13 0.59  Mac
+//normalP           0   0    0    0     0    0    0    0   1    Sing
+//icFractions      0.08 0.09 0.09 0.09 0.14 0.15 0.15 0.15 0.15  Mac
+//normal             0   0    0    0     0    0    0    0    0   Sing
+
+2. Series Fp[u, v, w, x, y, z]/u2x4+uvy4+v2z4
+//--------------------------------------------
+//characteristic p 2    3    5    7   11 
+//icFracP         0.07 0.22 9.67 143 12543 
+//normalP          0    0    5   42  1566
+//icFractions     1.16   *    *   *    *       *: > 6h
+//normal            0    0    0   0    0
+
+//3. Series Fp[u, v, w, x, y, z]/(u2xp+uvyp+v2zp) 
+//-----------------------------------------------
+//characteristic p  2    3    5    7    11   13  17 19 23 
+//icFracP          0.06 0.07 0.09 0.27 1.81 4.89 26 56 225 
+//normalP          0     0    0    0    1    2  6  10  27
+//icFractions      0.16 1.49 75.00 4009  *    *   *  *  * 
+//normal            0     0    2    836
+//### p=7 normal braucht 807 sec in:
+// ideal endid  = phi1(endid);      //### bottelneck'
+
+//1.
+int p = 2;  ring r = p,(u,v,x,y,z),dp; ideal i = x2v-y2u;
+//2.
+int p = 7; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2x4+uvy4+v2z4;
+//3.
+int p=23; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2*x^p+uv*y^p+v2*z^p;
+
+//IRREDUCIBLE EXAMPLES:
+//---------------------  
+//timing for MacBookPro 2.2GHz Intel Core 2 Duo, 4GB Ram
+//Sing. ix86Mac-darwin version 3-1-0 (3100-2008101314)  Oct 13 2008 14:46:59
+//if no time is given: < 1  sec
+
+//Apply:
+list nor = normal(i,"isPrim"); nor;
+list nor = normalP(i,"withRing","isPrim"); nor; 
+def R=nor[1][1]; setring R; norid; normap;
+setring r;
+norTest(i,nor);
+
+int tt = timer; 
+list nor = normalP(i,"withRing","isPrim"); nor; 
+timer-tt;
+int tt = timer; 
+list nor = normal(i,"isPrim");
+timer-tt;
+
+ring r=19,(x,y,u,v),dp;    //delta -1
+ideal i=x2v-y2u; 
+//norTest 2 sec
+
+ring r=2,(y,x2,x1),lp;     //delta -1
+ideal i=y^4+y^2*x2*x1+x2^3*x1^2+x2^2*x1^3;
+//### norid hat 1 Element nach interred
+
+ring r  = 11,(x,y,z),wp(2,1,2); //alles < 1 sec
+ideal i=z3 - xy4 + x2;          //not reduced, delta =0 ok
+ideal i=y4+x5+y2x;              //not reduced, delta -1
+//interred verkleinert norid
+
+ring r=3,(u,v,x,y,z),dp;   //delta -1
+ideal i=u2x3+uvy3+v2z3;
+
+ring r=3,(u,v,x,y,z),dp;   //delta -1
+ideal i=u2x4+uvy4+v2z4;
+//norTest(i,nor);  0 sec, norTest(i,nor) haengt!
+
+ring r=5,(u,v,x,y,z),dp;   //delta -1
+ideal i=u2x6+uvy6+v2z6;  
+//normalP 5sec, normalC 1sec
+//V5: norTest(i,nor); 45 sec bei normalP, V6 12 sec
+//28 sec bei normal
+
+ring r=5,(u,v,x,y,z),dp;   //delta -1
+ideal i=u2x5+uvy5+v2z5;
+//normalP 1sec, normalC 1 sec, 
+//norTest lange: minor(jacob(I),h,J) 193 (308)sec, haengt dann bei M = std(M); 
+//norTest(i,nor,2); verwenden!
+//Sing 3.0-4 orig  >9h! hŠngt bei Q = mstd(Q)[2];
+
+ring r=2,(y,x),wp(12,5);  //delta 3
+ideal i=y5+y2x4+y2x+yx2+x12; 
+//normalP 0 sec (Test 0 sec), normalC 2 sec (Test 2 sec)
+//normalC withGens (ohne interred) 0sec
+
+ring r=2,(y,x),dp;       //delta= 22
+ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9;
+//normalP 1sec, interred verkleinert norid betraechtlich
+//normalC haengt bei minor, ideal im loop wird zu gross ###
+//interred bei normalC vergroeesert string um Faktor 4000!
+//withGens haengt bei interred in loop 4 (> 10 h) oder 
+//(nach Ausschalten von interred) bei 
+//int delt=vdim(std(modulo(f,ideal(p)))); (>?h)
+
+//Leonard1: (1. Komponente von Leonard),  delta -1 
+ring r=2,(v,u,z,y,x),dp;
+ideal i = z3+zyx+y3x2+y2x3, uyx+z2,uz+z+y2x+yx2, u2+u+zy+zx,
+          v3+vux+vz2+vzyx+vzx+uz3+uz2y+z3+z2yx2;
+//normalP 5 sec (withRing 9 sec), norTest(i,nor,2); 45 sec
+//normalC 102sec, 99sec
+//### Zeit wird bei ideal Ann = quotient(SM[2],SL[1]); und bei
+// f  = quotient(p*J,J); verbraucht
+//withGens (ohne interred) 131sec, norTest(i,nor,2); 2min25sec
+//norTest(i,nor,2);  45 sec
+
+ ring r=2,(y,x),wp(25,21); //Leonard2, delta 232
+ ring r=2,(y,x),dp;
+ ideal i=
+ y^21+y^20*x +y^18*(x^3+x+1) +y^17*(x^3+1) +y^16*(x^4+x)
+ +y^15*(x^7+x^6+x^3+x+1) +y^14*x^7 +y^13*(x^8+x^7+x^6+x^4+x^3+1)
+ +y^12*(x^9+x^8+x^4+1) +y^11*(x^11+x^9+x^8+x^5+x^4+x^3+x^2)
+ +y^10*(x^12+x^9+x^8+x^7+x^5+x^3+x+1)
+ +y^9*(x^14+x^13+x^10+x^9+x^8+x^7+x^6+x^3+x^2+1)
+ +y^8*(x^13+x^9+x^8+x^6+x^4+x^3+x) +y^7*(x^16+x^15+x^13+x^12+x^11+x^7+x^3+x)
+ +y^6*(x^17+x^16+x^13+x^9+x^8+x) +y^5*(x^17+x^16+x^12+x^7+x^5+x^2+x+1)
+ +y^4*(x^19+x^16+x^15+x^12+x^6+x^5+x^3+1)
+ +y^3*(x^18+x^15+x^12+x^10+x^9+x^7+x^4+x)
+ +y^2*(x^22+x^21+x^20+x^18+x^13+x^12+x^9+x^8+x^7+x^5+x^4+x^3)
+ +y*(x^23+x^22+x^20+x^17+x^15+x^14+x^12+x^9)
+ +(x^25+x^23+x^19+x^17+x^15+x^13+x^11+x^5);
+//normalP: dp 2sec withRing 8sec, 
+//wp 4sec, withRing:51sec Zeit in lin = subst(lin, var(ii), vip); in elimpart ),
+//norTest(i,nor,2): haengt bei mstd(norid);
+//### normalC: (m. interred): haengt bei endid = interred(endid);
+//GEFIXTES INTERRED ABWARTEN. Dann interred aktivieren
+//interred(norid) haengt u. mst(norid) zu lange
+//(o.interred): haengt bei  haengt bei list SM = mstd(i); 
+//ideal in der Mitte zu gross 
+//i = Ideal (size 118, 13 var) fuer die neue Normalisierung
+
+REDUCIBLE EXAMPLES:
+------------------
+//Apply:
+int tt = timer; 
+list nor=normalP(i,"isPrim","withRing");
+timer-tt;
+
+list nor = normal(i); nor;
+list nor = normalC(i); nor;
+list nor = normalC(i, "withGens"); nor;
+list nor = normalP(i,"withRing"); nor;
+list nor = normalP(i); nor; 
+def R=nor[1][1]; setring R; norid; normap;
+
+//Leonhard 4 Komponenten, dim=2, delta: 0,0,0,-1
+ring r=2,(v,u,z,y,x),dp;      //lp zu lange
+ideal i=z3+zyx+y3x2+y2x3, uyx+z2, v3+vuyx+vux+vzyx+vzx+uy3x2+uy2x+zy3x+zy2x2;
+//normalP: 19 sec, withRing: 22 sec
+//normalC ohne (mit) interred: 112 (113)sec, equidim: 99sec
+//normalC 1. mal 111 sec, (2.mal) 450sec!! 3.mal 172 sec
+//(unterschiedlich lange primdec, mit Auswirkungen)
+//char 19: normalC: 15sec , withGens: 14sec (o.interr.)
+
+//----------------------Test Example for special cases -------------------
+int tt = timer; 
+list nor=normalP(i,"withRing");nor;
+//list nor=normalP(i,"withRing", "isPrim");nor;
+timer-tt;
+def R1 = nor[1][1]; setring R1;  norid; normap; interred(norid);
+setring r;
+
+int tt = timer; 
+list nor=normal(i,"isPrim");nor;
+timer-tt;
+
+ring r = 29,(x,y,z),dp;
+ideal i = x2y2,x2z2;       //Nicht equidimensional, equidim reduziert nicht, ok
+ideal i  = xyz*(z3-xy4);   //### interred(norid) verkuerzt
+//je 0 sec
+
+ideal j = x,y;
+ideal i = j*xy;
+equidim(i);            
+//hat eingebettete Komponente, equidim rechnet wie in Beschreibung (ok)
+    
+ring r  = 19,(x,y),dp;
+   ideal i = x3-y4;                   //delta = 3
+   ideal i = y*x*(x3-y4);             //delta = 11; 0,0,3
+   ideal i = (x2-y3)*(x3-y4);         //delta = 13; 1,3
+   ideal i = (x-y)*(x3+y2)*(x3-y4);   //delta = 23; 0,1,3     
+   ideal i = (x-1)*(x3+y2)*(x2-y3);   //delta = 16; 0,1,1
+   ideal i = (x-y^2)^2 - y*x^3;       //delta = 3
+   //singularities at not only at 0, hier rechnet equidim falsch
+
+// -------------------------- General Examples  ---------------------------//Huneke, irred., delta=2 (Version 3-0-4: < 1sec) 
+//Version 3-0-6 default: 1sec, mit gens 2sec, mit delta 5 sec
+//(prim,noFac):ca 7 Min, prim:ca 10 min(wg facstd)
+//
+// "equidim" < 1sec irred. 5sec
+// ring r=31991,(a,b,c,d,e),dp;
+ring r=2,(a,b,c,d,e),dp;                    //delta=2
 ideal i=
 5abcde-a5-b5-c5-d5-e5,
@@ -2458 +6127 @@
 a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
 b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
-
-
-//Vasconcelos (dauert laenger: 70 sec)
-ring r=32003,(x,y,z,w,t),dp;
-ideal i=
-x2+zw,
-y3+xwt,
-xw3+z3t+ywt2,
-y2w4-xy2z2t-w3t3;
-
-//Theo1
-ring r=32003,(x,y,z),wp(2,3,6);
+//normalC: char 2, 31991: 0 sec (isPrim); char 2, equidim: 7 sec
+//norTest(i,nor,2); 1sec
+//normalP char 2: 1sec (isPrim) 
+//size(norid); size(string(norid));21 1219 interred(norid): 21 1245 (0 sec)
+
+int tt = timer; 
+list nor=normalC(i);nor;
+timer-tt;
+
+list nor = normalP(i,"isPrim");
+
+//Vasconcelos irred., delta -1 (dauert laenger)
+//auf macbook pro = 20 sec mit alter Version, 
+//Sing 3-0-6:
+// Char 32003: "equidim" 30 sec, "noFac": 30sec
+//gens: nach 9 min abgebr (haengt in Lin = ideal(T*syzf);) !!!! Hans zu tun
+//Char 2: default (charp) 2 sec, normalC ca 30 sec
+//ring r=32003,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
+ring r=2,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
+ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3;
+//normalC: char 2: 22,  sec (mit und ohne isPrim)
+//normalP char 2: 0sec (isPrim)      o. interred
+//char 32003: ### haengt in ideal endid  = phi1(endid);
+
+//-------------------------------------------------------
+//kleine Beispiele:
+
+//Theo1, irred, delta=-1
+//normalC: 1sec, normalP: 3 sec
+ring r=32003,(x,y,z),wp(2,3,6); //dim 2,dim slocus 1
 ideal i=zy2-zx3-x6;
-
-//Theo1a (CohenMacaulay and regular in codimension 2)
-ring r=32003,(x,y,z,u),wp(2,3,6,6);
+//normalC: char 2,19,32003: 0  sec (isPrim)
+//normalP (isPrim) char 2,19: 0sec, char 29: 1sec
+
+//Theo1a, CohenMacaulay regular in codim 2, dim slocus=1, delta=0 
+//normalC: 0 sec, normalP: haegt in K=preimage(R,phi,L);
+ring r=32003,(x,y,z,u),dp;
 ideal i=zy2-zx3-x6+u2;
-
-
-//Theo2
-ring r=32003,(x,y,z),wp(3,4,12);
+//normalC: char 2,32003: 0  sec (isPrim)
+//normalP (isPrim) char 2: 0sec, char 19: haengt in K = preimage(Q,phi,L);
+
+//Theo2, irreduzibel, reduziert, < 1sec, delta -1
+ring r=0,(x,y,z),wp(3,4,12);
 ideal i=z*(y3-x4)+x8;
-
-//Theo2a
-ring r=32003,(T(1..4)),wp(3,4,12,17);
+//normalC: char 2,32003,0: 0  sec (isPrim)
+//normalP (isPrim) char 2: 0 1sec, char 19: 1sec char 29: 7 sec
+
+//Theo2a, reduiziert, 2-dim, dim_slocus=1, alte Version 3 sec,
+//normalP ca 30 sec, normalC ca 4sec, delta -1
+//ring r=32003,(T(1..4)),wp(3,4,12,17);
+//ring r=11,(T(1..4)),dp;
+ring r=11,(T(1..4)),wp(3,4,12,17);
 ideal i=
 T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
@@ -2488 +6184 @@
 T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
 T(1)^6*T(2)*T(3)+T(1)^2*T(2)^4*T(3)+T(1)^3*T(2)^2*T(4)-T(1)^2*T(2)*T(3)^2+T(4)^2;
-
-//Theo3
-ring r=32003,(x,y,z),wp(3,5,15);
+//normalC: char 2,32003: 0  sec (isPrim)
+//normalP (isPrim) char 2: 0sec, char 11 2se, char 19: 13sec 
+//norTest 48sec in char11
+//### interred verkuerzt
+//char 29: haengt in K = preimage(Q,phi,L);
+
+//Theo3, irred, 2-dim, 1-dim sing, < 1sec
+ring r=11,(x,y,z),wp(3,5,15);
 ideal i=z*(y3-x5)+x10;
-
-
-//Theo4
-ring r=32003,(x,y,z),dp;
-ideal i=(x-y)*(x-z)*(y-z);
-
-//Theo5
-ring r=32003,(x,y,z),wp(2,1,2);
-ideal i=z3-xy4;
+//normalC: char 2,0: 0  sec (withRing)
+//normalP (withRing) char 2,11: 0sec, char 19: 13sec norTest 12sec(char 11)
+
+//Theo4 reducible, delta (0,0,0) -1
+ring r=29,(x,y,z),dp;
+ideal i=(x-y)*(x-z)*(y-z); 
+//normalC: char 2,32003: 0  sec 
+//normalP char withRing 2, 29: 0sec, 6sec
 
 //Theo6
 ring r=32003,(x,y,z),dp;
 ideal i=x2y2+x2z2+y2z2;
-
-ring r=32003,(a,b,c,d,e,f),dp;
+//normalC: char 2,32003: 0  sec 
+//normalP char withRing 2, 29: 0sec, 4sec
+
+//Sturmfels, CM, 15 componenten, alle glatt
+ring r=0,(b,s,t,u,v,w,x,y,z),dp;
+ideal i= bv+su, bw+tu, sw+tv, by+sx, bz+tx, sz+ty,uy+vx,uz+wx,vz+wy,bvz;
+//normalC car 11, 0: 1sec, normalP 0 sec
+
+//riemenschneider, , dim 3, 5 Komp. delta (0,0,0,0,0), -1
+ring r=2,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
+ideal i=xz,vx,ux,su,qu,txy,stx,qtx,uv2z-uwz,uv3-uvw,puv2-puw;
+//alles 0 sec in char 2
+
+//4 Komponenten, alle glatt, 0sec
+ring r=11,(x,y,z,t),dp;
+ideal i=x2z+xzt,xyz,xy2-xyt,x2y+xyt;
+
+//dim 3, 2 Komponenten delta (-1,0), -1
+ring r=2,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
+ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;
+//alles 0 sec in char 2
+//---------------------------------------------------------
+int tt = timer; 
+list nor=normalP(i,"normalC","withRing");nor;
+timer-tt;
+
+//St_S/Y, 3 Komponenten, 2 glatt, 1 normal 
+//charp haengt (in char 20) in K=preimage(R,phi,L);
+//ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
+ring r=11,(b,s,t,u,v,w,x,y,z),dp;
+ideal i=wy-vz,vx-uy,tv-sw,su-bv,tuy-bvz;