source: git/Singular/LIB/normal.lib @ 0ae4ce

///////////////////////////////////////////////////////////////////////////////
// normal.lib
// algorithms for computing the normalization based on
// the criterion of Grauert/Remmert and ideas of De Jong & Vasconcelos
///////////////////////////////////////////////////////////////////////////////

version="$Id: normal.lib,v 1.27 2000-12-19 14:41:43 anne Exp $";
category="Commutative Algebra";
info="
LIBRARY:  normal.lib     PROCEDURES FOR NORMALIZATION
AUTHORS:  Gert-Martin Greuel, email: greuel@mathematik.uni-kl.de,
          Gerhard Pfister, email: 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";
///////////////////////////////////////////////////////////////////////////////
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
ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
         SBid = standard basis of id,
         J    = ideal of P containing the polynomial p,
         p    = nonzero divisor of R
COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure where
         R is the quotient ring of P modulo the standard basis SBid
RETURN:  a list of two objects
         _[1]: a polynomial ring, containing two ideals, 'endid' and 'endphi'
               s.t. _[1]/endid = Hom_R(J,J) and
               endphi describes the canonical map R -> Hom_R(J,J)
         _[2]: an integer which is 1 if phi is an isomorphism, 0 if not
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;
   intvec rw,rw1;
   list L;
   y = printlevel-voice+2;  // y=printlevel (default: y=0)
 def P = basering;
   ideal SBid, id, J = Li[1], Li[2], Li[3];
   poly p = Li[4];
   attrib(SBid,"isSB",1);
   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not

//---- set attributes for special cases where algorithm can be simplified -----
   if( homo==1 )
   {
      rw = ringweights(P);
   }
   if( typeof(attrib(id,"isPrim"))=="int" )
   {
      if(attrib(id,"isPrim")==1)  { isPr=1; }
   }
   if( typeof(attrib(id,"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";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);
//      export endid;
//      export endphi;
//      L = newR1;
      L=lastRing;
      L = insert(L,1,1);
      dbprint(y,"// case R = Hom(J,J)");
      if(y>=1)
      {
         "//   R=Hom(J,J)";
         "   ";
         lastRing;
         "   ";
         "//   the new ideal";
         endid;
         "   ";
         "//   the old ring";
         "   ";
         P;
         "   ";
         "//   the old ideal";
         "   ";
         setring P;
         id;
         "   ";
         setring lastRing;
         "//   the map";
         "   ";
         endphi;
         "   ";
         pause();
         newline;
      }
      setring P;
      return(L);
   }
   if(y>=1)
   {
      "// R is not equal to Hom(J,J), we have to try again";
      pause();
      newline;
   }
//---------- Hom(J,J) != R: create new ring and map form 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";
      " ";
      Lin;
      "   ";
   }
   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";
      "   ";
      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";
      "   ";
      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 and 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)
RETURN:  a list of rings (say nor), in each ring nor[i] are two ideals
         norid, 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)
NOTE:    to use the i-th ring type: def R=nor[i]; setring R;
         increasing printlevel displays more comments (default: printlevel=0)
