source: git/Singular/LIB/normal.lib @ a59d60

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