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

spielwiese
Last change on this file since 3c4dcc was 3c4dcc, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: DOS->UNIX and white space cleanup git-svn-id: file:///usr/local/Singular/svn/trunk@8073 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 69.3 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: normal.lib,v 1.44 2005-05-06 14:38:50 hannes 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
9MAIN PROCEDURES:
10 normal(I[,wd]);     computes the normalization of basering/I,
11                     resp. computes the normalization of basering/I and
12                     the delta invariant
13 HomJJ(L);           presentation of End_R(J) as affine ring, L a list
14 genus(I);           computes genus of the projective curve defined by I
15 primeClosure(L);    integral closure of R/p (p prime), L a list
16 closureFrac(L)      write poly in integral closure as element of Quot(R/p)
17
18AUXILIARY PROCEDURE:
19 deltaLoc(f,S);      (sum of) delta invariant(s) at conjugated singular points
20";
21
22LIB "general.lib";
23LIB "poly.lib";
24LIB "sing.lib";
25LIB "primdec.lib";
26LIB "elim.lib";
27LIB "presolve.lib";
28LIB "inout.lib";
29LIB "ring.lib";
30LIB "hnoether.lib";
31LIB "reesclos.lib";
32///////////////////////////////////////////////////////////////////////////////
33
34proc HomJJ (list Li)
35"USAGE:   HomJJ (Li);  Li = list: ideal SBid, ideal id, ideal J, poly p
36ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
37@*       SBid = standard basis of id,
38@*       J    = ideal of P containing the polynomial p,
39@*       p    = nonzero divisor of R
40COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
41         affine ring, together with the canonical map R --> Hom_R(J,J),
42         where R is the quotient ring of P modulo the standard basis SBid.
43RETURN:  a list l of two objects
44@format
45         l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
46               such that l[1]/endid = Hom_R(J,J) and
47               endphi describes the canonical map R -> Hom_R(J,J)
48         l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
49         l[3] : an integer, the contribution to delta
50@end format
51NOTE:    printlevel >=1: display comments (default: printlevel=0)
52EXAMPLE: example HomJJ;  shows an example
53"
54{
55//---------- initialisation ---------------------------------------------------
56   int isIso,isPr,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
57   intvec rw,rw1;
58   list L;
59   y = printlevel-voice+2;  // y=printlevel (default: y=0)
60   def P = basering;
61   ideal SBid, id, J = Li[1], Li[2], Li[3];
62   poly p = Li[4];
63   attrib(SBid,"isSB",1);
64   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not
65
66//---- set attributes for special cases where algorithm can be simplified -----
67   if( homo==1 )
68   {
69      rw = ringweights(P);
70   }
71   if( typeof(attrib(id,"isPrim"))=="int" )
72   {
73      if(attrib(id,"isPrim")==1)  { isPr=1; }
74   }
75   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
76   {
77      if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; }
78   }
79   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
80   {
81      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
82   }
83   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
84   {
85      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
86   }
87   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
88   {
89      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
90   }
91   if( typeof(attrib(id,"isEquidimensional"))=="int" )
92   {
93      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
94   }
95//-------------------------- go to quotient ring ------------------------------
96   qring R  = SBid;
97   ideal id = fetch(P,id);
98   ideal J  = fetch(P,J);
99   poly p   = fetch(P,p);
100   ideal f,rf,f2;
101   module syzf;
102//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
103   if ( y>=1 )
104   {
105     "// compute p*Hom(J,J) = p*J:J, p a non-zerodivisor";
106     "//   p is equal to:"; "";
107     p;
108     "";
109   }
110   f  = quotient(p*J,J);
111
112   if ( y>=1 )
113   { "// the module p*Hom(J,J) = p*J:J, p a non-zerodivisor";
114      "// p"; p;
115      "// f=p*J:J";f;
116   }
117   f2 = std(p);
118
119//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
120
121   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R
122   if ( size(rf) == 0 )
123   {
124      if ( homog(f) && find(ordstr(basering),"s")==0 )
125      {
126         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
127      }
128      else
129      {
130         ring newR1 = char(P),(X(1..nvars(P))),dp;
131      }
132      ideal endphi = maxideal(1);
133      ideal endid = fetch(P,id);
134      L=substpart(endid,endphi,homo,rw);
135      def lastRing=L[1];
136      setring lastRing;
137
138      attrib(endid,"onlySingularAtZero",oSAZ);
139      attrib(endid,"isCohenMacaulay",isCo);
140      attrib(endid,"isPrim",isPr);
141      attrib(endid,"isIsolatedSingularity",isIso);
142      attrib(endid,"isRegInCodim2",isRe);
143      attrib(endid,"isEqudimensional",isEq);
144      attrib(endid,"isCompleteIntersection",0);
145      attrib(endid,"isRad",0);
146      L=lastRing;
147      L = insert(L,1,1);
148      dbprint(y,"// case R = Hom(J,J)");
149      if(y>=1)
150      {
151         "//   R=Hom(J,J)";
152         "   ";
153         lastRing;
154         "   ";
155         "//   the new ideal";
156         endid;
157         "   ";
158         "//   the old ring";
159         "   ";
160         P;
161         "   ";
162         "//   the old ideal";
163         "   ";
164         setring P;
165         id;
166         "   ";
167         setring lastRing;
168         "//   the map";
169         "   ";
170         endphi;
171         "   ";
172         pause();
173         newline;
174      }
175      setring P;
176      L[3]=0;
177      return(L);
178   }
179   if(y>=1)
180   {
181      "// R is not equal to Hom(J,J), we have to try again";
182      pause();
183      newline;
184   }
185//---------- Hom(J,J) != R: create new ring and map from old ring -------------
186// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
187
188   f=mstd(f)[2];
189   ideal ann=quotient(f2,f);
190   int delt;
191   if(isIso&&isEq){delt=vdim(std(modulo(f,ideal(p))));}
192
193   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
194   q = size(f);
195
196   syzf = syz(f);
197
198   if ( homo==1 )
199   {
200      rw1 = rw,0;
201      for ( ii=2; ii<=q; ii++ )
202      {
203         rw  = rw, deg(f[ii])-deg(f[1]);
204         rw1 = rw1, deg(f[ii])-deg(f[1]);
205      }
206      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
207   }
208   else
209   {
210      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
211   }
212
213   map psi1 = P,maxideal(1);
214   ideal SBid = psi1(SBid);
215   attrib(SBid,"isSB",1);
216
217   qring newR = std(SBid);
218
219   map psi = R,ideal(X(1..nvars(R)));
220   ideal id = psi(id);
221   ideal f = psi(f);
222   module syzf = psi(syzf);
223   ideal pf,Lin,Quad,Q;
224   matrix T,A;
225   list L1;
226
227//---------- computation of Hom(J,J) as affine ring ---------------------------
228// determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
229// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course
230// the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
231// It is a fact, that the kernel is generated by the linear and the quadratic
232// relations
233
234   pf = f[1]*f;
235   T = matrix(ideal(T(1..q)),1,q);
236   Lin = ideal(T*syzf);
237   if(y>=1)
238   {
239      "// the ring structure of Hom(J,J) as R-algebra";
240      " ";
241      "//   the linear relations";
242      " ";
243      Lin;
244      "   ";
245   }
246
247   for (ii=2; ii<=q; ii++ )
248   {
249      for ( jj=2; jj<=ii; jj++ )
250      {
251         A = lift(pf,f[ii]*f[jj]);
252         Quad = Quad, ideal(T(jj)*T(ii) - T*A);          // quadratic relations
253      }
254   }
255
256   if(y>=1)
257   {
258      "//   the quadratic relations";
259      "   ";
260      interred(Quad);
261      pause();
262      newline;
263   }
264   Q = Lin,Quad;
265   Q = subst(Q,T(1),1);
266   Q=mstd(Q)[2];
267//---------- reduce number of variables by substitution, if possible ----------
268   if (homo==1)
269   {
270      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
271   }
272   else
273   {
274      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
275   }
276
277   ideal endid  = imap(newR,id),imap(newR,Q);
278   ideal endphi = ideal(X(1..nvars(R)));
279
280   L=substpart(endid,endphi,homo,rw);
281
282   def lastRing=L[1];
283   setring lastRing;
284
285   attrib(endid,"onlySingularAtZero",0);
286   map sigma=R,endphi;
287   ideal an=sigma(ann);
288   export(an);  //noetig?
289   ideal te=an,endid;
290   if(isIso&&(size(reduce(te,std(maxideal(1))))==0))
291   {
292      attrib(endid,"onlySingularAtZero",oSAZ);
293   }
294   kill te;
295   attrib(endid,"isCohenMacaulay",isCo);
296   attrib(endid,"isPrim",isPr);
297   attrib(endid,"isIsolatedSingularity",isIso);
298   attrib(endid,"isRegInCodim2",isRe);
299   attrib(endid,"isEquidimensional",isEq);
300   attrib(endid,"isCompleteIntersection",0);
301   attrib(endid,"isRad",0);
302   if(y>=1)
303   {
304      "//   the new ring after reduction of the number of variables";
305      "   ";
306      lastRing;
307      "   ";
308      "//   the new ideal";
309      "   ";
310      endid;
311      "   ";
312      "//   the old ring";
313      "   ";
314      P;
315      "   ";
316      "//   the old ideal";
317      "   ";
318      setring P;
319      id;
320      "   ";
321      setring lastRing;
322      "//   the map";
323      "   ";
324      endphi;
325      "   ";
326      pause();
327      newline;
328   }
329   L = lastRing;
330   L = insert(L,0,1);
331   L[3]=delt;
332   return(L);
333}
334example
335{"EXAMPLE:";  echo = 2;
336  ring r   = 0,(x,y),wp(2,3);
337  ideal id = y^2-x^3;
338  ideal J  = x,y;
339  poly p   = x;
340  list Li = std(id),id,J,p;
341  list L   = HomJJ(Li);
342  def end = L[1];    // defines ring L[1], containing ideals endid, endphi
343  setring end;       // makes end the basering
344  end;
345  endid;             // end/endid is isomorphic to End(r/id) as ring
346  map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
347  psi;
348}
349
350///////////////////////////////////////////////////////////////////////////////
351
352proc normal(ideal id, list #)
353"USAGE:   normal(i [,choose]);  i a radical ideal, choose empty, 1 ,\"g\" \"wd\"
354         if choose=1 the normalization of the associated primes is computed
355         (which is sometimes more efficient);
356         if @code{choose=\"wd\"} the delta invariant is computed
357         simultaneously; this may take much more time in the reducible case,
358         since the factorizing standard basis algorithm cannot be used.
359         if @code{choose=\"g\"} the generators of the integral closure
360         are computed.
361ASSUME:  The ideal must be radical, for non-radical ideals the output may
362         be wrong (i=radical(i); makes i radical)
363RETURN:  a list of rings, say nor and in case of @code{choose=\"wd\"} an
364         integer at the end of the list.
365         In case of @code{choose=\"g\"} it returns a list of polynomials
366         g_1,...,g_k+1 at the end of the list such that g_1/g_k+1,...
367         g_k/g_k+1 generate the integral closure.
368         Each ring @code{nor[i]} contains two ideals with given names
369         @code{norid} and @code{normap} such that@*
370         - the direct sum of the rings @code{nor[i]/norid} is the
371         normalization of basering/id;@*
372         - @code{normap} gives the normalization map from basering/id to
373         @code{nor[i]/norid} (for each i).
374NOTE:    to use the i-th ring type: @code{def R=nor[i]; setring R;}.
375@*       Increasing printlevel displays more comments (default: printlevel=0).
376@*       Not implemented for local or mixed orderings.
377@*       If the input ideal i is weighted homogeneous a weighted ordering may
378         be used (qhweight(i); computes weights).
379KEYWORDS: normalization; delta invariant.
380EXAMPLE: example normal; shows an example
381"
382{
383   int i,j,y,withdelta,withgens;
384   string sr;
385   list result,prim,keepresult;
386   y = printlevel-voice+2;
387   if(size(#)>0)
388   {
389     if(typeof(#[1])=="string")
390     {
391       if(#[1]=="g")
392       {
393          withgens=1;
394       }
395       else
396       {
397          withdelta=1;
398       }
399       kill #;
400       list #;
401     }
402   }
403   attrib(id,"isRadical",1);
404   if ( ord_test(basering) != 1)
405   {
406     "";
407     "// Not implemented for this ordering,";
408     "// please change to global ordering!";
409     return(result);
410   }
411   if(withgens)
412   {
413      list pr=minAssGTZ(id);
414      if(size(pr)>1){ERROR("ideal is not prime");}
415      def R=basering;
416      if(defined(ker)){kill ker;}
417      ideal ker=id;
418      export(ker);
419      list l=R;
420      l=primeClosure(l);
421      list gens=closureGenerators(l);
422      list resu=l[size(l)],gens;
423      dbprint(y+1,"
424// 'normal' created a list of one ring.
425// nor[2] is a list of elements of the basering such that
426// nor[2][i]/nor[2][size(nor[2])] generate the integral closure.
427// To see the rings, type (if the name of your list is nor):
428     show( nor);
429// To access the 1-st ring and map (similar for the others), type:
430     def R = nor[1]; setring R;  norid; normap;
431// R/norid is the 1-st ring of the normalization and
432// normap the map from the original basering to R/norid");
433
434      return(resu);
435   }
436   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
437   {
438      if(attrib(id,"isCompleteIntersection")==1)
439      {
440         attrib(id,"isCohenMacaulay",1);
441         attrib(id,"isEquidimensional",1);
442      }
443   }
444   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
445   {
446      if(attrib(id,"isCohenMacaulay")==1)
447      {
448         attrib(id,"isEquidimensional",1);
449      }
450   }
451   if( typeof(attrib(id,"isPrim"))=="int" )
452   {
453      if(attrib(id,"isPrim")==1)
454      {
455         attrib(id,"isEquidimensional",1);
456      }
457   }
458   if(size(#)==0)
459   {
460      if( typeof(attrib(id,"isEquidimensional"))=="int" )
461      {
462         if(attrib(id,"isEquidimensional")==1)
463         {
464            prim[1]=id;
465         }
466         else
467         {
468            prim=equidim(id);
469         }
470      }
471      else
472      {
473         prim=equidim(id);
474      }
475      if(y>=1)
476      {
477         "// we have ",size(prim),"equidimensional components";
478      }
479      if(withdelta &&(size(prim)>1))
480      {
481         withdelta=0;
482         "WARNING!  cannot compute delta,";
483         "the ideal is not equidimensional";
484      }
485      if(!withdelta)
486      {
487         option(redSB);
488         for(j=1;j<=size(prim);j++)
489         {
490            keepresult=keepresult+facstd(prim[j]);
491         }
492         prim=keepresult;
493         option(noredSB);
494      }
495   }
496   else
497   {
498      if( typeof(attrib(id,"isPrim"))=="int" )
499      {
500         if(attrib(id,"isPrim")==1)
501         {
502            prim[1]=id;
503         }
504         else
505         {
506            prim=minAssPrimes(id);
507         }
508      }
509      else
510      {
511         prim=minAssPrimes(id);
512      }
513      if(y>=1)
514      {
515         "// we have ",size(prim),"irreducible components";
516      }
517   }
518   for(i=1; i<=size(prim); i++)
519   {
520      if(y>=1)
521      {
522         "// we are in loop ",i;
523      }
524      attrib(prim[i],"isCohenMacaulay",0);
525      if(size(#)!=0)
526      {
527         attrib(prim[i],"isPrim",1);
528      }
529      else
530      {
531         attrib(prim[i],"isPrim",0);
532      }
533      attrib(prim[i],"isRegInCodim2",0);
534      attrib(prim[i],"isIsolatedSingularity",0);
535      attrib(prim[i],"isEquidimensional",1);
536      attrib(prim[i],"isCompleteIntersection",0);
537      attrib(prim[i],"onlySingularAtZero",0);
538      if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
539      {
540            if(attrib(id,"onlySingularAtZero")==1)
541             {attrib(prim[i],"onlySingularAtZero",1); }
542      }
543
544      if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
545      {
546            if(attrib(id,"isIsolatedSingularity")==1)
547             {attrib(prim[i],"isIsolatedSingularity",1); }
548      }
549
550      if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
551      {
552            if((attrib(id,"isIsolatedSingularity")==1)&&(size(#)==0))
553             {attrib(prim[i],"isIsolatedSingularity",1); }
554      }
555      keepresult=normalizationPrimes(prim[i],maxideal(1),0);
556
557      for(j=1;j<=size(keepresult)-1;j++)
558      {
559         result=insert(result,keepresult[j]);
560      }
561      sr = string(size(result));
562   }
563   if(withdelta)
564   {
565      sr = string(size(keepresult)-1);
566      result=keepresult;
567   }
568      dbprint(y+1,"
569// 'normal' created a list of "+sr+" ring(s).");
570      if(withdelta)
571      {
572        dbprint(y+1,"// nor["+sr+"+1] is the delta-invariant.");
573      }
574      dbprint(y+1,"// To see the rings, type (if the name of your list is nor):
575     show( nor);
576// To access the 1-st ring and map (similar for the others), type:
577     def R = nor[1]; setring R;  norid; normap;
578// R/norid is the 1-st ring of the normalization and
579// normap the map from the original basering to R/norid");
580
581      //kill endphi,endid;
582      return(result);
583}
584example
585{ "EXAMPLE:"; echo = 2;
586   ring r=32003,(x,y,z),wp(2,1,2);
587   ideal i=z3-xy4;
588   list nor=normal(i);
589   show(nor);
590   def r1=nor[1];
591   setring r1;
592   norid;
593   normap;
594
595   ring s=0,(x,y),dp;
596   ideal i=(x-y^2)^2 - y*x^3;
597   nor=normal(i,"wd");
598   //the delta-invariant
599   nor[size(nor)];
600   nor=normal(i,"g");
601   nor[size(nor)];
602}
603
604///////////////////////////////////////////////////////////////////////////////
605static proc normalizationPrimes(ideal i,ideal ihp,int delt, list #)
606"USAGE:   normalizationPrimes(i,ihp[,si]);  i equidimensional ideal, ihp map
607         (partial normalization),delta partial delta-invariant,
608          # ideal of singular locus
609RETURN:   a list of rings, say nor, and an integer, the delta-invariant
610          at the end of the list.
611          each ring nor[i] contains two ideals
612          with given names norid and normap such that
613           - the direct sum of the rings nor[i]/norid is
614             the normalization of basering/id;
615           - normap gives the normalization map from basering/id
616             to nor[i]/norid (for each i)
617EXAMPLE: example normalizationPrimes; shows an example
618"
619{
620
621   int y = printlevel-voice+2;  // y=printlevel (default: y=0)
622
623   if(y>=1)
624   {
625     "";
626     "// START a normalization loop with the ideal";  "";
627     i;  "";
628     basering;  "";
629     pause();
630     newline;
631   }
632
633   def BAS=basering;
634   list result,keepresult1,keepresult2;
635   ideal J,SB,MB;
636   int depth,lauf,prdim;
637   int ti=timer;
638
639   if(size(i)==0)
640   {
641      if(y>=1)
642      {
643          "// the ideal was the zero-ideal";
644      }
645         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
646                      +ordstr(basering)+");");
647         ideal norid=ideal(0);
648         ideal normap=fetch(BAS,ihp);
649         export norid;
650         export normap;
651         result=newR7;
652         result[size(result)+1]=delt;
653         setring BAS;
654         return(result);
655   }
656
657   if(y>=1)
658   {
659     "// SB-computation of the input ideal";
660   }
661   list SM=mstd(i);                //here the work starts
662   int dimSM =  dim(SM[1]);        //dimension of variety to normalize
663  // Case: Get an ideal containing a unit
664   if( dimSM == -1)
665   {  "";
666      "      // A unit ideal was found.";
667      "      // Stop with partial result computed so far";"";
668
669         MB=SM[2];
670         intvec rw;
671         list LL=substpart(MB,ihp,0,rw);
672         def newR6=LL[1];
673         setring newR6;
674         ideal norid=endid;
675         ideal normap=endphi;
676         kill endid,endphi;
677         export norid;
678         export normap;
679         result=newR6;
680         result[size(result)+1]=delt;
681         setring BAS;
682         return(result);
683   }
684
685   if(y>=1)
686   {
687      "//   the dimension is:"; "";
688      dimSM;"";
689   }
690   if(size(#)>0)
691   {
692      if(attrib(i,"onlySingularAtZero"))
693      {
694         list JM=maxideal(1),maxideal(1);
695         attrib(JM[1],"isSB",1);
696         attrib(JM[2],"isRad",1);
697      }
698      else
699      {
700         ideal te=#[1],SM[2];
701         list JM=mstd(te);
702         kill te;
703         if( typeof(attrib(#[1],"isRad"))!="int" )
704         {
705            attrib(JM[2],"isRad",0);
706         }
707      }
708   }
709
710   if(attrib(i,"isPrim")==1)
711   {
712      attrib(SM[2],"isPrim",1);
713   }
714   else
715   {
716      attrib(SM[2],"isPrim",0);
717   }
718   if(attrib(i,"isIsolatedSingularity")==1)
719   {
720      attrib(SM[2],"isIsolatedSingularity",1);
721   }
722   else
723   {
724      attrib(SM[2],"isIsolatedSingularity",0);
725   }
726   if(attrib(i,"isCohenMacaulay")==1)
727   {
728      attrib(SM[2],"isCohenMacaulay",1);
729   }
730   else
731   {
732      attrib(SM[2],"isCohenMacaulay",0);
733   }
734   if(attrib(i,"isRegInCodim2")==1)
735   {
736      attrib(SM[2],"isRegInCodim2",1);
737   }
738   else
739   {
740      attrib(SM[2],"isRegInCodim2",0);
741   }
742   if(attrib(i,"isEquidimensional")==1)
743   {
744      attrib(SM[2],"isEquidimensional",1);
745   }
746   else
747   {
748      attrib(SM[2],"isEquidimensional",0);
749   }
750   if(attrib(i,"isCompleteIntersection")==1)
751   {
752      attrib(SM[2],"isCompleteIntersection",1);
753   }
754   else
755   {
756      attrib(SM[2],"isCompleteIntersection",0);
757   }
758   if(attrib(i,"onlySingularAtZero")==1)
759   {
760      attrib(SM[2],"onlySingularAtZero",1);
761   }
762   else
763   {
764      attrib(SM[2],"onlySingularAtZero",0);
765   }
766   if((attrib(SM[2],"isIsolatedSingularity")==1)&&(homog(SM[2])==1))
767   {
768      attrib(SM[2],"onlySingularAtZero",1);
769   }
770
771   //the smooth case
772   if(size(#)>0)
773   {
774      if(dim(JM[1])==-1)
775      {
776         if(y>=1)
777         {
778            "// the ideal was smooth";
779         }
780         MB=SM[2];
781         intvec rw;;
782         list LL=substpart(MB,ihp,0,rw);
783         def newR6=LL[1];
784         setring newR6;
785         ideal norid=endid;
786         ideal normap=endphi;
787         kill endid,endphi;
788         export norid;
789         export normap;
790         result=newR6;
791         result[size(result)+1]=delt;
792         setring BAS;
793         return(result);
794     }
795   }
796
797   //the zero-dimensional case
798   if((dim(SM[1])==0)&&(homog(SM[2])==1))
799   {
800      if(y>=1)
801      {
802         "// the ideal was zero-dimensional and homogeneous";
803      }
804      MB=maxideal(1);
805      intvec rw;
806      list LL=substpart(MB,ihp,0,rw);
807      def newR5=LL[1];
808      setring newR5;
809      ideal norid=endid;
810      ideal normap=endphi;
811      kill endid,endphi;
812      export norid;
813      export normap;
814      result=newR5;
815      result[size(result)+1]=delt;
816      setring BAS;
817      return(result);
818   }
819
820   //the one-dimensional case
821   //in this case it is a line because
822   //it is irreducible and homogeneous
823   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
824        &&(homog(SM[2])==1))
825   {
826      if(y>=1)
827      {
828         "// the ideal defines a line";
829      }
830      MB=SM[2];
831      intvec rw;
832      list LL=substpart(MB,ihp,0,rw);
833      def newR4=LL[1];
834      setring newR4;
835      ideal norid=endid;
836      ideal normap=endphi;
837      kill endid,endphi;
838      export norid;
839      export normap;
840      result=newR4;
841      result[size(result)+1]=delt;
842      setring BAS;
843      return(result);
844   }
845
846   //the higher dimensional case
847   //we test first of all CohenMacaulay and
848   //complete intersection
849   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
850   {
851      //test for complete intersection
852      attrib(SM[2],"isCohenMacaulay",1);
853      attrib(SM[2],"isCompleteIntersection",1);
854      attrib(SM[2],"isEquidimensional",1);
855      if(y>=1)
856      {
857         "// the ideal is a complete intersection";
858      }
859   }
860   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(i)==1))
861   {
862      int tau=vdim(std(i+jacob(i)));
863      if(tau>=0)
864      {
865        execute("ring BASL="+charstr(BAS)+",("+varstr(BAS)+"),ds;");
866        ideal i=imap(BAS,i);
867        int tauloc=vdim(std(i+jacob(i)));
868        setring BAS;
869        attrib(SM[2],"onlySingularAtZero",(tau==tauloc));
870        kill BASL;
871      }
872   }
873   //compute the singular locus
874   if((attrib(SM[2],"onlySingularAtZero")==0)&&(size(#)==0))
875   {
876      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
877      if(y >=1 )
878      {
879         "// SB of singular locus will be computed";
880      }
881      ideal sin=J,SM[1];
882      list JM=mstd(sin);
883      attrib(JM[1],"isSB",1);
884
885      //JM[1] SB of singular locus, JM[2]=minbasis of singular locus
886      //SM[1] SB of the input ideal, SM[2] minbasis
887      if(y>=1)
888      {
889         "//   the dimension of the singular locus is:";"";
890         dim(JM[1]); "";
891      }
892
893      if(dim(JM[1])==-1)
894      {
895         if(y>=1)
896         {
897            "// the ideal is smooth";
898         }
899         MB=SM[2];
900         intvec rw;
901         list LL=substpart(MB,ihp,0,rw);
902         def newR3=LL[1];
903         setring newR3;
904         ideal norid=endid;
905         ideal normap=endphi;
906         kill endid,endphi;
907         export norid;
908         export normap;
909         result=newR3;
910         result[size(result)+1]=delt;
911         setring BAS;
912         return(result);
913      }
914      if(dim(JM[1])==0)
915      {
916         attrib(SM[2],"isIsolatedSingularity",1);
917         if(homog(SM[2])){attrib(SM[2],"onlySingularAtZero",1);}
918      }
919      if(dim(JM[1])<=dim(SM[1])-2)
920      {
921         attrib(SM[2],"isRegInCodim2",1);
922      }
923   }
924   else
925   {
926     if(size(#)==0)
927     {
928        list JM=maxideal(1),maxideal(1);
929
930        attrib(JM[1],"isSB",1);
931        if(dim(SM[1])>=2){attrib(SM[2],"isRegInCodim2",1);}
932     }
933   }
934   if((attrib(SM[2],"isRegInCodim2")==1)&&(attrib(SM[2],"isCohenMacaulay")==1))
935   {
936      if(y>=1)
937      {
938            "// the ideal was CohenMacaulay and regular in codimension 2";
939      }
940      MB=SM[2];
941      intvec rw;
942      list LL=substpart(MB,ihp,0,rw);
943      def newR6=LL[1];
944      setring newR6;
945      ideal norid=endid;
946      ideal normap=endphi;
947      kill endid,endphi;
948      export norid;
949      export normap;
950      result=newR6;
951      result[size(result)+1]=delt;
952      setring BAS;
953      return(result);
954   }
955
956   //if it is an isolated singularity only at 0 things are easier
957   //JM ideal of singular locus, SM ideal of variety
958   if(attrib(SM[2],"onlySingularAtZero"))
959   {
960      attrib(SM[2],"isIsolatedSingularity",1);
961      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
962      //vars not contained in ideal
963      ideal Ann=quotient(SM[2],SL[1]);
964      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
965      //qAnn=0 ==> the first var(=SL[1]) not contained in SM is a nzd of R/SM
966      if(size(qAnn)==0)
967      {
968         if(y>=1)
969         {
970            "";
971            "//   the ideal rad(J):";
972            "";
973            maxideal(1);
974            newline;
975         }
976         //again test for normality
977         //Hom(I,R)=R
978         list RR;
979         RR=SM[1],SM[2],maxideal(1),SL[1];
980
981         RR=HomJJ(RR,y);
982         if(RR[2]==0)
983         {
984            def newR=RR[1];
985            setring newR;
986            map psi=BAS,endphi;
987            list tluser=normalizationPrimes(endid,psi(ihp),delt+RR[3],an);
988            setring BAS;
989            return(tluser);
990         }
991         MB=SM[2];
992         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
993                      +ordstr(basering)+");");
994         ideal norid=fetch(BAS,MB);
995         ideal normap=fetch(BAS,ihp);
996         delt=delt+RR[3];
997         export norid;
998         export normap;
999         result=newR7;
1000         result[size(result)+1]=delt;
1001         setring BAS;
1002         return(result);
1003
1004       }
1005      //Now the case where qAnn!=0, i.e.SL[1] is a zero divisor of R/SM
1006      //and we have found a splitting: id and id1
1007      //id=Ann defines components of R/SM in the complement of V(SL[1])
1008      //id1 defines components of R/SM in the complement of V(id)
1009       else
1010       {
1011          ideal id=Ann;
1012          attrib(id,"isCohenMacaulay",0);
1013          attrib(id,"isPrim",0);
1014          attrib(id,"isIsolatedSingularity",1);
1015          attrib(id,"isRegInCodim2",0);
1016          attrib(id,"isCompleteIntersection",0);
1017          attrib(id,"isEquidimensional",0);
1018          attrib(id,"onlySingularAtZero",1);
1019
1020          ideal id1=quotient(SM[2],Ann);
1021          //ideal id=SL[1],SM[2];
1022          attrib(id1,"isCohenMacaulay",0);
1023          attrib(id1,"isPrim",0);
1024          attrib(id1,"isIsolatedSingularity",1);
1025          attrib(id1,"isRegInCodim2",0);
1026          attrib(id1,"isCompleteIntersection",0);
1027          attrib(id1,"isEquidimensional",0);
1028          attrib(id1,"onlySingularAtZero",1);
1029
1030          ideal t1=id,id1;
1031          int mul=vdim(std(t1));
1032          kill t1;
1033
1034          keepresult1=normalizationPrimes(id,ihp,0);
1035
1036          keepresult2=normalizationPrimes(id1,ihp,0);
1037
1038          delt=delt+mul+keepresult1[size(keepresult1)]
1039                         +keepresult1[size(keepresult1)];
1040
1041          for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
1042          {
1043             keepresult1=insert(keepresult1,keepresult2[lauf]);
1044          }
1045          keepresult1[size(keepresult1)]=delt;
1046          return(keepresult1);
1047       }
1048   }
1049
1050   //test for non-normality
1051   //Hom(I,I)<>R
1052   //we can use Hom(I,I) to continue
1053   ideal SL=simplify(reduce(JM[2],SM[1]),2);
1054   ideal Ann=quotient(SM[2],SL[1]);
1055   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
1056
1057   if(size(qAnn)==0)
1058   {
1059      list RR;
1060      list RS;
1061      //now we have to compute the radical
1062      if(y>=1)
1063      {
1064         "// radical computation of singular locus";
1065      }
1066      J=radical(JM[2]);   //the singular locus
1067      JM=mstd(J);
1068
1069      if((vdim(JM[1])==1)&&(size(reduce(J,std(maxideal(1))))==0))
1070      {
1071         attrib(SM[2],"onlySingularAtZero",1);
1072      }
1073      if(y>=1)
1074      {
1075        "//   radical is equal to:";"";
1076        J;
1077        "";
1078      }
1079      if(deg(SL[1])>deg(J[1]))
1080      {
1081         Ann=quotient(SM[2],J[1]);
1082         qAnn=simplify(reduce(Ann,SM[1]),2);
1083         if(size(qAnn)==0){SL[1]=J[1];}
1084      }
1085
1086      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
1087      RR=SM[1],SM[2],JM[2],SL[1];
1088
1089      //   evtl eine geeignete Potenz von JM?
1090     if(y>=1)
1091     {
1092        "// compute Hom(rad(J),rad(J))";
1093     }
1094
1095     RS=HomJJ(RR,y);
1096
1097      if(RS[2]==1)
1098      {
1099         def lastR=RS[1];
1100         setring lastR;
1101         map psi1=BAS,endphi;
1102         ideal norid=endid;
1103         ideal normap=psi1(ihp);
1104         kill endid,endphi;
1105         export norid;
1106         export normap;
1107         result=lastR;
1108         result[size(result)+1]=delt+RS[3];
1109         setring BAS;
1110         return(result);
1111      }
1112      int n=nvars(basering);
1113      ideal MJ=JM[2];
1114
1115      def newR=RS[1];
1116      setring newR;
1117
1118      map psi=BAS,endphi;
1119      list tluser=
1120             normalizationPrimes(endid,psi(ihp),delt+RS[3],psi(MJ));
1121      setring BAS;
1122      return(tluser);
1123   }
1124    // A component with singular locus the whole component found
1125   if( Ann == 1)
1126   {
1127      "// Input appeared not to be a radical ideal!";
1128      "// A (everywhere singular) component with ideal";
1129      "// equal to its Jacobian ideal was found";
1130      "// Procedure will stop with partial result computed so far";"";
1131
1132         MB=SM[2];
1133         intvec rw;
1134         list LL=substpart(MB,ihp,0,rw);
1135         def newR6=LL[1];
1136         setring newR6;
1137         ideal norid=endid;
1138         ideal normap=endphi;
1139         kill endid,endphi;
1140         export norid;
1141         export normap;
1142         result=newR6;
1143         result[size(result)+1]=delt;
1144         setring BAS;
1145         return(result);
1146   }
1147   else
1148   {
1149      int equi=attrib(SM[2],"isEquidimensional");
1150      int oSAZ=attrib(SM[2],"onlySingularAtZero");
1151      int isIs=attrib(SM[2],"isIsolatedSingularity");
1152
1153      ideal new1=Ann;
1154      ideal new2=quotient(SM[2],Ann);
1155      //ideal new2=SL[1],SM[2];
1156
1157      int mul;
1158      if(equi&&isIs)
1159      {
1160          ideal t2=new1,new2;
1161          mul=vdim(std(t2));
1162      }
1163      execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
1164                      +ordstr(basering)+");");
1165      if(y>=1)
1166      {
1167         "// zero-divisor found";
1168      }
1169      ideal vid=fetch(BAS,new1);
1170      ideal ihp=fetch(BAS,ihp);
1171      attrib(vid,"isCohenMacaulay",0);
1172      attrib(vid,"isPrim",0);
1173      attrib(vid,"isIsolatedSingularity",isIs);
1174      attrib(vid,"isRegInCodim2",0);
1175      attrib(vid,"onlySingularAtZero",oSAZ);
1176      attrib(vid,"isEquidimensional",equi);
1177      attrib(vid,"isCompleteIntersection",0);
1178
1179      keepresult1=normalizationPrimes(vid,ihp,0);
1180      int delta1=keepresult1[size(keepresult1)];
1181      setring BAS;
1182
1183      execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
1184                      +ordstr(basering)+");");
1185
1186      ideal vid=fetch(BAS,new2);
1187      ideal ihp=fetch(BAS,ihp);
1188      attrib(vid,"isCohenMacaulay",0);
1189      attrib(vid,"isPrim",0);
1190      attrib(vid,"isIsolatedSingularity",isIs);
1191      attrib(vid,"isRegInCodim2",0);
1192      attrib(vid,"isEquidimensional",equi);
1193      attrib(vid,"isCompleteIntersection",0);
1194      attrib(vid,"onlySingularAtZero",oSAZ);
1195
1196      keepresult2=normalizationPrimes(vid,ihp,0);
1197      int delta2=keepresult2[size(keepresult2)];
1198
1199      setring BAS;
1200
1201      for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
1202      {
1203         keepresult1=insert(keepresult1,keepresult2[lauf]);
1204      }
1205      keepresult1[size(keepresult1)]=delt+mul+delta1+delta2;
1206      return(keepresult1);
1207   }
1208}
1209example
1210{ "EXAMPLE:";echo = 2;
1211   // Huneke
1212   ring qr=31991,(a,b,c,d,e),dp;
1213   ideal i=
1214   5abcde-a5-b5-c5-d5-e5,
1215   ab3c+bc3d+a3be+cd3e+ade3,
1216   a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1217   abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1218   ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1219   a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1220   a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1221   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1222
1223   list pr=normalizationPrimes(i);
1224   def r1=pr[1];
1225   setring r1;
1226   norid;
1227   normap;
1228}
1229///////////////////////////////////////////////////////////////////////////////
1230proc substpart(ideal endid, ideal endphi, int homo, intvec rw)
1231
1232"//Repeated application of elimpart to endid, until no variables can be
1233//directy substituded. homo=1 if input is homogeneous, rw contains
1234//original weights, endphi (partial) normalization map";
1235
1236{
1237   def newRing=basering;
1238   int ii,jj;
1239   map phi = basering,maxideal(1);
1240
1241   list Le = elimpart(endid);
1242                                     //this proc and the next loop try to
1243   int q = size(Le[2]);              //substitute as many variables as possible
1244   intvec rw1 = 0;                   //indices of substituted variables
1245   rw1[nvars(basering)] = 0;
1246   rw1 = rw1+1;
1247
1248   while( size(Le[2]) != 0 )
1249   {
1250      endid = Le[1];
1251      map ps = newRing,Le[5];
1252
1253      phi = ps(phi);
1254      for(ii=1;ii<=size(Le[2])-1;ii++)
1255      {
1256         phi=phi(phi);
1257      }
1258      //eingefuegt wegen x2-y2z2+z3
1259      kill ps;
1260
1261      for( ii=1; ii<=size(rw1); ii++ )
1262      {
1263         if( Le[4][ii]==0 )
1264         {
1265            rw1[ii]=0;                             //look for substituted vars
1266         }
1267      }
1268      Le=elimpart(endid);
1269      q = q + size(Le[2]);
1270   }
1271   endphi = phi(endphi);
1272
1273//---------- return -----------------------------------------------------------
1274// in the homogeneous case put weights for the remaining vars correctly, i.e.
1275// delete from rw those weights for which the corresponding entry of rw1 is 0
1276
1277   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
1278   {
1279      jj=1;
1280      for( ii=2; ii<size(rw1); ii++)
1281      {
1282         jj++;
1283         if( rw1[ii]==0 )
1284         {
1285            rw=rw[1..jj-1],rw[jj+1..size(rw)];
1286            jj=jj-1;
1287         }
1288      }
1289      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
1290      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }
1291
1292      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),(a(rw),dp);
1293   }
1294   else
1295   {
1296      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
1297   }
1298
1299   ideal lastmap;
1300   q = 1;
1301   for(ii=1; ii<=size(rw1); ii++ )
1302   {
1303      if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; }
1304      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
1305   }
1306   map phi1 = newRing,lastmap;
1307   ideal endid  = phi1(endid);
1308   ideal endphi = phi1(endphi);
1309   export(endid);
1310   export(endphi);
1311   list L = lastRing;
1312   setring newRing;
1313   return(L);
1314}
1315/////////////////////////////////////////////////////////////////////////////
1316
1317proc genus(ideal I,list #)
1318"USAGE:   genus(I) or genus(i,1); I a 1-dimensional ideal
1319RETURN:  an integer, the geometric genus p_g = p_a - delta of the projective
1320         curve defined by I, where p_a is the arithmetic genus.
1321NOTE:    delta is the sum of all local delta-invariants of the singularities,
1322         i.e. dim(R'/R), R' the normalization of the local ring R of the
1323         singularity.
1324         genus(i,1) uses the normalization to compute delta. Usually this
1325         is slow but sometimes not.
1326EXAMPLE: example genus; shows an example
1327"
1328{
1329   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
1330
1331   def R0=basering;
1332   I=std(I);
1333   if(dim(I)!=1)
1334   {
1335      if(((homog(I))&&(dim(I)!=2))||(!homog(I)))
1336      {
1337        // ERROR("This is not a curve");
1338        if(w==1){"WARNING:This is not a curve";}
1339      }
1340   }
1341   list L=elimpart(I);
1342   if(size(L[2])!=0)
1343   {
1344      map psi=R0,L[5];
1345      I=std(psi(I));
1346   }
1347   if(size(I)==0)
1348   {
1349      return(0);
1350   }
1351   list N=findvars(I,0);
1352   if(size(N[1])==1)
1353   {
1354
1355      poly p=I[1];
1356     // if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");}
1357      return(-deg(p)+1);
1358   }
1359   if(size(N[1])<nvars(R0))
1360   {
1361      string newvar=string(N[1]);
1362      execute("ring R=("+charstr(R0)+"),("+newvar+"),dp;");
1363      ideal I =imap(R0,I);
1364      attrib(I,"isSB",1);
1365   }
1366   else
1367   {
1368      def R=basering;
1369   }
1370   if(dim(I)==2)
1371   {
1372      def newR=basering;
1373   }
1374   else
1375   {
1376      if(dim(I)==0)
1377      {
1378         execute("ring Rhelp=("+charstr(R0)+"),(@s,@t),dp;");
1379      }
1380      else
1381      {
1382         execute("ring Rhelp=("+charstr(R0)+"),(@s),dp;");
1383      }
1384      def newR=R+Rhelp;
1385      setring newR;
1386      ideal I=imap(R,I);
1387      I=homog(I,@s);
1388      attrib(I,"isSB",1);
1389   }
1390
1391   if((nvars(basering)<=3)&&(size(I)>1))
1392   {
1393       ERROR("This is not equidimensional");
1394   }
1395
1396   intvec hp=hilbPoly(I);
1397   int p_a=1-hp[1];
1398   int d=hp[2];
1399
1400   if(w>=1)
1401   {
1402      "";"The ideal of the projective curve:";"";I;"";
1403      "The coefficients of the Hilbert polynomial";hp;
1404      "arithmetic genus:";p_a;
1405      "degree:";d;"";
1406   }
1407
1408   intvec v = hilb(I,1);
1409   int i,o;
1410   if(nvars(basering)>3)
1411   {
1412      map phi=newR,maxideal(1);
1413      int de;
1414      ideal K,L1;
1415      matrix M;
1416      poly m=var(4);
1417      poly he;
1418      for(i=5;i<=nvars(basering);i++){m=m*var(i);}
1419      K=eliminate(I,m,v);
1420      if(size(K)==1){de=deg(K[1]);}
1421      m=var(1);
1422      for(i=2;i<=nvars(basering)-3;i++){m=m*var(i);}
1423      i=0;
1424      while(d!=de)
1425      {
1426         o=1;
1427         i++;
1428         K=phi(I);
1429         K=eliminate(K,m,v);
1430         if(size(K)==1){de=deg(K[1]);}
1431         if((i==5)&&(d!=de))
1432         {
1433            K=reduce(equidimMax(I),I);
1434            if(size(K)!=0){ERROR("This is not equidimensional");}
1435         }
1436         if(i==10)
1437         {
1438            J;K;
1439            ERROR("genus: did not find a good projection for to
1440                           the plain");
1441         }
1442         if(i<5)
1443         {
1444            M=sparsetriag(nvars(newR),nvars(newR),80-5*i,i);
1445         }
1446         else
1447         {
1448            if(i<8)
1449            {
1450               M=transpose(sparsetriag(nvars(newR),nvars(newR),80-5*i,i));
1451            }
1452            else
1453            {
1454               he=0;
1455               while(he==0)
1456               {
1457                  M=randommat(nvars(newR),nvars(newR),ideal(1),20);
1458                  he=det(M);
1459               }
1460            }
1461         }
1462         L1=M*transpose(maxideal(1));
1463         phi=newR,L1;
1464      }
1465      I=K;
1466   }
1467   poly p=I[1];
1468
1469   execute("ring S=("+charstr(R)+"),(x,y,t),dp;");
1470   ideal L=maxideal(1);
1471   execute("ring C=("+charstr(R)+"),(x,y),ds;");
1472   ideal I;
1473   execute("ring A=("+charstr(R)+"),(x,t),dp;");
1474   map phi=S,1,x,t;
1475   map psi=S,x,1,t;
1476   poly g,h;
1477   ideal I,I1;
1478   execute("ring B=("+charstr(R)+"),(x,t),ds;");
1479
1480   setring S;
1481   if(o)
1482   {
1483     for(i=1;i<=nvars(newR)-3;i++){L[i]=0;}
1484     L=L,maxideal(1);
1485   }
1486   map sigma=newR,L;
1487   poly F=sigma(p);
1488   if(w>=1){"the projected curve:";"";F;"";}
1489
1490   kill newR;
1491
1492   int genus=(d-1)*(d-2)/2;
1493   if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();}
1494
1495   int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing,
1496       tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf;
1497   list inv;
1498
1499   if(w>=1)
1500     {"";"analyse the singularities at oo";"";"singular locus at (1,x,0):";"";}
1501   setring A;
1502   g=phi(F);
1503   h=psi(F);
1504   I=g,jacob(g),var(2);
1505   I=std(I);
1506   if(deg(I[1])>0)
1507   {
1508      list qr=minAssGTZ(I);
1509      if(w>=1){qr;"";}
1510
1511      for(k=1;k<=size(qr);k++)
1512      {
1513         if(w>=1){ nsing=nsing+vdim(std(qr[k]));}
1514         inv=deltaLoc(g,qr[k]);
1515         deltainf=deltainf+inv[1];
1516         tauinf=tauinf+inv[2];
1517         l=vdim(std(qr[k]));
1518         if(inv[2]==l){nodesinf=nodesinf+l;}
1519         if(inv[2]==2*l){cuspsinf=cuspsinf+l;}
1520         nbranchinf=nbranchinf+inv[3];
1521      }
1522   }
1523   else
1524   {
1525     if(w>=1){"            the curve is smooth at (1,x,0)";"";}
1526   }
1527   if(w>=1){"singular locus at (0,1,0):";"";}
1528   inv=deltaLoc(h,maxideal(1));
1529   if((w>=1)&&(inv[2]!=0)){ nsing++;}
1530   deltainf=deltainf+inv[1];
1531   tauinf=tauinf+inv[2];
1532   if(inv[2]==1){nodesinf++;}
1533   if(inv[2]==2){cuspsinf++;}
1534
1535   if((w>=1)&&(inv[2]==0)){"            the curve is smooth at (0,1,0)";"";}
1536   if(inv[2]>0){nbranchinf=nbranchinf+inv[3];}
1537
1538   if(w>=1)
1539   {
1540      if(tauinf==0)
1541      {
1542        "            the curve is smooth at oo";"";
1543      }
1544      else
1545      {
1546         "number of singularities at oo:";nsing;
1547         "nodes at oo:";nodesinf;
1548         "cusps at oo:";cuspsinf;
1549         "branches at oo:";nbranchinf;
1550         "Tjurina number at oo:";tauinf;
1551         "delta at oo:";deltainf;
1552         "Milnor number at oo:";2*deltainf-nbranchinf+nsing;
1553         pause();
1554      }
1555      "singularities at (x,y,1):";"";
1556   }
1557   execute("ring newR=("+charstr(R)+"),(x,y),dp;");
1558   //the singularities at the affine part
1559   map sigma=S,var(1),var(2),1;
1560   ideal I=sigma(F);
1561
1562   if(size(#)!=0)
1563   {                              //uses the normalization to compute delta
1564      list nor=normal(I,"wd");
1565      delt=nor[size(nor)];
1566      genus=genus-delt-deltainf;
1567      setring R0;
1568      return(genus);
1569   }
1570
1571   ideal I1=jacob(I);
1572   matrix Hess[2][2]=jacob(I1);
1573   ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1
1574
1575   ideal radID=std(radical(ID));//the non-nodal locus
1576   if(w>=1){"the non-nodal locus:";"";radID;pause();"";}
1577   if(deg(radID[1])==0)
1578   {
1579     ideal IDsing=1;
1580   }
1581   else
1582   {
1583     ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID
1584   }
1585
1586   iglob=vdim(std(IDsing));
1587
1588   if(iglob!=0)//computation of the radical of IDsing
1589   {
1590      ideal radIDsing=reduce(IDsing,radID);
1591      if(size(radIDsing)==0)
1592      {
1593         radIDsing=radID;
1594         attrib(radIDsing,"isSB",1);
1595      }
1596      else
1597      {
1598         radIDsing=std(radical(IDsing));
1599      }
1600      iglob=vdim(radIDsing);
1601      if((w>=1)&&(iglob))
1602          {"the non-nodal-cuspidal locus:";radIDsing;pause();"";}
1603   }
1604   cusps=vdim(radID)-iglob;
1605   nsing=nsing+cusps;
1606
1607   if(iglob==0)
1608   {
1609      if(w>=1){"             there are only cusps and nodes";"";}
1610      tau=vdim(std(I+jacob(I)));
1611      tauinf=tauinf+tau;
1612      nodes=tau-2*cusps;
1613      delt=nodes+cusps;
1614      nbranch=2*tau-3*cusps;
1615      nsing=nsing+nodes;
1616   }
1617   else
1618   {
1619       if(w>=1){"the non-nodal-cuspidal singularities";"";}
1620       setring C;
1621       ideal I1=imap(newR,IDsing);
1622       iloc=vdim(std(I1));
1623       if(iglob==iloc)
1624       {
1625          if(w>=1){"only cusps and nodes outside (0,0,1)";}
1626          setring newR;
1627          tau=vdim(std(I+jacob(I)));
1628          tauinf=tauinf+tau;
1629          inv=deltaLoc(I[1],maxideal(1));
1630          delt=inv[1];
1631          tauloc=inv[2];
1632          nodes=tau-tauloc-2*cusps;
1633          nsing=nsing+nodes;
1634          nbranch=inv[3]+ 2*nodes+cusps;
1635          delt=delt+nodes+cusps;
1636          if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";}
1637        }
1638        else
1639        {
1640           setring newR;
1641           list pr=minAssGTZ(IDsing);
1642           if(w>=1){pr;}
1643
1644           for(k=1;k<=size(pr);k++)
1645           {
1646              if(w>=1){nsing=nsing+vdim(std(pr[k]));}
1647              inv=deltaLoc(I[1],pr[k]);
1648              delt=delt+inv[1];
1649              tausing=tausing+inv[2];
1650              nbranch=nbranch+inv[3];
1651           }
1652           tau=vdim(std(I+jacob(I)));
1653           tauinf=tauinf+tau;
1654           nodes=tau-tausing-2*cusps;
1655           nsing=nsing+nodes;
1656           delt=delt+nodes+cusps;
1657           nbranch=nbranch+2*nodes+cusps;
1658        }
1659   }
1660   genus=genus-delt-deltainf;
1661   if(w>=1)
1662   {
1663      "The projected plane curve has locally:";"";
1664      "singularities:";nsing;
1665      "branches:";nbranch+nbranchinf;
1666      "nodes:"; nodes+nodesinf;
1667      "cusps:";cusps+cuspsinf;
1668      "Tjurina number:";tauinf;
1669      "Milnor number:";2*(delt+deltainf)-nbranch-nbranchinf+nsing;
1670      "delta of the projected curve:";delt+deltainf;
1671      "delta of the curve:";p_a-genus;
1672      "genus:";genus;
1673      "====================================================";
1674      "";
1675   }
1676   setring R0;
1677   return(genus);
1678}
1679example
1680{ "EXAMPLE:"; echo = 2;
1681   ring r=0,(x,y),dp;
1682   ideal i=y^9 - x^2*(x - 1)^9;
1683   genus(i);
1684}
1685
1686///////////////////////////////////////////////////////////////////////////////
1687proc deltaLoc(poly f,ideal singL)
1688"USAGE:  deltaLoc(f,J);  f poly, J ideal
1689ASSUME: f is reduced bivariate polynomial; basering has exactly two variables;
1690        J is irreducible prime component of the singular locus of f (e.g., one
1691        entry of the output of @code{minAssGTZ(I);}, I = <f,jacob(f)>).
1692RETURN:  list L:
1693@texinfo
1694@table @asis
1695@item @code{L[1]}; int:
1696         the sum of (local) delta invariants of f at the (conjugated) singular
1697         points given by J.
1698@item @code{L[2]}; int:
1699         the sum of (local) Tjurina numbers of f at the (conjugated) singular
1700         points given by J.
1701@item @code{L[3]}; int:
1702         the sum of (local) number of branches of f at the (conjugated)
1703         singular points given by J.
1704@end table
1705@end texinfo
1706NOTE:    procedure makes use of @code{execute}; increasing printlevel displays
1707         more comments (default: printlevel=0).
1708SEE ALSO: delta, tjurina
1709KEYWORDS: delta invariant; Tjurina number
1710EXAMPLE: example deltaLoc;  shows an example
1711"
1712{
1713   option(redSB);
1714   def R=basering;
1715   execute("ring S=("+charstr(R)+"),(x,y),lp;");
1716   map phi=R,x,y;
1717   ideal singL=phi(singL);
1718   singL=simplify(std(singL),1);
1719   attrib(singL,"isSB",1);
1720   int d=vdim(singL);
1721   poly f=phi(f);
1722   int i;
1723   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
1724   if(d==1)
1725   {
1726      map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2];
1727      f=alpha(f);
1728      execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;");
1729      poly f=imap(S,f);
1730      ideal singL=imap(S,singL);
1731      if((w>=1)&&(ord(f)>=2))
1732      {
1733        "local analysis of the singularities";"";
1734        basering;
1735        singL;
1736        f;
1737        pause();
1738      }
1739   }
1740   else
1741   {
1742      poly p;
1743      poly c;
1744      map psi;
1745      number co;
1746
1747      while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1))
1748      {
1749         psi=S,x,y+random(-100,100)*x;
1750         singL=psi(singL);
1751         singL=std(singL);
1752          f=psi(f);
1753      }
1754
1755      if(deg(lead(singL[2]))==1)
1756      {
1757         p=singL[1];
1758         c=singL[2]-lead(singL[2]);
1759         co=leadcoef(singL[2]);
1760      }
1761      if(deg(lead(singL[1]))==1)
1762      {
1763         psi=S,y,x;
1764         f=psi(f);
1765         singL=psi(singL);
1766         p=singL[2];
1767         c=singL[1]-lead(singL[1]);
1768         co=leadcoef(singL[1]);
1769      }
1770
1771      execute("ring B=("+charstr(S)+"),a,dp;");
1772      map beta=S,a,a;
1773      poly p=beta(p);
1774
1775      execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;");
1776      number p=number(imap(B,p));
1777
1778      minpoly=p;
1779      map iota=S,a,a;
1780      number c=number(iota(c));
1781      number co=iota(co);
1782
1783      map alpha=S,x-c/co,y+a;
1784      poly f=alpha(f);
1785      f=cleardenom(f);
1786      if((w>=1)&&(ord(f)>=2))
1787      {
1788        "local analysis of the singularities";"";
1789        basering;
1790        alpha;
1791        f;
1792        pause();
1793        "";
1794      }
1795   }
1796   option(noredSB);
1797   ideal fstd=std(ideal(f)+jacob(f));
1798   poly hc=highcorner(fstd);
1799   int tau=vdim(fstd);
1800   int o=ord(f);
1801   int delt,nb;
1802
1803   if(tau==0)                 //smooth case
1804   {
1805      setring R;
1806      return(list(0,0,1));
1807   }
1808   if((char(basering)>=181)||(char(basering)==0))
1809   {
1810      if(o==2)                //A_k-singularity
1811      {
1812        if(w>=1){"A_k-singularity";"";}
1813         setring R;
1814         delt=(tau+1)/2;
1815         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
1816      }
1817      if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2)))
1818      {
1819        if(w>=1){"D_k- singularity";"";}
1820
1821         setring R;
1822         delt=(tau+2)/2;
1823         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
1824      }
1825
1826      int mu=vdim(std(jacob(f)));
1827      poly g=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon
1828
1829      list NP=newtonpoly(g);
1830      if(w>=1){"Newton-Polygon:";NP;"";}
1831      int s=size(NP);
1832
1833      if(is_NND(f,mu,NP))
1834      { // the Newton-polygon is non-degenerate
1835        // compute nb, the number of branches
1836        for(i=1;i<=s-1;i++)
1837        {
1838          nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]);
1839        }
1840        if(w>=1){"Newton-Polygon is non-degenerated";"";}
1841        return(list(d*(mu+nb-1)/2,d*tau,d*nb));
1842      }
1843
1844      if(w>=1){"Newton-Polygon is degenerated";"";}
1845
1846      // the following can certainly be made more efficient when replacing
1847      // 'reddevelop' (used only for computing number of branches) by
1848      // successive blowing-up + test if Newton polygon degenerate:
1849      if(s>2)    //  splitting of f
1850      {
1851         if(w>=1){"Newton polygon can be used for splitting";"";}
1852         intvec v=NP[1][2]-NP[2][2],NP[2][1];
1853         int de=w_deg(g,v);
1854         int st=w_deg(hc,v)+v[1]+v[2];
1855         poly f1=var(2)^NP[2][2];
1856         poly f2=jet(g,de,v)/var(2)^NP[2][2];
1857         poly h=g-f1*f2;
1858         de=w_deg(h,v);
1859         poly k;
1860         ideal wi=var(2)^NP[2][2],f2;
1861         matrix li;
1862         while(de<st)
1863         {
1864           k=jet(h,de,v);
1865           li=lift(wi,k);
1866           f1=f1+li[2,1];
1867           f2=f2+li[1,1];
1868           h=g-f1*f2;
1869           de=w_deg(h,v);
1870         }
1871         nb=deltaLoc(f1,maxideal(1))[3]+deltaLoc(f2,maxideal(1))[3];
1872         setring R;
1873         return(list(d*(mu+nb-1)/2,d*tau,d*nb));
1874      }
1875
1876      f=jet(f,deg(hc)+2);
1877      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
1878      ideal fac=factorize(f,1);
1879      if(size(fac)>1)
1880      {
1881         nb=0;
1882         for(i=1;i<=size(fac);i++)
1883         {
1884            nb=nb+deltaLoc(fac[i],maxideal(1))[3];
1885         }
1886         setring R;
1887         return(list(d*(mu+nb-1)/2,d*tau,d*nb));
1888      }
1889      list hne=reddevelop(f);
1890      nb=size(hne);
1891      setring R;
1892      kill HNEring;
1893      return(list(d*(mu+nb-1)/2,d*tau,d*nb));
1894   }
1895   else             //the case of small characteristic
1896   {
1897      f=jet(f,deg(hc)+2);
1898      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
1899      delt=delta(f);
1900      return(list(d*delt,d*tau,d));
1901   }
1902}
1903example
1904{ "EXAMPLE:"; echo = 2;
1905  ring r=0,(x,y),dp;
1906  poly f=(x2+y^2-1)^3 +27x2y2;
1907  ideal I=f,jacob(f);
1908  I=std(I);
1909  list qr=minAssGTZ(I);
1910  size(qr);
1911  // each component of the singular locus either describes a cusp or a pair
1912  // of conjugated nodes:
1913  deltaLoc(f,qr[1]);
1914  deltaLoc(f,qr[2]);
1915  deltaLoc(f,qr[3]);
1916  deltaLoc(f,qr[4]);
1917  deltaLoc(f,qr[5]);
1918  deltaLoc(f,qr[6]);
1919}
1920///////////////////////////////////////////////////////////////////////////////
1921// compute the weighted degree of p
1922static proc w_deg(poly p, intvec v)
1923{
1924   if(p==0){return(-1);}
1925   int d=0;
1926   while(jet(p,d,v)==0){d++;}
1927   d=(transpose(leadexp(jet(p,d,v)))*v)[1];
1928   return(d);
1929}
1930
1931//proc hilbPoly(ideal J)
1932//{
1933//   poly hp;
1934//   int i;
1935//   if(!attrib(J,"isSB")){J=std(J);}
1936//   intvec v = hilb(J,2);
1937//   for(i=1; i<=size(v); i++){ hp=hp+v[i]*(var(1)-i+2);}
1938//   return(hp);
1939//}
1940
1941
1942//////////////////////////////////////////////////////////////////////////////
1943
1944proc primeClosure (list L, list #)
1945"USAGE:    primeClosure(L [,c]); L a list of a ring containing a prime ideal
1946                                 ker, c an optional integer
1947RETURN:   a list L consisting of rings L[1],...,L[n] such that
1948          - L[1] is a copy of (not a reference to!) the input ring L[1]
1949          - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi
1950            such that
1951                    L[1]/ker --> ... --> L[n]/ker
1952            are injections given by the corresponding ideals phi, and L[n]/ker
1953            is the integral closure of L[1]/ker in its quotient field.
1954          - all rings L[i] contain a polynomial nzd such that elements of
1955            L[i]/ker are quotients of elements of L[i-1]/ker with denominator
1956            nzd via the injection phi.
1957NOTE:     - L is constructed by recursive calls of primeClosure itself.
1958          - c determines the choice of nzd:
1959               - c not given or equal to 0: first generator of the ideal SL,
1960                 the singular locus of Spec(L[i]/ker)
1961               - c<>0: the generator of SL with least number of monomials.
1962EXAMPLE:  example primeClosure; shows an example
1963"
1964{
1965// Start with a consistency check:
1966
1967  if (!(typeof(L[1])=="ring"))
1968    {
1969      "// Parameter must be a ring or a list containing a ring!";
1970      return(-1);
1971    }
1972
1973  int dblvl=printlevel-voice+2;
1974
1975
1976// Some auxiliary variables:
1977
1978  int n=size(L);
1979
1980// How to choose the non-zerodivisor later on?
1981
1982  int nzdoption=0;
1983  if (size(#)>0)
1984    {
1985      nzdoption=#[1];
1986    }
1987
1988// R0 is the ring to work with, if we are in step one, make a copy of the
1989// input ring, so that all objects are created in the copy, not in the original
1990// ring (therefore a copy, not a reference is defined).
1991
1992  if (n==1)
1993    {
1994      def R=L[1];
1995      def BAS=basering;
1996      setring R;
1997      if (!(typeof(ker)=="ideal"))
1998        {
1999          "// No ideal ker in the input ring!";
2000          return (-1);
2001        }
2002      ker=simplify(interred(ker),15);
2003      execute ("ring R0="+charstr(R)+",("+varstr(R)+"),("+ordstr(R)+");");
2004      ideal ker=fetch(R,ker);
2005
2006      // check whether we compute the normalization of the blow up of
2007      // an isolated singularity at the origin (checked in normalI)
2008
2009      if (typeof(attrib(L[1],"iso_sing_Rees"))=="int")
2010      {
2011        attrib(R0,"iso_sing_Rees",attrib(L[1],"iso_sing_Rees"));
2012      }
2013      L=R0;
2014    }
2015  else
2016    {
2017      def R0=L[n];
2018      setring R0;
2019    }
2020
2021
2022// In order to apply HomJJ from normal.lib, we need the radical of the singular
2023// locus of ker, J:=rad(ker):
2024
2025//  if (homog(ker)==1)
2026//    {
2027      list SM=mstd(ker);
2028//    }
2029/*  else
2030    {
2031      list SM=groebner(ker),ker;
2032    }*/
2033
2034// In the first iteration, we have to compute the singular locus "from
2035// scratch". In further iterations, we can fetch it from the previous one but
2036// have to compute its radical.
2037
2038  if (n==1)
2039    {
2040      if (typeof(attrib(R0,"iso_sing_Rees"))=="int")
2041      {
2042        ideal J;
2043        for (int s=1;s<=attrib(R0,"iso_sing_Rees");s++)
2044        {
2045          J=J,var(s);
2046        }
2047      }
2048      else
2049      {
2050        ideal J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
2051      }
2052      J=J+SM[2];
2053      if (homog(J)==1)
2054      {
2055        J=mstd(J)[2];
2056      }
2057      J=radical(interred(J));
2058    }
2059  else
2060    {
2061      J=radical(interred(J));
2062    }
2063
2064  // having computed the radical J of the ideal of the singular locus,
2065  // we now need to pick an element nzd of J
2066
2067  poly nzd=J[1];
2068  if (nzdoption)
2069  {
2070    for (int ii=2;ii<=ncols(J);ii++)
2071    {
2072      if ( (deg(nzd)>=deg(J[ii])) && (size(nzd)>size(J[ii])) )
2073      {
2074        nzd=J[ii];
2075      }
2076    }
2077  }
2078
2079  export nzd;
2080  list RR=SM[1],SM[2],J,nzd;
2081  list RS=HomJJ(RR);
2082
2083
2084//////////////////////////////////////////////////////////////////////////////
2085
2086
2087// If we've reached the integral closure (as determined by the result of
2088// HomJJ), then we are done, otherwise we have to prepare the next iteration.
2089
2090  if (RS[2]==1)               // we've reached the integral closure
2091    {
2092      kill(J);
2093      return(L);
2094    }
2095  else                        // prepare the next iteration
2096    {
2097      if (n==1)               // In the first iteration: keep only the data
2098        {                     // needed later on.
2099          kill (RR,SM);
2100
2101          export(ker);
2102        }
2103
2104      def R1=RS[1];           // The data of the next ring R1:
2105      setring R1;             // keep only what is necessary and kill
2106      ideal ker=endid;        // everything else.
2107      export(ker);
2108      ideal norid=endid;
2109      export(norid);
2110      kill(endid);
2111
2112      map phi=R0,endphi;      // fetch the singular locus
2113      ideal J=mstd(simplify(phi(J)+ker,4))[2];
2114      export(J);
2115      if(n>1)
2116      {
2117         ideal normap=phi(normap);
2118      }
2119      else
2120      {
2121         ideal normap=endphi;
2122      }
2123      export(normap);
2124      kill(phi);              // we save phi as ideal, not as map, so that
2125      ideal phi=endphi;       // we have more flexibility in the ring names
2126      kill(endphi);           // later on.
2127      export(phi);
2128      L=insert(L,R1,n);       // Add the new ring R1 and go on with the
2129      if (size(#)>0)
2130        {
2131          return (primeClosure(L,#));
2132        }
2133      else
2134        {
2135          return(primeClosure(L));         // next iteration.
2136        }
2137    }
2138}
2139example
2140{
2141  "EXAMPLE:"; echo=2;
2142  ring R=0,(x,y),dp;
2143  ideal I=x4,y4;
2144  def K=ReesAlgebra(I)[1];        // K contains ker such that K/ker=R[It]
2145  list L=primeClosure(K);
2146  def R(1)=L[1];                  // L[4] contains ker, L[4]/ker is the
2147  def R(4)=L[4];                  // integral closure of L[1]/ker
2148  setring R(1);
2149  R(1);
2150  ker;
2151  setring R(4);
2152  R(4);
2153  ker;
2154}
2155
2156
2157//////////////////////////////////////////////////////////////////////////////
2158//////////////////////////////////////////////////////////////////////////////
2159
2160proc closureRingtower(list L)
2161"USAGE:    closureRingtower(list L); L a list of rings
2162CREATE:   rings R(1),...,R(n) such that R(i)=L[i] for all i
2163EXAMPLE:  example closureRingtower; shows an example
2164"
2165{
2166  int n=size(L);
2167
2168  for (int i=1;i<=n;i++)
2169    {
2170      if (defined(R(i))) {
2171        string s="Fixed name R("+string(i)+") leads to conflict with existing "
2172              +"object having this name";
2173        ERROR(s);
2174      }
2175      def R(i)=L[i];
2176      export R(i);
2177    }
2178
2179  return();
2180}
2181example
2182{
2183  "EXAMPLE:"; echo=2;
2184  ring R=0,(x,y),dp;
2185  ideal I=x4,y4;
2186  list L=primeClosure(ReesAlgebra(I)[1]);
2187  closureRingtower(L);
2188  R(1);
2189  R(4);
2190  kill R(1),R(2),R(3),R(4);
2191}
2192
2193//////////////////////////////////////////////////////////////////////////////
2194//////////////////////////////////////////////////////////////////////////////
2195
2196proc closureFrac(list L)
2197"USAGE:    closureFrac (L); L a list as in the result of primeClosure, L[n]
2198                           containing an additional poly f
2199CREATE:   a list fraction of two elements of L[1], such that
2200          f=fraction[1]/fraction[2] via the injections phi L[i]-->L[i+1].
2201EXAMPLE:  example closureFrac; shows an example
2202"
2203{
2204// Define some auxiliary variables:
2205
2206  int n=size(L);
2207  int j,k,l;
2208  string mapstr;
2209
2210  for (int i=1;i<=n;i++) { def R(i)=L[i]; }
2211
2212// The quotient representing f is computed as in 'closureGenerators' with
2213// the differences that
2214//   - the loop is done twice: for the numerator and for the denominator;
2215//   - the result is stored in the list fraction and
2216//   - we have to make sure that no more objects of the rings R(i) survive.
2217
2218  for (j=1;j<=2;j++)
2219    {
2220      setring R(n);
2221      if (j==1)
2222        {
2223          poly p=f;
2224        }
2225      else
2226        {
2227          p=1;
2228        }
2229
2230      for (k=n;k>1;k--)
2231        {
2232
2233          if (j==1)
2234            {
2235              map phimap=R(k-1),phi;
2236            }
2237
2238          p=p*phimap(nzd);
2239
2240          if (j==2)
2241            {
2242              kill(phimap);
2243            }
2244
2245          if (j==1)
2246            {
2247              execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
2248                       Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
2249                       +"),dp("+string(ncols(phi))+"));");
2250              ideal phi = imap(R(k),phi);
2251              ideal J= imap (R(k),ker);
2252              for (l=1;l<=ncols(phi);l++)
2253                {
2254                  J=J+(Z(l)-phi[l]);
2255                }
2256              J=groebner(J);
2257              poly h=NF(imap(R(k),p),J);
2258            }
2259          else
2260            {
2261              setring S(k);
2262              h=NF(imap(R(k),p),J);
2263              setring R(k);
2264              kill(p);
2265            }
2266
2267          setring R(k-1);
2268
2269          if (j==1)
2270            {
2271              mapstr=" map backmap = S(k),";
2272              for (l=1;l<=nvars(R(k));l++)
2273                {
2274                  mapstr=mapstr+"0,";
2275                }
2276              execute (mapstr+"maxideal(1);");
2277
2278              poly p;
2279            }
2280          p=NF(backmap(h),std(ker));
2281          if (j==2)
2282            {
2283              kill(backmap);
2284            }
2285        }
2286
2287      if (j==1)
2288        {
2289          if (defined(fraction))
2290            {
2291              kill(fraction);
2292              list fraction=p;
2293            }
2294          else
2295            {
2296              list fraction=p;
2297            }
2298        }
2299      else
2300        {
2301          fraction=insert(fraction,p,1);
2302        }
2303    }
2304  export(fraction);
2305  return ();
2306}
2307example
2308{
2309  "EXAMPLE:"; echo=2;
2310  ring R=0,(x,y),dp;
2311  ideal ker=x2+y2;
2312  export ker;
2313  list L=primeClosure(R);          // We normalize R/ker
2314  for (int i=1;i<=size(L);i++) { def R(i)=L[i]; }
2315  setring R(2);
2316  kill(R);
2317  phi;                             // The map R(1)-->R(2)
2318  poly f=T(2);                     // We will get a representation of f
2319  export f;
2320  L[2]=R(2);
2321  closureFrac(L);
2322  setring R(1);
2323  kill R(2);
2324  fraction;                        // f=fraction[1]/fraction[2] via phi
2325  kill R(1);
2326}
2327
2328//////////////////////////////////////////////////////////////////////////////
2329//////////////////////////////////////////////////////////////////////////////
2330
2331static
2332proc closureGenerators(list L);    // called inside normalI
2333{
2334  // In the list L of rings R(1),...,R(n), compute representations of the ring
2335  // ring variables of the last ring R(n) as fractions of elements of R(1).
2336
2337  def Rees=basering;              // when called inside normalI, the Rees
2338                                  // Algebra is the current basering.
2339
2340  // First of all we need some variable declarations...
2341
2342  list preimages;
2343  int n=size(L);                  // the number of rings R(1)-->...-->R(n)
2344  int j,k,l;
2345  int length=nvars(L[n]);         // the number of variables of the last ring
2346  string mapstr;
2347
2348  for (int i=1;i<=n;i++) { def R(i)=L[i]; }
2349
2350  // For each variable (counter j) and for each intermediate ring (counter k):
2351  // Find a preimage of var_j*phi(nzd(k-1)) in R(k-1).
2352  // Finally, do the same for nzd.
2353
2354  for (j=1;j<=length+1;j++)
2355    {
2356      setring R(n);
2357
2358      if (j==1)
2359        {
2360          poly p;
2361        }
2362
2363      if (j<=nvars(R(n)))
2364        {
2365          p=var(j);
2366        }
2367      else
2368        {
2369          p=1;
2370        }
2371
2372      for (k=n;k>1;k--)
2373        {
2374
2375          if (j==1)
2376          {
2377            map phimap=R(k-1),phi;      // phimap is the map corresponding
2378          }                             // to phi
2379
2380          p=p*phimap(nzd);
2381
2382          // Compute the preimage of [p mod ker(k)] under phi in R(k-1):
2383          // As p is an element of Image(phi), there is a polynomial h such
2384          // that h is mapped to [p mod ker(k)], and h can be computed as the
2385          // normal form of p w.r.t. a Gröbner basis of
2386          // J(k):=<ker(k),Z(l)-phi(k)(l)> in R(k)[Z]=:S(k)
2387
2388          if (j==1)   // In the first iteration: Create S(k), fetch phi and
2389                      // ker(k) and construct the ideal J(k).
2390            {
2391              execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
2392                       Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
2393                       +"),dp("+string(ncols(phi))+"));");
2394              ideal phi=imap(R(k),phi);
2395              ideal J=imap (R(k),ker);
2396              for (l=1;l<=ncols(phi);l++)
2397                {
2398                  J=J+(Z(l)-phi[l]);
2399                }
2400              J=groebner(J);
2401              poly h=NF(imap(R(k),p),J);
2402            }
2403          else
2404            {
2405              setring S(k);
2406              h=NF(imap(R(k),p),J);
2407            }
2408
2409          setring R(k-1);
2410
2411          if (j==1)  // In the first iteration: Compute backmap:S(k)-->R(k-1)
2412            {
2413              mapstr=" map backmap = S(k),";
2414              for (l=1;l<=nvars(R(k));l++)
2415                {
2416                  mapstr=mapstr+"0,";
2417                }
2418              execute (mapstr+"maxideal(1);");
2419
2420              poly p;
2421            }
2422          p=NF(backmap(h),std(ker));
2423        }
2424
2425      // When down to R(1), store the result in the list preimages
2426
2427      preimages = insert(preimages,p,j-1);
2428    }
2429
2430  // R(1) was a copy of Rees, so we have to get back to the basering Rees from
2431  // the beginning and fetch the result (the list preimages) to this ring.
2432
2433  setring Rees;
2434  return (fetch(R(1),preimages));
2435}
2436///////////////////////////////////////////////////////////////////////////
2437
2438/*
2439                           Examples:
2440LIB"normal.lib";
2441//Huneke
2442ring qr=31991,(a,b,c,d,e),dp;
2443ideal i=
24445abcde-a5-b5-c5-d5-e5,
2445ab3c+bc3d+a3be+cd3e+ade3,
2446a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
2447abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
2448ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
2449a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
2450a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
2451b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
2452
2453
2454//Vasconcelos (dauert laenger: 70 sec)
2455ring r=32003,(x,y,z,w,t),dp;
2456ideal i=
2457x2+zw,
2458y3+xwt,
2459xw3+z3t+ywt2,
2460y2w4-xy2z2t-w3t3;
2461
2462//Theo1
2463ring r=32003,(x,y,z),wp(2,3,6);
2464ideal i=zy2-zx3-x6;
2465
2466//Theo1a (CohenMacaulay and regular in codimension 2)
2467ring r=32003,(x,y,z,u),wp(2,3,6,6);
2468ideal i=zy2-zx3-x6+u2;
2469
2470
2471//Theo2
2472ring r=32003,(x,y,z),wp(3,4,12);
2473ideal i=z*(y3-x4)+x8;
2474
2475//Theo2a
2476ring r=32003,(T(1..4)),wp(3,4,12,17);
2477ideal i=
2478T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
2479T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
2480T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
2481T(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;
2482
2483//Theo3
2484ring r=32003,(x,y,z),wp(3,5,15);
2485ideal i=z*(y3-x5)+x10;
2486
2487
2488//Theo4
2489ring r=32003,(x,y,z),dp;
2490ideal i=(x-y)*(x-z)*(y-z);
2491
2492//Theo5
2493ring r=32003,(x,y,z),wp(2,1,2);
2494ideal i=z3-xy4;
2495
2496//Theo6
2497ring r=32003,(x,y,z),dp;
2498ideal i=x2y2+x2z2+y2z2;
2499
2500ring r=32003,(a,b,c,d,e,f),dp;
2501ideal i=
2502bf,
2503af,
2504bd,
2505ad;
2506
2507//ist CM
2508//Sturmfels
2509ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
2510ideal i=
2511bv+su,
2512bw+tu,
2513sw+tv,
2514by+sx,
2515bz+tx,
2516sz+ty,
2517uy+vx,
2518uz+wx,
2519vz+wy,
2520bvz;
2521
2522//J S/Y
2523ring r=32003,(x,y,z,t),dp;
2524ideal i=
2525x2z+xzt,
2526xyz,
2527xy2-xyt,
2528x2y+xyt;
2529
2530//St_S/Y
2531ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
2532ideal i=
2533wy-vz,
2534vx-uy,
2535tv-sw,
2536su-bv,
2537tuy-bvz;
2538
2539//Horrocks:
2540ring r=32003,(a,b,c,d,e,f),dp;
2541ideal i=
2542adef-16000be2f+16001cef2,
2543ad2f+8002bdef+8001cdf2,
2544abdf-16000b2ef+16001bcf2,
2545a2df+8002abef+8001acf2,
2546ad2e-8000bde2-7999cdef,
2547acde-16000bce2+16001c2ef,
2548a2de-8000abe2-7999acef,
2549acd2+8002bcde+8001c2df,
2550abd2-8000b2de-7999bcdf,
2551a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
2552abde-8000b2e2-acdf-16001bcef-8001c2f2,
2553abcd-16000b2ce+16001bc2f,
2554a2cd+8002abce+8001ac2f,
2555a2bd-8000ab2e-7999abcf,
2556ab3f-3bdf3,
2557a2b2f-2adf3-16000bef3+16001cf4,
2558a3bf+4aef3,
2559ac3e-10668cde3,
2560a2c2e+10667ade3+16001be4+5334ce3f,
2561a3ce+10669ae3f,
2562bc3d+8001cd3e,
2563ac3d+8000bc3e+16001cd2e2+8001c4f,
2564b2c2d+16001ad4+4000bd3e+12001cd3f,
2565b2c2e-10668bc3f-10667cd2ef,
2566abc2e-cde2f,
2567b3cd-8000bd3f,
2568b3ce-10668b2c2f-10667bd2ef,
2569abc2f-cdef2,
2570a2bce-16000be3f+16001ce2f2,
2571ab3d-8000b4e-8001b3cf+16000bd2f2,
2572ab2cf-bdef2,
2573a2bcf-16000be2f2+16001cef3,
2574a4d-8000a3be+8001a3cf-2ae2f2;
2575
2576
2577ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
2578ideal k=
2579wy-vz,
2580vx-uy,
2581tv-sw,
2582su-bv,
2583tuy-bvz;
2584ideal j=x2y2+x2z2+y2z2;
2585ideal i=mstd(intersect(j,k))[2];
2586
2587
2588ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
2589ideal i=
2590wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
2591vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
2592tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
2593sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
2594tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;
2595
2596
2597//riemenschneider
2598ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
2599ideal i=
2600xz,
2601vx,
2602ux,
2603su,
2604qu,
2605txy,
2606stx,
2607qtx,
2608uv2z-uwz,
2609uv3-uvw,
2610puv2-puw;
2611
2612ring r=0,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
2613ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;
2614
2615//Yoshihiko Sakai
2616ring r=0,(x,y),dp;   //genus 0 4 nodes and 6 cusps
2617ideal i=(x2+y^2-1)^3 +27x2y2;
2618
2619ring r=0,(x,y),dp;   //genus 0
2620ideal i=(x-y^2)^2 - y*x^3;
2621
2622ring r=0,(x,y),dp;  //genus 4
2623ideal i=y3-x6+1;
2624
2625int m=9;           // q=9: genus 0
2626int p=2;
2627int q=9;//2,...,9
2628ring r=0,(x,y),dp;
2629ideal i=y^m - x^p*(x - 1)^q;
2630
2631ring r=0,(x,y),dp;  //genus 19
2632ideal 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;
2633
2634
2635ring r=0,(x,y),dp;  //genus 34
2636ideal i=y10+(-2494x2+474)*y8+(84366+2042158x4-660492)*y6
2637        +(128361096x4-47970216x2+6697080-761328152x6)*y4
2638        +(-12024807786x4-506101284x2+15052058268x6+202172841-3212x8)*y2
2639        +34263110700x4-228715574724x6+5431439286x2+201803238
2640        -9127158539954x10-3212722859346x8;
2641
2642//Rob Koelman
2643ring r=0,(x,y,z),dp;//genus 10  with 26 cusps
2644ideal i=
2645761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+
2646 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2-
2647 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4-
2648 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+
2649 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-
2650 z^10-474*z^8*y^2-84366*z^6*y^4;
2651
2652ring r=0,(x,y),dp;//genus 10  with 26 cusps
2653ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
2654-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
2655+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
2656-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;
2657
2658
2659ring r=0,(x,y),dp;   // genus 1  with 5 cusps
2660ideal i=57y5+516x4y-320x4+66y4-340x2y3+73y3+128x2-84x2y2-96x2y;
2661
2662//Mark van Hoeij
2663ring r=0,(x,y),dp;  //genus 19
2664ideal i=y20+y13x+x4y5+x3*(x+1)^2;
2665
2666ring r=0,(x,y),dp;  //genus 35
2667ideal i=y30+y13x+x4y5+x3*(x+1)^2;
2668
2669ring r=0,(x,y),dp;   //genus 55
2670ideal i=y40+y13x+x4y5+x3*(x+1)^2;
2671
2672ring r=0,(x,y),dp;  //genus 4
2673ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
2674
2675ring r=0,(y,z,w,u),dp; //genus -5
2676ideal i=y2+z2+w2+u2,w4-u4;
2677
2678ring r=0,(x,y,t),dp; //genus -5
2679ideal i=x8+8x7y+32x6y2+80x5y3+136x4y4+160x3y5+128x2y6+64xy7+16y8+4x6t2+24x5yt2+72x4y2t2+128x3y3t2+144x2y4t2+96xy5t2+32y6t2+14x4t4+56x3yt4+112x2y2t4+112xy3t4+40y4t4+20x2t6+40xyt6+8y2t6+9t8;
2680
2681ring r=0,(y,z,w,u),dp; //genus 9
2682ideal i=y2+z2+w2+u2,z4+w4+u4;
2683
2684ring r=0,(x,y,t),dp;
2685ideal i=
268625x8+200x7y+720x6y2+1520x5y3+2064x4y4+1856x3y5+1088x2y6+384xy7+64y8-12x6t2-72x5yt2-184x4y2t2-256x3y3t2-192x2y4t2-64xy5t2-2x4t4-8x3yt4+16xy3t4+16y4t4+4x2t6+8xyt6+8y2t6+t8;
2687
2688ring r=0,(x,y,t),dp;
2689ideal i=
269032761x8+786264x7y+8314416x6y2+50590224x5y3+193727376x4y4+478146240x3y5+742996800x2y6+664848000xy7+262440000y8+524176x7t+11007696x6yt+99772992x5y2t+505902240x4y3t+1549819008x3y4t+2868877440x2y5t+2971987200xy6t+1329696000y7t+3674308x6t2+66137544x5yt2+499561128x4y2t2+2026480896x3y3t2+4656222144x2y4t2+5746386240xy5t2+2976652800y6t2+14737840x5t3+221067600x4yt3+1335875904x3y2t3+4064449536x2y3t3+6226336512xy4t3+3842432640y5t3+36997422x4t4+443969064x3yt4+2012198112x2y2t4+4081745520xy3t4+3126751632y4t4+59524208x3t5+535717872x2yt5+1618766208xy2t5+1641991392y3t5+59938996x2t6+359633976xyt6+543382632y2t6+34539344xt7+103618032yt7+8720497t8;
2691
2692ring r=32003,(x,y,z,w,u),dp;
2693ideal i=x2+y2+z2+w2+u2,x3+y3+z3,z4+w4+u4;
2694
2695*/
2696
Note: See TracBrowser for help on using the repository browser.