source: git/Singular/LIB/normal.lib @ 9a3ece

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