source: git/Singular/LIB/normal.lib @ 3de58c

///////////////////////////////////////////////////////////////////////////////
version="$Id: normal.lib,v 1.35 2001-03-19 22:57:16 greuel Exp $";
category="Commutative Algebra";
info="
LIBRARY:  normal.lib     Normalization of Affine Rings
AUTHORS:  G.-M. Greuel,  greuel@mathematik.uni-kl.de,
@*        G. Pfister,    pfister@mathematik.uni-kl.de

PROCEDURES:
 normal(I);             computes the normalization of basering/I
 HomJJ(L);              presentation of End_R(J) as affine ring, L a list
";

LIB "general.lib";
LIB "sing.lib";
LIB "primdec.lib";
LIB "elim.lib";
LIB "presolve.lib";
LIB "inout.lib";
LIB "ring.lib";
///////////////////////////////////////////////////////////////////////////////
static proc isR_HomJR (list Li)
"USAGE:   isR_HomJR (Li);  Li = list: ideal SBid, ideal J, poly p
COMPUTE: module Hom_R(J,R) = R:J and compare with R
ASSUME:  R    = P/SBid,  P = basering
         SBid = standard basis of an ideal in P,
         J    = ideal in P containing the polynomial p,
         p    = nonzero divisor of R
RETURN:  1 if R = R:J, 0 if not
EXAMPLE: example isR_HomJR;  shows an example
"
{
   int n, ii;
 def P = basering;
   ideal SBid = Li[1];
   ideal J = Li[2];
   poly p = Li[3];
   attrib(SBid,"isSB",1);
   attrib(p,"isSB",1);
 qring R    = SBid;
   ideal J  = fetch(P,J);
   poly p   = fetch(P,p);
   ideal f  = quotient(p,J);
   ideal lp = std(p);
   n=1;
   for (ii=1; ii<=size(f); ii++ )
   {
      if ( reduce(f[ii],lp) != 0)
      { n = 0; break; }
   }
   return (n);
 //?spaeter hier einen Test ob Hom(I,R) = Hom(I,I)?
}
example
{"EXAMPLE:";  echo = 2;
  ring r   = 0,(x,y,z),dp;
  ideal id = y7-x5+z2;
  ideal J  = x3,y+z;
  poly p   = xy;
  list Li  = std(id),J,p;
  isR_HomJR (Li);

  ring s   = 0,(t,x,y),dp;
  ideal id = x2-y2*(y-t);
  ideal J  = jacob(id);
  poly p   = J[1];
  list Li  = std(id),J,p;
  isR_HomJR (Li);
}
///////////////////////////////////////////////////////////////////////////////

proc HomJJ (list Li)
"USAGE:   HomJJ (Li); Li list: ideal SBid, ideal id, ideal J, poly p, int count
ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
@*       SBid = standard basis of id,
@*       J    = ideal of P containing the polynomial p,
@*       p    = nonzero divisor of R
@*       count controls printlevel during recursive call
COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
         affine ring, together with the canonical map R --> Hom_R(J,J),
         where R is the quotient ring of P modulo the standard basis SBid.
RETURN:  a list l of two objects
@format
         l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
               such that l[1]/endid = Hom_R(J,J) and
               endphi describes the canonical map R -> Hom_R(J,J)
         l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
