source: git/Singular/LIB/normal.lib @ 82716e

spielwiese
Last change on this file since 82716e was 82716e, checked in by Hans Schönemann <hannes@…>, 26 years ago
*hannes: typos in the info-help-string git-svn-id: file:///usr/local/Singular/svn/trunk@1773 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 34.3 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.6 1998-05-14 18:45:10 Singular 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   ring  r = 0,(x,y,z),dp;
655
656}
657
658
659proc normalizationPrimes(ideal i,ideal ihp, list #)
660"USAGE:   normalizationPrimes(i);  i prime ideal
661RETURN:  a list of one ring L=R, in  R are two ideals
662         S,M such that R/M is the normalization
663         S is a standardbasis of M
664NOTE:    to use the ring: def r=L[1];setring r;
665EXAMPLE: example normalizationPrimes; shows an example
666"
667{
668   int y;
669   if(y==1)
670   {
671      "START one normalization loop with the ideal";
672      "        ";
673      i;
674      "        ";
675      basering;
676      "        ";
677      pause;
678   }
679
680   def BAS=basering;
681   list result,keepresult1,keepresult2;
682   ideal J,SB,MB,KK;
683   int depth,lauf,prdim;
684   int ti=timer;
685
686   if(size(i)==0)
687   {
688      if(y==1)
689      {
690          "the ideal was the zero-ideal";
691      }
692         execute "ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
693                      +ordstr(basering)+");";
694         ideal KK=ideal(0);
695         ideal PP=fetch(BAS,ihp);
696         export PP;
697         export KK;
698         result=newR7;
699         setring BAS;
700         return(result);
701
702   }
703
704   if(y==1)
705   {
706      "                                  ";
707      "                                  ";
708      " SB-computation of the input ideal";
709      "                                  ";
710   }
711   list SM=mstd(i);
712
713   if(y==1)
714   {
715      " the dimension is                 ";
716      dim(SM[1]);
717      "                                  ";
718   }
719
720   if(size(#)>0)
721   {
722      list JM=#[1],#[1];
723      attrib(JM[1],"isSB",1);
724      if( typeof(attrib(#[1],"isRad"))!="int" )
725      {
726         attrib(JM[2],"isRad",0);
727      }
728   }
729
730   if(attrib(i,"isPrim")==1)
731   {
732      attrib(SM[2],"isPrim",1);
733   }
734   else
735   {
736      attrib(SM[2],"isPrim",0);
737   }
738   if(attrib(i,"isIsolatedSingularity")==1)
739   {
740      attrib(SM[2],"isIsolatedSingularity",1);
741   }
742   else
743   {
744      attrib(SM[2],"isIsolatedSingularity",0);
745   }
746   if(attrib(i,"isCohenMacaulay")==1)
747   {
748      attrib(SM[2],"isCohenMacaulay",1);
749   }
750   else
751   {
752      attrib(SM[2],"isCohenMacaulay",0);
753   }
754   if(attrib(i,"isRegInCodim2")==1)
755   {
756      attrib(SM[2],"isRegInCodim2",1);
757   }
758   else
759   {
760      attrib(SM[2],"isRegInCodim2",0);
761   }
762    if(attrib(i,"isEquidimensional")==1)
763   {
764      attrib(SM[2],"isEquidimensional",1);
765   }
766   else
767   {
768      attrib(SM[2],"isEquidimensional",0);
769   }
770    if(attrib(i,"isCompleteIntersection")==1)
771   {
772      attrib(SM[2],"isCompleteIntersection",1);
773   }
774   else
775   {
776      attrib(SM[2],"isCompleteIntersection",0);
777   }
778
779   //the smooth case
780   if(size(#)>0)
781   {
782      if(dim(JM[1])==-1)
783      {
784         if(y==1)
785         {
786            "the ideal was smooth";
787         }
788         MB=SM[2];
789         execute "ring newR6="+charstr(basering)+",("+varstr(basering)+"),("
790                      +ordstr(basering)+");";
791         ideal KK=fetch(BAS,MB);
792         ideal PP=fetch(BAS,ihp);
793         export PP;
794         export KK;
795         result=newR6;
796         setring BAS;
797         return(result);
798     }
799   }
800
801   //the zero-dimensional case
802   if((dim(SM[1])==0)&&(homog(SM[2])==1))
803   {
804      if(y==1)
805      {
806         "the ideal was zero-dimensional and homogeneous";
807      }
808      MB=maxideal(1);
809      execute "ring newR5="+charstr(basering)+",("+varstr(basering)+"),("
810                      +ordstr(basering)+");";
811      ideal KK=fetch(BAS,MB);
812      ideal PP=fetch(BAS,ihp);
813      export PP;
814      export KK;
815      result=newR5;
816      setring BAS;
817      return(result);
818   }
819
820   //the one-dimensional case
821   //in this case it is a line because
822   //it is irreducible and homogeneous
823   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
824        &&(homog(SM[2])==1))
825   {
826      if(y==1)
827      {
828         "the ideal defines a line";
829      }
830      MB=SM[2];
831      execute "ring newR4="+charstr(basering)+",("+varstr(basering)+"),("
832                      +ordstr(basering)+");";
833      ideal KK=fetch(BAS,MB);
834      ideal PP=fetch(BAS,ihp);
835      export PP;
836      export KK;
837      result=newR4;
838      setring BAS;
839      return(result);
840   }
841
842   //the higher dimensional case
843   //we test first of all CohenMacaulay and
844   //complete intersection
845   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
846   {
847      //test for complete intersection
848      attrib(SM[2],"isCohenMacaulay",1);
849      attrib(SM[2],"isCompleteIntersection",1);
850      attrib(SM[2],"isEquidimensional",1);
851      if(y==1)
852      {
853         "it is a complete Intersection";
854      }
855   }
856//   if(attrib(i,"isIsolatedSingularity")==0)
857//   {
858//      list L=sres(SM[2],0);
859//      prdim=ncols(betti(L))-1;
860//      depth=nvars(basering)-prdim;
861//   }
862
863  // if(attrib(SM[2],"isCohenMacaulay")==0)
864  // {
865  //    test for CohenMacaulay
866  //    if((dim(SM[1]))==depth)
867  //    {
868  //    attrib(SM[2],"isCohenMacaulay",1);
869  //    "it is CohenMacaulay";
870  //    }
871  // }
872
873   //compute the singular locus+lower dimensional components
874   if(((attrib(SM[2],"isIsolatedSingularity")==0)||(homog(SM[2])==0))
875        &&(size(#)==0))
876   {
877
878      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
879      //ti=timer;
880      if(y==1)
881      {
882         "SB of singular locus will be computed";
883      }
884      ideal sin=J+SM[2];
885
886    //kills the embeded components
887//    if(homog(SM[2])==1)
888//    {
889//       sin=sat(sin,maxideal(1))[1];
890//    }
891
892      list JM=mstd(sin);
893      if(y==1)
894      {
895         "                                  ";
896         "the dimension of the singular locus is";
897         dim(JM[1]);
898         "                                  ";
899      }
900
901      attrib(JM[2],"isRad",0);
902      //   timer-ti;
903      attrib(JM[1],"isSB",1);
904      if(dim(JM[1])==-1)
905      {
906         if(y==1)
907         {
908            "it is smooth";
909         }
910         MB=SM[2];
911         execute "ring newR3="+charstr(basering)+",("+varstr(basering)+"),("
912                      +ordstr(basering)+");";
913         ideal KK=fetch(BAS,MB);
914         ideal PP=fetch(BAS,ihp);
915         export PP;
916         export KK;
917         result=newR3;
918         setring BAS;
919         return(result);
920      }
921      if(dim(JM[1])==0)
922      {
923         attrib(SM[2],"isIsolatedSingularity",1);
924      }
925      if(dim(JM[1])<=nvars(basering)-2)
926      {
927         attrib(SM[2],"isRegInCodim2",1);
928      }
929   }
930   else
931   {
932     if(size(#)==0)
933     {
934        list JM=maxideal(1),maxideal(1);
935        attrib(JM[1],"isSB",1);
936        attrib(SM[2],"isRegInCodim2",1);
937     }
938   }
939   //if it is an isolated singularity things are easier
940   if((dim(JM[1])==0)&&(homog(SM[2])==1))
941   {
942      attrib(SM[2],"isIsolatedSingularity",1);
943      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
944      ideal Ann=quotient(SM[2],SL[1]);
945      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
946
947      if(size(qAnn)==0)
948      {
949         if(y==1)
950         {
951            "                                  ";
952            "the ideal rad(J):                 ";
953            "                                  ";
954            maxideal(1);
955            "                                  ";
956            "                                  ";
957         }
958         //again test for normality
959         //Hom(I,R)=R
960         list RR;
961  //       list RR=SM[1],maxideal(1),SL[1];
962  //       "test with Hom(I,R)";
963  //       ti=timer;
964  //       if(isR_HomJR(RR)==1)
965  //       {
966  //         "it was normal";
967  //          timer-ti;
968  //          SB=SM[1];
969  //          SM=SM[2];
970  //          export SB,MB;
971  //          result=BAS;
972  //          return(result);
973  //       }
974  //       timer-ti;
975         RR=SM[1],SM[2],maxideal(1),SL[1];
976         ti=timer;
977         RR=HomJJ(RR,y);
978   //      timer-ti;
979         if(RR[2]==0)
980         {
981            def newR=RR[1];
982            setring newR;
983            map psi=BAS,endphi;
984         //   ti=timer;
985            list tluser=normalizationPrimes(endid,psi(ihp));
986
987        //    timer-ti;
988            setring BAS;
989            return(tluser);
990         }
991         MB=SM[2];
992         execute "ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
993                      +ordstr(basering)+");";
994         ideal KK=fetch(BAS,MB);
995         ideal PP=fetch(BAS,ihp);
996         export PP;
997         export KK;
998         result=newR7;
999         setring BAS;
1000         return(result);
1001
1002       }
1003       else
1004       {
1005          ideal id=qAnn+SM[2];
1006
1007          attrib(id,"isCohenMacaulay",0);
1008          attrib(id,"isPrim",0);
1009          attrib(id,"isIsolatedSingularity",1);
1010          attrib(id,"isRegInCodim2",0);
1011          attrib(id,"isCompleteIntersection",0);
1012          attrib(id,"isEquidimensional",0);
1013
1014          keepresult1=normalizationPrimes(id,ihp);
1015          ideal id1=quotient(SM[2],Ann)+SM[2];
1016//          evtl. qAnn nehmen
1017//          ideal id=SL[1]+SM[2];
1018
1019          attrib(id1,"isCohenMacaulay",0);
1020          attrib(id1,"isPrim",0);
1021          attrib(id1,"isIsolatedSingularity",1);
1022          attrib(id1,"isRegInCodim2",0);
1023          attrib(id1,"isCompleteIntersection",0);
1024          attrib(id1,"isEquidimensional",0);
1025
1026          keepresult2=normalizationPrimes(id1,ihp);
1027
1028          for(lauf=1;lauf<=size(keepresult2);lauf++)
1029          {
1030             keepresult1=insert(keepresult1,keepresult2[lauf]);
1031          }
1032          return(keepresult1);
1033       }
1034   }
1035
1036   //test for non-normality
1037   //Hom(I,I)<>R
1038   //we can use Hom(I,I) to continue
1039
1040   ideal SL=simplify(reduce(JM[2],SM[1]),2);
1041   ideal Ann=quotient(SM[2],SL[1]);
1042   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
1043
1044   if(size(qAnn)==0)
1045   {
1046      list RR;
1047 //      "test with Hom(I,I) before the radical computation";
1048        //again test for normality
1049        //Hom(I,R)=R
1050 //     list RR=SM[1],JM[2],SL[1];
1051 //     "test with Hom(I,R)";
1052 //     ti=timer;
1053 //     if(isR_HomJR(RR)==1)
1054 //     {
1055 //      "it was normal";
1056 //        timer-ti;
1057 //        SB=SM[1];
1058 //        SM=SM[2];
1059 //        export SB,MB;
1060 //        result=BAS;
1061 //        return(result);
1062 //     }
1063 //     timer-ti;
1064
1065 //     list  RR=SM[1],SM[2],JM[2],SL[1];
1066 //     ti=timer;
1067      list RS;
1068 //   list RS=HomJJ(RR);
1069 //   timer-ti;
1070 //     if(RS[2]==0)
1071 //     {
1072 //        "Hom(I,I) was sucessfull without radical";
1073 //        def newR=RS[1];
1074 //        setring newR;
1075 //        list tluser=normalizationPrimes(SM);
1076 //        setring BAS;
1077 //        return(tluser);
1078 //     }
1079
1080      //now we have to compute the radical
1081      if(y==1)
1082      {
1083         "radical computation";
1084      }
1085//      ti=timer;
1086
1087
1088      if((attrib(JM[2],"isRad")==0)&&(attrib(SM[2],"isEquidimensional")==0))
1089      {
1090           //J=radical(JM[2]);
1091          J=radical(SM[2]+ideal(SL[1]));
1092
1093          // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
1094      }
1095      if((attrib(JM[2],"isRad")==0)&&(attrib(SM[2],"isEquidimensional")==1))
1096      {
1097          ideal JJ=SM[2]+ideal(SL[1]);
1098         // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
1099          if(attrib(SM[2],"isCompleteIntersection")==0)
1100          {
1101             J=equiRadical(JM[2]);
1102             //J=equiRadical(JJ);
1103          }
1104          else
1105          {
1106             //J=radical(JM[2]);
1107             J=quotient(JJ,minor(jacob(JJ),size(JJ)));
1108          }
1109      }
1110//    timer-ti;
1111
1112      JM=J,J;
1113
1114      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
1115      // SL=simplify(reduce(J,SM[1]),2);
1116      // Ann=quotient(SM[2],SL[1]);
1117      // qAnn=simplify(reduce(Ann,SM[1]),2);
1118      // if(size(qAnn)!=0)
1119      // {
1120        // keepresult1=normalizationPrimes(qAnn+SM[2]);
1121        // keepresult2=normalizationPrimes(SL[1]+SM[2]);
1122        // for(lauf=1;lauf<=size(keepresult2);j++)
1123        // {
1124        //    keepresult1=insert(keepresult1,keepresult2[lauf]);
1125        // }
1126        // return(keepresult1);
1127      // }
1128      RR=SM[1],SM[2],JM[2],SL[1];
1129
1130//   evtl eine geeignete Potenz von JM?
1131     if(y==1)
1132     {
1133        "compute Hom(rad(J),rad(J)):";
1134     }
1135//   ti=timer;
1136
1137     RS=HomJJ(RR,y);
1138//   timer-ti;
1139
1140      if(RS[2]==1)
1141      {
1142         def lastR=RS[1];
1143         setring lastR;
1144         ideal KK=endid;
1145         ideal PP=fetch(BAS,ihp);
1146         export PP;
1147         export KK;
1148         setring BAS;
1149        // return(RS[1]);
1150         return(lastR);
1151      }
1152      int n=nvars(basering);
1153      ideal MJ=JM[2];
1154      def newR=RS[1];
1155      setring newR;
1156
1157      map psi=BAS,endphi;
1158      list tluser=
1159             normalizationPrimes(endid,psi(ihp),simplify(psi(MJ)+endid,4));
1160            // normalizationPrimes(endid);
1161      setring BAS;
1162      return(tluser);
1163   }
1164   else
1165   {
1166      int equi=attrib(SM[2],"isEquidimensional");
1167      ideal new1=qAnn+SM[2];
1168      execute "ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
1169                      +ordstr(basering)+");";
1170      if(y==1)
1171      {
1172         "zero-divisor found";
1173      }
1174      ideal vid=fetch(BAS,new1);
1175      ideal ihp=fetch(BAS,ihp);
1176      attrib(vid,"isCohenMacaulay",0);
1177      attrib(vid,"isPrim",0);
1178      attrib(vid,"isIsolatedSingularity",0);
1179      attrib(vid,"isRegInCodim2",0);
1180      if(equi==1)
1181      {
1182         attrib(vid,"isEquidimensional",1);
1183      }
1184      else
1185      {
1186         attrib(vid,"isEquidimensional",0);
1187      }
1188      attrib(vid,"isCompleteIntersection",0);
1189
1190      keepresult1=normalizationPrimes(vid,ihp);
1191
1192      setring BAS;
1193      ideal new2=quotient(SM[2],Ann)+SM[2];
1194// evtl. qAnn nehmen
1195//      ideal new2=SL[1]+SM[2];
1196      execute "ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
1197                      +ordstr(basering)+");";
1198
1199      ideal vid=fetch(BAS,new2);
1200      ideal ihp=fetch(BAS,ihp);
1201      attrib(vid,"isCohenMacaulay",0);
1202      attrib(vid,"isPrim",0);
1203      attrib(vid,"isIsolatedSingularity",0);
1204      attrib(vid,"isRegInCodim2",0);
1205      if(equi==1)
1206      {
1207         attrib(vid,"isEquidimensional",1);
1208      }
1209      else
1210      {
1211         attrib(vid,"isEquidimensional",0);
1212      }
1213      attrib(vid,"isCompleteIntersection",0);
1214
1215      keepresult2=normalizationPrimes(vid,ihp);
1216
1217      setring BAS;
1218      for(lauf=1;lauf<=size(keepresult2);lauf++)
1219      {
1220         keepresult1=insert(keepresult1,keepresult2[lauf]);
1221      }
1222      return(keepresult1);
1223   }
1224}
1225example
1226{ "EXAMPLE:";echo = 2;
1227   LIB"normal.lib";
1228   //Huneke
1229ring qr=31991,(a,b,c,d,e),dp;
1230ideal i=
12315abcde-a5-b5-c5-d5-e5,
1232ab3c+bc3d+a3be+cd3e+ade3,
1233a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1234abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1235ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1236a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1237a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1238b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1239
1240list pr=normal(i);
1241def r1=pr[1];
1242setring r1;
1243KK;
1244}
1245
1246
1247
1248/////////////////////////////////////////////////////////////////////////////
1249/*
1250LIB"normal.lib";
1251int aa=timer;list pr=normal(i);timer-aa;
1252
1253
1254
1255//Huneke
1256ring qr=31991,(a,b,c,d,e),dp;
1257ideal i=
12585abcde-a5-b5-c5-d5-e5,
1259ab3c+bc3d+a3be+cd3e+ade3,
1260a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1261abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1262ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1263a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1264a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1265b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1266
1267
1268//Vasconcelos
1269ring r=32003,(x,y,z,w,t),dp;
1270ideal i=
1271x2+zw,
1272y3+xwt,
1273xw3+z3t+ywt2,
1274y2w4-xy2z2t-w3t3;
1275
1276//Theo1
1277ring r=32003,(x,y,z),wp(2,3,6);
1278ideal i=zy2-zx3-x6;
1279
1280//Theo2
1281ring r=32003,(x,y,z),wp(3,4,12);
1282ideal i=z*(y3-x4)+x8;
1283
1284//Theo2a
1285ring r=32003,(T(1..4)),wp(3,4,12,17);
1286ideal i=
1287T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
1288T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
1289T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
1290T(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;
1291
1292//Theo3
1293ring r=32003,(x,y,z),wp(3,5,15);
1294ideal i=z*(y3-x5)+x10;
1295
1296
1297//Theo4
1298ring r=32003,(x,y,z),dp;
1299ideal i=(x-y)*(x-z)*(y-z);
1300
1301//Theo5
1302ring r=32003,(x,y,z),wp(2,1,2);
1303ideal i=z3-xy4;
1304
1305//Theo6
1306ring r=32003,(x,y,z),dp;
1307ideal i=x2y2+x2z2+y2z2;
1308
1309ring r=32003,(a,b,c,d,e,f),dp;
1310ideal i=
1311bf,
1312af,
1313bd,
1314ad;
1315
1316//Beispiel, wo vorher Primaerzerlegung schneller
1317//ist CM
1318//Sturmfels
1319ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1320ideal i=
1321bv+su,
1322bw+tu,
1323sw+tv,
1324by+sx,
1325bz+tx,
1326sz+ty,
1327uy+vx,
1328uz+wx,
1329vz+wy,
1330bvz;
1331
1332//J S/Y
1333ring r=32003,(x,y,z,t),dp;
1334ideal i=
1335x2z+xzt,
1336xyz,
1337xy2-xyt,
1338x2y+xyt;
1339
1340//St_S/Y
1341ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1342ideal i=
1343wy-vz,
1344vx-uy,
1345tv-sw,
1346su-bv,
1347tuy-bvz;
1348
1349//dauert laenger
1350//Horrocks:
1351ring r=32003,(a,b,c,d,e,f),dp;
1352ideal i=
1353adef-16000be2f+16001cef2,
1354ad2f+8002bdef+8001cdf2,
1355abdf-16000b2ef+16001bcf2,
1356a2df+8002abef+8001acf2,
1357ad2e-8000bde2-7999cdef,
1358acde-16000bce2+16001c2ef,
1359a2de-8000abe2-7999acef,
1360acd2+8002bcde+8001c2df,
1361abd2-8000b2de-7999bcdf,
1362a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
1363abde-8000b2e2-acdf-16001bcef-8001c2f2,
1364abcd-16000b2ce+16001bc2f,
1365a2cd+8002abce+8001ac2f,
1366a2bd-8000ab2e-7999abcf,
1367ab3f-3bdf3,
1368a2b2f-2adf3-16000bef3+16001cf4,
1369a3bf+4aef3,
1370ac3e-10668cde3,
1371a2c2e+10667ade3+16001be4+5334ce3f,
1372a3ce+10669ae3f,
1373bc3d+8001cd3e,
1374ac3d+8000bc3e+16001cd2e2+8001c4f,
1375b2c2d+16001ad4+4000bd3e+12001cd3f,
1376b2c2e-10668bc3f-10667cd2ef,
1377abc2e-cde2f,
1378b3cd-8000bd3f,
1379b3ce-10668b2c2f-10667bd2ef,
1380abc2f-cdef2,
1381a2bce-16000be3f+16001ce2f2,
1382ab3d-8000b4e-8001b3cf+16000bd2f2,
1383ab2cf-bdef2,
1384a2bcf-16000be2f2+16001cef3,
1385a4d-8000a3be+8001a3cf-2ae2f2;
1386
1387
1388ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1389
1390ideal k=
1391wy-vz,
1392vx-uy,
1393tv-sw,
1394su-bv,
1395tuy-bvz;
1396ideal j=x2y2+x2z2+y2z2;
1397ideal i=mstd(intersect(j,k))[2];
1398
1399//22
1400ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1401ideal i=
1402wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
1403vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
1404tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
1405sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
1406tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;
1407
1408
1409//riemenschneider
1410//33
1411//normal+primary 3
1412//primary 9
1413//radical 1
1414//minAssPrimes 2
1415ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
1416ideal i=
1417xz,
1418vx,
1419ux,
1420su,
1421qu,
1422txy,
1423stx,
1424qtx,
1425uv2z-uwz,
1426uv3-uvw,
1427puv2-puw;
1428
1429
1430*/
Note: See TracBrowser for help on using the repository browser.