source: git/Singular/LIB/normal.lib @ 67bd4c

spielwiese
Last change on this file since 67bd4c was 67bd4c, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes/pfister/greuel: update primdec.lib, add normal.lib git-svn-id: file:///usr/local/Singular/svn/trunk@884 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 31.2 KB
Line 
1// $Id: normal.lib,v 1.1 1997-11-05 18:16:55 Singular Exp $
2///////////////////////////////////////////////////////////////////////////////
3// normal.lib
4// algorithms for computing the normalization based on
5// the ideas of De Jong,Vasconcelos
6// written by Gert-Martin Greuel and Gerhard Pfister
7///////////////////////////////////////////////////////////////////////////////
8
9LIBRARY: normal.lib: PROCEDURE FOR NORMALIZATION (I)
10
11  normal(ideal I)
12  // computes a set of rings such that their product is the
13  // normalization of the reduced basering/I
14
15LIB "sing.lib";
16LIB "primdec.lib";
17LIB "elim.lib";
18LIB "presolve.lib";
19///////////////////////////////////////////////////////////////////////////////
20
21proc isR_HomJR (list Li)
22USAGE:   isR_HomJR (Li);  Li = list: ideal SBid, ideal J, poly p
23COMPUTE: module Hom_R(J,R) = R:J and compare with R
24ASSUME:  R    = P/SBid,  P = basering
25         SBid = standard basis of an ideal in P,
26         J    = ideal in P containing the polynomial p,
27         p    = nonzero divisor of R
28RETURN:  1 if R = R:J, 0 if not
29EXAMPLE: example isR_HomJR;  shows an example
30{
31   int n, ii;
32 def P = basering;
33   ideal SBid = Li[1];
34   ideal J = Li[2];
35   poly p = Li[3];
36   attrib(SBid,"isSB",1);
37   attrib(p,"isSB",1);
38 qring R    = SBid;
39   ideal J  = fetch(P,J);
40   poly p   = fetch(P,p);
41   ideal f  = quotient(p,J);
42   ideal lp = std(p);
43   n=1;
44   for (ii=1; ii<=size(f); ii++ )
45   {
46      if ( reduce(f[ii],lp) != 0)
47      { n = 0; break; }
48   }
49   return (n);
50 //?spaeter hier einen Test ob Hom(I,R) = Hom(I,I)?
51}
52example
53{"EXAMPLE:";  echo = 2;
54  ring r   = 0,(x,y,z),dp;
55  ideal id = y7-x5+z2;
56  ideal J  = x3,y+z;
57  poly p   = xy;
58  list Li  = std(id),J,p;
59  isR_HomJR (Li);
60
61  ring s   = 0,(t,x,y),dp;
62  ideal id = x2-y2*(y-t);
63  ideal J  = jacob(id);
64  poly p   = J[1];
65  list Li  = std(id),J,p;
66  isR_HomJR (Li);
67}
68///////////////////////////////////////////////////////////////////////////////
69
70proc extraweight (list # )
71USAGE:   extraweight (P); P=name of an existing ring (true name, not a string)
72RETURN:  intvec, size=nvars(P), consisting of the weights of the variables of P
73NOTE:    This is useful when enlarging P but keeping the weights of the old
74         variables
75EXAMPLE: example extraweight;  shows an example
76{
77   int ii,q,fi,fo,fia;
78   intvec rw,nw;
79   string os;
80   def P = #[1];
81   string osP = ordstr(P);
82   fo  = 1;
83//------------------------- find weights in ordstr(P) -------------------------
84   fi  = find(osP,"(",fo);
85   fia = find(osP,"a",fo)+find(osP,"w",fo)+find(osP,"W",fo);
86   while ( fia )
87   {
88      os = osP[fi+1,find(osP,")",fi)-fi-1];
89      if( find(os,",") )
90      {
91         execute "nw = "+os+";";
92         if( size(nw) > ii )
93         {
94             rw = rw,nw[ii+1..size(nw)];
95         }
96         else  {  ii = ii - size(nw); }
97
98         if( find(osP[1,fi],"a",fo) ) { ii = size(nw); }
99      }
100      else
101      {
102         execute "q = "+os+";";
103         if( q > ii )
104         {
105            nw = 0; nw[q-ii] = 0;
106            nw = nw + 1;          //creates an intvec 1,...,1 of length q-ii
107            rw = rw,nw;
108         }
109         else { ii = ii - q; }
110      }
111      fo  = fi+1;
112      fi  = find(osP,"(",fo);
113      fia = find(osP,"a",fo)+find(osP,"w",fo)+find(osP,"W",fo);
114   }
115//-------------- adjust weight vector to length = nvars(P)  -------------------
116   if( fo > 1 )
117   {                                            // case when weights were found
118      rw = rw[2..size(rw)];
119      if( size(rw) > nvars(P) )
120      {
121         rw = rw[1..nvars(P)];
122      }
123      if( size(rw) < nvars(P) )
124      {
125         nw=0; nw[nvars(P)-size(rw)]=0; nw=nw+1; rw=rw,nw;
126      }
127   }
128   else
129   {                                         // case when no weights were found
130      rw[nvars(P)]= 0; rw=rw+1;
131   }
132   return(rw);
133}
134example
135{"EXAMPLE:";  echo = 2;
136  ring r0 = 0,(x,y,z),dp;
137  extraweight(r0);
138  ring r1 = 0,x(1..5),(ds(3),wp(2,3));
139  extraweight(r1);
140  ring r2 = 0,x(1..5),(a(1,2,3,0),dp);
141  extraweight(r2);
142  ring r3 = 0,x(1..10),(a(1..5),dp(5),a(10..13),Wp(5..9));
143  extraweight(r3);
144  // an example for enlarging the ring:
145  intvec v = 6,2,3,4,5;
146  ring R = 0,x(1..10),(a(extraweight(r1),v),dp);
147  ordstr(R);
148}
149///////////////////////////////////////////////////////////////////////////////
150
151proc HomJJ (list Li,list #)
152USAGE:   HomJJ (Li);  Li = list: SBid,id,J,poly p
153ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
154         SBid = standard basis of id,
155         J    = ideal of P containing the polynomial p,
156         p    = nonzero divisor of R
157COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure where
158         R is the quotient ring of P modulo the standard basis SBid
159RETURN:  a list of two objects
160         _[1]: a polynomial ring, containing two ideals, 'endid' and 'endphi'
161               s.t. _[1]/endid = Hom_R(J,J) and
162               endphi describes the canonical map R -> Hom_R(J,J)
163         _[2]: an integer which is 1 if phi is an isomorphism, 0 if not
164EXAMPLE: example HomJJ;  shows an example
165{
166//---------- initialisation ---------------------------------------------------
167
168   int isIso,isPr,isCo,isRe,isEq,ii,jj,q,y;
169   intvec rw,rw1;
170   list L;
171   if( size(#) >=1 )
172   {
173      if( typeof(#[1]) == "int" ) { y = #[1]; }
174   }
175 def P = basering;
176   ideal SBid, id, J = Li[1], Li[2], Li[3];
177   poly p = Li[4];
178   attrib(SBid,"isSB",1);
179   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not
180
181//---- set attributes for special cases where algorithm can be simplified -----
182   if( homo==1 )
183   {
184      rw = extraweight(P);
185   }
186   if( typeof(attrib(id,"isPrim"))=="int" )
187   {
188      if(attrib(id,"isPrim")==1)  { isPr=1; }
189   }
190   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
191   {
192      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
193   }
194   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
195   {
196      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
197   }
198   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
199   {
200      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
201   }
202   if( typeof(attrib(id,"isEquidimensional"))=="int" )
203   {
204      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
205   }
206//-------------------------- go to quotient ring ------------------------------
207 qring R  = SBid;
208   ideal id = fetch(P,id);
209   ideal J  = fetch(P,J);
210   poly p   = fetch(P,p);
211   ideal f,rf,f2;
212   module syzf;
213
214//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
215   f  = quotient(p*J,J);
216   if(y==1)
217   {
218      "the module Hom(rad(J),rad(J)) presented by the values on";
219      "the non-zerodivisor";
220      "  ";
221      p;
222      " ";
223      f;
224   }
225   f2 = std(p);
226
227   if(isIso==0)
228   {
229     ideal f1=std(f);
230     attrib(f1,"isSB",1);
231    // if( codim(f1,f2) >= 0 )
232    // {
233    //  dbprint(printlevel-voice+3,"// dimension of non-normal locus is zero");
234    //    isIso=1;
235    // }
236  }
237//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
238   rf = interred(reduce(f,f2));       // represents p*Hom(J,J)/p*R = Hom(J,J)/R
239   if ( size(rf) == 0 )
240   {
241      if ( homog(f) && find(ordstr(basering),"s")==0 )
242      {
243         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
244      }
245      else
246      {
247         ring newR1 = char(P),(X(1..nvars(P))),dp;
248      }
249      ideal endphi = maxideal(1);
250      ideal endid = fetch(P,id);
251      attrib(endid,"isCohenMacaulay",isCo);
252      attrib(endid,"isPrim",isPr);
253      attrib(endid,"isIsolatedSingularity",isIso);
254      attrib(endid,"isRegInCodim2",isRe);
255      attrib(endid,"isEqudimensional",isEq);
256      attrib(endid,"isCompleteIntersection",0);
257      attrib(endid,"isRad",0);
258      export endid;
259      export endphi;
260      L = newR1;
261      L = insert(L,1,1);
262      dbprint(printlevel-voice+3,"// case R = Hom(J,J)");
263      if(y==1)
264      {
265         "R=Hom(rad(J),rad(J))";
266         "   ";
267         newR1;
268         "   ";
269         "the new ideal";
270         endid;
271         "   ";
272         "the old ring";
273         "   ";
274         P;
275         "   ";
276         "the old ideal";
277         "   ";
278         setring P;
279         id;
280         "   ";
281         setring newR1;
282         "the map";
283         "   ";
284         endphi;
285         "   ";
286         pause;
287      }
288      setring P;
289      return(L);
290   }
291   if(y==1)
292   {
293      "R is not equal to Hom(rad(J),rad(J)), we have to try again";
294      pause;
295   }
296//---------- Hom(J,j) != R: create new ring and map form old ring -------------
297// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
298
299   f = mstd(f)[2];                                   // minimizes f (partially)
300   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
301   q = size(f);
302   syzf = syz(f);
303
304   if ( homo==1 )
305   {
306      rw1 = rw,0;
307      for ( ii=2; ii<=q; ii++ )
308      {
309         rw  = rw, deg(f[ii])-deg(f[1]);
310         rw1 = rw1, deg(f[ii])-deg(f[1]);
311      }
312      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
313   }
314   else
315   {
316      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
317   }
318
319   map psi1 = P,maxideal(1);
320   ideal SBid = psi1(SBid);
321   attrib(SBid,"isSB",1);
322
323 qring newR = std(SBid);
324   map psi = R,ideal(X(1..nvars(R)));
325   ideal id = psi(id);
326   ideal f = psi(f);
327   module syzf = psi(syzf);
328   ideal pf,Lin,Quad,Q;
329   matrix T,A;
330   list L1;
331
332//---------- computation of Hom(J,J) as ring ----------------------------------
333// determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
334// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course
335// the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
336// It is a fact, that the kernel is generated by the linear and the quadratic
337// relations
338
339   pf = f[1]*f;
340   T = matrix(ideal(T(1..q)),1,q);
341   Lin = ideal(T*syzf);
342   if(y==1)
343   {
344      "the ring structure of Hom(rad(J),rad(J)) as R-algebra";
345      " ";
346      "the linear relations";
347      " ";
348      Lin;
349      "   ";
350   }
351   for (ii=2; ii<=q; ii++ )
352   {
353      for ( jj=2; jj<=ii; jj++ )
354      {
355         A = lift(pf,f[ii]*f[jj]);
356         Quad = Quad, ideal(T(jj)*T(ii) - T*A);          // quadratic relations
357      }
358   }
359   if(y==1)
360   {
361      "the quadratic relations";
362      "   ";
363      interred(Quad);
364      pause;
365   }
366   Q = Lin+Quad;
367   Q = subst(Q,T(1),1);
368   Q = interred(reduce(Q,std(0)));
369//---------- reduce number of variables by substitution, if possible ----------
370   if (homo==1)
371   {
372      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
373   }
374   else
375   {
376      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
377   }
378
379   ideal endid  = imap(newR,id)+imap(newR,Q);
380   ideal endphi = ideal(X(1..nvars(R)));
381
382   map phi = basering,maxideal(1);
383   list Le = elimpart(endid);
384           //this proc and the next loop try to
385   q = size(Le[2]);                 //substitute as many variables as possible
386   rw1 = 0;
387   rw1[nvars(basering)] = 0;
388   rw1 = rw1+1;
389
390   while( size(Le[2]) != 0 )
391   {
392      endid = Le[1];
393      map ps = newRing,Le[5];
394      phi = phi(ps);
395      kill ps;
396
397      for( ii=1; ii<=size(rw1); ii++ )
398      {
399         if( Le[4][ii]==0 )
400         {
401            rw1[ii]=0;                             //look for substituted vars
402         }
403      }
404      Le=elimpart(endid);
405      q = q + size(Le[2]);
406   }   
407   endphi = phi(endphi);
408
409//---------- return -----------------------------------------------------------
410// in the homogeneous case put weights for the remaining vars correctly, i.e.
411// delete from rw those weights for which the corresponding entry of rw1 is 0
412
413   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
414   {
415      jj=1;
416      for( ii=2; ii<size(rw1); ii++)
417      {
418         jj++;
419         if( rw1[ii]==0 )
420         {
421            rw=rw[1..jj-1],rw[jj+1..size(rw)];
422            jj=jj-1;
423         }
424      }
425      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
426      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }
427
428      ring lastRing = char(R),(T(1..nvars(newRing)-q)),(a(rw),dp);
429   }
430   else
431   {
432      ring lastRing = char(R),(T(1..nvars(newRing)-q)),dp;
433   }
434
435   ideal lastmap;
436   q = 1;
437   for(ii=1; ii<=size(rw1); ii++ )
438   {
439      if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; }
440      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
441   }
442   map phi = newRing,lastmap;
443   ideal endid  = phi(endid);
444   ideal endphi = phi(endphi);
445   attrib(endid,"isCohenMacaulay",isCo);
446   attrib(endid,"isPrim",isPr);
447   attrib(endid,"isIsolatedSingularity",isIso);
448   attrib(endid,"isRegInCodim2",isRe);
449   attrib(endid,"isEquidimensional",isEq);
450   attrib(endid,"isCompleteIntersection",0);
451   attrib(endid,"isRad",0);
452   export(endid);
453   export(endphi);
454   if(y==1)
455   {
456      "the new ring after reduction of the number of variables";
457      "   ";
458      lastRing;
459      "   ";
460      "the new ideal";
461      endid;
462      "   ";
463      "the old ring";
464      "   ";
465      P;
466      "   ";
467      "the old ideal";
468      "   ";
469      setring P;
470      id;
471      "   ";
472      setring lastRing;
473      "the map";
474      "   ";
475      endphi;
476      "   ";
477      pause;
478   }
479   L = lastRing;
480   L = insert(L,0,1);
481   return(L);
482}
483example
484{"EXAMPLE:";  echo = 2;
485  ring r   = 0,(x,y),wp(2,3);
486  ideal id = y^2-x^3;
487  ideal J  = x,y;
488  poly p   = x;
489  list Li = std(id),id,J,p;
490  list L   = HomJJ(Li);
491  def end = L[1];      // defines ring L[1], containing ideals endid and endphi
492  setring end;         // makes end the basering
493  end;
494  endid;               // end/endid is isomorphic to End(r/id) as ring
495  map psi = r,endphi;  // defines the canonical map r/id -> End(r/id)
496  psi;
497
498  ring r   = 32003,(x,y,z),dp;
499  ideal id = x2-xy-xz+yz;
500  ideal J =y-z,x-z;
501  poly p = x+16001y+16001z;
502  list Li = std(id),id,J,p;
503  list L   = HomJJ(Li,0);
504  def end = L[1];      // defines ring L[1], containing ideals endid and endphi
505  setring end;         // makes end the basering
506  end;
507  endid;               // end/endid is isomorphic to End(r/id) as ring
508
509}
510
511///////////////////////////////////////////////////////////////////////////////
512proc normal(ideal id, list #)
513USAGE:   normal(i,choose);  i ideal,choose empty or 1
514         in case you choose 1 the factorizing Buchberger algorithm
515         is not used which is sometimes more efficient
516RETURN:  a list of rings L=R1,...,Rk, in each Ri are two ideals
517         Si,Mi such that the product of Ri/Mi is the normalization
518         Si is a standardbasis of Mi
519NOTE:    to use the rings: def r=L[i];setring r;
520EXAMPLE: example normal; shows an example
521{
522   int i,j,y;
523   list result,prim,keepresult;
524
525   if(size(#)==0)
526   {
527      prim[1]=id;
528      if( typeof(attrib(id,"isEquidimensional"))=="int" )
529      {
530        if(attrib(id,"isEquidimensional")==1)
531        {
532           attrib(prim[1],"isEquidimensional",1);
533        }
534      }
535      else
536      {
537         attrib(prim[1],"isEquidimensional",0);
538      }
539      if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
540      {
541        if(attrib(id,"isCompleteIntersection")==1)
542        {
543           attrib(prim[1],"isCompleteIntersection",1);
544        }
545      }
546      else
547      {
548         attrib(prim[1],"isCompleteIntersection",0);
549      }
550
551      if( typeof(attrib(id,"isPrim"))=="int" )
552      {
553        if(attrib(id,"isPrim")==1)  {attrib(prim[1],"isPrim",1); }
554      }
555      else
556      {
557         attrib(prim[1],"isPrim",0);
558      }
559      if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
560      {
561         if(attrib(id,"isIsolatedSingularity")==1)
562             {attrib(prim[1],"isIsolatedSingularity",1); }
563      }
564      else
565      {
566         attrib(prim[1],"isIsolatedSingularity",0);
567      }
568      if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
569      {
570         if(attrib(id,"isCohenMacaulay")==1)
571           { attrib(prim[1],"isCohenMacaulay",1); }
572      }
573      else
574      {
575         attrib(prim[1],"isCohenMacaulay",0);
576      }
577      if( typeof(attrib(id,"isRegInCodim2"))=="int" )
578      {
579         if(attrib(id,"isRegInCodim2")==1)
580         { attrib(prim[1],"isRegInCodim2",1);}
581      }
582      else
583      {
584          attrib(prim[1],"isRegInCodim2",0);
585      }
586      return(normalizationPrimes(prim[1],maxideal(1)));
587   }
588   else
589   {
590      if(#[1]==0)
591      {
592         prim=minAssPrimes(id);
593      }
594      else
595      {
596         prim=minAssPrimes(id,1);
597      }
598
599      if(y==1)
600      {
601         "we have ";size(prim);"components";
602      }
603      for(i=1;i<=size(prim);i++)
604      {
605         if(y==1)
606         {
607            "we are in loop";i;
608         }
609         attrib(prim[i],"isCohenMacaulay",0);
610         attrib(prim[i],"isPrim",1);
611         attrib(prim[i],"isRegInCodim2",0);
612         attrib(prim[i],"isIsolatedSingularity",0);
613         attrib(prim[i],"isEquidimensional",0);
614         attrib(prim[i],"isCompleteIntersection",0);
615
616         if( typeof(attrib(id,"isEquidimensional"))=="int" )
617         {
618           if(attrib(id,"isEquidimensional")==1)
619           {
620              attrib(prim[i],"isEquidimensional",1);
621           }
622         }
623         else
624         {
625            attrib(prim[i],"isEquidimensional",0);
626         }
627         if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
628         {
629            if(attrib(id,"isIsolatedSingularity")==1)
630             {attrib(prim[i],"isIsolatedSingularity",1); }
631         }
632         else
633         {
634            attrib(prim[i],"isIsolatedSingularity",0);
635         }
636
637         keepresult=normalizationPrimes(prim[i],maxideal(1));
638         for(j=1;j<=size(keepresult);j++)
639         {
640            result=insert(result,keepresult[j]);
641         }
642      }
643      return(result);
644   }
645}
646example
647{ "EXAMPLE:"; echo = 2;
648   ring  r = 0,(x,y,z),dp;
649
650}
651
652
653proc normalizationPrimes(ideal i,ideal ihp, list #)
654USAGE:   normalizationPrimes(i);  i prime ideal
655RETURN:  a list of one ring L=R, in  R are two ideals
656         S,M such that R/M is the normalization
657         S is a standardbasis of M
658NOTE:    to use the ring: def r=L[1];setring r;
659EXAMPLE: example normalizationPrimes; shows an example
660{
661   int y;
662   if(y==1)
663   {
664      "START one normalization loop with the ideal";
665      "        ";
666      i;
667      "        ";
668      basering;
669      "        ";
670      pause;
671   }
672
673   def BAS=basering;
674   list result,keepresult1,keepresult2;
675   ideal J,SB,MB,KK;
676   int depth,lauf,prdim;
677   int ti=timer;
678
679   if(size(i)==0)
680   {
681      if(y==1)
682      {
683          "the ideal was the zero-ideal";
684      }
685         execute "ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
686                      +ordstr(basering)+");";
687         ideal KK=ideal(0);
688         ideal PP=fetch(BAS,ihp);
689         export PP;
690         export KK;
691         result=newR7;
692         setring BAS;
693         return(result);
694
695   }
696
697   if(y==1)
698   {
699      "                                  ";
700      "                                  ";
701      " SB-computation of the input ideal";
702      "                                  ";
703   }
704   list SM=mstd(i);
705
706   if(y==1)
707   {
708      " the dimension is                 ";
709      dim(SM[1]);
710      "                                  ";
711   }
712
713   if(size(#)>0)
714   {
715      list JM=#[1],#[1];
716      attrib(JM[1],"isSB",1);
717      if( typeof(attrib(#[1],"isRad"))!="int" )
718      {
719         attrib(JM[2],"isRad",0);
720      }
721   }
722
723   if(attrib(i,"isPrim")==1)
724   {
725      attrib(SM[2],"isPrim",1);
726   }
727   else
728   {
729      attrib(SM[2],"isPrim",0);
730   }
731   if(attrib(i,"isIsolatedSingularity")==1)
732   {
733      attrib(SM[2],"isIsolatedSingularity",1);
734   }
735   else
736   {
737      attrib(SM[2],"isIsolatedSingularity",0);
738   }
739   if(attrib(i,"isCohenMacaulay")==1)
740   {
741      attrib(SM[2],"isCohenMacaulay",1);
742   }
743   else
744   {
745      attrib(SM[2],"isCohenMacaulay",0);
746   }
747   if(attrib(i,"isRegInCodim2")==1)
748   {
749      attrib(SM[2],"isRegInCodim2",1);
750   }
751   else
752   {
753      attrib(SM[2],"isRegInCodim2",0);
754   }
755    if(attrib(i,"isEquidimensional")==1)
756   {
757      attrib(SM[2],"isEquidimensional",1);
758   }
759   else
760   {
761      attrib(SM[2],"isEquidimensional",0);
762   }
763    if(attrib(i,"isCompleteIntersection")==1)
764   {
765      attrib(SM[2],"isCompleteIntersection",1);
766   }
767   else
768   {
769      attrib(SM[2],"isCompleteIntersection",0);
770   }
771
772   //the smooth case
773   if(size(#)>0)
774   {
775      if(dim(JM[1])==-1)
776      {
777         if(y==1)
778         {
779            "the ideal was smooth";
780         }
781         MB=SM[2];
782         execute "ring newR6="+charstr(basering)+",("+varstr(basering)+"),("
783                      +ordstr(basering)+");";
784         ideal KK=fetch(BAS,MB);
785         ideal PP=fetch(BAS,ihp);
786         export PP;
787         export KK;
788         result=newR6;
789         setring BAS;
790         return(result);
791     }
792   }
793
794   //the zero-dimensional case
795   if((dim(SM[1])==0)&&(homog(SM[2])==1))
796   {
797      if(y==1)
798      {
799         "the ideal was zero-dimensional and homogeneous";
800      }
801      MB=maxideal(1);
802      execute "ring newR5="+charstr(basering)+",("+varstr(basering)+"),("
803                      +ordstr(basering)+");";
804      ideal KK=fetch(BAS,MB);
805      ideal PP=fetch(BAS,ihp);
806      export PP;
807      export KK;
808      result=newR5;
809      setring BAS;
810      return(result);
811   }
812
813   //the one-dimensional case
814   //in this case it is a line because
815   //it is irreducible and homogeneous
816   if((dim(SM[1])==1)&&(attrib(SM[2],"isPrim")==1)
817        &&(homog(SM[2])==1))
818   {
819      if(y==1)
820      {
821         "the ideal defines a line";
822      }
823      MB=SM[2];
824      execute "ring newR4="+charstr(basering)+",("+varstr(basering)+"),("
825                      +ordstr(basering)+");";
826      ideal KK=fetch(BAS,MB);
827      ideal PP=fetch(BAS,ihp);
828      export PP;
829      export KK;
830      result=newR4;
831      setring BAS;
832      return(result);
833   }
834
835   //the higher dimensional case
836   //we test first of all CohenMacaulay and
837   //complete intersection
838   if(((size(SM[2])+dim(SM[1]))==nvars(basering))&&(homog(SM[2])==1))
839   {
840      //test for complete intersection
841      attrib(SM[2],"isCohenMacaulay",1);
842      attrib(SM[2],"isCompleteIntersection",1);
843      attrib(SM[2],"isEquidimensional",1);
844      if(y==1)
845      {
846         "it is a complete Intersection";
847      }
848   }
849//   if(attrib(i,"isIsolatedSingularity")==0)
850//   {
851//      list L=sres(SM[2],0);
852//      prdim=ncols(betti(L))-1;
853//      depth=nvars(basering)-prdim;
854//   }
855
856  // if(attrib(SM[2],"isCohenMacaulay")==0)
857  // {
858  //    test for CohenMacaulay
859  //    if((dim(SM[1]))==depth)
860  //    {
861  //    attrib(SM[2],"isCohenMacaulay",1);
862  //    "it is CohenMacaulay";
863  //    }
864  // }
865
866   //compute the singular locus+lower dimensional components
867   if(((attrib(SM[2],"isIsolatedSingularity")==0)||(homog(SM[2])==0))
868        &&(size(#)==0))
869   {
870
871      J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
872      //ti=timer;
873      if(y==1)
874      {
875         "SB of singular locus will be computed";
876      }
877      ideal sin=J+SM[2];
878
879    //kills the embeded components
880//    if(homog(SM[2])==1)
881//    {
882//       sin=sat(sin,maxideal(1))[1];
883//    }
884
885      list JM=mstd(sin);
886      if(y==1)
887      {
888         "                                  ";
889         "the dimension of the singular locus is";
890         dim(JM[1]);
891         "                                  ";
892      }
893
894      attrib(JM[2],"isRad",0);
895      //   timer-ti;
896      attrib(JM[1],"isSB",1);
897      if(dim(JM[1])==-1)
898      {
899         if(y==1)
900         {
901            "it is smooth";
902         }
903         MB=SM[2];
904         execute "ring newR3="+charstr(basering)+",("+varstr(basering)+"),("
905                      +ordstr(basering)+");";
906         ideal KK=fetch(BAS,MB);
907         ideal PP=fetch(BAS,ihp);
908         export PP;
909         export KK;
910         result=newR3;
911         setring BAS;
912         return(result);
913      }
914      if(dim(JM[1])==0)
915      {
916         attrib(SM[2],"isIsolatedSingularity",1);
917      }
918      if(dim(JM[1])<=nvars(basering)-2)
919      {
920         attrib(SM[2],"isRegInCodim2",1);
921      }
922   }
923   else
924   {
925     if(size(#)==0)
926     {
927        list JM=maxideal(1),maxideal(1);
928        attrib(JM[1],"isSB",1);
929        attrib(SM[2],"isRegInCodim2",1);
930     }
931   }
932   //if it is an isolated singularity things are easier
933   if((dim(JM[1])==0)&&(homog(SM[2])==1))
934   {
935      attrib(SM[2],"isIsolatedSingularity",1);
936      ideal SL=simplify(reduce(maxideal(1),SM[1]),2);
937      ideal Ann=quotient(SM[2],SL[1]);
938      ideal qAnn=simplify(reduce(Ann,SM[1]),2);
939
940      if(size(qAnn)==0)
941      {
942         if(y==1)
943         {
944            "                                  ";
945            "the ideal rad(J):                 ";
946            "                                  ";
947            maxideal(1);
948            "                                  ";
949            "                                  ";
950         }
951         //again test for normality
952         //Hom(I,R)=R
953         list RR;
954  //       list RR=SM[1],maxideal(1),SL[1];
955  //       "test with Hom(I,R)";
956  //       ti=timer;
957  //       if(isR_HomJR(RR)==1)
958  //       {
959  //         "it was normal";
960  //          timer-ti;
961  //          SB=SM[1];
962  //          SM=SM[2];
963  //          export SB,MB;
964  //          result=BAS;
965  //          return(result);
966  //       }
967  //       timer-ti;
968         RR=SM[1],SM[2],maxideal(1),SL[1];
969         ti=timer;
970         RR=HomJJ(RR,y);
971   //      timer-ti;
972         if(RR[2]==0)
973         {
974            def newR=RR[1];
975            setring newR;
976            map psi=BAS,endphi;
977         //   ti=timer;
978            list tluser=normalizationPrimes(endid,psi(ihp));
979
980        //    timer-ti;
981            setring BAS;
982            return(tluser);
983         }
984         MB=SM[2];
985         execute "ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
986                      +ordstr(basering)+");";
987         ideal KK=fetch(BAS,MB);
988         ideal PP=fetch(BAS,ihp);
989         export PP;
990         export KK;
991         result=newR7;
992         setring BAS;
993         return(result);
994
995       }
996       else
997       {       
998          ideal id=qAnn+SM[2];
999
1000          attrib(id,"isCohenMacaulay",0);
1001          attrib(id,"isPrim",0);
1002          attrib(id,"isIsolatedSingularity",1);
1003          attrib(id,"isRegInCodim2",0);
1004          attrib(id,"isCompleteIntersection",0);
1005          attrib(id,"isEquidimensional",0);
1006
1007          keepresult1=normalizationPrimes(id,ihp);
1008          ideal id1=quotient(SM[2],Ann)+SM[2];
1009//          evtl. qAnn nehmen
1010//          ideal id=SL[1]+SM[2];
1011
1012          attrib(id1,"isCohenMacaulay",0);
1013          attrib(id1,"isPrim",0);
1014          attrib(id1,"isIsolatedSingularity",1);
1015          attrib(id1,"isRegInCodim2",0);
1016          attrib(id1,"isCompleteIntersection",0);
1017          attrib(id1,"isEquidimensional",0);
1018
1019          keepresult2=normalizationPrimes(id1,ihp);
1020
1021          for(lauf=1;lauf<=size(keepresult2);lauf++)
1022          {
1023             keepresult1=insert(keepresult1,keepresult2[lauf]);
1024          }
1025          return(keepresult1);
1026       }
1027   }
1028
1029   //test for non-normality
1030   //Hom(I,I)<>R
1031   //we can use Hom(I,I) to continue
1032
1033   ideal SL=simplify(reduce(JM[2],SM[1]),2);
1034   ideal Ann=quotient(SM[2],SL[1]);
1035   ideal qAnn=simplify(reduce(Ann,SM[1]),2);
1036
1037   if(size(qAnn)==0)
1038   {
1039      list RR;
1040 //      "test with Hom(I,I) before the radical computation";
1041        //again test for normality
1042        //Hom(I,R)=R
1043 //     list RR=SM[1],JM[2],SL[1];
1044 //     "test with Hom(I,R)";
1045 //     ti=timer;
1046 //     if(isR_HomJR(RR)==1)
1047 //     {
1048 //      "it was normal";
1049 //        timer-ti;
1050 //        SB=SM[1];
1051 //        SM=SM[2];
1052 //        export SB,MB;
1053 //        result=BAS;
1054 //        return(result);
1055 //     }
1056 //     timer-ti;
1057
1058 //     list  RR=SM[1],SM[2],JM[2],SL[1];
1059 //     ti=timer;
1060      list RS;
1061 //   list RS=HomJJ(RR);
1062 //   timer-ti;
1063 //     if(RS[2]==0)
1064 //     {
1065 //        "Hom(I,I) was sucessfull without radical";
1066 //        def newR=RS[1];
1067 //        setring newR;
1068 //        list tluser=normalizationPrimes(SM);
1069 //        setring BAS;
1070 //        return(tluser);
1071 //     }
1072
1073      //now we have to compute the radical
1074      if(y==1)
1075      {
1076         "radical computation";
1077      }
1078//      ti=timer;
1079
1080
1081      if((attrib(JM[2],"isRad")==0)&&(attrib(SM[2],"isEquidimensional")==0))
1082      {
1083           //J=radical(JM[2]);
1084          J=radical(SM[2]+ideal(SL[1]));
1085
1086          // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
1087      }
1088      if((attrib(JM[2],"isRad")==0)&&(attrib(SM[2],"isEquidimensional")==1))
1089      {
1090          ideal JJ=SM[2]+ideal(SL[1]);
1091         // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal
1092          if(attrib(SM[2],"isCompleteIntersection")==0)
1093          {
1094             J=equiRadical(JM[2]);
1095             //J=equiRadical(JJ);
1096          }
1097          else
1098          {
1099             //J=radical(JM[2]);
1100             J=quotient(JJ,minor(jacob(JJ),size(JJ)));
1101          }
1102      }
1103//    timer-ti;
1104
1105      JM=J,J;
1106
1107      //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen
1108      // SL=simplify(reduce(J,SM[1]),2);
1109      // Ann=quotient(SM[2],SL[1]);
1110      // qAnn=simplify(reduce(Ann,SM[1]),2);
1111      // if(size(qAnn)!=0)
1112      // {
1113        // keepresult1=normalizationPrimes(qAnn+SM[2]);
1114        // keepresult2=normalizationPrimes(SL[1]+SM[2]);
1115        // for(lauf=1;lauf<=size(keepresult2);j++)
1116        // {
1117        //    keepresult1=insert(keepresult1,keepresult2[lauf]);
1118        // }
1119        // return(keepresult1);
1120      // }
1121      RR=SM[1],SM[2],JM[2],SL[1];
1122
1123//   evtl eine geeignete Potenz von JM?
1124     if(y==1)
1125     {
1126        "compute Hom(rad(J),rad(J)):";
1127     }
1128//   ti=timer;
1129
1130     RS=HomJJ(RR,y);
1131//   timer-ti;
1132
1133      if(RS[2]==1)
1134      {
1135         def lastR=RS[1];
1136         setring lastR;
1137         ideal KK=endid;
1138         ideal PP=fetch(BAS,ihp);
1139         export PP;
1140         export KK;
1141         setring BAS;
1142        // return(RS[1]);
1143         return(lastR);
1144      }
1145      int n=nvars(basering);
1146      ideal MJ=JM[2];
1147      def newR=RS[1];
1148      setring newR;
1149
1150      map psi=BAS,endphi;
1151      list tluser=
1152             normalizationPrimes(endid,psi(ihp),simplify(psi(MJ)+endid,4));
1153            // normalizationPrimes(endid);
1154      setring BAS;
1155      return(tluser);
1156   }
1157   else
1158   {
1159      int equi=attrib(SM[2],"isEquidimensional");
1160      ideal new1=qAnn+SM[2];
1161      execute "ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
1162                      +ordstr(basering)+");";
1163      if(y==1)
1164      {
1165         "zero-divisor found";
1166      }
1167      ideal vid=fetch(BAS,new1);
1168      ideal ihp=fetch(BAS,ihp);
1169      attrib(vid,"isCohenMacaulay",0);
1170      attrib(vid,"isPrim",0);
1171      attrib(vid,"isIsolatedSingularity",0);
1172      attrib(vid,"isRegInCodim2",0);
1173      if(equi==1)
1174      {
1175         attrib(vid,"isEquidimensional",1);
1176      }
1177      else
1178      {
1179         attrib(vid,"isEquidimensional",0);
1180      }
1181      attrib(vid,"isCompleteIntersection",0);
1182
1183      keepresult1=normalizationPrimes(vid,ihp);
1184
1185      setring BAS;
1186      ideal new2=quotient(SM[2],Ann)+SM[2];
1187// evtl. qAnn nehmen
1188//      ideal new2=SL[1]+SM[2];
1189      execute "ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
1190                      +ordstr(basering)+");";
1191
1192      ideal vid=fetch(BAS,new2);
1193      ideal ihp=fetch(BAS,ihp);
1194      attrib(vid,"isCohenMacaulay",0);
1195      attrib(vid,"isPrim",0);
1196      attrib(vid,"isIsolatedSingularity",0);
1197      attrib(vid,"isRegInCodim2",0);
1198      if(equi==1)
1199      {
1200         attrib(vid,"isEquidimensional",1);
1201      }
1202      else
1203      {
1204         attrib(vid,"isEquidimensional",0);
1205      }
1206      attrib(vid,"isCompleteIntersection",0);
1207
1208      keepresult2=normalizationPrimes(vid,ihp);
1209
1210      setring BAS;
1211      for(lauf=1;lauf<=size(keepresult2);lauf++)
1212      {
1213         keepresult1=insert(keepresult1,keepresult2[lauf]);
1214      }
1215      return(keepresult1);
1216   }
1217}
1218example
1219{ "EXAMPLE:";echo = 2;
1220   //Huneke
1221  ring qr=31991,(a,b,c,d,e),dp;
1222  ideal i=
1223  5abcde-a5-b5-c5-d5-e5,
1224  ab3c+bc3d+a3be+cd3e+ade3,
1225  a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1226  abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1227  ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1228  a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1229  a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1230  b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1231 
1232  list pr=normal(i);
1233  def r1=pr[1];
1234  setring r1;
1235  KK;
1236}
Note: See TracBrowser for help on using the repository browser.