@end format
NOTE:    printlevel >=1: display comments (default: printlevel=0)
EXAMPLE: example HomJJ;  shows an example
"
{
//---------- initialisation ---------------------------------------------------

   int isIso,isPr,isCo,isRe,isEq,ii,jj,q,y,count;
   intvec rw,rw1;
   list L;
   if(size(Li)>=5)
   {
     count = Li[5];
   }
   y = printlevel-voice+count;  
 def P = basering;
   ideal SBid, id, J = Li[1], Li[2], Li[3];
   poly p = Li[4];
   attrib(SBid,"isSB",1);
   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not

//---- set attributes for special cases where algorithm can be simplified -----
   if( homo==1 )
   {
      rw = ringweights(P);
   }
   if( typeof(attrib(id,"isPrim"))=="int" )
   {
      if(attrib(id,"isPrim")==1)  { isPr=1; }
   }
   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
   {
      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
   }
   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
   {
      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
   }
   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
   {
      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
   }
   if( typeof(attrib(id,"isEquidimensional"))=="int" )
   {
      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
   }
//-------------------------- go to quotient ring ------------------------------
 qring R  = SBid;
   ideal id = fetch(P,id);
   ideal J  = fetch(P,J);
   poly p   = fetch(P,p);
   ideal f,rf,f2;
   module syzf;

//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
   if ( y>=1 )
   {
     "// compute p*Hom(J,J) = p*J:J, p a non-zerodivisor";
     "//   p is equal to:"; "";
     p;
     "";
   }
   f  = quotient(p*J,J);
   if ( y>=1 )
   { "// the module p*Hom(J,J) = p*J:J, p a non-zerodivisor";
      "// p"; p;
      "// f=p*J:J";
      if( y>=2 ) { f; }      
   }
   f2 = std(p);

   if(isIso==0)
   {
     ideal f1=std(f);
     attrib(f1,"isSB",1);
    // if( codim(f1,f2) >= 0 )
    // {
    //  dbprint(printlevel-voice+3,"// dimension of non-normal locus is zero");
    //    isIso=1;
    // }
  }
//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------

   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R
   if ( size(rf) == 0 )
   {
      if ( homog(f) && find(ordstr(basering),"s")==0 )
      {
         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
      }
      else
      {
         ring newR1 = char(P),(X(1..nvars(P))),dp;
      }
      ideal endphi = maxideal(1);
      ideal endid = fetch(P,id);
      L=substpart(endid,endphi,homo,rw);
      def lastRing=L[1];
      setring lastRing;

      attrib(endid,"isCohenMacaulay",isCo);
      attrib(endid,"isPrim",isPr);
      attrib(endid,"isIsolatedSingularity",isIso);
      attrib(endid,"isRegInCodim2",isRe);
      attrib(endid,"isEqudimensional",isEq);
      attrib(endid,"isCompleteIntersection",0);
      attrib(endid,"isRad",0);
      L=lastRing;
      L = insert(L,1,1);
      if(y>=1)
      {
         "// case R = Hom(J,J), we are ready with this component";
         "   ";
       if( y>=2 ) 
       {
        lastRing;
         "   ";
         "//   the new ideal";
         if( y>=2 ) { endid; } 
         "   ";
         "//   the old ring";
         "   ";
         P;
         "   ";
         "//   the old ideal";
         "   ";
         setring P;
         id;
         "   ";
         setring lastRing;
         "//   the map";
         "   ";
         endphi;
         "   ";
         pause();
         newline;
       }
      }
      setring P;
      return(L);
   }
   if(y>=1)
   {
      "// R is not equal to Hom(J,J), we have to try again";
    if( y>=2 )
    {  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 = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
   q = size(f);
   syzf = syz(f);

   if ( homo==1 )
   {
      rw1 = rw,0;
      for ( ii=2; ii<=q; ii++ )
      {
         rw  = rw, deg(f[ii])-deg(f[1]);
         rw1 = rw1, deg(f[ii])-deg(f[1]);
      }
      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
   }
   else
   {
      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
   }

   map psi1 = P,maxideal(1);
   ideal SBid = psi1(SBid);
   attrib(SBid,"isSB",1);

 qring newR = std(SBid);
   map psi = R,ideal(X(1..nvars(R)));
   ideal id = psi(id);
   ideal f = psi(f);
   module syzf = psi(syzf);
   ideal pf,Lin,Quad,Q;
   matrix T,A;
   list L1;

//---------- computation of Hom(J,J) as affine ring ---------------------------
// determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course
// the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
// It is a fact, that the kernel is generated by the linear and the quadratic
// relations

   pf = f[1]*f;
   T = matrix(ideal(T(1..q)),1,q);
   Lin = ideal(T*syzf);
   if(y>=1)
   {
      "// the ring structure of Hom(J,J) as R-algebra";
      " ";
      "//   the linear relations";
      " ";
      if( y>=2 )
      {  Lin;      
         pause();
         "";
      }
   }
   for (ii=2; ii<=q; ii++ )
   {
      for ( jj=2; jj<=ii; jj++ )
      {
         A = lift(pf,f[ii]*f[jj]);
         Quad = Quad, ideal(T(jj)*T(ii) - T*A);          // quadratic relations
      }
   }
   if(y>=1)
   {
      "//   the quadratic relations";
        if( y>=2 )
      {      
          "   ";
         interred(Quad);
         pause();
         newline;
      }
   }
   Q = Lin+Quad;
   Q = subst(Q,T(1),1);
   Q = interred(reduce(Q,std(0)));
//---------- reduce number of variables by substitution, if possible ----------
   if (homo==1)
   {
      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
   }
   else
   {
      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
   }

   ideal endid  = imap(newR,id)+imap(newR,Q);
   ideal endphi = ideal(X(1..nvars(R)));

   L=substpart(endid,endphi,homo,rw);
   def lastRing=L[1];
   setring lastRing;
   attrib(endid,"isCohenMacaulay",isCo);
   attrib(endid,"isPrim",isPr);
   attrib(endid,"isIsolatedSingularity",isIso);
   attrib(endid,"isRegInCodim2",isRe);
   attrib(endid,"isEquidimensional",isEq);
   attrib(endid,"isCompleteIntersection",0);
   attrib(endid,"isRad",0);
  // export(endid);
  // export(endphi);
   if(y>=1)
   {
      "//   the new ring after reduction of the number of variables";
      show(lastRing);
      pause(); "
      ";
    if(y >=2)
    {
      lastRing;
      "   ";
      "//   the new ideal";
      "   ";
      endid;
      "   ";
      "//   the old ring";
      "   ";
      P;
      "   ";
      "//   the old ideal";
      "   ";
      setring P;
      id;
      "   ";
      setring lastRing;
      "//   the map";
      "   ";
      endphi;
      "   ";
      pause();
      newline;
    }
   }
   L = lastRing;
   L = insert(L,0,1);
   return(L);
}
example
{"EXAMPLE:";  echo = 2;
  ring r   = 0,(x,y),wp(2,3);
  ideal id = y^2-x^3;
  ideal J  = x,y;
  poly p   = x;
  list Li = std(id),id,J,p;
  list L   = HomJJ(Li);
  def end = L[1];    // defines ring L[1], containing ideals endid, endphi
  setring end;       // makes end the basering
  end;
  endid;             // end/endid is isomorphic to End(r/id) as ring
  map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
  psi;
}

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

proc normal(ideal id, list #)
"USAGE:   normal(i [,choose]);  i a radical ideal, choose empty or 1
         if choose=1 the normalization of the associated primes is computed
         (which is sometimes more efficient)
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:
@format
         each ring nor[i] contains two ideals
         with given names norid and normap such that
           - the direct sum of the rings nor[i]/norid is
             the normalization of basering/id;
           - normap gives the normalization map from basering/id
             to nor[i]/norid (for each i)
@end format
NOTE:    to use the i-th ring type: def R=nor[i]; setring R;.
@*       Not implemented for local or mixed orderings and quotient rings.
@*       If the input ideal i is weighted homogeneous a weighted ordering may
         be used (qhweight(i); computes weights).
@*       printlevel = 1: count normalization loops (default: printlevel=0)
@*       printlevel > 1: protocoll of normalization steps
@*       printlevel > 2: protocoll of all normalization steps, pauses
@*       printlevel > 3: display some (>4 all) intermediate results, pauses
EXAMPLE: example normal; shows an example
"
{
   int i,j,y;
   string sr;
   list result,prim,keepresult;
   y = printlevel-voice+2;

   attrib(id,"isRadical",1);
   if ( ord_test(basering) != 1)
   {
     "";
     "// Not implemented for this ordering,";
     "// please change to global ordering!";
     return(result);
   }
   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
   {
      if(attrib(id,"isCompleteIntersection")==1)
      {
         attrib(id,"isCohenMacaulay",1);
         attrib(id,"isEquidimensional",1);
      }
   }
   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
   {
      if(attrib(id,"isCohenMacaulay")==1)
      {
         attrib(id,"isEquidimensional",1);
      }
   }
   if( typeof(attrib(id,"isPrim"))=="int" )
   {
      if(attrib(id,"isPrim")==1)
      {
         attrib(id,"isEquidimensional",1);
      }
   }
   if(size(#)==0)
   {
      if( typeof(attrib(id,"isEquidimensional"))=="int" )
      {
         if(attrib(id,"isEquidimensional")==1)
         {
            prim[1]=id;
         }
         else
         {
            prim=equidim(id);
         }
      }
      else
      {
         prim=equidim(id);
      }
      if(y>=1)
      {
        "";
         "// we have ",size(prim),"equidimensional component(s)";
      }
   }
   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 component(s)";
      }
   }
   for(i=1; i<=size(prim); i++)
   {
      if(y>=1)
      {
         "";
         "// BEGIN with equidimensional/irreducible component",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);

      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);j++)
      {
         result=insert(result,keepresult[j]);
      }
      sr = string(size(result));
   }
      dbprint(y+1,"
// 'normal' created a list of "+sr+" ring(s).
// To see the rings, type (if the name of your list is nor):
     show(nor);
// To access the 1-st ring and map (similar for the others), type:
     def R = nor[1]; setring R;  norid; normap;
// R/norid is the 1-st ring of the normalization and
// normap the map from the original basering to R/norid");

      //kill endphi,endid;
      return(result);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),wp(2,1,2);
   ideal i=z3-xy4;
   list nor=normal(i);
   show(nor);
   def r1=nor[1];
   setring r1;
   norid;
   normap;
}

///////////////////////////////////////////////////////////////////////////////
static proc normalizationPrimes(ideal i,ideal ihp, int count, list #)
"USAGE:   normalizationPrimes(i,ihp[,si,countt]);  i prime ideal, ihp map
         (partial normalization), si ideal of singular locus, 
         count = integer to count the number of normalization loops
RETURN:  a list of one ring L=R, in  R are two ideals
         S,M such that R/M is the normalization of original basering
         S is a standardbasis of M
NOTE:    to use the ring: def r=L[1];setring r;
         printlevel = 1: count normalization loops (default: printlevel=0)
         printlevel > 1: protocoll of normalization steps
         printlevel > 2: protocoll of all normalization steps, pauses
         printlevel > 3: display some (>4 all) intermediate results, pauses
EXAMPLE: example normalizationPrimes; shows an example
"
{
   count = count+1;   
   int y = printlevel-voice+count+1;  // y=printlevel (default: y=0)
   if(y>=0)
   {
     "// START normalization LOOP ",count;
   }
   if( y>=3)
   {
     "// with ideal";  "";
     i;  "";
     basering;  "";
     pause();
     newline;
   }


   def BAS=basering;
   list result,keepresult1,keepresult2;
   ideal J,SB,MB;
   int depth,lauf,prdim;
   int ti=timer;

   if(size(i)==0)
   {
      if(y>=1)
      {
          "// the ideal was the zero-ideal";
      }
         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");
         ideal norid=ideal(0);
         ideal normap=fetch(BAS,ihp);
         export norid;
         export normap;
         result=newR7;
         setring BAS;
         return(result);
   }

   if(y>=1)
   {
   "// SB-computation of ideal of size",size(i),"in ring with",nvars(basering),"variables";
   }
// -------------- SB-computation of the input ideal ----------------------
   list SM=mstd(i);                //here the work starts, SM[1] a SB of i
                                   //SM[2] a smaller set of generators of i
   int dimSM =  dim(SM[1]);        //dimension of i, V(i)=variety to normalize
  // Case: Get an ideal containing a unit
   if( dimSM == -1)
   {  "";
      "      // A unit ideal was found.";
      "      // Stop with partial result computed so far";"";

         MB=SM[2];
         intvec rw;
         list LL=substpart(MB,ihp,0,rw);
         def newR6=LL[1];
         setring newR6;
         ideal norid=endid;
         ideal normap=endphi;
         kill endid,endphi;
         export norid;
         export normap;
         result=newR6;
         setring BAS;
         return(result);
   }

   if(y>=1)
   {
      "//   the dimension is:"; "";
      dimSM;"";
   }

   if(size(#)>0)                   //ideal of sing locus is given in #[1]
   {
      list JM=mstd(#[1]);
      if( typeof(attrib(#[1],"isRad"))!="int" )
      {
         attrib(JM[2],"isRad",0);
      }
   }

// -------- transfer attributes from ideal i to SM[2], SM = mstd(i) --------
   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);
   }

//************* case 0: check and prepare special cases ****************
//no computation for the normalization in case 0

//-------------------------- 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;
         setring BAS;
         return(result);
     }
   }
   // recall: SM = mstd(i), i = ideal to start with in normaliztion loop
   //         JM = mstd(#[1]), #[1]= ideal of singular locus of i
   //              #[1] is given after the first normalization loop

//---------------- the homogeneous zero-dimensional case ----------------
   if((dim(SM[1])==0)&&(homog(SM[2])==1))
   {
      if(y>=1)
      {
         "// the ideal was zero-dimensional and homogeneous";
      }
      MB=maxideal(1);
      intvec rw;
      list LL=substpart(MB,ihp,0,rw);
      def newR5=LL[1];
      setring newR5;
      ideal norid=endid;
      ideal normap=endphi;
      kill endid,endphi;
      export norid;
      export normap;
      result=newR5;
      setring BAS;
      return(result);
   }

//------------- the homogeneous one-dimensional case -------------------
//in this case i defines a line because it is irreducible and homogeneous
   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
        &&(homog(SM[2])==1))
   {
      if(y>=1)
      {
         "// the ideal defines a line";
      }
      MB=SM[2];
      intvec rw;
      list LL=substpart(MB,ihp,0,rw);
      def newR4=LL[1];
      setring newR4;
      ideal norid=endid;
      ideal normap=endphi;
      kill endid,endphi;
      export norid;
      export normap;
      result=newR4;
      setring BAS;
      return(result);
   }
//----------------------- the higher dimensional case ----------------------
   //we test first of all CohenMacaulay and
   //complete intersection
   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
   {
      //test for complete intersection
      attrib(SM[2],"isCohenMacaulay",1);
      attrib(SM[2],"isCompleteIntersection",1);
      attrib(SM[2],"isEquidimensional",1);
      if(y>=1)
      {
         "// the ideal is a complete intersection";
      }
   }

   //compute the singular locus+lower dimensional components

//------- case: not(isolated singularity and homogeneous) -------------------
   if((attrib(SM[2],"isIsolatedSingularity")==0) && (size(#)==0))
   {
//---------- computation of ideal of singular locus -------------------------
      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
      //ti=timer;
      if(y >=1 )
      {
         "// SB of singular locus will be computed ";
         if(y >=2 )
         {
           "//in ring:";
           show(basering);
         }
      }
      
      J = simplify(J,16);  //this makes ist for huge J much faster
      ideal sin=J,SM[2]; 
      list JM=mstd(sin);
      
      //JM[1] SB of singular locus, JM[2] minbasis of singular locus
      //SM[1] SB of ideal in normalisation loop, SM[2] minbasis

      if(y>=1)
      {
         "//   the dimension of the singular locus is:";"";
         dim(JM[1]); "";
      }
      //   timer-ti;
      attrib(JM[1],"isSB",1);
      if(dim(JM[1])==-1)
      {
         if(y>=1)
         {
            "// the ideal is smooth";
         }
         MB=SM[2];
         intvec rw;
         list LL=substpart(MB,ihp,0,rw);
         def newR3=LL[1];
         setring newR3;
         ideal norid=endid;
         ideal normap=endphi;
         kill endid,endphi;
         export norid;
         export normap;
         result=newR3;
         setring BAS;
         return(result);
      }
      if(dim(JM[1])==0 && (homog(SM[2])==1))
      {
         attrib(SM[2],"isIsolatedSingularity",1);
      }
      if(dim(JM[1])<=dim(SM[1])-2)
      {
         attrib(SM[2],"isRegInCodim2",1);
      }
   }
   
//------------ case: isolated singularity and homogeneous -----------------
   if(attrib(SM[2],"isIsolatedSingularity")==1)   
   {
     if(size(#)==0)
     {
        if(defined(JM)==voice)
        {  JM=maxideal(1),maxideal(1);
        }
        else
        {  list JM=maxideal(1),maxideal(1);
        
        }
        attrib(JM[1],"isSB",1);
     }
   }

//------------ case: Cohen Macaulay and regular in codim 2 ----------------- 
//in this case we are ready, the ring is normal  
   if((attrib(SM[2],"isRegInCodim2")==1)&&(attrib(SM[2],"isCohenMacaulay")==1))
   {
      if(y>=1)
      {
            "// the ideal was CohenMacaulay and regular in codimension 2";
      }
      MB=SM[2];
      intvec rw;
      list LL=substpart(MB,ihp,0,rw);
      def newR6=LL[1];
      setring newR6;
      ideal norid=endid;
      ideal normap=endphi;
      kill endid,endphi;
      export norid;
      export normap;
      result=newR6;
      setring BAS;
      return(result);
   }

//************* case 1: isolated singularity and homogeneous ****************

   // recall: SM = mstd(i), i = ideal to start with in normaliztion loop
   //         JM = mstd(#[1]), #[1]= ideal of singular locus of i
   //              #[1] is given after the first normalization loop
   //if SM is an isolated singularity and homogeneous, we know that this 
   //persists in the following normalization loops and things are easier
   //since the radical of the singular locus is the maximal ideal
   //JM ideal of singular locus, SM ideal of variety

   if(attrib(SM[2],"isIsolatedSingularity")==1)
   {
      if(y>=1)
      {
         "// CASE 1: unique isolated singularity at 0";
         "// radial of singular locus is the maximal ideal";
      }
      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
      //SL = 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
  //--------------- case: a nonzero-divisor was found ---------------
      if(size(qAnn)==0)
      {
         if(y>=1)
         {
            "";
            "//   a nonzero-divisor was found";
            "//   the ideal rad(J):";
            "";
            maxideal(1);
            newline;
            "// TEST for normality: R=Hom(J,J)?";
            "";
         }
    //test for normality:
         //?spaeter: test for normality, Hom(I,R)==R?
         list RR;
         RR=SM[1],SM[2],maxideal(1),SL[1],count;
         //  ti=timer;
         RR=HomJJ(RR);
         //  timer-ti;
         if(RR[2]==0) 
    //the ring is not normal, start a new normalization loop
         {
            def newR=RR[1];
            setring newR;
            map psi=BAS,endphi;
         //  ti=timer;
            list tluser=normalizationPrimes(endid,psi(ihp),count);
         //  timer-ti;
            setring BAS;
            return(tluser);
         }
    //the ring is normal, prepare output and go home
         MB=SM[2];
         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");
         ideal norid=fetch(BAS,MB);
         ideal normap=fetch(BAS,ihp);
         export norid;
         export normap;
         result=newR7;
         setring BAS;
         return(result);

       }
  //--------------- case: a zero-divisor was found ---------------
      //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=qAnn+SM[2] defines components of R/SM in the complement of V(SL[1])
      //id1 defines components of R/SM in the complement of V(id)
      //?????instead of id1 we can take SL[1]+Ann+SM[2]???????????
       else
       {
         if(y>=0)
         {
            "//   a zero-divisor was found, we have SPLITTING of ideals";
           if(y >=2 )
           {
             "// in ring:";
             show(basering);
            }
            "";
         }
         
          ideal id=qAnn+SM[2];

          attrib(id,"isCohenMacaulay",0);
          attrib(id,"isPrim",0);
          attrib(id,"isIsolatedSingularity",1);
          attrib(id,"isRegInCodim2",0);
          attrib(id,"isCompleteIntersection",0);
          attrib(id,"isEquidimensional",0);

          if(y>=0)
          {
    "//   start a normalization loop with 1st split ideal (size",size(id),"in",nvars(basering),"vars)";
          }
          keepresult1=normalizationPrimes(id,ihp,count);
          ideal id1=quotient(SM[2],Ann)+SM[2];
//          evtl. qAnn statt Ann nehmen
//          ideal id=SL[1]+SM[2];

          attrib(id1,"isCohenMacaulay",0);
          attrib(id1,"isPrim",0);
          attrib(id1,"isIsolatedSingularity",1);
          attrib(id1,"isRegInCodim2",0);
          attrib(id1,"isCompleteIntersection",0);
          attrib(id1,"isEquidimensional",0);

          if(y>=0)
          {
              "//   start a normalization loop with 2nd split ideal (size",size(id),"in",nvars(basering),"vars)";
          }
          keepresult2=normalizationPrimes(id1,ihp,count);

          for(lauf=1;lauf<=size(keepresult2);lauf++)
          {
             keepresult1=insert(keepresult1,keepresult2[lauf]);
          }
          return(keepresult1);
       }
   }

//********** case 2: no unique isolated singularity at 0 *************

   //in this case the radical must be computed
   // recall: SM = mstd(i), i = ideal to start with in normaliztion loop
   //         JM = mstd(#[1]), #[1]= ideal of singular locus of i
   //              #[1] is given after the first normalization loop
   //test for normality:  Hom(J,J)=R
   //test for non-normality: Hom(J,J)!=R, we can use Hom(J,J) to continue

      if(y>=1)
      {
         "// CASE 2: no unique isolated singularity";
         "// radical has to be computed";
      }

   ideal SL=simplify(reduce(JM[2],SM[1]),2);
   //SL = elements of ideal of singular locus not contained in ideal i
   ideal Ann=quotient(SM[2],SL[1]);
   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
   //qAnn=0 ==> SL[1] is a nonzero-divisor of R/SM

  //--------------- case: a nonzero-divisor was found ---------------
   if(size(qAnn)==0)
   {
      list RR;
      list RS;
    //now we have to compute the radical
      if(y>=1)
      {
        "";
         "//   a nonzero-divisor was found";
         "//   radical computation of ideal of singular locus";
      }
      
    //-------------- computation of the radical J --------------------
    //We have at least two possibilities:
    //J=radical(JM[2]), the radical of the full singular locus, or
    //J=radical(SM[2]+ideal(SL[1])), JM[2] contains SM[2]+ideal(SL[1])
      //the first is usually better!

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

      J=radical(JM[2]);
      // alternativ:   J=radical(SM[2]+ideal(SL[1])); 

      if(y>=2)
      {
          "// the radical is equal to:";        
          J;
          newline;
      }
      if(y>=1)
      {        
        "// TEST for normality: R=Hom(J,J)?";
        "";
      }

      JM=J,J;              //J = new ideal for singular locus
      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen

    //test for normality:
      RR=SM[1],SM[2],JM[2],SL[1],count;
      RS=HomJJ(RR);
      
    //--- the ring is normal, prepare output and go home
      if(RS[2]==1)
      {
         def lastR=RS[1];
         setring lastR;
         map psi1=BAS,endphi;
         ideal norid=endid;
         ideal normap=psi1(ihp);
         kill endid,endphi;
         export norid;
         export normap;
         setring BAS;
         return(lastR);
      }

    //--- the ring is not normal, start a new normalization loop       
      int n=nvars(basering);
      ideal MJ=JM[2];

      def newR=RS[1];
      setring newR;

      map psi=BAS,endphi;
      list tluser=
           normalizationPrimes(endid,psi(ihp),count,simplify(psi(MJ)+endid,4));
      setring BAS;
      return(tluser);
   }
  //--------------- case: a zero-divisor was found ---------------
  //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=qAnn+SM[2] defines components of R/SM in the complement of V(SL[1])
  //id1 defines components of R/SM in the complement of V(id)
      //?????instead of id1 we can take SL[1]+Ann+SM[2]???????????

   // A component with singular locus the whole component found:
   if( Ann == 1)
   {
      "// Input appeared not to be a radical ideal!";
      "// A (everywhere singular) component with ideal";
      "// equal to its Jacobian ideal was found";
      "// Procedure will stop with partial result computed so far";"";

         MB=SM[2];
         intvec rw;
         list LL=substpart(MB,ihp,0,rw);
         def newR6=LL[1];
         setring newR6;
         ideal norid=endid;
         ideal normap=endphi;
         kill endid,endphi;
         export norid;
         export normap;
         result=newR6;
         setring BAS;
         return(result);
   }
   // The general case with splitting of ring, i.e. qAnn!=0
   else
   {
      int equi=attrib(SM[2],"isEquidimensional");
      ideal new1=qAnn+SM[2];
      execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");
      if(y>=0)
      {
            "//   a zero-divisor was found, we have SPLITTING of ideals";
            "";
      }
      ideal vid=fetch(BAS,new1);
      ideal ihp=fetch(BAS,ihp);
      attrib(vid,"isCohenMacaulay",0);
      attrib(vid,"isPrim",0);
      attrib(vid,"isIsolatedSingularity",0);
      attrib(vid,"isRegInCodim2",0);
      if(equi==1)
      {
         attrib(vid,"isEquidimensional",1);
      }
      else
      {
         attrib(vid,"isEquidimensional",0);
      }
      attrib(vid,"isCompleteIntersection",0);

      if(y>=0)
      {
            "// start a normalization loop with 1st split ideal (size",size(vid),"in",nvars(basering),"vars)";
      }

      keepresult1=normalizationPrimes(vid,ihp,count);

      setring BAS;
      ideal new2=quotient(SM[2],Ann)+SM[2];
// evtl. qAnn nehmen
      execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
                      +ordstr(basering)+");");

      ideal vid=fetch(BAS,new2);
      ideal ihp=fetch(BAS,ihp);
      attrib(vid,"isCohenMacaulay",0);
      attrib(vid,"isPrim",0);
      attrib(vid,"isIsolatedSingularity",0);
      attrib(vid,"isRegInCodim2",0);
      if(equi==1)
      {
         attrib(vid,"isEquidimensional",1);
      }
      else
      {
         attrib(vid,"isEquidimensional",0);
      }
      attrib(vid,"isCompleteIntersection",0);

      if(y>=0)
      {
            "// start a normalization loop with 2nd split ideal (size",size(vid),"in",nvars(basering),"vars)";
      }

      keepresult2=normalizationPrimes(vid,ihp,count);

      setring BAS;
      for(lauf=1;lauf<=size(keepresult2);lauf++)
      {
         keepresult1=insert(keepresult1,keepresult2[lauf]);
      }
      return(keepresult1);
   }
}
example
{ "EXAMPLE:";echo = 2;
   // Huneke
   ring qr=31991,(a,b,c,d,e),dp;
   ideal i=
   5abcde-a5-b5-c5-d5-e5,
   ab3c+bc3d+a3be+cd3e+ade3,
   a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
   abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
   ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
   a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
   a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;

   list pr=normalizationPrimes(i,maxideal(1),1);
   def r1=pr[1];
   setring r1;
   norid;
   normap;
}
///////////////////////////////////////////////////////////////////////////////
static proc substpart(ideal endid, ideal endphi, int homo, intvec rw)

"//Repeated application of elimpart to endid, until no variables can be
//directy substituded. homo=1 if input is homogeneous, rw contains
//original weights, endphi (partial) normalization map";

{
   def newRing=basering;
   int ii,jj;
   map phi = basering,maxideal(1);

   //endid=diagon(endid);

   list Le = elimpart(endid);
                                     //this proc and the next loop try to
   int q = size(Le[2]);              //substitute as many variables as possible
   intvec rw1 = 0;                   //indices of substituted variables
   rw1[nvars(basering)] = 0;
   rw1 = rw1+1;

   while( size(Le[2]) != 0 )
   {
      endid = Le[1];
      map ps = newRing,Le[5];

      phi = ps(phi);
      for(ii=1;ii<=size(Le[2])-1;ii++)
      {
         phi=phi(phi);
      }
      //eingefuegt wegen x2-y2z2+z3
      kill ps;

      for( ii=1; ii<=size(rw1); ii++ )
      {
         if( Le[4][ii]==0 )
         {
            rw1[ii]=0;                             //look for substituted vars
         }
      }
      Le=elimpart(endid);
      q = q + size(Le[2]);
   }
   endphi = phi(endphi);

//---------- return -----------------------------------------------------------
// in the homogeneous case put weights for the remaining vars correctly, i.e.
// delete from rw those weights for which the corresponding entry of rw1 is 0

   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
   {
      jj=1;
      for( ii=2; ii<size(rw1); ii++)
      {
         jj++;
         if( rw1[ii]==0 )
         {
            rw=rw[1..jj-1],rw[jj+1..size(rw)];
            jj=jj-1;
         }
      }
      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }

      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),(a(rw),dp);
   }
   else
   {
      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
   }

   ideal lastmap;
   q = 1;
   for(ii=1; ii<=size(rw1); ii++ )
   {
      if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; }
      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
   }
   map phi1 = newRing,lastmap;
   ideal endid  = phi1(endid);
   ideal endphi = phi1(endphi);
   export(endid);
   export(endphi);
   list L = lastRing;
   setring newRing;
   return(L);
}
///////////////////////////////////////////////////////////////////////////////
static
proc diagon(ideal i)
{
   matrix m;
   intvec iv = option(get);
   option(redSB);
   ideal j=liftstd(jet(i,1),m);
   option(set,iv);
   return(ideal(matrix(i)*m));
}
/////////////////////////////////////////////////////////////////////////////
/*
Aenderungen:
1. normal kommentiert, bei Berechnung de singulaeren Ortes ein simplify(J,16)
eingefuehrt, um bei riesigen Minorenzahlen, das Ideal zu verkleinern (bis 
Faktor 10 Beschleunigung). 
Protokoll mit printlevel so gesteuert, dass es bei rekursivem Aufruf korrekt 
arbeitet.
list nor = normal(i);     //mit equidim Zerlegung
list nor = normal(i,1);   //mit prim Zerlegung
Zeiten au sony_pumuckel (P2, 500)
                           Examples:
LIB"normal.lib";
//1. Huneke, 1 Komponente
//prim: 2 sec equidim:1sec
ring qr=31991,(a,b,c,d,e),dp;
ideal i=
5abcde-a5-b5-c5-d5-e5,
ab3c+bc3d+a3be+cd3e+ade3,
a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;


//2. Vasconcelos (dauert laenger: 83 sec auf sony_pumuckel)
ring r=32003,(x,y,z,w,t),dp;
ideal i=
x2+zw,
y3+xwt,
xw3+z3t+ywt2,
y2w4-xy2z2t-w3t3;

//2a. Vasconcelos verkleinert
//prim:2Komp, 2 Ringe, 16 sec (manchmal lange, haengt am Faktorisierer)
//equidim: 1 Komp, 7 loops, 2 ringe, 12sec
ring r=32003,(x,y,z,w,t),dp;
ideal i=
x+zw,
y3+xwt,
xw3+z3t+ywt2;

//3. GM1
// irreducible, 13 normalization loops, 1 Ring 
//2sec mit prim, 1 sec mit equidim
ring r=32003,(x,y,z,u),dp;
ideal i = x2+y3,z2+u3,y2+z3;

//3a. GM1
ring r=32003,(x,y),dp;
ideal i = intersect(y3+x2,y2+x3);  // beide 0 sec
ideal i = intersect(y3+x2,y2+x3,y2+x2+x3);
//prim 0sec,
//equidim sehr lange (Relationen in HomJJ zu gross)

//4. GM2
//i nicht reduziert, equidim bricht ab (0 sec)
//prim: 11 irred comp, 11 loops, (1sec)
ring r=32003,(x,y,z,u),dp;
ideal i = x2+y3,z2+u3,y2+z3,x2+z3;

//5. GM3,  radical von GM2
//prim: 11 irred comp, 11 loops, 11 Ringe, (1sec)
//equidim: 1 equidim comp, 1 loop, 1 Ring, (0sec)
ring r=32003,(x,y,z,u),dp;
ideal i =yu+u,yz+z,y2+y,xy+x,x2+y,u3+z2,z3-y;

//GM4
//equidim: 2 equidim comp, 3 Ringe, (0sec);
//prim: 3 Komp, (0sec)
ring r=32003,(x,y,z,u),dp;
ideal i1 = x2+y3,u,z;
ideal i2 = u2+z3,x,y;
ideal i3 = x2+y2+z2+u4;
ideal i = intersect(i1,i2,i3);

//GM4a  Hier dauert prim laenger! (## facstd)
//equidim: 3 equidim Komp, 4 Ringe(0sec)
//prim 13 Komp (63 sec), wegen facstd 
//##ev minAssGTZ(i,1) verwenden (ohne facstd, ist noch fehlerhaft)
ring r=32003,(x,y,z,u),dp;
ideal i1 = x2+y3,u,z;
ideal i2 = u2+z3,x,y;
ideal i3 = x2+y2+z2+u4;
ideal i4 =yu+u,yz+z,y2+y,xy+x,x2+y,u3+z2,z3-y;
ideal i = intersect(i1,i2,i3,i4);

//GM5
//equidim: 4 Komp,0sec, prim 13 Komp, 1sec
ring r=32003,(x,y,z,u,v),dp;
ideal i1 = x2+y3,v;                            //2dim
ideal i2 = u2+z3,x,y;                          //3dim
ideal i3 = x2+y2+z2+u4;                        //4dim
ideal i4 =yu+u,yz+z,y2+y,xy+x,x2+y,u3+z2,z3-y; //1dim
ideal i = intersect(i1,i2,i3,i4);

//cyclic 5
//equidim: 1 Komp 0sec, prim: 20 Komp, 3 sec
ring r=32003,(x,y,z,u,v),dp;
ideal i = 
x+y+z+u+v,
xy+yz+zu+xv+uv,
xyz+yzu+xyv+xuv+zuv,
xyzu+xyzv+xyuv+xzuv+yzuv,
xyzuv-1;

// cyclic 5 hat Normalisierung (1 embim weniger)
///equidim: 1 Komp 0sec, prim: 20 Komp, 2 sec
ring r=32003,(x,y,z,u),dp;
ideal i = 
x2+xz-yz+2xu+yu+u2,
xy2-xyz+y2z-y2u+xzu+yzu+z2u-xu2-2yu2+zu2-u3,
xyz2+xyzu+y2zu-xz2u+yz2u-z3u-xyu2-xzu2-2z2u2+xu3+yu3-zu3+u4,
2xyzu2+y2zu2+2yz2u2-xyu3-2xzu3-yzu3-z2u3+xu4+yu4-2zu4+u5-1;

//cyclic(6)
//equidim: 1Komp in 5 vars 1sec
//prim: 90 (!) Ringe, 12 sec

//Theo1
ring r=32003,(x,y,z),wp(2,3,6);
ideal i=zy2-zx3-x6;

//Theo1a (CohenMacaulay and regular in codimension 2)
ring r=32003,(x,y,z,u),wp(2,3,6,6);
ideal i=zy2-zx3-x6+u2;


//Theo2
ring r=32003,(x,y,z),wp(3,4,12);
ideal i=z*(y3-x4)+x8;

//Theo2a
ring r=32003,(T(1..4)),wp(3,4,12,17);
ideal i=
T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
T(1)^6*T(2)*T(3)+T(1)^2*T(2)^4*T(3)+T(1)^3*T(2)^2*T(4)-T(1)^2*T(2)*T(3)^2+T(4)^2;

//Theo3
ring r=32003,(x,y,z),wp(3,5,15);
ideal i=z*(y3-x5)+x10;


//Theo4
ring r=32003,(x,y,z),dp;
ideal i=(x-y)*(x-z)*(y-z);

//Theo5
ring r=32003,(x,y,z),wp(2,1,2);
ideal i=z3-xy4;

//Theo6
ring r=32003,(x,y,z),dp;
ideal i=x2y2+x2z2+y2z2;

ring r=32003,(a,b,c,d,e,f),dp;
ideal i=
bf,
af,
bd,
ad;

//Sturmfels, wo vorher Prim schneller (2 sec,sonst 860 sec)
//ist CM
//prim: 15 loops, 15 Komp, 1 sec, 
//equidim:15 Ringe, 93 sec mit simplify(J,16), 
//ohne simlify(J,16) 860sec?, 
//andere simplify sind z.T. viel langsamer
ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ideal i=
bv+su,
bw+tu,
sw+tv,
by+sx,
bz+tx,
sz+ty,
uy+vx,
uz+wx,
vz+wy,
bvz;

//J S/Y
ring r=32003,(x,y,z,t),dp;
ideal i=
x2z+xzt,
xyz,
xy2-xyt,
x2y+xyt;

//St_S/Y
ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ideal i=
wy-vz,
vx-uy,
tv-sw,
su-bv,
tuy-bvz;

//Horrocks:
//CHAR 32003:mit prim 1 sec, equidim: 115 sec, beide 6 Ringe 
//           Singulaere Ort hat zu Beginn > 106 000 Erzeuger!!
//char 31991: prim 1sec, 8 Ringe, equidim: 25 Ringe(!), 162 sec, 
//            nicht reduziert!
//            Singulaere Ort hat zu Beginn > 28 000 Erzeuger!!
//            i=radical(i) -> 8 Ringe, 1sec (radical <1 sec)
//Horrocks:
ring r=31991,(a,b,c,d,e,f),dp;
ideal i=
adef-16000be2f+16001cef2,
ad2f+8002bdef+8001cdf2,
abdf-16000b2ef+16001bcf2,
a2df+8002abef+8001acf2,
ad2e-8000bde2-7999cdef,
acde-16000bce2+16001c2ef,
a2de-8000abe2-7999acef,
acd2+8002bcde+8001c2df,
abd2-8000b2de-7999bcdf,
a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
abde-8000b2e2-acdf-16001bcef-8001c2f2,
abcd-16000b2ce+16001bc2f,
a2cd+8002abce+8001ac2f,
a2bd-8000ab2e-7999abcf,
ab3f-3bdf3,
a2b2f-2adf3-16000bef3+16001cf4,
a3bf+4aef3,
ac3e-10668cde3,
a2c2e+10667ade3+16001be4+5334ce3f,
a3ce+10669ae3f,
bc3d+8001cd3e,
ac3d+8000bc3e+16001cd2e2+8001c4f,
b2c2d+16001ad4+4000bd3e+12001cd3f,
b2c2e-10668bc3f-10667cd2ef,
abc2e-cde2f,
b3cd-8000bd3f,
b3ce-10668b2c2f-10667bd2ef,
abc2f-cdef2,
a2bce-16000be3f+16001ce2f2,
ab3d-8000b4e-8001b3cf+16000bd2f2,
ab2cf-bdef2,
a2bcf-16000be2f2+16001cef3,
a4d-8000a3be+8001a3cf-2ae2f2;


ring r=32003,(b,s,t,u,v,w,x,y,z),dp;   //3sec
ideal k=
wy-vz,
vx-uy,
tv-sw,
su-bv,
tuy-bvz;
ideal j=x2y2+x2z2+y2z2;
ideal i=mstd(intersect(j,k))[2];

//22,
// neu, prim: 3 sec, equidim 1 sec, je 4 Ringe
ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
ideal i=
wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;


//riemenschneider, 5 Komponenten
//33(alte Zeiten), normal+primary 3, primary 9, radical 1, minAssGTZ; 2
//neu: prim 0sec, equi 1 sec, je 5 Ringe
ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
ideal i=
xz,
vx,
ux,
su,
qu,
txy,
stx,
qtx,
uv2z-uwz,
uv3-uvw,
puv2-puw;

ring r=0,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;   //0sec
*/