source: git/Singular/LIB/normal.lib @ 7497ef3

spielwiese
Last change on this file since 7497ef3 was 7497ef3, checked in by Gerhard Pfister <pfister@…>, 22 years ago
Bug bei genus gefixed git-svn-id: file:///usr/local/Singular/svn/trunk@5652 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 48.1 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: normal.lib,v 1.37 2001-10-23 10:26:24 pfister Exp $";
3category="Commutative Algebra";
4info="
5LIBRARY:  normal.lib     Normalization of Affine Rings
6AUTHORS:  G.-M. Greuel,  greuel@mathematik.uni-kl.de,
7@*        G. Pfister,    pfister@mathematik.uni-kl.de
8
9PROCEDURES:
10 normal(I[,wd]);        computes the normalization of basering/I
11                        resp. computes the normalization of basering/I and
12                        the delta-invariante
13 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
16";
17
18LIB "general.lib";
19LIB "sing.lib";
20LIB "primdec.lib";
21LIB "elim.lib";
22LIB "presolve.lib";
23LIB "inout.lib";
24LIB "ring.lib";
25LIB "hnoether.lib";
26///////////////////////////////////////////////////////////////////////////////
27
28proc HomJJ (list Li)
29"USAGE:   HomJJ (Li);  Li = list: ideal SBid, ideal id, ideal J, poly p
30ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
31@*       SBid = standard basis of id,
32@*       J    = ideal of P containing the polynomial p,
33@*       p    = nonzero divisor of R
34COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
35         affine ring, together with the canonical map R --> Hom_R(J,J),
36         where R is the quotient ring of P modulo the standard basis SBid.
37RETURN:  a list l of two objects
38@format
39         l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
40               such that l[1]/endid = Hom_R(J,J) and
41               endphi describes the canonical map R -> Hom_R(J,J)
42         l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
43         l[3] : an integer, the contribution to delta
44@end format
45NOTE:    printlevel >=1: display comments (default: printlevel=0)
46EXAMPLE: example HomJJ;  shows an example
47"
48{
49//---------- initialisation ---------------------------------------------------
50
51   int isIso,isPr,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
52   intvec rw,rw1;
53   list L;
54   y = printlevel-voice+2;  // y=printlevel (default: y=0)
55   def P = basering;
56   ideal SBid, id, J = Li[1], Li[2], Li[3];
57   poly p = Li[4];
58   attrib(SBid,"isSB",1);
59   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not
60
61//---- set attributes for special cases where algorithm can be simplified -----
62   if( homo==1 )
63   {
64      rw = ringweights(P);
65   }
66   if( typeof(attrib(id,"isPrim"))=="int" )
67   {
68      if(attrib(id,"isPrim")==1)  { isPr=1; }
69   }
70   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
71   {
72      if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; }
73   }
74   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
75   {
76      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
77   }
78   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
79   {
80      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
81   }
82   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
83   {
84      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
85   }
86   if( typeof(attrib(id,"isEquidimensional"))=="int" )
87   {
88      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
89   }
90//-------------------------- go to quotient ring ------------------------------
91   qring R  = SBid;
92   ideal id = fetch(P,id);
93   ideal J  = fetch(P,J);
94   poly p   = fetch(P,p);
95   ideal f,rf,f2;
96   module syzf;
97//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
98   if ( y>=1 )
99   {
100     "// compute p*Hom(J,J) = p*J:J, p a non-zerodivisor";
101     "//   p is equal to:"; "";
102     p;
103     "";
104   }
105   f  = quotient(p*J,J);
106
107   if ( y>=1 )
108   { "// the module p*Hom(J,J) = p*J:J, p a non-zerodivisor";
109      "// p"; p;
110      "// f=p*J:J";f;
111   }
112   f2 = std(p);
113 
114//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
115
116   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R
117
118   if ( size(rf) == 0 )
119   {
120      if ( homog(f) && find(ordstr(basering),"s")==0 )
121      {
122         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
123      }
124      else
125      {
126         ring newR1 = char(P),(X(1..nvars(P))),dp;
127      }
128      ideal endphi = maxideal(1);
129      ideal endid = fetch(P,id);
130      L=substpart(endid,endphi,homo,rw);
131      def lastRing=L[1];
132      setring lastRing;
133
134      attrib(endid,"onlySingularAtZero",oSAZ);
135      attrib(endid,"isCohenMacaulay",isCo);
136      attrib(endid,"isPrim",isPr);
137      attrib(endid,"isIsolatedSingularity",isIso);
138      attrib(endid,"isRegInCodim2",isRe);
139      attrib(endid,"isEqudimensional",isEq);
140      attrib(endid,"isCompleteIntersection",0);
141      attrib(endid,"isRad",0);
142      L=lastRing;
143      L = insert(L,1,1);
144      dbprint(y,"// case R = Hom(J,J)");
145      if(y>=1)
146      {
147         "//   R=Hom(J,J)";
148         "   ";
149         lastRing;
150         "   ";
151         "//   the new ideal";
152         endid;
153         "   ";
154         "//   the old ring";
155         "   ";
156         P;
157         "   ";
158         "//   the old ideal";
159         "   ";
160         setring P;
161         id;
162         "   ";
163         setring lastRing;
164         "//   the map";
165         "   ";
166         endphi;
167         "   ";
168         pause();
169         newline;
170      }
171      setring P;
172      L[3]=0;
173      return(L);
174   }
175   if(y>=1)
176   {
177      "// R is not equal to Hom(J,J), we have to try again";
178      pause();
179      newline;
180   }
181//---------- Hom(J,J) != R: create new ring and map from old ring -------------
182// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
183
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
189   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
190   q = size(f);
191
192   syzf = syz(f);
193
194   if ( homo==1 )
195   {
196      rw1 = rw,0;
197      for ( ii=2; ii<=q; ii++ )
198      {
199         rw  = rw, deg(f[ii])-deg(f[1]);
200         rw1 = rw1, deg(f[ii])-deg(f[1]);
201      }
202      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
203   }
204   else
205   {
206      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
207   }
208
209   map psi1 = P,maxideal(1);
210   ideal SBid = psi1(SBid);
211   attrib(SBid,"isSB",1);
212
213   qring newR = std(SBid);
214
215   map psi = R,ideal(X(1..nvars(R)));
216   ideal id = psi(id);
217   ideal f = psi(f);
218   module syzf = psi(syzf);
219   ideal pf,Lin,Quad,Q;
220   matrix T,A;
221   list L1;
222
223//---------- computation of Hom(J,J) as affine ring ---------------------------
224// determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
225// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course
226// the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
227// It is a fact, that the kernel is generated by the linear and the quadratic
228// relations
229
230   pf = f[1]*f;
231   T = matrix(ideal(T(1..q)),1,q);
232   Lin = ideal(T*syzf);
233   if(y>=1)
234   {
235      "// the ring structure of Hom(J,J) as R-algebra";
236      " ";
237      "//   the linear relations";
238      " ";
239      Lin;
240      "   ";
241   }
242
243   for (ii=2; ii<=q; ii++ )
244   {
245      for ( jj=2; jj<=ii; jj++ )
246      {
247         A = lift(pf,f[ii]*f[jj]);
248         Quad = Quad, ideal(T(jj)*T(ii) - T*A);          // quadratic relations
249      }
250   }
251
252   if(y>=1)
253   {
254      "//   the quadratic relations";
255      "   ";
256      interred(Quad);
257      pause();
258      newline;
259   }
260   Q = Lin,Quad;
261   Q = subst(Q,T(1),1);
262   Q=mstd(Q)[2];
263//---------- reduce number of variables by substitution, if possible ----------
264   if (homo==1)
265   {
266      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
267   }
268   else
269   {
270      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
271   }
272
273   ideal endid  = imap(newR,id),imap(newR,Q);
274   ideal endphi = ideal(X(1..nvars(R)));
275
276   L=substpart(endid,endphi,homo,rw);
277
278   def lastRing=L[1];
279   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;
291   attrib(endid,"isCohenMacaulay",isCo);
292   attrib(endid,"isPrim",isPr);
293   attrib(endid,"isIsolatedSingularity",isIso);
294   attrib(endid,"isRegInCodim2",isRe);
295   attrib(endid,"isEquidimensional",isEq);
296   attrib(endid,"isCompleteIntersection",0);
297   attrib(endid,"isRad",0);
298   if(y>=1)
299   {
300      "//   the new ring after reduction of the number of variables";
301      "   ";
302      lastRing;
303      "   ";
304      "//   the new ideal";
305      "   ";
306      endid;
307      "   ";
308      "//   the old ring";
309      "   ";
310      P;
311      "   ";
312      "//   the old ideal";
313      "   ";
314      setring P;
315      id;
316      "   ";
317      setring lastRing;
318      "//   the map";
319      "   ";
320      endphi;
321      "   ";
322      pause();
323      newline;
324   }
325   L = lastRing;
326   L = insert(L,0,1);
327   L[3]=delta;
328   return(L);
329}
330example
331{"EXAMPLE:";  echo = 2;
332  ring r   = 0,(x,y),wp(2,3);
333  ideal id = y^2-x^3;
334  ideal J  = x,y;
335  poly p   = x;
336  list Li = std(id),id,J,p;
337  list L   = HomJJ(Li);
338  def end = L[1];    // defines ring L[1], containing ideals endid, endphi
339  setring end;       // makes end the basering
340  end;
341  endid;             // end/endid is isomorphic to End(r/id) as ring
342  map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
343  psi;
344}
345
346///////////////////////////////////////////////////////////////////////////////
347
348proc normal(ideal id, list #)
349"USAGE:   normal(i [,choose]);  i a radical ideal, choose empty, 1 or "wd"
350         if choose=1 the normalization of the associated primes is computed
351         (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.
355ASSUME:  The ideal must be radical, for non radical ideals the output may
356         be wrong (i=radical(i); makes i radical)
357RETURN:   a list of rings, say nor and in case of choose="wd" an integer
358@format
359         at the end of the list.
360         each ring nor[i] contains two ideals
361         with given names norid and normap such that
362           - the direct sum of the rings nor[i]/norid is
363             the normalization of basering/id;
364           - normap gives the normalization map from basering/id
365             to nor[i]/norid (for each i)
366@end format
367NOTE:    to use the i-th ring type: def R=nor[i]; setring R;.
368@*       Increasing printlevel displays more comments (default: printlevel=0)
369@*       Not implemented for local or mixed orderings.
370@*       If the input ideal i is weighted homogeneous a weighted ordering may
371         be used (qhweight(i); computes weights).
372EXAMPLE: example normal; shows an example
373"
374{
375   int i,j,y,withdelta;
376   string sr;
377   list result,prim,keepresult;
378   y = printlevel-voice+2;
379   if(size(#)>0)
380   {
381     if(typeof(#[1])=="string")
382     {
383       kill #;
384       list #;
385       withdelta=1;
386     }
387   }
388   attrib(id,"isRadical",1);
389   if ( ord_test(basering) != 1)
390   {
391     "";
392     "// Not implemented for this ordering,";
393     "// please change to global ordering!";
394     return(result);
395   }
396   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
397   {
398      if(attrib(id,"isCompleteIntersection")==1)
399      {
400         attrib(id,"isCohenMacaulay",1);
401         attrib(id,"isEquidimensional",1);
402      }
403   }
404   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
405   {
406      if(attrib(id,"isCohenMacaulay")==1)
407      {
408         attrib(id,"isEquidimensional",1);
409      }
410   }
411   if( typeof(attrib(id,"isPrim"))=="int" )
412   {
413      if(attrib(id,"isPrim")==1)
414      {
415         attrib(id,"isEquidimensional",1);
416      }
417   }
418   if(size(#)==0)
419   {
420      if( typeof(attrib(id,"isEquidimensional"))=="int" )
421      {
422         if(attrib(id,"isEquidimensional")==1)
423         {
424            prim[1]=id;
425         }
426         else
427         {
428            prim=equidim(id);
429         }
430      }
431      else
432      {
433         prim=equidim(id);
434      }
435      if(y>=1)
436      {
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);
454      }
455   }
456   else
457   {
458      if( typeof(attrib(id,"isPrim"))=="int" )
459      {
460         if(attrib(id,"isPrim")==1)
461         {
462            prim[1]=id;
463         }
464         else
465         {
466            prim=minAssPrimes(id);
467         }
468      }
469      else
470      {
471         prim=minAssPrimes(id);
472      }
473      if(y>=1)
474      {
475         "// we have ",size(prim),"irreducible components";
476      }
477   }
478   for(i=1; i<=size(prim); i++)
479   {
480      if(y>=1)
481      {
482         "// we are in loop ",i;
483      }
484      attrib(prim[i],"isCohenMacaulay",0);
485      if(size(#)!=0)
486      {
487         attrib(prim[i],"isPrim",1);
488      }
489      else
490      {
491         attrib(prim[i],"isPrim",0);
492      }
493      attrib(prim[i],"isRegInCodim2",0);
494      attrib(prim[i],"isIsolatedSingularity",0);
495      attrib(prim[i],"isEquidimensional",1);
496      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      }
503
504      if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
505      {
506            if(attrib(id,"isIsolatedSingularity")==1)
507             {attrib(prim[i],"isIsolatedSingularity",1); }
508      }
509
510      if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
511      {
512            if((attrib(id,"isIsolatedSingularity")==1)&&(size(#)==0))
513             {attrib(prim[i],"isIsolatedSingularity",1); }
514      }
515      keepresult=normalizationPrimes(prim[i],maxideal(1),0);
516
517      for(j=1;j<=size(keepresult)-1;j++)
518      {
519         result=insert(result,keepresult[j]);
520      }
521      sr = string(size(result));     
522   }
523   if(withdelta)
524   {
525      sr = string(size(keepresult)-1);
526      result=keepresult;
527   }
528      dbprint(y+1,"
529// 'normal' created a list of "+sr+" ring(s).
530// nor["+sr+"+1] is the delta-invariant in case of choose=wd.
531// To see the rings, type (if the name of your list is nor):
532     show( nor);
533// To access the 1-st ring and map (similar for the others), type:
534     def R = nor[1]; setring R;  norid; normap;
535// R/norid is the 1-st ring of the normalization and
536// normap the map from the original basering to R/norid");
537
538      //kill endphi,endid;
539      return(result);
540}
541example
542{ "EXAMPLE:"; echo = 2;
543   ring r=32003,(x,y,z),wp(2,1,2);
544   ideal i=z3-xy4;
545   list nor=normal(i);
546   show(nor);
547   def r1=nor[1];
548   setring r1;
549   norid;
550   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)];
557}
558
559///////////////////////////////////////////////////////////////////////////////
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)
572EXAMPLE: example normalizationPrimes; shows an example
573"
574{
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";  "";
582     i;  "";
583     basering;  "";
584     pause();
585     newline;
586   }
587
588   def BAS=basering;
589   list result,keepresult1,keepresult2;
590   ideal J,SB,MB;
591   int depth,lauf,prdim;
592   int ti=timer;
593
594   if(size(i)==0)
595   {
596      if(y>=1)
597      {
598          "// the ideal was the zero-ideal";
599      }
600         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
601                      +ordstr(basering)+");");
602         ideal norid=ideal(0);
603         ideal normap=fetch(BAS,ihp);
604         export norid;
605         export normap;
606         result=newR7;
607         result[size(result)+1]=delta;
608         setring BAS;
609         return(result);
610   }
611
612   if(y>=1)
613   {
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
618  // Case: Get an ideal containing a unit
619   if( dimSM == -1)
620   {  "";
621      "      // A unit ideal was found.";
622      "      // Stop with partial result computed so far";"";
623
624         MB=SM[2];
625         intvec rw;
626         list LL=substpart(MB,ihp,0,rw);
627         def newR6=LL[1];
628         setring newR6;
629         ideal norid=endid;
630         ideal normap=endphi;
631         kill endid,endphi;
632         export norid;
633         export normap;
634         result=newR6;
635         result[size(result)+1]=delta;
636         setring BAS;
637         return(result);
638   }
639
640   if(y>=1)
641   {
642      "//   the dimension is:"; "";
643      dimSM;"";
644   }
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
665   if(attrib(i,"isPrim")==1)
666   {
667      attrib(SM[2],"isPrim",1);
668   }
669   else
670   {
671      attrib(SM[2],"isPrim",0);
672   }
673   if(attrib(i,"isIsolatedSingularity")==1)
674   {
675      attrib(SM[2],"isIsolatedSingularity",1);
676   }
677   else
678   {
679      attrib(SM[2],"isIsolatedSingularity",0);
680   }
681   if(attrib(i,"isCohenMacaulay")==1)
682   {
683      attrib(SM[2],"isCohenMacaulay",1);
684   }
685   else
686   {
687      attrib(SM[2],"isCohenMacaulay",0);
688   }
689   if(attrib(i,"isRegInCodim2")==1)
690   {
691      attrib(SM[2],"isRegInCodim2",1);
692   }
693   else
694   {
695      attrib(SM[2],"isRegInCodim2",0);
696   }
697   if(attrib(i,"isEquidimensional")==1)
698   {
699      attrib(SM[2],"isEquidimensional",1);
700   }
701   else
702   {
703      attrib(SM[2],"isEquidimensional",0);
704   }
705   if(attrib(i,"isCompleteIntersection")==1)
706   {
707      attrib(SM[2],"isCompleteIntersection",1);
708   }
709   else
710   {
711      attrib(SM[2],"isCompleteIntersection",0);
712   }
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
727   if(size(#)>0)
728   {
729      if(dim(JM[1])==-1)
730      {
731         if(y>=1)
732         {
733            "// the ideal was smooth";
734         }
735         MB=SM[2];
736         intvec rw;;
737         list LL=substpart(MB,ihp,0,rw);
738         def newR6=LL[1];
739         setring newR6;
740         ideal norid=endid;
741         ideal normap=endphi;
742         kill endid,endphi;
743         export norid;
744         export normap;
745         result=newR6;
746         result[size(result)+1]=delta;
747         setring BAS;
748         return(result);
749     }
750   }
751
752   //the zero-dimensional case
753   if((dim(SM[1])==0)&&(homog(SM[2])==1))
754   {
755      if(y>=1)
756      {
757         "// the ideal was zero-dimensional and homogeneous";
758      }
759      MB=maxideal(1);
760      intvec rw;
761      list LL=substpart(MB,ihp,0,rw);
762      def newR5=LL[1];
763      setring newR5;
764      ideal norid=endid;
765      ideal normap=endphi;
766      kill endid,endphi;
767      export norid;
768      export normap;
769      result=newR5;
770      result[size(result)+1]=delta;
771      setring BAS;
772      return(result);
773   }
774
775   //the one-dimensional case
776   //in this case it is a line because
777   //it is irreducible and homogeneous
778   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
779        &&(homog(SM[2])==1))
780   {
781      if(y>=1)
782      {
783         "// the ideal defines a line";
784      }
785      MB=SM[2];
786      intvec rw;
787      list LL=substpart(MB,ihp,0,rw);
788      def newR4=LL[1];
789      setring newR4;
790      ideal norid=endid;
791      ideal normap=endphi;
792      kill endid,endphi;
793      export norid;
794      export normap;
795      result=newR4;
796      result[size(result)+1]=delta;
797      setring BAS;
798      return(result);
799   }
800
801   //the higher dimensional case
802   //we test first of all CohenMacaulay and
803   //complete intersection
804   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
805   {
806      //test for complete intersection
807      attrib(SM[2],"isCohenMacaulay",1);
808      attrib(SM[2],"isCompleteIntersection",1);
809      attrib(SM[2],"isEquidimensional",1);
810      if(y>=1)
811      {
812         "// the ideal is a complete intersection";
813      }
814   }
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   {
832      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
833      if(y >=1 )
834      {
835         "// SB of singular locus will be computed";
836      }
837      ideal sin=J,SM[1];
838      list JM=mstd(sin);
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
843      if(y>=1)
844      {
845         "//   the dimension of the singular locus is:";"";
846         dim(JM[1]); "";
847      }
848
849      if(dim(JM[1])==-1)
850      {
851         if(y>=1)
852         {
853            "// the ideal is smooth";
854         }
855         MB=SM[2];
856         intvec rw;
857         list LL=substpart(MB,ihp,0,rw);
858         def newR3=LL[1];
859         setring newR3;
860         ideal norid=endid;
861         ideal normap=endphi;
862         kill endid,endphi;
863         export norid;
864         export normap;
865         result=newR3;
866         result[size(result)+1]=delta;
867         setring BAS;
868         return(result);
869      }
870      if(dim(JM[1])==0)
871      {
872         attrib(SM[2],"isIsolatedSingularity",1);
873         if(homog(SM[2])){attrib(SM[2],"onlySingularAtZero",1);}
874      }
875      if(dim(JM[1])<=dim(SM[1])-2)
876      {
877         attrib(SM[2],"isRegInCodim2",1);
878      }
879   }
880   else
881   {
882     if(size(#)==0)
883     {
884        list JM=maxideal(1),maxideal(1);
885
886        attrib(JM[1],"isSB",1);
887        if(dim(SM[1])>=2){attrib(SM[2],"isRegInCodim2",1);}
888     }
889   }
890   if((attrib(SM[2],"isRegInCodim2")==1)&&(attrib(SM[2],"isCohenMacaulay")==1))
891   {
892      if(y>=1)
893      {
894            "// the ideal was CohenMacaulay and regular in codimension 2";
895      }
896      MB=SM[2];
897      intvec rw;
898      list LL=substpart(MB,ihp,0,rw);
899      def newR6=LL[1];
900      setring newR6;
901      ideal norid=endid;
902      ideal normap=endphi;
903      kill endid,endphi;
904      export norid;
905      export normap;
906      result=newR6;
907      result[size(result)+1]=delta;
908      setring BAS;
909      return(result);
910   }
911
912   //if it is an isolated singularity only at 0 things are easier
913   //JM ideal of singular locus, SM ideal of variety
914   if(attrib(SM[2],"onlySingularAtZero"))         
915   {
916      attrib(SM[2],"isIsolatedSingularity",1);
917      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
918      //vars not contained in ideal
919      ideal Ann=quotient(SM[2],SL[1]);
920      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
921      //qAnn=0 ==> the first var(=SL[1]) not contained in SM is a nzd of R/SM
922      if(size(qAnn)==0)
923      {
924         if(y>=1)
925         {
926            "";
927            "//   the ideal rad(J):";
928            "";
929            maxideal(1);
930            newline;
931         }
932         //again test for normality
933         //Hom(I,R)=R
934         list RR;
935         RR=SM[1],SM[2],maxideal(1),SL[1];
936
937         RR=HomJJ(RR,y);
938         if(RR[2]==0)
939         {
940            def newR=RR[1];
941            setring newR;
942            map psi=BAS,endphi;
943            list tluser=normalizationPrimes(endid,psi(ihp),delta+RR[3],an);
944            setring BAS;
945            return(tluser);
946         }
947         MB=SM[2];
948         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
949                      +ordstr(basering)+");");
950         ideal norid=fetch(BAS,MB);
951         ideal normap=fetch(BAS,ihp);
952         delta=delta+RR[3];
953         export norid;
954         export normap;
955         result=newR7;
956         result[size(result)+1]=delta;
957         setring BAS;
958         return(result);
959
960       }
961      //Now the case where qAnn!=0, i.e.SL[1] is a zero divisor of R/SM
962      //and we have found a splitting: id and id1
963      //id=Ann defines components of R/SM in the complement of V(SL[1])
964      //id1 defines components of R/SM in the complement of V(id)
965       else
966       {
967          ideal id=Ann;
968          attrib(id,"isCohenMacaulay",0);
969          attrib(id,"isPrim",0);
970          attrib(id,"isIsolatedSingularity",1);
971          attrib(id,"isRegInCodim2",0);
972          attrib(id,"isCompleteIntersection",0);
973          attrib(id,"isEquidimensional",0);
974          attrib(id,"onlySingularAtZero",1);
975
976          ideal id1=quotient(SM[2],Ann);
977          //ideal id=SL[1],SM[2];
978          attrib(id1,"isCohenMacaulay",0);
979          attrib(id1,"isPrim",0);
980          attrib(id1,"isIsolatedSingularity",1);
981          attrib(id1,"isRegInCodim2",0);
982          attrib(id1,"isCompleteIntersection",0);
983          attrib(id1,"isEquidimensional",0);
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++)
998          {
999             keepresult1=insert(keepresult1,keepresult2[lauf]);
1000          }
1001          keepresult1[size(keepresult1)]=delta;
1002          return(keepresult1);
1003       }
1004   }
1005
1006   //test for non-normality
1007   //Hom(I,I)<>R
1008   //we can use Hom(I,I) to continue
1009   ideal SL=simplify(reduce(JM[2],SM[1]),2);
1010   ideal Ann=quotient(SM[2],SL[1]);
1011   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
1012
1013   if(size(qAnn)==0)
1014   {
1015      list RR;
1016      list RS;
1017      //now we have to compute the radical
1018      if(y>=1)
1019      {
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;
1034        "";
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
1044      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
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
1055      if(RS[2]==1)
1056      {
1057         def lastR=RS[1];
1058         setring lastR;
1059         map psi1=BAS,endphi;
1060         ideal norid=endid;
1061         ideal normap=psi1(ihp);
1062         kill endid,endphi;
1063         export norid;
1064         export normap;
1065         result=lastR;
1066         result[size(result)+1]=delta+RS[3];
1067         setring BAS;
1068         return(result);
1069      }
1070      int n=nvars(basering);
1071      ideal MJ=JM[2];
1072
1073      def newR=RS[1];
1074      setring newR;
1075
1076      map psi=BAS,endphi;
1077      list tluser=
1078             normalizationPrimes(endid,psi(ihp),delta+RS[3],psi(MJ));
1079      setring BAS;
1080      return(tluser);
1081   }
1082    // A component with singular locus the whole component found
1083   if( Ann == 1)
1084   {
1085      "// Input appeared not to be a radical ideal!";
1086      "// A (everywhere singular) component with ideal";
1087      "// equal to its Jacobian ideal was found";
1088      "// Procedure will stop with partial result computed so far";"";
1089
1090         MB=SM[2];
1091         intvec rw;
1092         list LL=substpart(MB,ihp,0,rw);
1093         def newR6=LL[1];
1094         setring newR6;
1095         ideal norid=endid;
1096         ideal normap=endphi;
1097         kill endid,endphi;
1098         export norid;
1099         export normap;
1100         result=newR6;
1101         result[size(result)+1]=delta;
1102         setring BAS;
1103         return(result);
1104   }
1105   else
1106   {
1107      int equi=attrib(SM[2],"isEquidimensional");
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      }
1121      execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
1122                      +ordstr(basering)+");");
1123      if(y>=1)
1124      {
1125         "// zero-divisor found";
1126      }
1127      ideal vid=fetch(BAS,new1);
1128      ideal ihp=fetch(BAS,ihp);
1129      attrib(vid,"isCohenMacaulay",0);
1130      attrib(vid,"isPrim",0);
1131      attrib(vid,"isIsolatedSingularity",isIs);
1132      attrib(vid,"isRegInCodim2",0);
1133      attrib(vid,"onlySingularAtZero",oSAZ);
1134      attrib(vid,"isEquidimensional",equi);     
1135      attrib(vid,"isCompleteIntersection",0);
1136
1137      keepresult1=normalizationPrimes(vid,ihp,0);
1138      int delta1=keepresult1[size(keepresult1)];
1139      setring BAS;
1140
1141      execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
1142                      +ordstr(basering)+");");
1143
1144      ideal vid=fetch(BAS,new2);
1145      ideal ihp=fetch(BAS,ihp);
1146      attrib(vid,"isCohenMacaulay",0);
1147      attrib(vid,"isPrim",0);
1148      attrib(vid,"isIsolatedSingularity",isIs);
1149      attrib(vid,"isRegInCodim2",0);
1150      attrib(vid,"isEquidimensional",equi);     
1151      attrib(vid,"isCompleteIntersection",0);
1152      attrib(vid,"onlySingularAtZero",oSAZ);
1153
1154      keepresult2=normalizationPrimes(vid,ihp,0);
1155      int delta2=keepresult2[size(keepresult2)];
1156
1157      setring BAS;
1158     
1159      for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
1160      {
1161         keepresult1=insert(keepresult1,keepresult2[lauf]);
1162      }
1163      keepresult1[size(keepresult1)]=delta+mul+delta1+delta2;
1164      return(keepresult1);
1165   }
1166}
1167example
1168{ "EXAMPLE:";echo = 2;
1169   // Huneke
1170   ring qr=31991,(a,b,c,d,e),dp;
1171   ideal i=
1172   5abcde-a5-b5-c5-d5-e5,
1173   ab3c+bc3d+a3be+cd3e+ade3,
1174   a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1175   abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1176   ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1177   a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1178   a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1179   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1180
1181   list pr=normalizationPrimes(i);
1182   def r1=pr[1];
1183   setring r1;
1184   norid;
1185   normap;
1186}
1187///////////////////////////////////////////////////////////////////////////////
1188static proc substpart(ideal endid, ideal endphi, int homo, intvec rw)
1189
1190"//Repeated application of elimpart to endid, until no variables can be
1191//directy substituded. homo=1 if input is homogeneous, rw contains
1192//original weights, endphi (partial) normalization map";
1193
1194{
1195   def newRing=basering;
1196   int ii,jj;
1197   map phi = basering,maxideal(1);
1198
1199   list Le = elimpart(endid);
1200                                     //this proc and the next loop try to
1201   int q = size(Le[2]);              //substitute as many variables as possible
1202   intvec rw1 = 0;                   //indices of substituted variables
1203   rw1[nvars(basering)] = 0;
1204   rw1 = rw1+1;
1205
1206   while( size(Le[2]) != 0 )
1207   {
1208      endid = Le[1];
1209      map ps = newRing,Le[5];
1210
1211      phi = ps(phi);
1212      for(ii=1;ii<=size(Le[2])-1;ii++)
1213      {
1214         phi=phi(phi);
1215      }
1216      //eingefuegt wegen x2-y2z2+z3
1217      kill ps;
1218
1219      for( ii=1; ii<=size(rw1); ii++ )
1220      {
1221         if( Le[4][ii]==0 )
1222         {
1223            rw1[ii]=0;                             //look for substituted vars
1224         }
1225      }
1226      Le=elimpart(endid);
1227      q = q + size(Le[2]);
1228   }
1229   endphi = phi(endphi);
1230
1231//---------- return -----------------------------------------------------------
1232// in the homogeneous case put weights for the remaining vars correctly, i.e.
1233// delete from rw those weights for which the corresponding entry of rw1 is 0
1234
1235   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
1236   {
1237      jj=1;
1238      for( ii=2; ii<size(rw1); ii++)
1239      {
1240         jj++;
1241         if( rw1[ii]==0 )
1242         {
1243            rw=rw[1..jj-1],rw[jj+1..size(rw)];
1244            jj=jj-1;
1245         }
1246      }
1247      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
1248      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }
1249
1250      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),(a(rw),dp);
1251   }
1252   else
1253   {
1254      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
1255   }
1256
1257   ideal lastmap;
1258   q = 1;
1259   for(ii=1; ii<=size(rw1); ii++ )
1260   {
1261      if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; }
1262      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
1263   }
1264   map phi1 = newRing,lastmap;
1265   ideal endid  = phi1(endid);
1266   ideal endphi = phi1(endphi);
1267   export(endid);
1268   export(endphi);
1269   list L = lastRing;
1270   setring newRing;
1271   return(L);
1272}
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"
1286{
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          setring newR;
1431          tau=vdim(std(I+jacob(I)));
1432          inv=deltaLoc(I[1],maxideal(1));
1433          delta=inv[1];
1434          tauloc=inv[2];
1435          nodes=tau-tauloc-2*cusps;
1436          nbranch=inv[3]+ 2*nodes+cusps;         
1437          delta=delta+nodes+cusps;
1438        }
1439        else
1440        {
1441           setring newR;
1442           list pr=minAssGTZ(IDsing);                   
1443           for(k=1;k<=size(pr);k++)
1444           {
1445              inv=deltaLoc(I[1],pr[k]);
1446              delta=delta+inv[1];
1447              tausing=tausing+inv[2];
1448              nbranch=nbranch+inv[3];
1449           }
1450           tau=vdim(std(I+jacob(I)));
1451
1452           nodes=tau-tausing-2*cusps;
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);
1469}
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
1558      int s=size(NP);
1559      int nN=-NP[1][2]-NP[s][1]+1;  // computation of the Newton-number
1560      intmat m[2][2];
1561      for(i=1;i<=s-1;i++)
1562      {
1563         m=NP[i+1],NP[i];
1564         nN=nN+det(m);
1565      }
1566
1567      if(mu==nN)                   // the Newton-polygon is non-degenerate
1568      {                            // compute nb, the number of branches
1569        for(i=1;i<=s-1;i++)
1570        {
1571          nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]);
1572        }
1573        return(list(d*(mu+nb-1)/2,d*tau,d*nb));
1574      }
1575
1576   //da reddevelop nur benutzt wird, um die Anzahl der Zweige zu bestimmen
1577   //kann man das sicher schneller machen:
1578   //die Aufblasung durchfuehren und stets testen, ob das Newton-polyeder
1579   //nicht ausgeartet ist.
1580
1581      if(s>2)    //  splitting of f
1582      {
1583         intvec v=NP[1][2]-NP[2][2],NP[2][1];
1584         int de=w_deg(g,v);
1585         int st=w_deg(hc,v)+v[1]+v[2];
1586         poly f1=var(2)^NP[2][2];
1587         poly f2=jet(g,de,v)/var(2)^NP[2][2];
1588         poly h=g-f1*f2;
1589         de=w_deg(h,v);
1590         poly k;
1591         ideal wi=var(2)^NP[2][2],f2;
1592         matrix li;
1593         while(de<st)
1594         {
1595           k=jet(h,de,v);
1596           li=lift(wi,k);
1597           f1=f1+li[2,1];
1598           f2=f2+li[1,1];
1599           h=g-f1*f2;
1600           de=w_deg(h,v);
1601         }
1602         nb=deltaLoc(f1,maxideal(1))[3]+deltaLoc(f2,maxideal(1))[3];
1603         return(list(d*(mu+nb-1)/2,d*tau,d*nb));
1604      }
1605
1606      f=jet(f,deg(hc)+2);
1607      list hne=reddevelop(f);
1608      nb=size(hne);
1609      setring R;
1610      kill HNEring;       
1611      return(list(d*(mu+nb-1)/2,d*tau,d*nb));
1612   }
1613   else             //the case of small characteristic
1614   {
1615      f=jet(f,deg(hc)+2);
1616      list hne=reddevelop(f);
1617      nb=size(hne);
1618      if(nb==1)
1619      {
1620         delta=invariants(hne[1])[5]/2;
1621         setring R;
1622         kill HNEring;
1623         return(list(d*delta,d*tau,d));
1624      }
1625      setring R;
1626      kill HNEring;
1627      //delta direkt aus reddevelop zurueckgeben
1628      ERROR("the case of small characteristic is not fully implemented yet");
1629   }
1630}
1631
1632proc w_deg(poly p, intvec v)
1633{
1634   if(p==0){return(-1);}
1635   int d=0;
1636   while(jet(p,d,v)==0){d++;}
1637   d=(transpose(leadexp(jet(p,d,v)))*v)[1];
1638   return(d); 
1639}
1640
1641proc newton (poly f)
1642{
1643   def R1=basering;
1644   execute("ring R2=("+charstr(R1)+"),("+varstr(R1)+"),ls;");
1645   poly f=imap(R1,f);
1646   intvec A=(0,ord(subst(f,var(1),0)));
1647   intvec B=(ord(subst(f,var(2),0)),0);
1648   intvec C,H; list L;
1649   int abbruch,i;
1650   poly hilf;
1651   L[1]=A;
1652   f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
1653   map xytausch=R2,var(2),var(1);
1654   for (i=2; f!=0; i++)
1655   {
1656      abbruch=0;
1657      while (abbruch==0)
1658      {
1659         C=leadexp(f);         
1660         if(jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
1661         {
1662            abbruch=1;
1663         }       
1664         else
1665         {
1666            f=jet(f,-C[1]-1,intvec(-1,0));
1667         }
1668     }
1669     hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));
1670     H=leadexp(xytausch(hilf));
1671     A=H[2],H[1];
1672     L[i]=A;
1673     f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]-A[1]));
1674   }
1675   L[i]=B;
1676   setring R1;
1677   return(L);
1678}
1679
1680///////////////////////////////////////////////////////////////////////////
1681
1682/*
1683                           Examples:
1684LIB"normal.lib";
1685//Huneke
1686ring qr=31991,(a,b,c,d,e),dp;
1687ideal i=
16885abcde-a5-b5-c5-d5-e5,
1689ab3c+bc3d+a3be+cd3e+ade3,
1690a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1691abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1692ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1693a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1694a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1695b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1696
1697
1698//Vasconcelos (dauert laenger: 70 sec)
1699ring r=32003,(x,y,z,w,t),dp;
1700ideal i=
1701x2+zw,
1702y3+xwt,
1703xw3+z3t+ywt2,
1704y2w4-xy2z2t-w3t3;
1705
1706//Theo1
1707ring r=32003,(x,y,z),wp(2,3,6);
1708ideal i=zy2-zx3-x6;
1709
1710//Theo1a (CohenMacaulay and regular in codimension 2)
1711ring r=32003,(x,y,z,u),wp(2,3,6,6);
1712ideal i=zy2-zx3-x6+u2;
1713
1714
1715//Theo2
1716ring r=32003,(x,y,z),wp(3,4,12);
1717ideal i=z*(y3-x4)+x8;
1718
1719//Theo2a
1720ring r=32003,(T(1..4)),wp(3,4,12,17);
1721ideal i=
1722T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
1723T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
1724T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
1725T(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;
1726
1727//Theo3
1728ring r=32003,(x,y,z),wp(3,5,15);
1729ideal i=z*(y3-x5)+x10;
1730
1731
1732//Theo4
1733ring r=32003,(x,y,z),dp;
1734ideal i=(x-y)*(x-z)*(y-z);
1735
1736//Theo5
1737ring r=32003,(x,y,z),wp(2,1,2);
1738ideal i=z3-xy4;
1739
1740//Theo6
1741ring r=32003,(x,y,z),dp;
1742ideal i=x2y2+x2z2+y2z2;
1743
1744ring r=32003,(a,b,c,d,e,f),dp;
1745ideal i=
1746bf,
1747af,
1748bd,
1749ad;
1750
1751//ist CM
1752//Sturmfels
1753ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1754ideal i=
1755bv+su,
1756bw+tu,
1757sw+tv,
1758by+sx,
1759bz+tx,
1760sz+ty,
1761uy+vx,
1762uz+wx,
1763vz+wy,
1764bvz;
1765
1766//J S/Y
1767ring r=32003,(x,y,z,t),dp;
1768ideal i=
1769x2z+xzt,
1770xyz,
1771xy2-xyt,
1772x2y+xyt;
1773
1774//St_S/Y
1775ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1776ideal i=
1777wy-vz,
1778vx-uy,
1779tv-sw,
1780su-bv,
1781tuy-bvz;
1782
1783//Horrocks:
1784ring r=32003,(a,b,c,d,e,f),dp;
1785ideal i=
1786adef-16000be2f+16001cef2,
1787ad2f+8002bdef+8001cdf2,
1788abdf-16000b2ef+16001bcf2,
1789a2df+8002abef+8001acf2,
1790ad2e-8000bde2-7999cdef,
1791acde-16000bce2+16001c2ef,
1792a2de-8000abe2-7999acef,
1793acd2+8002bcde+8001c2df,
1794abd2-8000b2de-7999bcdf,
1795a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
1796abde-8000b2e2-acdf-16001bcef-8001c2f2,
1797abcd-16000b2ce+16001bc2f,
1798a2cd+8002abce+8001ac2f,
1799a2bd-8000ab2e-7999abcf,
1800ab3f-3bdf3,
1801a2b2f-2adf3-16000bef3+16001cf4,
1802a3bf+4aef3,
1803ac3e-10668cde3,
1804a2c2e+10667ade3+16001be4+5334ce3f,
1805a3ce+10669ae3f,
1806bc3d+8001cd3e,
1807ac3d+8000bc3e+16001cd2e2+8001c4f,
1808b2c2d+16001ad4+4000bd3e+12001cd3f,
1809b2c2e-10668bc3f-10667cd2ef,
1810abc2e-cde2f,
1811b3cd-8000bd3f,
1812b3ce-10668b2c2f-10667bd2ef,
1813abc2f-cdef2,
1814a2bce-16000be3f+16001ce2f2,
1815ab3d-8000b4e-8001b3cf+16000bd2f2,
1816ab2cf-bdef2,
1817a2bcf-16000be2f2+16001cef3,
1818a4d-8000a3be+8001a3cf-2ae2f2;
1819
1820
1821ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1822ideal k=
1823wy-vz,
1824vx-uy,
1825tv-sw,
1826su-bv,
1827tuy-bvz;
1828ideal j=x2y2+x2z2+y2z2;
1829ideal i=mstd(intersect(j,k))[2];
1830
1831
1832ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1833ideal i=
1834wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
1835vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
1836tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
1837sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
1838tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;
1839
1840
1841//riemenschneider 
1842ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
1843ideal i=
1844xz,
1845vx,
1846ux,
1847su,
1848qu,
1849txy,
1850stx,
1851qtx,
1852uv2z-uwz,
1853uv3-uvw,
1854puv2-puw;
1855
1856ring r=0,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
1857ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;
1858
1859//Yoshihiko Sakai
1860ring r=0,(x,y),dp;   //genus 0 4 nodes and 6 cusps   
1861ideal i=(x2+y^2-1)^3 +27x2y2;
1862
1863ring r=0,(x,y),dp;   //genus 0
1864ideal i=(x-y^2)^2 - y*x^3;
1865
1866ring r=0,(x,y),dp;  //genus 4
1867ideal i=y3-x6+1;
1868
1869int m=9;           // q=9: genus 0
1870int p=2;
1871int q=9;//2,...,9
1872ring r=0,(x,y),dp; 
1873ideal i=y^m - x^p*(x - 1)^q;
1874
1875ring r=0,(x,y),dp;  //genus 19
1876ideal 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;
1877
1878
1879ring r=0,(x,y),dp;  //genus 34
1880ideal i=y10+(-2494x2+474)*y8+(84366+2042158x4-660492)*y6
1881        +(128361096x4-47970216x2+6697080-761328152x6)*y4
1882        +(-12024807786x4-506101284x2+15052058268x6+202172841-3212x8)*y2
1883        +34263110700x4-228715574724x6+5431439286x2+201803238
1884        -9127158539954x10-3212722859346x8;
1885
1886//Rob Koelman
1887ring r=0,(x,y,z),dp;//genus 10  with 26 cusps
1888ideal i=
1889761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+
1890 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2-
1891 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4-
1892 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+
1893 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-
1894 z^10-474*z^8*y^2-84366*z^6*y^4;
1895
1896ring r=0,(x,y),dp;//genus 10  with 26 cusps
1897ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
1898-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
1899+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
1900-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;
1901
1902
1903ring r=0,(x,y),dp;   // genuss 1  with 5 cusps
1904ideal i=57y5+516x4y-320x4+66y4-340x2y3+73y3+128x2-84x2y2-96x2y;;
1905
1906//Mark van Hoeij
1907ring r=0,(x,y),dp;  //genus 19
1908ideal i=y20+y13x+x4y5+x3*(x+1)^2;
1909
1910ring r=0,(x,y),dp;  //genus 35
1911ideal i=y30+y13x+x4y5+x3*(x+1)^2;
1912
1913ring r=0,(x,y),dp;   //genus 55
1914ideal i=y40+y13x+x4y5+x3*(x+1)^2;
1915
1916ring r=0,(x,y),dp;  //genus 4
1917ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
1918
1919*/
1920
Note: See TracBrowser for help on using the repository browser.