source: git/Singular/LIB/normal.lib @ 72c1b9

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