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

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