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

spielwiese
Last change on this file since 0ae4ce was 0ae4ce, checked in by Anne Frühbis-Krüger <anne@…>, 23 years ago
• Property mode set to `100644`
File size: 35.0 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2// normal.lib
3// algorithms for computing the normalization based on
4// the criterion of Grauert/Remmert and ideas of De Jong & Vasconcelos
5///////////////////////////////////////////////////////////////////////////////
6
7version="\$Id: normal.lib,v 1.27 2000-12-19 14:41:43 anne Exp \$";
8category="Commutative Algebra";
9info="
10LIBRARY:  normal.lib     PROCEDURES FOR NORMALIZATION
11AUTHORS:  Gert-Martin Greuel, email: greuel@mathematik.uni-kl.de,
12          Gerhard Pfister, email: pfister@mathematik.uni-kl.de
13
14PROCEDURES:
15 normal(I);             computes the normalization of basering/I
16 HomJJ(L);              presentation of End_R(J) as affine ring, L a list
17";
18
19LIB "general.lib";
20LIB "sing.lib";
21LIB "primdec.lib";
22LIB "elim.lib";
23LIB "presolve.lib";
24LIB "inout.lib";
25///////////////////////////////////////////////////////////////////////////////
26static
27proc isR_HomJR (list Li)
28"USAGE:   isR_HomJR (Li);  Li = list: ideal SBid, ideal J, poly p
29COMPUTE: module Hom_R(J,R) = R:J and compare with R
30ASSUME:  R    = P/SBid,  P = basering
31         SBid = standard basis of an ideal in P,
32         J    = ideal in P containing the polynomial p,
33         p    = nonzero divisor of R
34RETURN:  1 if R = R:J, 0 if not
35EXAMPLE: example isR_HomJR;  shows an example
36"
37{
38   int n, ii;
39 def P = basering;
40   ideal SBid = Li[1];
41   ideal J = Li[2];
42   poly p = Li[3];
43   attrib(SBid,"isSB",1);
44   attrib(p,"isSB",1);
45 qring R    = SBid;
46   ideal J  = fetch(P,J);
47   poly p   = fetch(P,p);
48   ideal f  = quotient(p,J);
49   ideal lp = std(p);
50   n=1;
51   for (ii=1; ii<=size(f); ii++ )
52   {
53      if ( reduce(f[ii],lp) != 0)
54      { n = 0; break; }
55   }
56   return (n);
57 //?spaeter hier einen Test ob Hom(I,R) = Hom(I,I)?
58}
59example
60{"EXAMPLE:";  echo = 2;
61  ring r   = 0,(x,y,z),dp;
62  ideal id = y7-x5+z2;
63  ideal J  = x3,y+z;
64  poly p   = xy;
65  list Li  = std(id),J,p;
66  isR_HomJR (Li);
67
68  ring s   = 0,(t,x,y),dp;
69  ideal id = x2-y2*(y-t);
70  ideal J  = jacob(id);
71  poly p   = J[1];
72  list Li  = std(id),J,p;
73  isR_HomJR (Li);
74}
75///////////////////////////////////////////////////////////////////////////////
76
77proc HomJJ (list Li)
78"USAGE:   HomJJ (Li);  Li = list: ideal SBid, ideal id, ideal J, poly p
79ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
80         SBid = standard basis of id,
81         J    = ideal of P containing the polynomial p,
82         p    = nonzero divisor of R
83COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure where
84         R is the quotient ring of P modulo the standard basis SBid
85RETURN:  a list of two objects
86         _[1]: a polynomial ring, containing two ideals, 'endid' and 'endphi'
87               s.t. _[1]/endid = Hom_R(J,J) and
88               endphi describes the canonical map R -> Hom_R(J,J)
89         _[2]: an integer which is 1 if phi is an isomorphism, 0 if not
90NOTE:    printlevel >=1: display comments (default: printlevel=0)
91EXAMPLE: example HomJJ;  shows an example
92"
93{
94//---------- initialisation ---------------------------------------------------
95
96   int isIso,isPr,isCo,isRe,isEq,ii,jj,q,y;
97   intvec rw,rw1;
98   list L;
99   y = printlevel-voice+2;  // y=printlevel (default: y=0)
100 def P = basering;
101   ideal SBid, id, J = Li[1], Li[2], Li[3];
102   poly p = Li[4];
103   attrib(SBid,"isSB",1);
104   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not
105
106//---- set attributes for special cases where algorithm can be simplified -----
107   if( homo==1 )
108   {
109      rw = ringweights(P);
110   }
111   if( typeof(attrib(id,"isPrim"))=="int" )
112   {
113      if(attrib(id,"isPrim")==1)  { isPr=1; }
114   }
115   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
116   {
117      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
118   }
119   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
120   {
121      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
122   }
123   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
124   {
125      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
126   }
127   if( typeof(attrib(id,"isEquidimensional"))=="int" )
128   {
129      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
130   }
131//-------------------------- go to quotient ring ------------------------------
132 qring R  = SBid;
133   ideal id = fetch(P,id);
134   ideal J  = fetch(P,J);
135   poly p   = fetch(P,p);
136   ideal f,rf,f2;
137   module syzf;
138
139//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
140   if ( y>=1 )
141   {
142     "// compute p*Hom(J,J) = p*J:J, p a non-zerodivisor";
143     "//   p is equal to:"; "";
144     p;
145     "";
146   }
147   f  = quotient(p*J,J);
148   if ( y>=1 )
149   { "// the module p*Hom(J,J) = p*J:J, p a non-zerodivisor";
150      "// p"; p;
151      "// f=p*J:J";f;
152   }
153   f2 = std(p);
154
155   if(isIso==0)
156   {
157     ideal f1=std(f);
158     attrib(f1,"isSB",1);
159    // if( codim(f1,f2) >= 0 )
160    // {
161    //  dbprint(printlevel-voice+3,"// dimension of non-normal locus is zero");
162    //    isIso=1;
163    // }
164  }
165//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
166
167   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R
168   if ( size(rf) == 0 )
169   {
170      if ( homog(f) && find(ordstr(basering),"s")==0 )
171      {
172         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
173      }
174      else
175      {
176         ring newR1 = char(P),(X(1..nvars(P))),dp;
177      }
178      ideal endphi = maxideal(1);
179      ideal endid = fetch(P,id);
180      L=substpart(endid,endphi,homo,rw);
181      def lastRing=L[1];
182      setring lastRing;
183
184      attrib(endid,"isCohenMacaulay",isCo);
185      attrib(endid,"isPrim",isPr);
186      attrib(endid,"isIsolatedSingularity",isIso);
187      attrib(endid,"isRegInCodim2",isRe);
188      attrib(endid,"isEqudimensional",isEq);
189      attrib(endid,"isCompleteIntersection",0);
191//      export endid;
192//      export endphi;
193//      L = newR1;
194      L=lastRing;
195      L = insert(L,1,1);
196      dbprint(y,"// case R = Hom(J,J)");
197      if(y>=1)
198      {
199         "//   R=Hom(J,J)";
200         "   ";
201         lastRing;
202         "   ";
203         "//   the new ideal";
204         endid;
205         "   ";
206         "//   the old ring";
207         "   ";
208         P;
209         "   ";
210         "//   the old ideal";
211         "   ";
212         setring P;
213         id;
214         "   ";
215         setring lastRing;
216         "//   the map";
217         "   ";
218         endphi;
219         "   ";
220         pause();
221         newline;
222      }
223      setring P;
224      return(L);
225   }
226   if(y>=1)
227   {
228      "// R is not equal to Hom(J,J), we have to try again";
229      pause();
230      newline;
231   }
232//---------- Hom(J,J) != R: create new ring and map form old ring -------------
233// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
234
235   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
236   q = size(f);
237   syzf = syz(f);
238
239   if ( homo==1 )
240   {
241      rw1 = rw,0;
242      for ( ii=2; ii<=q; ii++ )
243      {
244         rw  = rw, deg(f[ii])-deg(f[1]);
245         rw1 = rw1, deg(f[ii])-deg(f[1]);
246      }
247      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
248   }
249   else
250   {
251      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
252   }
253
254   map psi1 = P,maxideal(1);
255   ideal SBid = psi1(SBid);
256   attrib(SBid,"isSB",1);
257
258 qring newR = std(SBid);
259   map psi = R,ideal(X(1..nvars(R)));
260   ideal id = psi(id);
261   ideal f = psi(f);
262   module syzf = psi(syzf);
264   matrix T,A;
265   list L1;
266
267//---------- computation of Hom(J,J) as affine ring ---------------------------
268// determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
269// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course
270// the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
271// It is a fact, that the kernel is generated by the linear and the quadratic
272// relations
273
274   pf = f[1]*f;
275   T = matrix(ideal(T(1..q)),1,q);
276   Lin = ideal(T*syzf);
277   if(y>=1)
278   {
279      "// the ring structure of Hom(J,J) as R-algebra";
280      " ";
281      "//   the linear relations";
282      " ";
283      Lin;
284      "   ";
285   }
286   for (ii=2; ii<=q; ii++ )
287   {
288      for ( jj=2; jj<=ii; jj++ )
289      {
290         A = lift(pf,f[ii]*f[jj]);
292      }
293   }
294   if(y>=1)
295   {
297      "   ";
299      pause();
300      newline;
301   }
303   Q = subst(Q,T(1),1);
304   Q = interred(reduce(Q,std(0)));
305//---------- reduce number of variables by substitution, if possible ----------
306   if (homo==1)
307   {
308      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
309   }
310   else
311   {
312      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
313   }
314
315   ideal endid  = imap(newR,id)+imap(newR,Q);
316   ideal endphi = ideal(X(1..nvars(R)));
317
318   L=substpart(endid,endphi,homo,rw);
319   def lastRing=L[1];
320   setring lastRing;
321   attrib(endid,"isCohenMacaulay",isCo);
322   attrib(endid,"isPrim",isPr);
323   attrib(endid,"isIsolatedSingularity",isIso);
324   attrib(endid,"isRegInCodim2",isRe);
325   attrib(endid,"isEquidimensional",isEq);
326   attrib(endid,"isCompleteIntersection",0);
328  // export(endid);
329  // export(endphi);
330   if(y>=1)
331   {
332      "//   the new ring after reduction of the number of variables";
333      "   ";
334      lastRing;
335      "   ";
336      "//   the new ideal";
337      "   ";
338      endid;
339      "   ";
340      "//   the old ring";
341      "   ";
342      P;
343      "   ";
344      "//   the old ideal";
345      "   ";
346      setring P;
347      id;
348      "   ";
349      setring lastRing;
350      "//   the map";
351      "   ";
352      endphi;
353      "   ";
354      pause();
355      newline;
356   }
357   L = lastRing;
358   L = insert(L,0,1);
359   return(L);
360}
361example
362{"EXAMPLE:";  echo = 2;
363  ring r   = 0,(x,y),wp(2,3);
364  ideal id = y^2-x^3;
365  ideal J  = x,y;
366  poly p   = x;
367  list Li = std(id),id,J,p;
368  list L   = HomJJ(Li);
369  def end = L[1];    // defines ring L[1], containing ideals endid and endphi
370  setring end;       // makes end the basering
371  end;
372  endid;             // end/endid is isomorphic to End(r/id) as ring
373  map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
374  psi;
375}
376
377///////////////////////////////////////////////////////////////////////////////
378proc normal(ideal id, list #)
379"USAGE:   normal(i [,choose]);  i a radical ideal, choose empty or 1
380         if choose=1 the normalization of the associated primes is computed
381         (which is sometimes more efficient)
382RETURN:  a list of rings (say nor), in each ring nor[i] are two ideals
383         norid, normap such that the direct sum of the rings nor[i]/norid is
384         the normalization of basering/id; normap gives the normalization map
385         from basering/id to nor[i]/norid (for each i)
386NOTE:    to use the i-th ring type: def R=nor[i]; setring R;
387         increasing printlevel displays more comments (default: printlevel=0)
388COMMENT: Not implemented for local or mixed orderings.
389         If the input ideal i is weighted homogeneous a weighted ordering may
390         be used (qhweight(i); computes weights).
391CAUTION: The proc does not check whether the input is radical, for non radical
393EXAMPLE: example normal; shows an example
394"
395{
396   int i,j,y;
397   string sr;
398   list result,prim,keepresult;
399   y = printlevel-voice+2;
400
402   if ( find(ordstr(basering),"s")+find(ordstr(basering),"M") != 0)
403   {
404     "";
405     "// Not implemented for this ordering,";
406     "// please change to global ordering!";
407     return(result);
408   }
409   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
410   {
411      if(attrib(id,"isCompleteIntersection")==1)
412      {
413         attrib(id,"isCohenMacaulay",1);
414         attrib(id,"isEquidimensional",1);
415      }
416   }
417   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
418   {
419      if(attrib(id,"isCohenMacaulay")==1)
420      {
421         attrib(id,"isEquidimensional",1);
422      }
423   }
424   if( typeof(attrib(id,"isPrim"))=="int" )
425   {
426      if(attrib(id,"isPrim")==1)
427      {
428         attrib(id,"isEquidimensional",1);
429      }
430   }
431   if(size(#)==0)
432   {
433      if( typeof(attrib(id,"isEquidimensional"))=="int" )
434      {
435         if(attrib(id,"isEquidimensional")==1)
436         {
437            prim[1]=id;
438         }
439         else
440         {
441            prim=equidim(id);
442         }
443      }
444      else
445      {
446         prim=equidim(id);
447      }
448      if(y>=1)
449      {
450         "// we have ",size(prim),"equidimensional components";
451      }
452   }
453   else
454   {
455      if( typeof(attrib(id,"isPrim"))=="int" )
456      {
457         if(attrib(id,"isPrim")==1)
458         {
459            prim[1]=id;
460         }
461         else
462         {
463            prim=minAssPrimes(id);
464         }
465      }
466      else
467      {
468         prim=minAssPrimes(id);
469      }
470      if(y>=1)
471      {
472         "// we have ",size(prim),"irreducible components";
473      }
474   }
475   for(i=1; i<=size(prim); i++)
476   {
477      if(y>=1)
478      {
479         "// we are in loop ",i;
480      }
481      attrib(prim[i],"isCohenMacaulay",0);
482      if(size(#)!=0)
483      {
484         attrib(prim[i],"isPrim",1);
485      }
486      else
487      {
488         attrib(prim[i],"isPrim",0);
489      }
490      attrib(prim[i],"isRegInCodim2",0);
491      attrib(prim[i],"isIsolatedSingularity",0);
492      attrib(prim[i],"isEquidimensional",1);
493      attrib(prim[i],"isCompleteIntersection",0);
494
495      if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
496      {
497            if(attrib(id,"isIsolatedSingularity")==1)
498             {attrib(prim[i],"isIsolatedSingularity",1); }
499      }
500
501      if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
502      {
503            if((attrib(id,"isIsolatedSingularity")==1)&&(size(#)==0))
504             {attrib(prim[i],"isIsolatedSingularity",1); }
505      }
506      keepresult=normalizationPrimes(prim[i],maxideal(1));
507      for(j=1;j<=size(keepresult);j++)
508      {
509         result=insert(result,keepresult[j]);
510      }
511      sr = string(size(result));
512   }
513      dbprint(y+1,"
514// 'normal' created a list of "+sr+" ring(s).
515// To see the rings, type (if the name of your list is nor):
516     show( nor);
517// To access the 1-st ring and map (and similair for the others), type:
518     def R = nor[1]; setring R;  norid; normap;
519// R/norid is the 1-st ring of the normalization and
520// normap the map from the original basering to R/norid");
521
522      //kill endphi,endid;
523      return(result);
524}
525example
526{ "EXAMPLE:"; echo = 2;
527   ring r=32003,(x,y,z),wp(2,1,2);
528   ideal i=z3-xy4;
529   list nor=normal(i);
530   show(nor);
531   def r1=nor[1];
532   setring r1;
533   norid;
534   normap;
535}
536
537///////////////////////////////////////////////////////////////////////////////
538static
539proc normalizationPrimes(ideal i,ideal ihp, list #)
540"USAGE:   normalizationPrimes(i,ihp[,si]);  i prime ideal, ihp map
541         (partial normalization), si SB of singular locus
542RETURN:  a list of one ring L=R, in  R are two ideals
543         S,M such that R/M is the normalization
544         S is a standardbasis of M
545NOTE:    to use the ring: def r=L[1];setring r;
546         printlevel >= voice+1: display comments (default: printlevel=0)
547EXAMPLE: example normalizationPrimes; shows an example
548"
549{
550   int y = printlevel-voice+2;  // y=printlevel (default: y=0)
551
552   if(y>=1)
553   {
554     "";
555     "// START a normalization loop with the ideal";  "";
556     i;  "";
557     basering;  "";
558     pause();
559     newline;
560   }
561
562   def BAS=basering;
563   list result,keepresult1,keepresult2;
564   ideal J,SB,MB;
565   int depth,lauf,prdim;
566   int ti=timer;
567
568   if(size(i)==0)
569   {
570      if(y>=1)
571      {
572          "// the ideal was the zero-ideal";
573      }
574         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
575                      +ordstr(basering)+");");
576         ideal norid=ideal(0);
577         ideal normap=fetch(BAS,ihp);
578         export norid;
579         export normap;
580         result=newR7;
581         setring BAS;
582         return(result);
583   }
584
585   if(y>=1)
586   {
587     "// SB-computation of the input ideal";
588   }
589
590   list SM=mstd(i);                //here the work starts
591   int dimSM =  dim(SM[1]);        //dimension of variety to normalize
592  // Case: Get an ideal containing a unit
593   if( dimSM == -1)
594   {  "";
595      "      // A unit ideal was found.";
596      "      // Stop with partial result computed so far";"";
597
598         MB=SM[2];
599         intvec rw;
600         list LL=substpart(MB,ihp,0,rw);
601         def newR6=LL[1];
602         setring newR6;
603         ideal norid=endid;
604         ideal normap=endphi;
605         kill endid,endphi;
606         export norid;
607         export normap;
608         result=newR6;
609         setring BAS;
610         return(result);
611   }
612
613   if(y>=1)
614   {
615      "//   the dimension is:"; "";
616      dimSM;"";
617   }
618
619   if(size(#)>0)
620   {
621      list JM=mstd(#[1]);
623      {
625      }
626   }
627
628   if(attrib(i,"isPrim")==1)
629   {
630      attrib(SM[2],"isPrim",1);
631   }
632   else
633   {
634      attrib(SM[2],"isPrim",0);
635   }
636   if(attrib(i,"isIsolatedSingularity")==1)
637   {
638      attrib(SM[2],"isIsolatedSingularity",1);
639   }
640   else
641   {
642      attrib(SM[2],"isIsolatedSingularity",0);
643   }
644   if(attrib(i,"isCohenMacaulay")==1)
645   {
646      attrib(SM[2],"isCohenMacaulay",1);
647   }
648   else
649   {
650      attrib(SM[2],"isCohenMacaulay",0);
651   }
652   if(attrib(i,"isRegInCodim2")==1)
653   {
654      attrib(SM[2],"isRegInCodim2",1);
655   }
656   else
657   {
658      attrib(SM[2],"isRegInCodim2",0);
659   }
660   if(attrib(i,"isEquidimensional")==1)
661   {
662      attrib(SM[2],"isEquidimensional",1);
663   }
664   else
665   {
666      attrib(SM[2],"isEquidimensional",0);
667   }
668    if(attrib(i,"isCompleteIntersection")==1)
669   {
670      attrib(SM[2],"isCompleteIntersection",1);
671   }
672   else
673   {
674      attrib(SM[2],"isCompleteIntersection",0);
675   }
676
677   //the smooth case
678   if(size(#)>0)
679   {
680      if(dim(JM[1])==-1)
681      {
682         if(y>=1)
683         {
684            "// the ideal was smooth";
685         }
686         MB=SM[2];
687         intvec rw;
688         list LL=substpart(MB,ihp,0,rw);
689         def newR6=LL[1];
690         setring newR6;
691         ideal norid=endid;
692         ideal normap=endphi;
693         kill endid,endphi;
694         export norid;
695         export normap;
696         result=newR6;
697         setring BAS;
698         return(result);
699     }
700   }
701
702   //the zero-dimensional case
703   if((dim(SM[1])==0)&&(homog(SM[2])==1))
704   {
705      if(y>=1)
706      {
707         "// the ideal was zero-dimensional and homogeneous";
708      }
709      MB=maxideal(1);
710      intvec rw;
711      list LL=substpart(MB,ihp,0,rw);
712      def newR5=LL[1];
713      setring newR5;
714      ideal norid=endid;
715      ideal normap=endphi;
716      kill endid,endphi;
717      export norid;
718      export normap;
719      result=newR5;
720      setring BAS;
721      return(result);
722   }
723
724   //the one-dimensional case
725   //in this case it is a line because
726   //it is irreducible and homogeneous
727   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
728        &&(homog(SM[2])==1))
729   {
730      if(y>=1)
731      {
732         "// the ideal defines a line";
733      }
734      MB=SM[2];
735      intvec rw;
736      list LL=substpart(MB,ihp,0,rw);
737      def newR4=LL[1];
738      setring newR4;
739      ideal norid=endid;
740      ideal normap=endphi;
741      kill endid,endphi;
742      export norid;
743      export normap;
744      result=newR4;
745      setring BAS;
746      return(result);
747   }
748
749   //the higher dimensional case
750   //we test first of all CohenMacaulay and
751   //complete intersection
752   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
753   {
754      //test for complete intersection
755      attrib(SM[2],"isCohenMacaulay",1);
756      attrib(SM[2],"isCompleteIntersection",1);
757      attrib(SM[2],"isEquidimensional",1);
758      if(y>=1)
759      {
760         "// the ideal is a complete intersection";
761      }
762   }
763
764   //compute the singular locus+lower dimensional components
765   if(((attrib(SM[2],"isIsolatedSingularity")==0)||(homog(SM[2])==0))
766        &&(size(#)==0))
767   {
768/*
769write (":a normal-fehler" ,
770         "basering:",string(basering),"nvars:", nvars(basering),
771       "dim(SM[1]):",dim(SM[1]),"ncols(jacob(SM[2]))",ncols(jacob(SM[2])),
772       "SM:", SM);
773
774     pause();
775*/
776      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
777      //ti=timer;
778      if(y >=1 )
779      {
780         "// SB of singular locus will be computed";
781      }
782      ideal sin=J+SM[2];
783
784    //kills the embeded components
785
786      list JM=mstd(sin);
787      //JM[1] SB os singular locus, JM[2]=minbasis of singular locus
788      //SM[1] SB of irreducible component, SM[2] minbasis
789      if(y>=1)
790      {
791         "//   the dimension of the singular locus is:";"";
792         dim(JM[1]); "";
793      }
795      //   timer-ti;
796      attrib(JM[1],"isSB",1);
797      if(dim(JM[1])==-1)
798      {
799         if(y>=1)
800         {
801            "// the ideal is smooth";
802         }
803         MB=SM[2];
804         intvec rw;
805         list LL=substpart(MB,ihp,0,rw);
806         def newR3=LL[1];
807         setring newR3;
808         ideal norid=endid;
809         ideal normap=endphi;
810         kill endid,endphi;
811         export norid;
812         export normap;
813         result=newR3;
814         setring BAS;
815         return(result);
816      }
817      if(dim(JM[1])==0)
818      {
819         attrib(SM[2],"isIsolatedSingularity",1);
820      }
821      if(dim(JM[1])<=dim(SM[1])-2)
822      {
823         attrib(SM[2],"isRegInCodim2",1);
824      }
825   }
826   else
827   {
828     if(size(#)==0)
829     {
830        list JM=maxideal(1),maxideal(1);
831        attrib(JM[1],"isSB",1);
832        attrib(SM[2],"isRegInCodim2",1);
833     }
834   }
835   if((attrib(SM[2],"isRegInCodim2")==1)&&(attrib(SM[2],"isCohenMacaulay")==1))
836   {
837      if(y>=1)
838      {
839            "// the ideal was CohenMacaulay and regular in codimension 2";
840      }
841      MB=SM[2];
842      intvec rw;
843      list LL=substpart(MB,ihp,0,rw);
844      def newR6=LL[1];
845      setring newR6;
846      ideal norid=endid;
847      ideal normap=endphi;
848      kill endid,endphi;
849      export norid;
850      export normap;
851      result=newR6;
852      setring BAS;
853      return(result);
854   }
855   //if it is an isolated singularity things are easier
856   //JM ideal of singular locus, SM ideal of variety
857   if((dim(JM[1])==0)&&(homog(SM[2])==1))            //isolated sing. case
858   {
859      attrib(SM[2],"isIsolatedSingularity",1);
860      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
861           //vars not contained in ideal
862      ideal Ann=quotient(SM[2],SL[1]);
863      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
864
865      //qAnn=0 ==> the first var(=SL[1]) not contained in SM is a nzd of R/SM
866      if(size(qAnn)==0)
867      {
868         if(y>=1)
869         {
870            "";
872            "";
873            maxideal(1);
874            newline;
875         }
876         //again test for normality
877         //Hom(I,R)=R
878         list RR;
879         RR=SM[1],SM[2],maxideal(1),SL[1];
880         ti=timer;
881         RR=HomJJ(RR,y);
882         if(RR[2]==0)
883         {
884            def newR=RR[1];
885            setring newR;
886            map psi=BAS,endphi;
887         //   ti=timer;
888            list tluser=normalizationPrimes(endid,psi(ihp));
889
890        //    timer-ti;
891            setring BAS;
892            return(tluser);
893         }
894         MB=SM[2];
895         execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
896                      +ordstr(basering)+");");
897         ideal norid=fetch(BAS,MB);
898         ideal normap=fetch(BAS,ihp);
899         export norid;
900         export normap;
901         result=newR7;
902         // the following 2 lines don't work : nor is not defined
903         //def R = nor[1]; setring R;     //make the 1-st ring the basering
904         //norid; normap;                 //data of the normalization)
905         setring BAS;
906         return(result);
907
908       }
909      //Now the case where qAnn!=0, i.e.SL[1] is a zero divisor of R/SM
910      //and we have found a splitting: id and id1
911      //id=qAnn+SM[2] defines components of R/SM in the complement of V(SL[1])
912      //id1 defines components of R/SM in the complement of V(id)
913      //?????instead of id1 we can take SL[1]+Ann+SM[2]???????????
914       else
915       {
916          ideal id=qAnn+SM[2];
917
918          attrib(id,"isCohenMacaulay",0);
919          attrib(id,"isPrim",0);
920          attrib(id,"isIsolatedSingularity",1);
921          attrib(id,"isRegInCodim2",0);
922          attrib(id,"isCompleteIntersection",0);
923          attrib(id,"isEquidimensional",0);
924
925          keepresult1=normalizationPrimes(id,ihp);
926          ideal id1=quotient(SM[2],Ann)+SM[2];
927//          evtl. qAnn statt Ann nehmen
928//          ideal id=SL[1]+SM[2];
929
930          attrib(id1,"isCohenMacaulay",0);
931          attrib(id1,"isPrim",0);
932          attrib(id1,"isIsolatedSingularity",1);
933          attrib(id1,"isRegInCodim2",0);
934          attrib(id1,"isCompleteIntersection",0);
935          attrib(id1,"isEquidimensional",0);
936
937          keepresult2=normalizationPrimes(id1,ihp);
938
939          for(lauf=1;lauf<=size(keepresult2);lauf++)
940          {
941             keepresult1=insert(keepresult1,keepresult2[lauf]);
942          }
943          return(keepresult1);
944       }
945   }
946
947   //test for non-normality
948   //Hom(I,I)<>R
949   //we can use Hom(I,I) to continue
950
951   ideal SL=simplify(reduce(JM[2],SM[1]),2);
952   ideal Ann=quotient(SM[2],SL[1]);
953   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
954
955   if(size(qAnn)==0)
956   {
957      list RR;
958      list RS;
959      //now we have to compute the radical
960      if(y>=1)
961      {
962         "// radical computation of singular locus";
963      }
964
966      {
969
970          // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
971      }
973      {
974          ideal JJ=SM[2]+ideal(SL[1]);
975         // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
976          if(attrib(SM[2],"isCompleteIntersection")==0)
977          {
980          }
981          else
982          {
984             J=quotient(JJ,minor(jacob(JJ),size(JJ)));
985          }
986      }
987      if(y>=1)
988      {
989        "//   radical is equal to:";"";
990        J;
991        "";
992      }
993
994      JM=J,J;
995
996      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
997      RR=SM[1],SM[2],JM[2],SL[1];
998
999      //   evtl eine geeignete Potenz von JM?
1000     if(y>=1)
1001     {
1003     }
1004
1005     RS=HomJJ(RR,y);
1006
1007      if(RS[2]==1)
1008      {
1009         def lastR=RS[1];
1010         setring lastR;
1011         map psi1=BAS,endphi;
1012         ideal norid=endid;
1013         ideal normap=psi1(ihp);
1014         kill endid,endphi;
1015         export norid;
1016         export normap;
1017         setring BAS;
1018         return(lastR);
1019      }
1020      int n=nvars(basering);
1021      ideal MJ=JM[2];
1022
1023      def newR=RS[1];
1024      setring newR;
1025
1026      map psi=BAS,endphi;
1027      list tluser=
1028             normalizationPrimes(endid,psi(ihp),simplify(psi(MJ)+endid,4));
1029      setring BAS;
1030      return(tluser);
1031   }
1032    // A component with singular locus the whole component found
1033   if( Ann == 1)
1034   {
1035      "// Input appeared not to be a radical ideal!";
1036      "// A (everywhere singular) component with ideal";
1037      "// equal to its Jacobian ideal was found";
1038      "// Procedure will stop with partial result computed so far";"";
1039
1040         MB=SM[2];
1041         intvec rw;
1042         list LL=substpart(MB,ihp,0,rw);
1043         def newR6=LL[1];
1044         setring newR6;
1045         ideal norid=endid;
1046         ideal normap=endphi;
1047         kill endid,endphi;
1048         export norid;
1049         export normap;
1050         result=newR6;
1051         setring BAS;
1052         return(result);
1053   }
1054   else
1055   {
1056      int equi=attrib(SM[2],"isEquidimensional");
1057      ideal new1=qAnn+SM[2];
1058      execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
1059                      +ordstr(basering)+");");
1060      if(y>=1)
1061      {
1062         "// zero-divisor found";
1063      }
1064      ideal vid=fetch(BAS,new1);
1065      ideal ihp=fetch(BAS,ihp);
1066      attrib(vid,"isCohenMacaulay",0);
1067      attrib(vid,"isPrim",0);
1068      attrib(vid,"isIsolatedSingularity",0);
1069      attrib(vid,"isRegInCodim2",0);
1070      if(equi==1)
1071      {
1072         attrib(vid,"isEquidimensional",1);
1073      }
1074      else
1075      {
1076         attrib(vid,"isEquidimensional",0);
1077      }
1078      attrib(vid,"isCompleteIntersection",0);
1079
1080      keepresult1=normalizationPrimes(vid,ihp);
1081
1082      setring BAS;
1083      ideal new2=quotient(SM[2],Ann)+SM[2];
1084// evtl. qAnn nehmen
1085      execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
1086                      +ordstr(basering)+");");
1087
1088      ideal vid=fetch(BAS,new2);
1089      ideal ihp=fetch(BAS,ihp);
1090      attrib(vid,"isCohenMacaulay",0);
1091      attrib(vid,"isPrim",0);
1092      attrib(vid,"isIsolatedSingularity",0);
1093      attrib(vid,"isRegInCodim2",0);
1094      if(equi==1)
1095      {
1096         attrib(vid,"isEquidimensional",1);
1097      }
1098      else
1099      {
1100         attrib(vid,"isEquidimensional",0);
1101      }
1102      attrib(vid,"isCompleteIntersection",0);
1103
1104      keepresult2=normalizationPrimes(vid,ihp);
1105
1106      setring BAS;
1107      for(lauf=1;lauf<=size(keepresult2);lauf++)
1108      {
1109         keepresult1=insert(keepresult1,keepresult2[lauf]);
1110      }
1111      return(keepresult1);
1112   }
1113}
1114example
1115{ "EXAMPLE:";echo = 2;
1116   //Huneke
1117   ring qr=31991,(a,b,c,d,e),dp;
1118   ideal i=
1119   5abcde-a5-b5-c5-d5-e5,
1121   a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1122   abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1124   a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1126   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1127
1128   list pr=normalizationPrimes(i);
1129   def r1=pr[1];
1130   setring r1;
1131   norid;
1132   normap;
1133}
1134///////////////////////////////////////////////////////////////////////////////
1135static
1136proc substpart(ideal endid, ideal endphi, int homo, intvec rw)
1137
1138"//Repeated application of elimpart to endid, until no variables can be
1139//directy substituded. homo=1 if input is homogeneous, rw contains
1140//original weights, endphi (partial) normalization map";
1141
1142{
1143   def newRing=basering;
1144   int ii,jj;
1145   map phi = basering,maxideal(1);
1146
1147   //endid=diagon(endid);
1148
1149   list Le = elimpart(endid);
1150                                     //this proc and the next loop try to
1151   int q = size(Le[2]);              //substitute as many variables as possible
1152   intvec rw1 = 0;                   //indices of substituted variables
1153   rw1[nvars(basering)] = 0;
1154   rw1 = rw1+1;
1155
1156   while( size(Le[2]) != 0 )
1157   {
1158      endid = Le[1];
1159      map ps = newRing,Le[5];
1160
1161      phi = ps(phi);
1162      for(ii=1;ii<=size(Le[2])-1;ii++)
1163      {
1164         phi=phi(phi);
1165      }
1166      //eingefuegt wegen x2-y2z2+z3
1167      kill ps;
1168
1169      for( ii=1; ii<=size(rw1); ii++ )
1170      {
1171         if( Le[4][ii]==0 )
1172         {
1173            rw1[ii]=0;                             //look for substituted vars
1174         }
1175      }
1176      Le=elimpart(endid);
1177      q = q + size(Le[2]);
1178   }
1179   endphi = phi(endphi);
1180
1181//---------- return -----------------------------------------------------------
1182// in the homogeneous case put weights for the remaining vars correctly, i.e.
1183// delete from rw those weights for which the corresponding entry of rw1 is 0
1184
1185   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
1186   {
1187      jj=1;
1188      for( ii=2; ii<size(rw1); ii++)
1189      {
1190         jj++;
1191         if( rw1[ii]==0 )
1192         {
1193            rw=rw[1..jj-1],rw[jj+1..size(rw)];
1194            jj=jj-1;
1195         }
1196      }
1197      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
1198      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }
1199
1200      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),(a(rw),dp);
1201   }
1202   else
1203   {
1204      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
1205   }
1206
1207   ideal lastmap;
1208   q = 1;
1209   for(ii=1; ii<=size(rw1); ii++ )
1210   {
1211      if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; }
1212      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
1213   }
1214   map phi1 = newRing,lastmap;
1215   ideal endid  = phi1(endid);
1216   ideal endphi = phi1(endphi);
1217   export(endid);
1218   export(endphi);
1219   list L = lastRing;
1220   setring newRing;
1221   return(L);
1222}
1223///////////////////////////////////////////////////////////////////////////////
1224static
1225proc diagon(ideal i)
1226{
1227   matrix m;
1228   intvec iv = option(get);
1229   option(redSB);
1230   ideal j=liftstd(jet(i,1),m);
1231   option(set,iv);
1232   return(ideal(matrix(i)*m));
1233}
1234/////////////////////////////////////////////////////////////////////////////
1235/*
1236                           Examples:
1237LIB"normal.lib";
1238//Huneke
1239ring qr=31991,(a,b,c,d,e),dp;
1240ideal i=
12415abcde-a5-b5-c5-d5-e5,
1243a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1244abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1246a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1248b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1249
1250
1251//Vasconcelos
1252ring r=32003,(x,y,z,w,t),dp;
1253ideal i=
1254x2+zw,
1255y3+xwt,
1256xw3+z3t+ywt2,
1257y2w4-xy2z2t-w3t3;
1258
1259//Theo1
1260ring r=32003,(x,y,z),wp(2,3,6);
1261ideal i=zy2-zx3-x6;
1262
1263//Theo1a (CohenMacaulay and regular in codimension 2)
1264ring r=32003,(x,y,z,u),wp(2,3,6,6);
1265ideal i=zy2-zx3-x6+u2;
1266
1267
1268//Theo2
1269ring r=32003,(x,y,z),wp(3,4,12);
1270ideal i=z*(y3-x4)+x8;
1271
1272//Theo2a
1273ring r=32003,(T(1..4)),wp(3,4,12,17);
1274ideal i=
1275T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
1276T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
1277T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
1278T(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;
1279
1280//Theo3
1281ring r=32003,(x,y,z),wp(3,5,15);
1282ideal i=z*(y3-x5)+x10;
1283
1284
1285//Theo4
1286ring r=32003,(x,y,z),dp;
1287ideal i=(x-y)*(x-z)*(y-z);
1288
1289//Theo5
1290ring r=32003,(x,y,z),wp(2,1,2);
1291ideal i=z3-xy4;
1292
1293//Theo6
1294ring r=32003,(x,y,z),dp;
1295ideal i=x2y2+x2z2+y2z2;
1296
1297ring r=32003,(a,b,c,d,e,f),dp;
1298ideal i=
1299bf,
1300af,
1301bd,
1303
1304//Beispiel, wo vorher Primaerzerlegung schneller
1305//ist CM
1306//Sturmfels
1307ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1308ideal i=
1309bv+su,
1310bw+tu,
1311sw+tv,
1312by+sx,
1313bz+tx,
1314sz+ty,
1315uy+vx,
1316uz+wx,
1317vz+wy,
1318bvz;
1319
1320//J S/Y
1321ring r=32003,(x,y,z,t),dp;
1322ideal i=
1323x2z+xzt,
1324xyz,
1325xy2-xyt,
1326x2y+xyt;
1327
1328//St_S/Y
1329ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1330ideal i=
1331wy-vz,
1332vx-uy,
1333tv-sw,
1334su-bv,
1335tuy-bvz;
1336
1337//dauert laenger
1338//Horrocks:
1339ring r=32003,(a,b,c,d,e,f),dp;
1340ideal i=
1343abdf-16000b2ef+16001bcf2,
1344a2df+8002abef+8001acf2,
1346acde-16000bce2+16001c2ef,
1347a2de-8000abe2-7999acef,
1348acd2+8002bcde+8001c2df,
1349abd2-8000b2de-7999bcdf,
1350a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
1351abde-8000b2e2-acdf-16001bcef-8001c2f2,
1352abcd-16000b2ce+16001bc2f,
1353a2cd+8002abce+8001ac2f,
1354a2bd-8000ab2e-7999abcf,
1355ab3f-3bdf3,
1357a3bf+4aef3,
1358ac3e-10668cde3,
1360a3ce+10669ae3f,
1361bc3d+8001cd3e,
1362ac3d+8000bc3e+16001cd2e2+8001c4f,
1364b2c2e-10668bc3f-10667cd2ef,
1365abc2e-cde2f,
1366b3cd-8000bd3f,
1367b3ce-10668b2c2f-10667bd2ef,
1368abc2f-cdef2,
1369a2bce-16000be3f+16001ce2f2,
1370ab3d-8000b4e-8001b3cf+16000bd2f2,
1371ab2cf-bdef2,
1372a2bcf-16000be2f2+16001cef3,
1373a4d-8000a3be+8001a3cf-2ae2f2;
1374
1375
1376ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1377
1378ideal k=
1379wy-vz,
1380vx-uy,
1381tv-sw,
1382su-bv,
1383tuy-bvz;
1384ideal j=x2y2+x2z2+y2z2;
1385ideal i=mstd(intersect(j,k))[2];
1386
1387//22
1388ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1389ideal i=
1390wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
1391vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
1392tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
1393sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
1394tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;
1395
1396
1397//riemenschneider
1398//33
1399//normal+primary 3
1400//primary 9
1402//minAssPrimes 2
1403ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
1404ideal i=
1405xz,
1406vx,
1407ux,
1408su,
1409qu,
1410txy,
1411stx,
1412qtx,
1413uv2z-uwz,
1414uv3-uvw,
1415puv2-puw;
1416
1417ring r=0,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
1418ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;
1419
1420
1421
1422*/
1423
Note: See TracBrowser for help on using the repository browser.