Changeset 50cbdc in git for Singular/LIB/normal.lib


Ignore:
Timestamp:
Aug 27, 2001, 4:48:02 PM (23 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
2567b5a6cb7109be5a83e53eb94abb1c38fb9945
Parents:
3de58c9ca0aeaafdf5cb29f986967bffa405b542
Message:
*hannes: merge-2-0-2


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/normal.lib

    r3de58c r50cbdc  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: normal.lib,v 1.35 2001-03-19 22:57:16 greuel Exp $";
     2version="$Id: normal.lib,v 1.36 2001-08-27 14:47:56 Singular Exp $";
    33category="Commutative Algebra";
    44info="
    … …  
    88
    99PROCEDURES:
    10  normal(I);             computes the normalization of basering/I
     10 normal(I[,wd]);        computes the normalization of basering/I
     11                        resp. computes the normalization of basering/I and
     12                        the delta-invariante
    1113 HomJJ(L);              presentation of End_R(J) as affine ring, L a list
     14 genus(I);              computes the genus of the projective curve defined
     15                        by I
    1216";
    1317
    … …  
    1923LIB "inout.lib";
    2024LIB "ring.lib";
     25LIB "hnoether.lib";
    2126///////////////////////////////////////////////////////////////////////////////
    22 static proc isR_HomJR (list Li)
    23 "USAGE:   isR_HomJR (Li);  Li = list: ideal SBid, ideal J, poly p
    24 COMPUTE: module Hom_R(J,R) = R:J and compare with R
    25 ASSUME:  R    = P/SBid,  P = basering
    26          SBid = standard basis of an ideal in P,
    27          J    = ideal in P containing the polynomial p,
    28          p    = nonzero divisor of R
    29 RETURN:  1 if R = R:J, 0 if not
    30 EXAMPLE: example isR_HomJR;  shows an example
    31 "
    32 {
    33    int n, ii;
    34  def P = basering;
    35    ideal SBid = Li[1];
    36    ideal J = Li[2];
    37    poly p = Li[3];
    38    attrib(SBid,"isSB",1);
    39    attrib(p,"isSB",1);
    40  qring R    = SBid;
    41    ideal J  = fetch(P,J);
    42    poly p   = fetch(P,p);
    43    ideal f  = quotient(p,J);
    44    ideal lp = std(p);
    45    n=1;
    46    for (ii=1; ii<=size(f); ii++ )
    47    {
    48       if ( reduce(f[ii],lp) != 0)
    49       { n = 0; break; }
    50    }
    51    return (n);
    52  //?spaeter hier einen Test ob Hom(I,R) = Hom(I,I)?
    53 }
    54 example
    55 {"EXAMPLE:";  echo = 2;
    56   ring r   = 0,(x,y,z),dp;
    57   ideal id = y7-x5+z2;
    58   ideal J  = x3,y+z;
    59   poly p   = xy;
    60   list Li  = std(id),J,p;
    61   isR_HomJR (Li);
    62 
    63   ring s   = 0,(t,x,y),dp;
    64   ideal id = x2-y2*(y-t);
    65   ideal J  = jacob(id);
    66   poly p   = J[1];
    67   list Li  = std(id),J,p;
    68   isR_HomJR (Li);
    69 }
    70 ///////////////////////////////////////////////////////////////////////////////
    7127
    7228proc HomJJ (list Li)
    73 "USAGE:   HomJJ (Li); Li list: ideal SBid, ideal id, ideal J, poly p, int count
     29"USAGE:   HomJJ (Li);  Li = list: ideal SBid, ideal id, ideal J, poly p
    7430ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
    7531@*       SBid = standard basis of id,
    7632@*       J    = ideal of P containing the polynomial p,
    7733@*       p    = nonzero divisor of R
    78 @*       count controls printlevel during recursive call
    7934COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
    8035         affine ring, together with the canonical map R --> Hom_R(J,J),
    … …  
    8641               endphi describes the canonical map R -> Hom_R(J,J)
    8742         l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
     43         l[3] : an integer, the contribution to delta
    8844@end format
    8945NOTE:    printlevel >=1: display comments (default: printlevel=0)
    … …  
    9349//---------- initialisation ---------------------------------------------------
    9450
    95    int isIso,isPr,isCo,isRe,isEq,ii,jj,q,y,count;
     51   int isIso,isPr,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
    9652   intvec rw,rw1;
    9753   list L;
    98    if(size(Li)>=5)
    99    {
    100      count = Li[5];
    101    }
    102    y = printlevel-voice+count; 
    103  def P = basering;
     54   y = printlevel-voice+2;  // y=printlevel (default: y=0)
     55   def P = basering;
    10456   ideal SBid, id, J = Li[1], Li[2], Li[3];
    10557   poly p = Li[4];
    … …  
    11668      if(attrib(id,"isPrim")==1)  { isPr=1; }
    11769   }
     70   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
     71   {
     72      if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; }
     73   }
    11874   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
    11975   {
    … …  
    13389   }
    13490//-------------------------- go to quotient ring ------------------------------
    135  qring R  = SBid;
     91   qring R  = SBid;
    13692   ideal id = fetch(P,id);
    13793   ideal J  = fetch(P,J);
    … …  
    13995   ideal f,rf,f2;
    14096   module syzf;
    141 
    14297//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
    14398   if ( y>=1 )
    … …  
    149104   }
    150105   f  = quotient(p*J,J);
     106
    151107   if ( y>=1 )
    152108   { "// the module p*Hom(J,J) = p*J:J, p a non-zerodivisor";
    153109      "// p"; p;
    154       "// f=p*J:J";
    155       if( y>=2 ) { f; }     
     110      "// f=p*J:J";f;
    156111   }
    157112   f2 = std(p);
    158 
    159    if(isIso==0)
    160    {
    161      ideal f1=std(f);
    162      attrib(f1,"isSB",1);
    163     // if( codim(f1,f2) >= 0 )
    164     // {
    165     //  dbprint(printlevel-voice+3,"// dimension of non-normal locus is zero");
    166     //    isIso=1;
    167     // }
    168   }
     113 
    169114//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
    170115
    171116   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R
     117
    172118   if ( size(rf) == 0 )
    173119   {
    … …  
    186132      setring lastRing;
    187133
     134      attrib(endid,"onlySingularAtZero",oSAZ);
    188135      attrib(endid,"isCohenMacaulay",isCo);
    189136      attrib(endid,"isPrim",isPr);
    … …  
    195142      L=lastRing;
    196143      L = insert(L,1,1);
     144      dbprint(y,"// case R = Hom(J,J)");
    197145      if(y>=1)
    198146      {
    199          "// case R = Hom(J,J), we are ready with this component";
     147         "//   R=Hom(J,J)";
    200148         "   ";
    201        if( y>=2 )
    202        {
    203         lastRing;
     149         lastRing;
    204150         "   ";
    205151         "//   the new ideal";
    206          if( y>=2 ) { endid; }
     152         endid;
    207153         "   ";
    208154         "//   the old ring";
    … …  
    222168         pause();
    223169         newline;
    224        }
    225170      }
    226171      setring P;
     172      L[3]=0;
    227173      return(L);
    228174   }
    … …  
    230176   {
    231177      "// R is not equal to Hom(J,J), we have to try again";
    232     if( y>=2 )
    233     {  pause();
    234        newline;
    235     }   
     178      pause();
     179      newline;
    236180   }
    237181//---------- Hom(J,J) != R: create new ring and map from old ring -------------
    238182// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
    239183
     184   f=mstd(f)[2];
     185   ideal ann=quotient(f2,f);
     186   int delta;
     187   if(isIso&&isEq){delta=vdim(std(modulo(f,ideal(p))));}
     188
    240189   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
    241190   q = size(f);
     191
    242192   syzf = syz(f);
    243193
    … …  
    261211   attrib(SBid,"isSB",1);
    262212
    263  qring newR = std(SBid);
     213   qring newR = std(SBid);
     214
    264215   map psi = R,ideal(X(1..nvars(R)));
    265216   ideal id = psi(id);
    … …  
    286237      "//   the linear relations";
    287238      " ";
    288       if( y>=2 )
    289       {  Lin;     
    290          pause();
    291          "";
    292       }
    293    }
     239      Lin;
     240      "   ";
     241   }
     242
    294243   for (ii=2; ii<=q; ii++ )
    295244   {
    … …  
    300249      }
    301250   }
     251
    302252   if(y>=1)
    303253   {
    304254      "//   the quadratic relations";
    305         if( y>=2 )
    306       {     
    307           "   ";
    308          interred(Quad);
    309          pause();
    310          newline;
    311       }
    312    }
    313    Q = Lin+Quad;
     255      "   ";
     256      interred(Quad);
     257      pause();
     258      newline;
     259   }
     260   Q = Lin,Quad;
    314261   Q = subst(Q,T(1),1);
    315    Q = interred(reduce(Q,std(0)));
     262   Q=mstd(Q)[2];
    316263//---------- reduce number of variables by substitution, if possible ----------
    317264   if (homo==1)
    … …  
    324271   }
    325272
    326    ideal endid  = imap(newR,id)+imap(newR,Q);
     273   ideal endid  = imap(newR,id),imap(newR,Q);
    327274   ideal endphi = ideal(X(1..nvars(R)));
    328275
    329276   L=substpart(endid,endphi,homo,rw);
     277
    330278   def lastRing=L[1];
    331279   setring lastRing;
     280
     281   attrib(endid,"onlySingularAtZero",0);
     282   map sigma=R,endphi;
     283   ideal an=sigma(ann);
     284   export(an);  //noetig?
     285   ideal te=an,endid;
     286   if(isIso&&(size(reduce(te,std(maxideal(1))))==0))
     287   {
     288      attrib(endid,"onlySingularAtZero",oSAZ);
     289   }
     290   kill te;
    332291   attrib(endid,"isCohenMacaulay",isCo);
    333292   attrib(endid,"isPrim",isPr);
    … …  
    337296   attrib(endid,"isCompleteIntersection",0);
    338297   attrib(endid,"isRad",0);
    339   // export(endid);
    340   // export(endphi);
    341298   if(y>=1)
    342299   {
    343300      "//   the new ring after reduction of the number of variables";
    344       show(lastRing);
    345       pause(); "
    346       ";
    347     if(y >=2)
    348     {
     301      "   ";
    349302      lastRing;
    350303      "   ";
    … …  
    369322      pause();
    370323      newline;
    371     }
    372324   }
    373325   L = lastRing;
    374326   L = insert(L,0,1);
     327   L[3]=delta;
    375328   return(L);
    376329}
    … …  
    394347
    395348proc normal(ideal id, list #)
    396 "USAGE:   normal(i [,choose]);  i a radical ideal, choose empty or 1
     349"USAGE:   normal(i [,choose]);  i a radical ideal, choose empty, 1 or "wd"
    397350         if choose=1 the normalization of the associated primes is computed
    398351         (which is sometimes more efficient)
     352         if choose="wd" the delta-invariant is computed simultaneously
     353         this may take much more time in the reducible case because the
     354         factorizing standardbasis algorithm cannot be used.
    399355ASSUME:  The ideal must be radical, for non radical ideals the output may
    400356         be wrong (i=radical(i); makes i radical)
    401 RETURN:   a list of rings, say nor:
     357RETURN:   a list of rings, say nor and in case of choose="wd" an integer
    402358@format
     359         at the end of the list.
    403360         each ring nor[i] contains two ideals
    404361         with given names norid and normap such that
    … …  
    409366@end format
    410367NOTE:    to use the i-th ring type: def R=nor[i]; setring R;.
    411 @*       Not implemented for local or mixed orderings and quotient rings.
     368@*       Increasing printlevel displays more comments (default: printlevel=0)
     369@*       Not implemented for local or mixed orderings.
    412370@*       If the input ideal i is weighted homogeneous a weighted ordering may
    413371         be used (qhweight(i); computes weights).
    414 @*       printlevel = 1: count normalization loops (default: printlevel=0)
    415 @*       printlevel > 1: protocoll of normalization steps
    416 @*       printlevel > 2: protocoll of all normalization steps, pauses
    417 @*       printlevel > 3: display some (>4 all) intermediate results, pauses
    418372EXAMPLE: example normal; shows an example
    419373"
    420374{
    421    int i,j,y;
     375   int i,j,y,withdelta;
    422376   string sr;
    423377   list result,prim,keepresult;
    424378   y = printlevel-voice+2;
    425 
     379   if(size(#)>0)
     380   {
     381     if(typeof(#[1])=="string")
     382     {
     383       kill #;
     384       list #;
     385       withdelta=1;
     386     }
     387   }
    426388   attrib(id,"isRadical",1);
    427389   if ( ord_test(basering) != 1)
    … …  
    473435      if(y>=1)
    474436      {
    475         "";
    476          "// we have ",size(prim),"equidimensional component(s)";
     437         "// we have ",size(prim),"equidimensional components";
     438      }
     439      if(withdelta &&(size(prim)>1))
     440      {
     441         withdelta=0;
     442         "WARNING!  cannot compute delta,";
     443         "the ideal is not equidimensional";
     444      }
     445      if(!withdelta)
     446      {
     447         option(redSB);
     448         for(j=1;j<=size(prim);j++)
     449         {
     450            keepresult=keepresult+facstd(prim[j]);
     451         }
     452         prim=keepresult;
     453         option(noredSB);
    477454      }
    478455   }
    … …  
    487464         else
    488465         {
    489             prim=minAssGTZ(id);
     466            prim=minAssPrimes(id);
    490467         }
    491468      }
    492469      else
    493470      {
    494          prim=minAssGTZ(id);
     471         prim=minAssPrimes(id);
    495472      }
    496473      if(y>=1)
    497474      {
    498          "";
    499          "// we have ",size(prim),"irreducible component(s)";
     475         "// we have ",size(prim),"irreducible components";
    500476      }
    501477   }
    … …  
    504480      if(y>=1)
    505481      {
    506          "";
    507          "// BEGIN with equidimensional/irreducible component",i;
     482         "// we are in loop ",i;
    508483      }
    509484      attrib(prim[i],"isCohenMacaulay",0);
    … …  
    520495      attrib(prim[i],"isEquidimensional",1);
    521496      attrib(prim[i],"isCompleteIntersection",0);
     497      attrib(prim[i],"onlySingularAtZero",0);
     498      if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
     499      {
     500            if(attrib(id,"onlySingularAtZero")==1)
     501             {attrib(prim[i],"onlySingularAtZero",1); }
     502      }
    522503
    523504      if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
    … …  
    533514      }
    534515      keepresult=normalizationPrimes(prim[i],maxideal(1),0);
    535       for(j=1;j<=size(keepresult);j++)
     516
     517      for(j=1;j<=size(keepresult)-1;j++)
    536518      {
    537519         result=insert(result,keepresult[j]);
    538520      }
    539       sr = string(size(result));
     521      sr = string(size(result));     
     522   }
     523   if(withdelta)
     524   {
     525      sr = string(size(keepresult)-1);
     526      result=keepresult;
    540527   }
    541528      dbprint(y+1,"
    542529// 'normal' created a list of "+sr+" ring(s).
     530// nor["+sr+"+1] is the delta-invariant in case of choose=wd.
    543531// To see the rings, type (if the name of your list is nor):
    544      show(nor);
     532     show( nor);
    545533// To access the 1-st ring and map (similar for the others), type:
    546534     def R = nor[1]; setring R;  norid; normap;
    … …  
    561549   norid;
    562550   normap;
     551
     552   ring s=0,(x,y),dp;
     553   ideal i=(x-y^2)^2 - y*x^3;
     554   nor=normal(i,"wd");
     555   //the delta-invariant
     556   nor[size(nor)];
    563557}
    564558
    565559///////////////////////////////////////////////////////////////////////////////
    566 static proc normalizationPrimes(ideal i,ideal ihp, int count, list #)
    567 "USAGE:   normalizationPrimes(i,ihp[,si,countt]);  i prime ideal, ihp map
    568          (partial normalization), si ideal of singular locus,
    569          count = integer to count the number of normalization loops
    570 RETURN:  a list of one ring L=R, in  R are two ideals
    571          S,M such that R/M is the normalization of original basering
    572          S is a standardbasis of M
    573 NOTE:    to use the ring: def r=L[1];setring r;
    574          printlevel = 1: count normalization loops (default: printlevel=0)
    575          printlevel > 1: protocoll of normalization steps
    576          printlevel > 2: protocoll of all normalization steps, pauses
    577          printlevel > 3: display some (>4 all) intermediate results, pauses
     560static proc normalizationPrimes(ideal i,ideal ihp,int delta, list #)
     561"USAGE:   normalizationPrimes(i,ihp[,si]);  i equidimensional ideal, ihp map
     562         (partial normalization),delta partial delta-invariant,
     563          # ideal of singular locus
     564RETURN:   a list of rings, say nor, and an integer, the delta-invariant
     565          at the end of the list.
     566          each ring nor[i] contains two ideals
     567          with given names norid and normap such that
     568           - the direct sum of the rings nor[i]/norid is
     569             the normalization of basering/id;
     570           - normap gives the normalization map from basering/id
     571             to nor[i]/norid (for each i)
    578572EXAMPLE: example normalizationPrimes; shows an example
    579573"
    580574{
    581    count = count+1;   
    582    int y = printlevel-voice+count+1;  // y=printlevel (default: y=0)
    583    if(y>=0)
    584    {
    585      "// START normalization LOOP ",count;
    586    }
    587    if( y>=3)
    588    {
    589      "// with ideal";  "";
     575
     576   int y = printlevel-voice+2;  // y=printlevel (default: y=0)
     577
     578   if(y>=1)
     579   {
     580     "";
     581     "// START a normalization loop with the ideal";  "";
    590582     i;  "";
    591583     basering;  "";
    … …  
    593585     newline;
    594586   }
    595 
    596587
    597588   def BAS=basering;
    … …  
    614605         export normap;
    615606         result=newR7;
     607         result[size(result)+1]=delta;
    616608         setring BAS;
    617609         return(result);
    … …  
    620612   if(y>=1)
    621613   {
    622    "// SB-computation of ideal of size",size(i),"in ring with",nvars(basering),"variables";
    623    }
    624 // -------------- SB-computation of the input ideal ----------------------
    625    list SM=mstd(i);                //here the work starts, SM[1] a SB of i
    626                                    //SM[2] a smaller set of generators of i
    627    int dimSM =  dim(SM[1]);        //dimension of i, V(i)=variety to normalize
     614     "// SB-computation of the input ideal";
     615   }
     616   list SM=mstd(i);                //here the work starts
     617   int dimSM =  dim(SM[1]);        //dimension of variety to normalize
    628618  // Case: Get an ideal containing a unit
    629619   if( dimSM == -1)
    … …  
    643633         export normap;
    644634         result=newR6;
     635         result[size(result)+1]=delta;
    645636         setring BAS;
    646637         return(result);
    … …  
    652643      dimSM;"";
    653644   }
    654 
    655    if(size(#)>0)                   //ideal of sing locus is given in #[1]
    656    {
    657       list JM=mstd(#[1]);
    658       if( typeof(attrib(#[1],"isRad"))!="int" )
    659       {
    660          attrib(JM[2],"isRad",0);
    661       }
    662    }
    663 
    664 // -------- transfer attributes from ideal i to SM[2], SM = mstd(i) --------
     645   if(size(#)>0)
     646   {
     647      if(attrib(i,"onlySingularAtZero"))
     648      {
     649         list JM=maxideal(1),maxideal(1);
     650         attrib(JM[1],"isSB",1);
     651         attrib(JM[2],"isRad",1);
     652      }
     653      else
     654      {
     655         ideal te=#[1],SM[2];
     656         list JM=mstd(te);
     657         kill te;
     658         if( typeof(attrib(#[1],"isRad"))!="int" )
     659         {
     660            attrib(JM[2],"isRad",0);
     661         }
     662      }
     663   }
     664
    665665   if(attrib(i,"isPrim")==1)
    666666   {
    … …  
    703703      attrib(SM[2],"isEquidimensional",0);
    704704   }
    705     if(attrib(i,"isCompleteIntersection")==1)
     705   if(attrib(i,"isCompleteIntersection")==1)
    706706   {
    707707      attrib(SM[2],"isCompleteIntersection",1);
    … …  
    711711      attrib(SM[2],"isCompleteIntersection",0);
    712712   }
    713 
    714 //************* case 0: check and prepare special cases ****************
    715 //no computation for the normalization in case 0
    716 
    717 //-------------------------- the smooth case ----------------------------
     713   if(attrib(i,"onlySingularAtZero")==1)
     714   {
     715      attrib(SM[2],"onlySingularAtZero",1);
     716   }
     717   else
     718   {
     719      attrib(SM[2],"onlySingularAtZero",0);
     720   }
     721   if((attrib(SM[2],"isIsolatedSingularity")==1)&&(homog(SM[2])==1))
     722   {
     723      attrib(SM[2],"onlySingularAtZero",1);
     724   }
     725
     726   //the smooth case
    718727   if(size(#)>0)
    719728   {
    … …  
    725734         }
    726735         MB=SM[2];
    727          intvec rw;
     736         intvec rw;;
    728737         list LL=substpart(MB,ihp,0,rw);
    729738         def newR6=LL[1];
    … …  
    735744         export normap;
    736745         result=newR6;
     746         result[size(result)+1]=delta;
    737747         setring BAS;
    738748         return(result);
    739749     }
    740750   }
    741    // recall: SM = mstd(i), i = ideal to start with in normaliztion loop
    742    //         JM = mstd(#[1]), #[1]= ideal of singular locus of i
    743    //              #[1] is given after the first normalization loop
    744 
    745 //---------------- the homogeneous zero-dimensional case ----------------
     751
     752   //the zero-dimensional case
    746753   if((dim(SM[1])==0)&&(homog(SM[2])==1))
    747754   {
    … …  
    761768      export normap;
    762769      result=newR5;
     770      result[size(result)+1]=delta;
    763771      setring BAS;
    764772      return(result);
    765773   }
    766774
    767 //------------- the homogeneous one-dimensional case -------------------
    768 //in this case i defines a line because it is irreducible and homogeneous
     775   //the one-dimensional case
     776   //in this case it is a line because
     777   //it is irreducible and homogeneous
    769778   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
    770779        &&(homog(SM[2])==1))
    … …  
    785794      export normap;
    786795      result=newR4;
     796      result[size(result)+1]=delta;
    787797      setring BAS;
    788798      return(result);
    789799   }
    790 //----------------------- the higher dimensional case ----------------------
     800
     801   //the higher dimensional case
    791802   //we test first of all CohenMacaulay and
    792803   //complete intersection
    … …  
    802813      }
    803814   }
    804 
    805    //compute the singular locus+lower dimensional components
    806 
    807 //------- case: not(isolated singularity and homogeneous) -------------------
    808    if((attrib(SM[2],"isIsolatedSingularity")==0) && (size(#)==0))
    809    {
    810 //---------- computation of ideal of singular locus -------------------------
     815   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(i)==1))
     816   {
     817      int tau=vdim(std(i+jacob(i)));
     818      if(tau>=0)
     819      {
     820        execute("ring BASL="+charstr(BAS)+",("+varstr(BAS)+"),ds;");
     821        ideal i=imap(BAS,i);
     822        int tauloc=vdim(std(i+jacob(i)));
     823        setring BAS;
     824        attrib(SM[2],"onlySingularAtZero",(tau==tauloc));
     825        kill BASL;
     826      }
     827   }
     828
     829   //compute the singular locus
     830   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(#)==0))
     831   {
    811832      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
    812       //ti=timer;
    813833      if(y >=1 )
    814834      {
    815          "// SB of singular locus will be computed ";
    816          if(y >=2 )
    817          {
    818            "//in ring:";
    819            show(basering);
    820          }
    821       }
    822      
    823       J = simplify(J,16);  //this makes ist for huge J much faster
    824       ideal sin=J,SM[2];
     835         "// SB of singular locus will be computed";
     836      }
     837      ideal sin=J,SM[1];
    825838      list JM=mstd(sin);
    826      
    827       //JM[1] SB of singular locus, JM[2] minbasis of singular locus
    828       //SM[1] SB of ideal in normalisation loop, SM[2] minbasis
    829 
     839      attrib(JM[1],"isSB",1);
     840
     841      //JM[1] SB of singular locus, JM[2]=minbasis of singular locus
     842      //SM[1] SB of the input ideal, SM[2] minbasis
    830843      if(y>=1)
    831844      {
    … …  
    833846         dim(JM[1]); "";
    834847      }
    835       //   timer-ti;
    836       attrib(JM[1],"isSB",1);
     848
    837849      if(dim(JM[1])==-1)
    838850      {
    … …  
    852864         export normap;
    853865         result=newR3;
     866         result[size(result)+1]=delta;
    854867         setring BAS;
    855868         return(result);
    856869      }
    857       if(dim(JM[1])==0 && (homog(SM[2])==1))
     870      if(dim(JM[1])==0)
    858871      {
    859872         attrib(SM[2],"isIsolatedSingularity",1);
     873         if(homog(SM[2])){attrib(SM[2],"onlySingularAtZero",1);}
    860874      }
    861875      if(dim(JM[1])<=dim(SM[1])-2)
    … …  
    864878      }
    865879   }
    866    
    867 //------------ case: isolated singularity and homogeneous -----------------
    868    if(attrib(SM[2],"isIsolatedSingularity")==1)   
     880   else
    869881   {
    870882     if(size(#)==0)
    871883     {
    872         if(defined(JM)==voice)
    873         {  JM=maxideal(1),maxideal(1);
    874         }
    875         else
    876         {  list JM=maxideal(1),maxideal(1);
    877        
    878         }
     884        list JM=maxideal(1),maxideal(1);
     885
    879886        attrib(JM[1],"isSB",1);
     887        if(dim(SM[1])>=2){attrib(SM[2],"isRegInCodim2",1);}
    880888     }
    881889   }
    882 
    883 //------------ case: Cohen Macaulay and regular in codim 2 -----------------
    884 //in this case we are ready, the ring is normal 
    885890   if((attrib(SM[2],"isRegInCodim2")==1)&&(attrib(SM[2],"isCohenMacaulay")==1))
    886891   {
    … …  
    900905      export normap;
    901906      result=newR6;
     907      result[size(result)+1]=delta;
    902908      setring BAS;
    903909      return(result);
    904910   }
    905911
    906 //************* case 1: isolated singularity and homogeneous ****************
    907 
    908    // recall: SM = mstd(i), i = ideal to start with in normaliztion loop
    909    //         JM = mstd(#[1]), #[1]= ideal of singular locus of i
    910    //              #[1] is given after the first normalization loop
    911    //if SM is an isolated singularity and homogeneous, we know that this
    912    //persists in the following normalization loops and things are easier
    913    //since the radical of the singular locus is the maximal ideal
     912   //if it is an isolated singularity only at 0 things are easier
    914913   //JM ideal of singular locus, SM ideal of variety
    915 
    916    if(attrib(SM[2],"isIsolatedSingularity")==1)
    917    {
    918       if(y>=1)
    919       {
    920          "// CASE 1: unique isolated singularity at 0";
    921          "// radial of singular locus is the maximal ideal";
    922       }
     914   if(attrib(SM[2],"onlySingularAtZero"))         
     915   {
     916      attrib(SM[2],"isIsolatedSingularity",1);
    923917      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
    924       //SL = vars not contained in ideal
     918      //vars not contained in ideal
    925919      ideal Ann=quotient(SM[2],SL[1]);
    926       ideal qAnn=simplify(reduce(Ann,SM[1]),2);     
     920      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
    927921      //qAnn=0 ==> the first var(=SL[1]) not contained in SM is a nzd of R/SM
    928   //--------------- case: a nonzero-divisor was found ---------------
    929922      if(size(qAnn)==0)
    930923      {
    … …  
    932925         {
    933926            "";
    934             "//   a nonzero-divisor was found";
    935927            "//   the ideal rad(J):";
    936928            "";
    937929            maxideal(1);
    938930            newline;
    939             "// TEST for normality: R=Hom(J,J)?";
    940             "";
    941931         }
    942     //test for normality:
    943          //?spaeter: test for normality, Hom(I,R)==R?
     932         //again test for normality
     933         //Hom(I,R)=R
    944934         list RR;
    945          RR=SM[1],SM[2],maxideal(1),SL[1],count;
    946          //  ti=timer;
    947          RR=HomJJ(RR);
    948          //  timer-ti;
    949          if(RR[2]==0)
    950     //the ring is not normal, start a new normalization loop
     935         RR=SM[1],SM[2],maxideal(1),SL[1];
     936
     937         RR=HomJJ(RR,y);
     938         if(RR[2]==0)
    951939         {
    952940            def newR=RR[1];
    953941            setring newR;
    954942            map psi=BAS,endphi;
    955          //  ti=timer;
    956             list tluser=normalizationPrimes(endid,psi(ihp),count);
    957          //  timer-ti;
     943            list tluser=normalizationPrimes(endid,psi(ihp),delta+RR[3],an);
    958944            setring BAS;
    959945            return(tluser);
    960946         }
    961     //the ring is normal, prepare output and go home
    962947         MB=SM[2];
    963948         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
    … …  
    965950         ideal norid=fetch(BAS,MB);
    966951         ideal normap=fetch(BAS,ihp);
     952         delta=delta+RR[3];
    967953         export norid;
    968954         export normap;
    969955         result=newR7;
     956         result[size(result)+1]=delta;
    970957         setring BAS;
    971958         return(result);
    972959
    973960       }
    974   //--------------- case: a zero-divisor was found ---------------
    975       //Now the case where qAnn!=0, i.e.SL[1] is a zero-divisor of R/SM
     961      //Now the case where qAnn!=0, i.e.SL[1] is a zero divisor of R/SM
    976962      //and we have found a splitting: id and id1
    977       //id=qAnn+SM[2] defines components of R/SM in the complement of V(SL[1])
     963      //id=Ann defines components of R/SM in the complement of V(SL[1])
    978964      //id1 defines components of R/SM in the complement of V(id)
    979       //?????instead of id1 we can take SL[1]+Ann+SM[2]???????????
    980965       else
    981966       {
    982          if(y>=0)
    983          {
    984             "//   a zero-divisor was found, we have SPLITTING of ideals";
    985            if(y >=2 )
    986            {
    987              "// in ring:";
    988              show(basering);
    989             }
    990             "";
    991          }
    992          
    993           ideal id=qAnn+SM[2];
    994 
     967          ideal id=Ann;
    995968          attrib(id,"isCohenMacaulay",0);
    996969          attrib(id,"isPrim",0);
    … …  
    999972          attrib(id,"isCompleteIntersection",0);
    1000973          attrib(id,"isEquidimensional",0);
    1001 
    1002           if(y>=0)
    1003           {
    1004     "//   start a normalization loop with 1st split ideal (size",size(id),"in",nvars(basering),"vars)";
    1005           }
    1006           keepresult1=normalizationPrimes(id,ihp,count);
    1007           ideal id1=quotient(SM[2],Ann)+SM[2];
    1008 //          evtl. qAnn statt Ann nehmen
    1009 //          ideal id=SL[1]+SM[2];
    1010 
     974          attrib(id,"onlySingularAtZero",1);
     975
     976          ideal id1=quotient(SM[2],Ann);
     977          //ideal id=SL[1],SM[2];
    1011978          attrib(id1,"isCohenMacaulay",0);
    1012979          attrib(id1,"isPrim",0);
    … …  
    1015982          attrib(id1,"isCompleteIntersection",0);
    1016983          attrib(id1,"isEquidimensional",0);
    1017 
    1018           if(y>=0)
    1019           {
    1020               "//   start a normalization loop with 2nd split ideal (size",size(id),"in",nvars(basering),"vars)";
    1021           }
    1022           keepresult2=normalizationPrimes(id1,ihp,count);
    1023 
    1024           for(lauf=1;lauf<=size(keepresult2);lauf++)
     984          attrib(id1,"onlySingularAtZero",1);
     985         
     986          ideal t1=id,id1;
     987          int mul=vdim(std(t1));
     988          kill t1;
     989
     990          keepresult1=normalizationPrimes(id,ihp,0);
     991
     992          keepresult2=normalizationPrimes(id1,ihp,0);
     993           
     994          delta=delta+mul+keepresult1[size(keepresult1)]
     995                         +keepresult1[size(keepresult1)];
     996
     997          for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
    1025998          {
    1026999             keepresult1=insert(keepresult1,keepresult2[lauf]);
    10271000          }
     1001          keepresult1[size(keepresult1)]=delta;
    10281002          return(keepresult1);
    10291003       }
    10301004   }
    10311005
    1032 //********** case 2: no unique isolated singularity at 0 *************
    1033 
    1034    //in this case the radical must be computed
    1035    // recall: SM = mstd(i), i = ideal to start with in normaliztion loop
    1036    //         JM = mstd(#[1]), #[1]= ideal of singular locus of i
    1037    //              #[1] is given after the first normalization loop
    1038    //test for normality:  Hom(J,J)=R
    1039    //test for non-normality: Hom(J,J)!=R, we can use Hom(J,J) to continue
    1040 
    1041       if(y>=1)
    1042       {
    1043          "// CASE 2: no unique isolated singularity";
    1044          "// radical has to be computed";
    1045       }
    1046 
     1006   //test for non-normality
     1007   //Hom(I,I)<>R
     1008   //we can use Hom(I,I) to continue
    10471009   ideal SL=simplify(reduce(JM[2],SM[1]),2);
    1048    //SL = elements of ideal of singular locus not contained in ideal i
    10491010   ideal Ann=quotient(SM[2],SL[1]);
    10501011   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
    1051    //qAnn=0 ==> SL[1] is a nonzero-divisor of R/SM
    1052 
    1053   //--------------- case: a nonzero-divisor was found ---------------
     1012
    10541013   if(size(qAnn)==0)
    10551014   {
    10561015      list RR;
    10571016      list RS;
    1058     //now we have to compute the radical
     1017      //now we have to compute the radical
    10591018      if(y>=1)
    10601019      {
     1020         "// radical computation of singular locus";
     1021      }
     1022
     1023      J=radical(JM[2]);   //the singular locus
     1024      JM=mstd(J);
     1025     
     1026      if((vdim(JM[1])==1)&&(size(reduce(J,std(maxideal(1))))==0))
     1027      {
     1028         attrib(SM[2],"onlySingularAtZero",1);
     1029      }
     1030      if(y>=1)
     1031      {
     1032        "//   radical is equal to:";"";
     1033        J;
    10611034        "";
    1062          "//   a nonzero-divisor was found";
    1063          "//   radical computation of ideal of singular locus";
    1064       }
    1065      
    1066     //-------------- computation of the radical J --------------------
    1067     //We have at least two possibilities:
    1068     //J=radical(JM[2]), the radical of the full singular locus, or
    1069     //J=radical(SM[2]+ideal(SL[1])), JM[2] contains SM[2]+ideal(SL[1])
    1070       //the first is usually better!
    1071 
    1072       if(y>=1)
    1073       {
    1074         "// compute radical J of ideal of singular locus";"";
    1075       }
    1076 
    1077       J=radical(JM[2]);
    1078       // alternativ:   J=radical(SM[2]+ideal(SL[1]));
    1079 
    1080       if(y>=2)
    1081       {
    1082           "// the radical is equal to:";       
    1083           J;
    1084           newline;
    1085       }
    1086       if(y>=1)
    1087       {       
    1088         "// TEST for normality: R=Hom(J,J)?";
    1089         "";
    1090       }
    1091 
    1092       JM=J,J;              //J = new ideal for singular locus
     1035      }
     1036
     1037      if(deg(SL[1])>deg(J[1]))
     1038      {
     1039         Ann=quotient(SM[2],J[1]);
     1040         qAnn=simplify(reduce(Ann,SM[1]),2);
     1041         if(size(qAnn)==0){SL[1]=J[1];}
     1042      }
     1043
    10931044      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
    1094 
    1095     //test for normality:
    1096       RR=SM[1],SM[2],JM[2],SL[1],count;
    1097       RS=HomJJ(RR);
    1098      
    1099     //--- the ring is normal, prepare output and go home
     1045      RR=SM[1],SM[2],JM[2],SL[1];
     1046
     1047      //   evtl eine geeignete Potenz von JM?
     1048     if(y>=1)
     1049     {
     1050        "// compute Hom(rad(J),rad(J))";
     1051     }
     1052
     1053     RS=HomJJ(RR,y);
     1054
    11001055      if(RS[2]==1)
    11011056      {
    … …  
    11081063         export norid;
    11091064         export normap;
     1065         result=lastR;
     1066         result[size(result)+1]=delta+RS[3];
    11101067         setring BAS;
    1111          return(lastR);
    1112       }
    1113 
    1114     //--- the ring is not normal, start a new normalization loop       
     1068         return(result);
     1069      }
    11151070      int n=nvars(basering);
    11161071      ideal MJ=JM[2];
    … …  
    11211076      map psi=BAS,endphi;
    11221077      list tluser=
    1123            normalizationPrimes(endid,psi(ihp),count,simplify(psi(MJ)+endid,4));
     1078             normalizationPrimes(endid,psi(ihp),delta+RS[3],psi(MJ));
    11241079      setring BAS;
    11251080      return(tluser);
    11261081   }
    1127   //--------------- case: a zero-divisor was found ---------------
    1128   //Now the case where qAnn!=0, i.e.SL[1] is a zero-divisor of R/SM
    1129   //and we have found a splitting: id and id1
    1130   //id=qAnn+SM[2] defines components of R/SM in the complement of V(SL[1])
    1131   //id1 defines components of R/SM in the complement of V(id)
    1132       //?????instead of id1 we can take SL[1]+Ann+SM[2]???????????
    1133 
    1134    // A component with singular locus the whole component found:
     1082    // A component with singular locus the whole component found
    11351083   if( Ann == 1)
    11361084   {
    … …  
    11511099         export normap;
    11521100         result=newR6;
     1101         result[size(result)+1]=delta;
    11531102         setring BAS;
    11541103         return(result);
    11551104   }
    1156    // The general case with splitting of ring, i.e. qAnn!=0
    11571105   else
    11581106   {
    11591107      int equi=attrib(SM[2],"isEquidimensional");
    1160       ideal new1=qAnn+SM[2];
     1108      int oSAZ=attrib(SM[2],"onlySingularAtZero");
     1109      int isIs=attrib(SM[2],"isIsolatedSingularity");
     1110
     1111      ideal new1=Ann;
     1112      ideal new2=quotient(SM[2],Ann);
     1113      //ideal new2=SL[1],SM[2];
     1114
     1115      int mul;
     1116      if(equi&&isIs)
     1117      {
     1118          ideal t2=new1,new2;
     1119          mul=vdim(std(t2));
     1120      }
    11611121      execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
    11621122                      +ordstr(basering)+");");
    1163       if(y>=0)
    1164       {
    1165             "//   a zero-divisor was found, we have SPLITTING of ideals";
    1166             "";
     1123      if(y>=1)
     1124      {
     1125         "// zero-divisor found";
    11671126      }
    11681127      ideal vid=fetch(BAS,new1);
    … …  
    11701129      attrib(vid,"isCohenMacaulay",0);
    11711130      attrib(vid,"isPrim",0);
    1172       attrib(vid,"isIsolatedSingularity",0);
     1131      attrib(vid,"isIsolatedSingularity",isIs);
    11731132      attrib(vid,"isRegInCodim2",0);
    1174       if(equi==1)
    1175       {
    1176          attrib(vid,"isEquidimensional",1);
    1177       }
    1178       else
    1179       {
    1180          attrib(vid,"isEquidimensional",0);
    1181       }
     1133      attrib(vid,"onlySingularAtZero",oSAZ);
     1134      attrib(vid,"isEquidimensional",equi);     
    11821135      attrib(vid,"isCompleteIntersection",0);
    11831136
    1184       if(y>=0)
    1185       {
    1186             "// start a normalization loop with 1st split ideal (size",size(vid),"in",nvars(basering),"vars)";
    1187       }
    1188 
    1189       keepresult1=normalizationPrimes(vid,ihp,count);
    1190 
     1137      keepresult1=normalizationPrimes(vid,ihp,0);
     1138      int delta1=keepresult1[size(keepresult1)];
    11911139      setring BAS;
    1192       ideal new2=quotient(SM[2],Ann)+SM[2];
    1193 // evtl. qAnn nehmen
     1140
    11941141      execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
    11951142                      +ordstr(basering)+");");
    … …  
    11991146      attrib(vid,"isCohenMacaulay",0);
    12001147      attrib(vid,"isPrim",0);
    1201       attrib(vid,"isIsolatedSingularity",0);
     1148      attrib(vid,"isIsolatedSingularity",isIs);
    12021149      attrib(vid,"isRegInCodim2",0);
    1203       if(equi==1)
    1204       {
    1205          attrib(vid,"isEquidimensional",1);
    1206       }
    1207       else
    1208       {
    1209          attrib(vid,"isEquidimensional",0);
    1210       }
     1150      attrib(vid,"isEquidimensional",equi);     
    12111151      attrib(vid,"isCompleteIntersection",0);
    1212 
    1213       if(y>=0)
    1214       {
    1215             "// start a normalization loop with 2nd split ideal (size",size(vid),"in",nvars(basering),"vars)";
    1216       }
    1217 
    1218       keepresult2=normalizationPrimes(vid,ihp,count);
     1152      attrib(vid,"onlySingularAtZero",oSAZ);
     1153
     1154      keepresult2=normalizationPrimes(vid,ihp,0);
     1155      int delta2=keepresult2[size(keepresult2)];
    12191156
    12201157      setring BAS;
    1221       for(lauf=1;lauf<=size(keepresult2);lauf++)
     1158     
     1159      for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
    12221160      {
    12231161         keepresult1=insert(keepresult1,keepresult2[lauf]);
    12241162      }
     1163      keepresult1[size(keepresult1)]=delta+mul+delta1+delta2;
    12251164      return(keepresult1);
    12261165   }
    … …  
    12401179   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
    12411180
    1242    list pr=normalizationPrimes(i,maxideal(1),1);
     1181   list pr=normalizationPrimes(i);
    12431182   def r1=pr[1];
    12441183   setring r1;
    … …  
    12571196   int ii,jj;
    12581197   map phi = basering,maxideal(1);
    1259 
    1260    //endid=diagon(endid);
    12611198
    12621199   list Le = elimpart(endid);
    … …  
    13341271   return(L);
    13351272}
    1336 ///////////////////////////////////////////////////////////////////////////////
    1337 static
    1338 proc diagon(ideal i)
     1273/////////////////////////////////////////////////////////////////////////////
     1274
     1275proc genus(ideal K,list #)
     1276"USAGE:   genus(I) or genus(i,1); I a 1-dimensional ideal
     1277RETURN:  an integer, the geometric genus p_g = p_a - delta of the projective
     1278         curve defined by I, where p_a is the arithmetic genus.
     1279NOTE:    delta is the sum of all local delta-invariants of the singularities,
     1280         i.e. dim(R'/R), R' the normalization of the local ring R of the
     1281         singularity.
     1282         genus(i,1) uses the normalization to compute delta. Usually this
     1283         is slow but sometimes not.
     1284EXAMPLE: example genus; shows an example
     1285"
    13391286{
    1340    matrix m;
    1341    intvec iv = option(get);
    1342    option(redSB);
    1343    ideal j=liftstd(jet(i,1),m);
    1344    option(set,iv);
    1345    return(ideal(matrix(i)*m));
     1287   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
     1288
     1289   def R=basering;
     1290   K=std(K);
     1291 
     1292   if(nvars(R)==3)
     1293   {
     1294      if((dim(K)!=2)||(!homog(K))||(size(K)>1)){ERROR("Input is not a curve");}
     1295      execute("ring newR=("+charstr(R)+"),(x,y),dp;"); 
     1296      map kappa=R,x,y,1;
     1297      ideal K=kappa(K); 
     1298   }
     1299   if((nvars(R)>3)||(size(K)>1))
     1300   {  // hier geeignet projezieren
     1301      ERROR("This case is not implemented yet");
     1302   }
     1303   if(nvars(R)==2)
     1304   {
     1305      execute("ring newR=("+charstr(R)+"),(x,y),dp;"); 
     1306      map kappa=R,x,y;
     1307      ideal K=kappa(K); 
     1308   }
     1309   
     1310   // assume now that R is a ring with two variables 
     1311   poly p=K[1];   
     1312   ideal I;
     1313   if(homog(p))
     1314   {
     1315      if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");}
     1316      return(-deg(p)+1);
     1317   }
     1318   execute("ring S=("+charstr(R)+"),(x,y,t),dp;");
     1319   execute("ring C=("+charstr(R)+"),(x,y),ds;");
     1320   ideal I;
     1321   execute("ring A=("+charstr(R)+"),(x,t),dp;");
     1322   map phi=S,1,x,t;
     1323   map psi=S,x,1,t;
     1324   poly g,h;
     1325   ideal I,I1;
     1326   execute("ring B=("+charstr(R)+"),(x,t),ds;");
     1327   
     1328   setring S;
     1329   poly F=imap(newR,p);
     1330   F=homog(F,t);
     1331   int d=deg(F);
     1332   int genus=(d-1)*(d-2)/2;
     1333   
     1334//   if(w>=1){"test for smoothness";}
     1335//   if(dim(std(jacob(F)))==0)  //the smooth case
     1336//   {
     1337//      setring R;
     1338//      return(genus);
     1339//   }
     1340 
     1341   int delta,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,
     1342       tauloc,tausing,k,rat,nbranchinf,nbranch,nodes;
     1343   list inv;
     1344
     1345   if(w>=1){"singularities at oo";}
     1346   setring A;       
     1347   g=phi(F);
     1348   h=psi(F);
     1349   I=g,jacob(g),var(2);
     1350   I=std(I);
     1351   if(deg(I[1])>0)
     1352   {
     1353      list qr=minAssGTZ(I);
     1354      for(k=1;k<=size(qr);k++)
     1355      {
     1356         inv=deltaLoc(g,qr[k]);
     1357         deltainf=deltainf+inv[1];
     1358         tauinf=tauinf+inv[2];
     1359         nbranchinf=nbranchinf+inv[3];
     1360      }
     1361   }
     1362   inv=deltaLoc(h,maxideal(1));
     1363   deltainf=deltainf+inv[1];
     1364   tauinf=tauinf+inv[2];
     1365   if(inv[2]>0){nbranchinf=nbranchinf+inv[3];}
     1366 
     1367   if(w>=1)
     1368   {
     1369      "branches at oo:";nbranchinf;
     1370      "tau at oo:";tauinf;
     1371      "delta at oo:";deltainf;
     1372      "singularities not at oo";
     1373   }
     1374                         
     1375   setring newR;         //the singularities at the affine part
     1376   map sigma=S,var(1),var(2),1;
     1377   I=sigma(F);
     1378
     1379   if((size(#)!=0)||((char(basering)<181)&&(char(basering)!=0)))     
     1380   {                              //uses the normalization to compute delta
     1381      list nor=normal(I,"wd");
     1382      delta=nor[size(nor)];
     1383      genus=genus-delta-deltainf;
     1384      setring R;
     1385      return(genus);
     1386   }
     1387
     1388   ideal I1=jacob(I);
     1389   matrix Hess[2][2]=jacob(I1);
     1390   ideal ID=I+I1+ideal(det(Hess));
     1391
     1392   if(w>=1){"the cusps and nodes";}
     1393   
     1394   ideal radID=std(radical(ID));
     1395   ideal IDsing=minor(jacob(ID),2)+radID;
     1396   iglob=vdim(std(IDsing));
     1397
     1398   if(iglob!=0)
     1399   {
     1400      ideal radIDsing=reduce(IDsing,radID);
     1401      if(size(radIDsing)==0)
     1402      {
     1403         radIDsing=radID;
     1404         attrib(radIDsing,"isSB",1);
     1405      }
     1406      else
     1407      {
     1408         radIDsing=std(radical(IDsing));
     1409      }
     1410      iglob=vdim(radIDsing);
     1411   }
     1412   cusps=vdim(radID)-iglob;
     1413       
     1414   if(w>=1){"the other singularities";}
     1415 
     1416   if(iglob==0)   //only cusps and double points
     1417   {
     1418      tau=vdim(std(I+jacob(I)));
     1419      nodes=tau-2*cusps;
     1420      delta=nodes+cusps;
     1421      nbranch=2*tau-3*cusps;   
     1422   }
     1423   else
     1424   {
     1425       setring C;                     
     1426       ideal I1=imap(newR,IDsing);
     1427       iloc=vdim(std(I1));
     1428       if(iglob==iloc)  // only cusps and nodes outside 0
     1429       { 
     1430          I=imap(newR,I);
     1431          tauloc=vdim(std(I+jacob(I)));
     1432          setring newR;
     1433          inv=deltaLoc(I[1],maxideal(1));
     1434          delta=inv[1];
     1435          tau=inv[2];
     1436          nodes=tau-tauloc-2*cusps;
     1437          nbranch=inv[3]+ 2*nodes+cusps;         
     1438          delta=delta+nodes+cusps;
     1439        }
     1440        else
     1441        {
     1442           setring newR;
     1443           list pr=minAssGTZ(IDsing);                   
     1444           for(k=1;k<=size(pr);k++)
     1445           {
     1446              inv=deltaLoc(I[1],pr[k]);
     1447              delta=delta+inv[1];
     1448              tausing=tausing+inv[2];
     1449              nbranch=nbranch+inv[3];
     1450           }
     1451           nodes=tau-tauloc-2*cusps;
     1452           tau=vdim(std(I+jacob(I)));
     1453           delta=delta+nodes+cusps; 
     1454           nbranch=nbranch+2*nodes+cusps;           
     1455        }
     1456   }
     1457
     1458   if(w>=1)
     1459   {
     1460      "branches :";nbranch;
     1461      "nodes:"; nodes;
     1462      "cusps:";cusps;
     1463      "tau :";tau;
     1464      "delta:";delta;
     1465   }
     1466   genus=genus-delta-deltainf;
     1467   setring R;
     1468   return(genus);
    13461469}
    1347 /////////////////////////////////////////////////////////////////////////////
     1470example
     1471{ "EXAMPLE:"; echo = 2;
     1472   ring r=0,(x,y),dp;
     1473   ideal i=y^9 - x^2*(x - 1)^9;
     1474   genus(i);
     1475}
     1476///////////////////////////////////////////////////////////////////////////
     1477proc deltaLoc(poly f,ideal singL)
     1478{
     1479   def R=basering;
     1480   execute("ring S=("+charstr(R)+"),(x,y),lp;");
     1481   map phi=R,x,y;
     1482   ideal singL=phi(singL);
     1483   singL=std(singL);
     1484   int d=vdim(singL);
     1485   poly f=phi(f);
     1486   int i;
     1487
     1488   if(d==1)
     1489   {
     1490      map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2];
     1491      f=alpha(f);
     1492      execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;");
     1493      poly f=imap(S,f);
     1494   }
     1495   else
     1496   {
     1497      poly p;
     1498      poly c;
     1499      map psi;
     1500      while((deg(singL[1])>1)&&(deg(singL[2])>1))
     1501      {
     1502         psi=S,x,y+random(-100,100)*x;
     1503         singL=psi(singL);
     1504         singL=std(singL);
     1505      }
     1506      if(deg(singL[2])==1){p=singL[1];c=singL[2][2];}
     1507      if(deg(singL[1])==1)
     1508      {
     1509         psi=S,y,x;
     1510         f=psi(f);
     1511         singL=psi(singL);
     1512         p=singL[2];
     1513         c=singL[1][2];
     1514      }
     1515      execute("ring B=("+charstr(S)+"),a,dp;");
     1516      map beta=S,a,a;
     1517      poly p=beta(p);
     1518      execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;");
     1519      number p=number(imap(B,p));
     1520      minpoly=p;
     1521      number c=number(imap(S,c));
     1522      map alpha=S,x-c,y+a;
     1523     
     1524      poly f=alpha(f);
     1525      f=cleardenom(f);
     1526   }
     1527   ideal fstd=std(ideal(f)+jacob(f));
     1528   poly hc=highcorner(fstd);
     1529   int tau=vdim(fstd);     
     1530   int o=ord(f);
     1531   int delta,nb;
     1532
     1533   if(tau==0)                 //smooth case
     1534   {
     1535      setring R;
     1536      return(list(0,0,1));
     1537   }
     1538   if((char(basering)>=181)||(char(basering)==0))
     1539   {
     1540      if(o==2)                //A_k-singularity
     1541      {
     1542         setring R;
     1543         delta=(tau+1)/2;
     1544         return(list(d*delta,d*tau,d*(2*delta-tau+1)));   
     1545      }
     1546      if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2)))
     1547      {//D_k- singularity
     1548         setring R;
     1549         delta=(tau+2)/2;
     1550         return(list(d*delta,d*tau,d*(2*delta-tau+1)));
     1551      }
     1552
     1553      int mu=vdim(std(jacob(f)));
     1554      poly g=f+var(1)^mu+var(2)^mu;  //to obtain a convinient Newton-polygon
     1555 
     1556      list NP=newton(g);
     1557      int s=size(NP);
     1558      int nN=-NP[1][2]-NP[s][1]+1;  // computation of the Newton-number
     1559      intmat m[2][2];
     1560      for(i=1;i<=s-1;i++)
     1561      {
     1562         m=NP[i+1],NP[i];
     1563         nN=nN+det(m);
     1564      }
     1565
     1566      if(mu==nN)                   // the Newton-polygon is non-degenerate
     1567      {                            // compute nb, the number of branches
     1568        for(i=1;i<=s-1;i++)
     1569        {
     1570          nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]);
     1571        }
     1572        return(list(d*(mu+nb-1)/2,d*tau,d*nb));
     1573      }
     1574
     1575   //da reddevelop nur benutzt wird, um die Anzahl der Zweige zu bestimmen
     1576   //kann man das sicher schneller machen:
     1577   //die Aufblasung durchfuehren und stets testen, ob das Newton-polyeder
     1578   //nicht ausgeartet ist.
     1579
     1580      if(s>2)    //  splitting of f
     1581      {
     1582         intvec v=NP[1][2]-NP[2][2],NP[2][1];
     1583         int de=w_deg(g,v);
     1584         int st=w_deg(hc,v)+v[1]+v[2];
     1585         poly f1=var(2)^NP[2][2];
     1586         poly f2=jet(g,de,v)/var(2)^NP[2][2];
     1587         poly h=g-f1*f2;
     1588         de=w_deg(h,v);
     1589         poly k;
     1590         ideal wi=var(2)^NP[2][2],f2;
     1591         matrix li;
     1592         while(de<st)
     1593         {
     1594           k=jet(h,de,v);
     1595           li=lift(wi,k);
     1596           f1=f1+li[2,1];
     1597           f2=f2+li[1,1];
     1598           h=g-f1*f2;
     1599           de=w_deg(h,v);
     1600         }
     1601         nb=deltaLoc(f1,maxideal(1))[3]+deltaLoc(f2,maxideal(1))[3];
     1602         return(list(d*(mu+nb-1)/2,d*tau,d*nb));
     1603      }
     1604
     1605      f=jet(f,deg(hc)+2);
     1606      list hne=reddevelop(f);
     1607      nb=size(hne);
     1608      setring R;
     1609      kill HNEring;       
     1610      return(list(d*(mu+nb-1)/2,d*tau,d*nb));
     1611   }
     1612   else             //the case of small characteristic
     1613   {
     1614      f=jet(f,deg(hc)+2);
     1615      list hne=reddevelop(f);
     1616      nb=size(hne);
     1617      if(nb==1)
     1618      {
     1619         delta=invariants(hne[1])[5]/2;
     1620         setring R;
     1621         kill HNEring;
     1622         return(list(d*delta,d*tau,d));
     1623      }
     1624      setring R;
     1625      kill HNEring;
     1626      //delta direkt aus reddevelop zurueckgeben
     1627      ERROR("the case of small characteristic is not fully implemented yet");
     1628   }
     1629}
     1630
     1631proc w_deg(poly p, intvec v)
     1632{
     1633   if(p==0){return(-1);}
     1634   int d=0;
     1635   while(jet(p,d,v)==0){d++;}
     1636   d=(transpose(leadexp(jet(p,d,v)))*v)[1];
     1637   return(d); 
     1638}
     1639
     1640proc newton (poly f)
     1641{
     1642   def R1=basering;
     1643   execute("ring R2=("+charstr(R1)+"),("+varstr(R1)+"),ls;");
     1644   poly f=imap(R1,f);
     1645   intvec A=(0,ord(subst(f,var(1),0)));
     1646   intvec B=(ord(subst(f,var(2),0)),0);
     1647   intvec C,H; list L;
     1648   int abbruch,i;
     1649   poly hilf;
     1650   L[1]=A;
     1651   f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
     1652   map xytausch=R2,var(2),var(1);
     1653   for (i=2; f!=0; i++)
     1654   {
     1655      abbruch=0;
     1656      while (abbruch==0)
     1657      {
     1658         C=leadexp(f);         
     1659         if(jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
     1660         {
     1661            abbruch=1;
     1662         }       
     1663         else
     1664         {
     1665            f=jet(f,-C[1]-1,intvec(-1,0));
     1666         }
     1667     }
     1668     hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));
     1669     H=leadexp(xytausch(hilf));
     1670     A=H[2],H[1];
     1671     L[i]=A;
     1672     f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]-A[1]));
     1673   }
     1674   L[i]=B;
     1675   setring R1;
     1676   return(L);
     1677}
     1678
     1679///////////////////////////////////////////////////////////////////////////
     1680
    13481681/*
    1349 Aenderungen:
    1350 1. normal kommentiert, bei Berechnung de singulaeren Ortes ein simplify(J,16)
    1351 eingefuehrt, um bei riesigen Minorenzahlen, das Ideal zu verkleinern (bis
    1352 Faktor 10 Beschleunigung).
    1353 Protokoll mit printlevel so gesteuert, dass es bei rekursivem Aufruf korrekt
    1354 arbeitet.
    1355 list nor = normal(i);     //mit equidim Zerlegung
    1356 list nor = normal(i,1);   //mit prim Zerlegung
    1357 Zeiten au sony_pumuckel (P2, 500)
    13581682                           Examples:
    13591683LIB"normal.lib";
    1360 //1. Huneke, 1 Komponente
    1361 //prim: 2 sec equidim:1sec
     1684//Huneke
    13621685ring qr=31991,(a,b,c,d,e),dp;
    13631686ideal i=
    … …  
    13721695
    13731696
    1374 //2. Vasconcelos (dauert laenger: 83 sec auf sony_pumuckel)
     1697//Vasconcelos (dauert laenger: 70 sec)
    13751698ring r=32003,(x,y,z,w,t),dp;
    13761699ideal i=
    … …  
    13791702xw3+z3t+ywt2,
    13801703y2w4-xy2z2t-w3t3;
    1381 
    1382 //2a. Vasconcelos verkleinert
    1383 //prim:2Komp, 2 Ringe, 16 sec (manchmal lange, haengt am Faktorisierer)
    1384 //equidim: 1 Komp, 7 loops, 2 ringe, 12sec
    1385 ring r=32003,(x,y,z,w,t),dp;
    1386 ideal i=
    1387 x+zw,
    1388 y3+xwt,
    1389 xw3+z3t+ywt2;
    1390 
    1391 //3. GM1
    1392 // irreducible, 13 normalization loops, 1 Ring
    1393 //2sec mit prim, 1 sec mit equidim
    1394 ring r=32003,(x,y,z,u),dp;
    1395 ideal i = x2+y3,z2+u3,y2+z3;
    1396 
    1397 //3a. GM1
    1398 ring r=32003,(x,y),dp;
    1399 ideal i = intersect(y3+x2,y2+x3);  // beide 0 sec
    1400 ideal i = intersect(y3+x2,y2+x3,y2+x2+x3);
    1401 //prim 0sec,
    1402 //equidim sehr lange (Relationen in HomJJ zu gross)
    1403 
    1404 //4. GM2
    1405 //i nicht reduziert, equidim bricht ab (0 sec)
    1406 //prim: 11 irred comp, 11 loops, (1sec)
    1407 ring r=32003,(x,y,z,u),dp;
    1408 ideal i = x2+y3,z2+u3,y2+z3,x2+z3;
    1409 
    1410 //5. GM3,  radical von GM2
    1411 //prim: 11 irred comp, 11 loops, 11 Ringe, (1sec)
    1412 //equidim: 1 equidim comp, 1 loop, 1 Ring, (0sec)
    1413 ring r=32003,(x,y,z,u),dp;
    1414 ideal i =yu+u,yz+z,y2+y,xy+x,x2+y,u3+z2,z3-y;
    1415 
    1416 //GM4
    1417 //equidim: 2 equidim comp, 3 Ringe, (0sec);
    1418 //prim: 3 Komp, (0sec)
    1419 ring r=32003,(x,y,z,u),dp;
    1420 ideal i1 = x2+y3,u,z;
    1421 ideal i2 = u2+z3,x,y;
    1422 ideal i3 = x2+y2+z2+u4;
    1423 ideal i = intersect(i1,i2,i3);
    1424 
    1425 //GM4a  Hier dauert prim laenger! (## facstd)
    1426 //equidim: 3 equidim Komp, 4 Ringe(0sec)
    1427 //prim 13 Komp (63 sec), wegen facstd
    1428 //##ev minAssGTZ(i,1) verwenden (ohne facstd, ist noch fehlerhaft)
    1429 ring r=32003,(x,y,z,u),dp;
    1430 ideal i1 = x2+y3,u,z;
    1431 ideal i2 = u2+z3,x,y;
    1432 ideal i3 = x2+y2+z2+u4;
    1433 ideal i4 =yu+u,yz+z,y2+y,xy+x,x2+y,u3+z2,z3-y;
    1434 ideal i = intersect(i1,i2,i3,i4);
    1435 
    1436 //GM5
    1437 //equidim: 4 Komp,0sec, prim 13 Komp, 1sec
    1438 ring r=32003,(x,y,z,u,v),dp;
    1439 ideal i1 = x2+y3,v;                            //2dim
    1440 ideal i2 = u2+z3,x,y;                          //3dim
    1441 ideal i3 = x2+y2+z2+u4;                        //4dim
    1442 ideal i4 =yu+u,yz+z,y2+y,xy+x,x2+y,u3+z2,z3-y; //1dim
    1443 ideal i = intersect(i1,i2,i3,i4);
    1444 
    1445 //cyclic 5
    1446 //equidim: 1 Komp 0sec, prim: 20 Komp, 3 sec
    1447 ring r=32003,(x,y,z,u,v),dp;
    1448 ideal i =
    1449 x+y+z+u+v,
    1450 xy+yz+zu+xv+uv,
    1451 xyz+yzu+xyv+xuv+zuv,
    1452 xyzu+xyzv+xyuv+xzuv+yzuv,
    1453 xyzuv-1;
    1454 
    1455 // cyclic 5 hat Normalisierung (1 embim weniger)
    1456 ///equidim: 1 Komp 0sec, prim: 20 Komp, 2 sec
    1457 ring r=32003,(x,y,z,u),dp;
    1458 ideal i =
    1459 x2+xz-yz+2xu+yu+u2,
    1460 xy2-xyz+y2z-y2u+xzu+yzu+z2u-xu2-2yu2+zu2-u3,
    1461 xyz2+xyzu+y2zu-xz2u+yz2u-z3u-xyu2-xzu2-2z2u2+xu3+yu3-zu3+u4,
    1462 2xyzu2+y2zu2+2yz2u2-xyu3-2xzu3-yzu3-z2u3+xu4+yu4-2zu4+u5-1;
    1463 
    1464 //cyclic(6)
    1465 //equidim: 1Komp in 5 vars 1sec
    1466 //prim: 90 (!) Ringe, 12 sec
    14671704
    14681705//Theo1
    … …  
    15111748ad;
    15121749
    1513 //Sturmfels, wo vorher Prim schneller (2 sec,sonst 860 sec)
    15141750//ist CM
    1515 //prim: 15 loops, 15 Komp, 1 sec,
    1516 //equidim:15 Ringe, 93 sec mit simplify(J,16),
    1517 //ohne simlify(J,16) 860sec?,
    1518 //andere simplify sind z.T. viel langsamer
     1751//Sturmfels
    15191752ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
    15201753ideal i=
    … …  
    15481781
    15491782//Horrocks:
    1550 //CHAR 32003:mit prim 1 sec, equidim: 115 sec, beide 6 Ringe
    1551 //           Singulaere Ort hat zu Beginn > 106 000 Erzeuger!!
    1552 //char 31991: prim 1sec, 8 Ringe, equidim: 25 Ringe(!), 162 sec,
    1553 //            nicht reduziert!
    1554 //            Singulaere Ort hat zu Beginn > 28 000 Erzeuger!!
    1555 //            i=radical(i) -> 8 Ringe, 1sec (radical <1 sec)
    1556 //Horrocks:
    1557 ring r=31991,(a,b,c,d,e,f),dp;
     1783ring r=32003,(a,b,c,d,e,f),dp;
    15581784ideal i=
    15591785adef-16000be2f+16001cef2,
    … …  
    15921818
    15931819
    1594 ring r=32003,(b,s,t,u,v,w,x,y,z),dp;   //3sec
     1820ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
    15951821ideal k=
    15961822wy-vz,
    … …  
    16021828ideal i=mstd(intersect(j,k))[2];
    16031829
    1604 //22,
    1605 // neu, prim: 3 sec, equidim 1 sec, je 4 Ringe
     1830
    16061831ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
    16071832ideal i=
    … …  
    16131838
    16141839
    1615 //riemenschneider, 5 Komponenten
    1616 //33(alte Zeiten), normal+primary 3, primary 9, radical 1, minAssGTZ; 2
    1617 //neu: prim 0sec, equi 1 sec, je 5 Ringe
     1840//riemenschneider 
    16181841ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
    16191842ideal i=
    … …  
    16311854
    16321855ring r=0,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
    1633 ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;   //0sec
     1856ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;
     1857
     1858//Yoshihiko Sakai
     1859ring r=0,(x,y),dp;   //genus 0 4 nodes and 6 cusps   
     1860ideal i=(x2+y^2-1)^3 +27x2y2;
     1861
     1862ring r=0,(x,y),dp;   //genus 0
     1863ideal i=(x-y^2)^2 - y*x^3;
     1864
     1865ring r=0,(x,y),dp;  //genus 4
     1866ideal i=y3-x6+1;
     1867
     1868int m=9;           // q=9: genus 0
     1869int p=2;
     1870int q=9;//2,...,9
     1871ring r=0,(x,y),dp; 
     1872ideal i=y^m - x^p*(x - 1)^q;
     1873
     1874ring r=0,(x,y),dp;  //genus 19
     1875ideal i=55*x^8+66*y^2*x^9+837*x^2*y^6-75*y^4*x^2-70*y^6-97*y^7*x^2;
     1876
     1877
     1878ring r=0,(x,y),dp;  //genus 34
     1879ideal i=y10+(-2494x2+474)*y8+(84366+2042158x4-660492)*y6
     1880        +(128361096x4-47970216x2+6697080-761328152x6)*y4
     1881        +(-12024807786x4-506101284x2+15052058268x6+202172841-3212x8)*y2
     1882        +34263110700x4-228715574724x6+5431439286x2+201803238
     1883        -9127158539954x10-3212722859346x8;
     1884
     1885//Rob Koelman
     1886ring r=0,(x,y,z),dp;//genus 10  with 26 cusps
     1887ideal i=
     1888761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+
     1889 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2-
     1890 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4-
     1891 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+
     1892 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-
     1893 z^10-474*z^8*y^2-84366*z^6*y^4;
     1894
     1895ring r=0,(x,y),dp;//genus 10  with 26 cusps
     1896ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
     1897-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
     1898+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
     1899-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;
     1900
     1901
     1902ring r=0,(x,y),dp;   // genuss 1  with 5 cusps
     1903ideal i=57y5+516x4y-320x4+66y4-340x2y3+73y3+128x2-84x2y2-96x2y;;
     1904
     1905//Mark van Hoeij
     1906ring r=0,(x,y),dp;  //genus 19
     1907ideal i=y20+y13x+x4y5+x3*(x+1)^2;
     1908
     1909ring r=0,(x,y),dp;  //genus 35
     1910ideal i=y30+y13x+x4y5+x3*(x+1)^2;
     1911
     1912ring r=0,(x,y),dp;   //genus 55
     1913ideal i=y40+y13x+x4y5+x3*(x+1)^2;
     1914
     1915ring r=0,(x,y),dp;  //genus 55
     1916ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
     1917
    16341918*/
    16351919
Note: See TracChangeset for help on using the changeset viewer.