source: git/Singular/LIB/normal.lib @ 75089b

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