COMMENT: 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).
CAUTION: The proc does not check whether the input is radical, for non radical
         ideals the output may be wrong (i=radical(i); makes i radical)
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 ( find(ordstr(basering),"s")+find(ordstr(basering),"M") != 0)
   {
     "";
     "// 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 components";
      }
   }
   else
   {  
      if( typeof(attrib(id,"isPrim"))=="int" )
      {
         if(attrib(id,"isPrim")==1)
         {
            prim[1]=id;
         }
         else
         {   
            prim=minAssPrimes(id);
         }
      }
      else
      {
         prim=minAssPrimes(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);

      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));
      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 (and similair 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, list #)
"USAGE:   normalizationPrimes(i,ihp[,si]);  i prime ideal, ihp map 
         (partial normalization), si SB of singular locus
RETURN:  a list of one ring L=R, in  R are two ideals
         S,M such that R/M is the normalization
         S is a standardbasis of M
NOTE:    to use the ring: def r=L[1];setring r;
         printlevel >= voice+1: display comments (default: printlevel=0)
EXAMPLE: example normalizationPrimes; shows an example
"
{
   int y = printlevel-voice+2;  // y=printlevel (default: y=0)
  
   if(y>=1)
   {
     "";
     "// START a normalization loop with the 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 the input ideal";
   }
   
   list SM=mstd(i);                //here the work starts
   int dimSM =  dim(SM[1]);        //dimension of 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)
   {
      list JM=mstd(#[1]);
      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);
   }

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

   //the 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 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))
   {
      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
   if(((attrib(SM[2],"isIsolatedSingularity")==0)||(homog(SM[2])==0))
        &&(size(#)==0))
   {
/*
write (":a normal-fehler" ,
         "basering:",string(basering),"nvars:", nvars(basering),
       "dim(SM[1]):",dim(SM[1]),"ncols(jacob(SM[2]))",ncols(jacob(SM[2])),
       "SM:", SM);
 
     pause();
*/     
      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
      //ti=timer;
      if(y >=1 )
      {
         "// SB of singular locus will be computed";
      }
      ideal sin=J+SM[2];

    //kills the embeded components

      list JM=mstd(sin);
      //JM[1] SB os singular locus, JM[2]=minbasis of singular locus
      //SM[1] SB of irreducible component, SM[2] minbasis
      if(y>=1)
      {
         "//   the dimension of the singular locus is:";"";
         dim(JM[1]); "";
      }
      attrib(JM[2],"isRad",0);
      //   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)
      {
         attrib(SM[2],"isIsolatedSingularity",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);
        attrib(SM[2],"isRegInCodim2",1);
     }
   }
   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);          
   }
   //if it is an isolated singularity things are easier
   //JM ideal of singular locus, SM ideal of variety
   if((dim(JM[1])==0)&&(homog(SM[2])==1))            //isolated sing. case
   {
      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)
      {
         if(y>=1)
         {
            "";
            "//   the ideal rad(J):";
            "";
            maxideal(1);
            newline;
         }
         //again test for normality
         //Hom(I,R)=R
         list RR;
         RR=SM[1],SM[2],maxideal(1),SL[1];
         ti=timer;
         RR=HomJJ(RR,y);
         if(RR[2]==0)
         {
            def newR=RR[1];
            setring newR;
            map psi=BAS,endphi;
         //   ti=timer;
            list tluser=normalizationPrimes(endid,psi(ihp));

        //    timer-ti;
            setring BAS;
            return(tluser);
         }
         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;
         // the following 2 lines don't work : nor is not defined
         //def R = nor[1]; setring R;     //make the 1-st ring the basering
         //norid; normap;                 //data of the normalization)
         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=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
       {
          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);

          keepresult1=normalizationPrimes(id,ihp);
          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);

          keepresult2=normalizationPrimes(id1,ihp);

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

   //test for non-normality
   //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)
   {
      list RR;
      list RS;
      //now we have to compute the radical
      if(y>=1)
      {
         "// radical computation of singular locus";
      }

      if((attrib(JM[2],"isRad")==0)&&(attrib(SM[2],"isEquidimensional")==0))
      {
           //J=radical(JM[2]);
          J=radical(SM[2]+ideal(SL[1]));

          // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
      }
      if((attrib(JM[2],"isRad")==0)&&(attrib(SM[2],"isEquidimensional")==1))
      {
          ideal JJ=SM[2]+ideal(SL[1]);
         // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
          if(attrib(SM[2],"isCompleteIntersection")==0)
          {
            J=equiRadical(JM[2]);
             //J=equiRadical(JJ);
          }
          else
          {
             //J=radical(JM[2]);
             J=quotient(JJ,minor(jacob(JJ),size(JJ)));
          }
      }
      if(y>=1)
      {
        "//   radical is equal to:";"";
        J;
        "";
      }

      JM=J,J;

      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
      RR=SM[1],SM[2],JM[2],SL[1];

      //   evtl eine geeignete Potenz von JM?
     if(y>=1)
     {
        "// compute Hom(rad(J),rad(J))";
     }

     RS=HomJJ(RR,y);

      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);
      }
      int n=nvars(basering);
      ideal MJ=JM[2];

      def newR=RS[1];
      setring newR;

      map psi=BAS,endphi;
      list tluser=
             normalizationPrimes(endid,psi(ihp),simplify(psi(MJ)+endid,4));
      setring BAS;
      return(tluser);
   }
    // 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);
   }
   else
   {
      int equi=attrib(SM[2],"isEquidimensional");
      ideal new1=qAnn+SM[2];
      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);
      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);

      keepresult1=normalizationPrimes(vid,ihp);

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

      keepresult2=normalizationPrimes(vid,ihp);

      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);
   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));
}
/////////////////////////////////////////////////////////////////////////////
/*
                           Examples:
LIB"normal.lib";
//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;


//Vasconcelos
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);
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;

//Beispiel, wo vorher Primaerzerlegung schneller
//ist CM
//Sturmfels
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;

//dauert laenger
//Horrocks:
ring r=32003,(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;

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
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
//33
//normal+primary 3
//primary 9
//radical 1
//minAssPrimes 2
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;



*/