source: git/Singular/LIB/ehv.lib @ f2aefad

spielwiese
Last change on this file since f2aefad was f2aefad, checked in by Andreas Steenpass <steenpass@…>, 5 years ago
fix replacements of execute() by create_ring()
  • Property mode set to 100644
File size: 41.0 KB
Line 
1////////////////////////////////////////////////////////////////////
2version="version ehv.lib 4.1.2.0 Feb_2019 "; // $Id$
3category="Commutative Algebra";
4
5info="
6LIBRARY: EHV.lib      PROCEDURES FOR PRIMARY DECOMPOSITION OF IDEALS
7AUTHORS: Kai Dehmann, dehmann@mathematik.uni-kl.de;
8
9OVERVIEW:
10  Algorithms for primary decomposition and radical-computation
11  based on the ideas of Eisenbud, Huneke, and Vasconcelos.
12
13PROCEDURES:
14  equiMaxEHV(I);                      equidimensional part of I
15  removeComponent(I,e);               intersection of the primary components
16                                      of I of dimension >= e
17  AssOfDim(I,e);                      an ideal such that the associated primes
18                                      are exactly the associated primes of I
19                                      having dimension e
20  equiRadEHV(I [,Strategy]);          equidimensional radical of I
21  radEHV(I [,Strategy]);              radical of I
22  IntAssOfDim1(I,e);                  intersection of the associated primes of I
23                                      having dimension e
24  IntAssOfDim2(I,e);                  another way of computing the intersection
25                                      of the associated primes of I
26                                      having dimension e
27  decompEHV(I);                       decomposition of a zero-dimensional
28                                      radical ideal I
29  AssEHV(I [,Strategy]);              associated primes of I
30  minAssEHV(I [,Strategy]);           minimal associated primes of I
31  localize(I,P,l);                    the contraction of the ideal generated by I
32                                      in the localization w.r.t P
33  componentEHV(I,P,L [,Strategy]);    a P-primary component for I
34  primdecEHV(I [,Strategy]);          a minimal primary decomposition of I
35  primDecsAreEquivalent(L, K);        procedure for comparing the output of
36                                      primary decomposition algorithms (checks
37                                      if the computed associated primes coincide)
38";
39
40LIB "ring.lib";
41LIB "general.lib";
42LIB "elim.lib";
43LIB "poly.lib";
44LIB "random.lib";
45LIB "inout.lib";
46LIB "matrix.lib";
47LIB "algebra.lib";
48LIB "normal.lib";
49
50
51/////////////////////////////////////////////////////////////////////
52//                                                                 //
53//             G E N E R A L   A L G O R I T H M S                 //
54//                                                                 //
55/////////////////////////////////////////////////////////////////////
56
57/////////////////////////////////////////////////////////////////////
58static proc AnnExtEHV(int n,list re)
59"USAGE:   AnnExtEHV(n,re); n integer, re resolution
60RETURN:  ideal, the annihilator of Ext^n(R/I,R) with given
61         resolution re of I"
62{
63  if(printlevel > 2){"Entering AnnExtEHV.";}
64
65  if(n < 0)
66    {
67      ideal ann = ideal(1);
68      if(printlevel > 2){"Leaving AnnExtEHV.";}
69      return(ann);
70    }
71  int l = size(re);
72
73  if(n < l)
74    {
75      matrix f = transpose(re[n+1]);
76      if(n == 0)
77        {
78          matrix g = 0*gen(ncols(f));
79        }
80      else
81        {
82          matrix g = transpose(re[n]);
83        }
84      module k = syz(f);
85      ideal ann = quotient1(g,k);
86      if(printlevel > 2){"Leaving AnnExtEHV.";}
87      return(ann);
88    }
89
90  if(n == l)
91    {
92      ideal ann = Ann(transpose(re[n]));
93      if(printlevel > 2){"Leaving AnnExtEHV.";}
94      return(ann);
95    }
96
97  ideal ann = ideal(1);
98  if(printlevel > 2){"Leaving AnnExtEHV.";}
99  return(ann);
100}
101
102
103/////////////////////////////////////////////////////////////////////
104//static
105proc isSubset(ideal I,ideal J)
106"USAGE:   isSubset(I,J); I, J ideals
107RETURN:  integer, 1 if I is a subset of J and 0 otherwise
108"
109{
110  if ( attrib(J,"isSB") == 0) { J = groebner(J); }
111  int s = ncols(I);
112  for(int i=1; i<=s; i++)
113    {
114      if(reduce(I[i],J)!=0)
115        {
116          return(0);
117        }
118    }
119  return(1);
120}
121
122
123/////////////////////////////////////////////////////////////////////
124//                                                                 //
125//         T H E   E Q U I D I M E N S I O N A L   P A R T         //
126//                                                                 //
127/////////////////////////////////////////////////////////////////////
128
129/////////////////////////////////////////////////////////////////////
130proc equiMaxEHV(ideal I)
131"USAGE:   equiMaxEHV(I); I ideal
132RETURN:  ideal, the equidimensional part of I.
133NOTE:    Uses algorithm of Eisenbud, Huneke, and Vasconcelos.
134EXAMPLE: example equiMaxEHV; shows an example
135"
136{
137  if(printlevel > 2){"Entering equiMaxEHV.";}
138  if(attrib(basering,"global")==0)
139    {
140      ERROR("// Not implemented for this ordering, please change to global ordering.");
141    }
142  ideal J = groebner(I);
143  int cod = nvars(basering)-dim(J);
144
145  //If I is the entire ring...
146  if(cod > nvars(basering))
147    {
148      //...then return the ideal generated by 1.
149      return(ideal(1));
150    }
151
152  //Compute a resolution of I.
153  if(printlevel > 2){"Computing resolution.";}
154  if(homog(I)==1)
155    {
156      list re = sres(J,cod+1);
157      re = minres(re);
158    }
159  else
160    {
161      list re = mres(I,cod+1);
162    }
163  if(printlevel > 2){"Finished computing resolution.";}
164
165  //Compute the annihilator of the cod-th EXT-module.
166  ideal ann = AnnExtEHV(cod,re);
167  attrib(ann,"isEquidimensional",1);
168  if(printlevel > 2){"Leaving equiMaxEHV.";}
169  return(ann);
170}
171
172example
173{
174  "EXAMPLE:";
175  echo = 2;
176  ring  r = 0,(x,y,z),dp;
177  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
178  equiMaxEHV(I);
179}
180
181
182/////////////////////////////////////////////////////////////////////
183proc removeComponent(ideal I, int e)
184"USAGE:   removeComponent(I,e); I ideal, e integer
185RETURN:  ideal, the intersection of the primary components
186         of I of dimension >= e
187EXAMPLE: example removeComponent; shows an example"
188{
189  if(attrib(basering,"global")==0)
190    {
191      ERROR("// Not implemented for this ordering, please change to global ordering.");
192    }
193
194  ideal J = groebner(I);
195
196  //Compute a resolution of I
197  if(homog(I)==1)
198    {
199      list re = sres(J,0);
200      re = minres(re);
201    }
202  else
203    {
204      list re = mres(I,0);
205    }
206
207  int f = nvars(basering);
208  int cod;
209  ideal ann;
210  int g = nvars(basering) - e;
211  while(f > g)
212    {
213      ann = AnnExtEHV(f,re);
214      cod = nvars(basering) - dim(groebner(ann));
215      if( cod == f )
216        {
217          I = quotient(I,ann);
218        }
219      f = f-1;
220    }
221  return(I);
222}
223
224example
225{
226  "EXAMPLE:";
227  echo = 2;
228  ring  r = 0,(x,y,z),dp;
229  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
230  removeComponent(I,1);
231}
232
233
234/////////////////////////////////////////////////////////////////////
235proc AssOfDim(ideal I, int e)
236"USAGE:   AssOfDim(I,e); I ideal, e integer
237RETURN:  ideal, such that the associated primes are exactly
238         the associated primes of I having dimension e
239EXAMPLE: example AssOfDim; shows an example"
240{
241  if(attrib(basering,"global")==0)
242    {
243      ERROR("// Not implemented for this ordering, please change to global ordering.");
244    }
245  int g = nvars(basering) - e;
246
247  //Compute a resolution of I.
248  ideal J = std(I);
249  if(homog(I)==1)
250    {
251      list re = sres(J,g+1);
252      re = minres(re);
253    }
254  else
255    {
256      list re = mres(I,g+1);
257    }
258
259  ideal ann = AnnExtEHV(g,re);
260  int cod = nvars(basering) - dim(std(ann));
261
262  //If the codimension of I_g:=Ann(Ext^g(R/I,R)) equals g...
263  if(cod == g)
264    {
265      //...then return the equidimensional part of I_g...
266      ann = equiMaxEHV(ann);
267      attrib(ann,"isEquidimensional",1);
268      return(ann);
269    }
270  //...otherwise...
271  else
272    {
273      //...I has no associated primes of dimension e.
274      return(ideal(1));
275    }
276}
277
278example
279{
280  "EXAMPLE:";
281  echo = 2;
282  ring  r = 0,(x,y,z),dp;
283  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
284  AssOfDim(I,1);
285}
286
287/////////////////////////////////////////////////////////////////////
288//                                                                 //
289//                      T H E   R A D I C A L                      //
290//                                                                 //
291/////////////////////////////////////////////////////////////////////
292
293/////////////////////////////////////////////////////////////////////
294static proc aJacob(ideal I, int a)
295"USAGE:   aJacob(I,a); I ideal, a integer
296RETURN:  ideal, the ath-Jacobian ideal of I"
297{
298  matrix M = jacob(I);
299  int n = nvars(basering);
300  if(n-a <= 0)
301    {
302      return(ideal(1));
303    }
304  if(n-a > nrows(M) or n-a > ncols(M))
305    {
306      return(ideal(0));
307    }
308  ideal J = minor(M,n-a);
309  return(J);
310}
311
312
313/////////////////////////////////////////////////////////////////////
314proc equiRadEHV(ideal I, list #)
315"USAGE:   equiRadEHV(I [,Strategy]); I ideal, Strategy list
316RETURN:  ideal, the equidimensional radical of I,
317         i.e. the intersection of the minimal associated primes of I
318         having the same dimension as I
319NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos,
320         Works only in characteristic 0 or p large.
321         The (optional) second argument determines the strategy used:
322         Strategy[1]        > strategy for the equidimensional part
323                     = 0    : uses equiMaxEHV
324                     = 1    : uses equidimMax
325         Strategy[2]        > strategy for the radical
326                     = 0    : combination of strategy 1 and 2
327                     = 1    : computation of the radical just with the
328                              help of regular sequences
329                     = 2    : does not try to find a regular sequence
330         If no second argument is given then Strategy=(0,0,0) is used.
331EXAMPLE: example equiRadEHV; shows an example"
332{
333  if(printlevel > 2){"Entering equiRadEHV.";}
334  if(attrib(basering,"global")==0)
335    {
336      ERROR("// Not implemented for this ordering, please change to global ordering.");
337    }
338  if((char(basering)<100)&&(char(basering)!=0))
339    {
340      "WARNING: The characteristic is too small, the result may be wrong";
341    }
342
343  //Define the Strategy to be used.
344  if(size(#) > 0)
345    {
346      if(#[1]!=1)
347        {
348          int equStr = 0;
349        }
350      else
351        {
352          int equStr = 1;
353        }
354      if(size(#) > 1)
355        {
356          if(#[2]!=1 and #[2]!=2)
357            {
358              int strategy = 0;
359            }
360          else
361            {
362              int strategy = #[2];
363            }
364        }
365      else
366        {
367          int strategy = 0;
368        }
369    }
370  else
371    {
372      int equStr = 0;
373      int strategy = 0;
374    }
375
376  ideal J,I0,radI0,L,radI1,I2,radI2;
377  int l,n;
378  intvec op = option(get);
379  matrix M;
380
381  option(redSB);
382  list m = mstd(I);
383  option(set,op);
384
385  int d = dim(m[1]);
386  if(d==-1)
387    {
388      return(ideal(1));
389    }
390
391  if(strategy != 2)
392    {
393      /////////////////////////////////////////////
394      // Computing the equidimensional radical   //
395      // via regular sequenves                   //
396      /////////////////////////////////////////////
397
398      if(printlevel > 2){"Trying to find a regular sequence.";}
399      int cod = nvars(basering)-d;
400
401      //Complete intersection case:
402      if(cod==size(m[2]))
403        {
404          J = aJacob(m[2],d);
405          if(printlevel > 2){"Leaving equiRadEHV.";}
406          return(quotient(m[2],J));
407        }
408
409      //First codim elements of I are a complete intersection:
410      for(l=1; l<=cod; l++)
411        {
412          I0[l] = m[2][l];
413        }
414      n = dim(groebner(I0))+cod-nvars(basering);
415
416      //Last codim elements of I are a complete intersection:
417      if(n!=0)
418              {
419          for(l=1; l<=cod; l++)
420            {
421              I0[l] = m[2][size(m[2])-l+1];
422            }
423          n = dim(groebner(I0))+cod-nvars(basering);
424        }
425
426      //Taking a generic linear combination of the input:
427      if(n!=0)
428              {
429          M = transpose(sparsetriag(size(m[2]),cod,95,1));
430          I0 = ideal(M*transpose(m[2]));
431          n = dim(groebner(I0))+cod-nvars(basering);
432        }
433
434      //Taking a more generic linear combination of the input:
435      if(n!=0)
436        {
437          while(strategy == 1 and n!=0)
438            {
439              M = transpose(sparsetriag(size(m[2]),cod,0,100));
440              I0 = ideal(M*transpose(m[2]));
441              n = dim(groebner(I0))+cod-nvars(basering);
442            }
443        }
444
445      if(n==0)
446        {
447          J = aJacob(I0,d);
448          if(printlevel > 2){"1st quotient.";}
449          radI0 = quotient(I0,J);
450          if(printlevel > 2){"2nd quotient.";}
451          L = quotient(radI0,m[2]);
452          if(printlevel > 2){"3rd quotient.";}
453          radI1 = quotient(radI0,L);
454          attrib(radI1,"isEquidimensional",1);
455          attrib(radI1,"isRadical",1);
456          if(printlevel > 2){"Leaving equiRadEHV.";}
457          return(radI1);
458        }
459    }
460
461  ////////////////////////////////////////////////////
462  // Computing the equidimensional radical directly //
463  ////////////////////////////////////////////////////
464
465  if(printlevel > 2){"Computing the equidimensional radical directly";}
466
467  //Compute the equidimensional part depending on the chosen strategy
468  if(equStr == 0)
469    {
470      I = equiMaxEHV(I);
471    }
472  if(equStr == 1)
473    {
474      I = equidimMax(I);
475    }
476  int a = nvars(basering)-1;
477
478  while(a > d)
479    {
480      if(printlevel > 2){"While-Loop: "+string(a);}
481      J = aJacob(I,a);
482      while(dim(groebner(J+I))==d)
483        {
484          if(printlevel > 2){"Quotient-Computation.";}
485          I = quotient(I,J);
486          if(printlevel > 2){"Computing the a-th Jacobian";}
487            J = aJacob(I,a);
488        }
489      a = a-1;
490    }
491  if(printlevel > 2){"We left While-Loop.";}
492  if(printlevel > 2){"Computing the a-th Jacobian";}
493  J = aJacob(I,d);
494  if(printlevel > 2){"Quotient-Computation.";}
495  I = quotient(I,J);
496  attrib(I,"isEquidimensional",1);
497  attrib(I,"isRadical",1);
498  if(printlevel > 2){"Leaving equiRadEHV.";}
499  return(I);
500}
501
502example
503{
504  "EXAMPLE:";
505  echo = 2;
506  ring  r = 0,(x,y,z),dp;
507  poly  p = z2+1;
508  poly  q = z3+2;
509  ideal i = p*q^2,y-z2;
510  ideal pr= equiRadEHV(i);
511  pr;
512}
513
514/////////////////////////////////////////////////////////////////////
515proc radEHV(ideal I, list #)
516"USAGE:   radEHV(I [,Strategy]); ideal I, Strategy list
517RETURN:  ideal, the radical of I
518NOTE:    uses the algorithm of Eisenbud/Huneke/Vasconcelos
519         Works only in characteristic 0 or p large.
520         The (optional) second argument determines the strategy used:
521         Strategy[1]        > strategy for the equidimensional part
522                     = 0    : uses equiMaxEHV
523                     = 1    : uses equidimMax
524         Strategy[2]        > strategy for the radical
525                     = 0    : combination of strategy 1 and 2
526                     = 1    : computation of the radical just with the
527                              help of regular sequences
528                     = 2    : does not try to find a regular sequence
529         Strategy[3]        > strategy for the computation of ideal quotients
530                     = n    : uses quot(.,.,n) for the ideal quotient computations
531         If no second argument is given then Strategy=(0,0,0) is used.
532EXAMPLE: example radEHV; shows an example"
533{
534  if(printlevel > 2){"Entering radEHV.";}
535  if(attrib(basering,"global")==0)
536    {
537      ERROR("// Not implemented for this ordering, please change to global ordering.");
538    }
539  if (size(I) == 0)  { return( ideal(0) ); }
540  //Compute the equidimensional radical J of I.
541  ideal J = equiRadEHV(I,#);
542
543  //If I is the entire ring...
544  if(deg(J[1]) <= 0)
545    {
546      //...then return the ideal generated by 1...
547      return(ideal(1));
548    }
549
550  //...else remove the maximal dimensional components and
551  //compute the radical K of the lower dimensional components ...
552  ideal K = radEHV(sat(I,J)[1],#);
553
554  //..and intersect it with J.
555  K = intersect(J,K);
556  attrib(K,"isRadical",1);
557  return(K);
558}
559
560example
561{
562  "EXAMPLE:";
563  echo = 2;
564  ring  r = 0,(x,y,z),dp;
565  poly  p = z2+1;
566  poly  q = z3+2;
567  ideal i = p*q^2,y-z2;
568  ideal pr= radical(i);
569  pr;
570}
571
572/////////////////////////////////////////////////////////////////////
573proc IntAssOfDim1(ideal I, int e)
574"USAGE:   IntAssOfDim1(I,e); I idea, e integer
575RETURN:  ideal, the intersection of the associated primes of I having dimension e
576EXAMPLE: example IntAssOfDim1; shows an example"
577{
578  if(attrib(basering,"global")==0)
579    {
580      ERROR("// Not implemented for this ordering, please change to global ordering.");
581    }
582  int g = nvars(basering) - e;
583
584  //Compute a resolution of I.
585  ideal J = groebner(I);
586  if(homog(I)==1)
587    {
588      list re = sres(J,g+1);
589      re = minres(re);
590    }
591  else
592    {
593      list re = mres(I,g+1);
594    }
595
596  ideal ann = AnnExtEHV(g,re);
597  int cod = nvars(basering) - dim(groebner(ann));
598  //If the codimension of I_g:=Ann(Ext^g(R/I,I)) equals g...
599  if(cod == g)
600    {
601      //...then return the equidimensional radical of I_g...
602      ann = equiRadEHV(ann);
603      attrib(ann,"isEquidimensional",1);
604      attrib(ann,"isRadical",1);
605      return(ann);
606    }
607  else
608    {
609      //...else I has no associated primes of dimension e.
610      return(ideal(1));
611    }
612}
613
614example
615{
616  "EXAMPLE:";
617  echo = 2;
618  ring  r = 0,(x,y,z),dp;
619  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
620  IntAssOfDim1(I,1);
621}
622
623/////////////////////////////////////////////////////////////////////
624proc IntAssOfDim2(ideal I, int e)
625"USAGE:   IntAssOfDim2(I,e); I ideal, e integer
626RETURN:  ideal, the intersection of the associated primes of I having dimension e
627EXAMPLE: example IntAssOfDim2; shows an example"
628{
629  if(attrib(basering,"global")==0)
630    {
631      ERROR("// Not implemented for this ordering, please change to global ordering.");
632    }
633  ideal I1 = removeComponent(I,e);
634  ideal I2 = removeComponent(I,e+1);
635  ideal b = quotient(I1,I2);
636  b = equiRadEHV(b);
637  return(b);
638}
639
640example
641{
642  "EXAMPLE:";
643  echo = 2;
644  ring  r = 0,(x,y,z),dp;
645  ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
646  IntAssOfDim2(I,1);
647}
648
649
650/////////////////////////////////////////////////////////////////////
651//                                                                 //
652//   Z E R O - D I M E N S I O N A L   D E C O M P O S I T I O N   //
653//                                                                 //
654/////////////////////////////////////////////////////////////////////
655
656/////////////////////////////////////////////////////////////////////
657proc decompEHV(ideal I)
658"USAGE:   decompEHV(I); I zero-dimensional radical ideal
659RETURN:  list, the associated primes of I
660EXAMPLE: example decompEHV; shows an example"
661{
662  if(printlevel > 2){"Entering decompEHV.";}
663  if(attrib(basering,"global")==0)
664    {
665      ERROR("// Not implemented for this ordering, please change to global ordering.");
666    }
667
668  int e,m,vd,nfact;
669  list l,k,L;
670  poly f,h;
671  ideal J;
672  def base = basering;
673
674  if( attrib(I,"isSB")!=1 )
675    {
676      I = groebner(I);
677    }
678
679  while(1)
680    {
681      //Choose a random polynomial f from R.
682      e = random(0,100);
683      f = sparsepoly(e);
684
685      //Check if f lies not in I.
686      if(reduce(f,I)!=0)
687        {
688          J = quotient1(I,f);
689
690          //If f is a zerodivisor modulo I...
691          if( isSubset(J,I) == 0 )
692            {
693              //...then use recursion...
694              if(printlevel > 2){"We found a zerodivisor -- recursion";}
695              l = decompEHV(J) + decompEHV(I+f);
696              return(l);
697            }
698          if(printlevel > 2){"We found a non-zero-divisor.";}
699
700          //...else compute the vectorspace dimension vd of I and...
701          vd = vdim(I);
702
703          //...compute m minimal such that 1,f,f^2,...,f^m are linearly dependent.
704          qring Q = I;
705          poly g = fetch(base,f);
706          k = algDependent(g);
707          def R = k[2];
708          setring R;
709          if(size(ker)!=1)
710            {
711              //Calculate a generator for ker.
712              ker = mstd(ker)[2];
713            }
714          poly g = ker[1];
715          m = deg(g);
716
717          //If m and vd coincide...
718          if(m==vd)
719            {
720              if(printlevel > 2){"We have a good candidate.";}
721              //...then factorize g.
722              if(printlevel > 2){"Factorizing.";}
723              L = factorize(g,2);    //returns non-constant factors and multiplicities
724              nfact = size(L[1]);
725
726              //If g is irreducible...
727              if(nfact==1 and L[2][1]==1)
728                {
729                  if(printlevel > 2){"The element is irreducible.";}
730                  setring base;
731                  //..then I is a maximal ideal...
732                  l[1] = I;
733                  kill R,Q;
734                  return(l);
735                }
736              //...else...
737              else
738                {
739
740                  if(printlevel > 2){"The element is not irreducible -- recursion.";}
741                  //...take a non-trivial factor g1 of g ...
742                  poly g1 = L[1][1];
743                  //..and insert f in g1...
744                  ring newR = create_ring(ringlist(R)[1], "(y(1))", "(" + ordstr(R) + ")", "no_minpoly");
745                  poly g2 = imap(R,g1);
746                  setring base;
747                  h = fetch(newR,g2);
748                  h = subst(h,var(1),f);
749                  //...and use recursion...
750                  kill R,Q,newR;
751                  l = l + decompEHV(quotient1(I,h));
752                  l = l + decompEHV(I+h);
753                  return(l);
754                }
755            }
756          setring base;
757          kill R,Q;
758          }
759    }
760}
761
762example
763{
764  "EXAMPLE:";
765  echo = 2;
766  ring r = 32003,(x,y),dp;
767  ideal i = x2+y2-10,x2+xy+2y2-16;
768  decompEHV(i);
769}
770
771
772/////////////////////////////////////////////////////////////////////
773//                                                                 //
774//           A S S O C I A T E D   P R I M E S                     //
775//                                                                 //
776/////////////////////////////////////////////////////////////////////
777
778
779/////////////////////////////////////////////////////////////////////
780static proc idempotent(ideal I)
781"USAGE:   idempotent(I); ideal I (weighted) homogeneous radical ideal,
782         I intersected K[x(1),...,x(k)] zero-dimensional
783         where deg(x(i))=0 for all i <= k and deg(x(i))>0 for all i>k.
784RETURN:  a list of ideals I(1),...,I(t) such that
785         K[x]/I = K[x]/I(1) x ... x K[x]/I(t)"
786{
787  if(printlevel > 2){"Entering idempotent.";}
788  if(attrib(basering,"global")==0)
789    {
790      ERROR("// Not implemented for this ordering, please change to global ordering.");
791    }
792  int n = nvars(basering);
793  poly f,g;
794  string j;
795  int i,k,splits;
796  list l;
797
798  //Collect all variables of degree 0.
799  for(i=1; i<=n; i++)
800    {
801      if(deg(var(i)) > 0)
802        {
803          f = f*var(i);
804        }
805      else
806        {
807          splits = 1;
808          j = j +  string(var(i)) + "," ;
809        }
810    }
811  //If there are no variables of degree 0
812  //then there are no idempotents and we are done...
813  if(splits == 0)
814    {
815      l[1]=I;
816      return(l);
817    }
818
819  //...else compute J = I intersected K[x(1),...,x(k)]...
820  ideal J = eliminate(I,f);
821  def base = basering;
822  j = j[1,size(j)-1];
823  j = "ring @r = (" + charstr(basering) + "),(" + j + "),(" + ordstr(basering) + ");";
824  execute(j);
825  ideal J = imap(base,J);
826
827  //...and compute the associated primes of the zero-dimensional ideal J.
828  list L = decompEHV(J);
829  int s = size(L);
830  ideal K,Z;
831  poly g;
832  //For each associated prime ideal P_i of J...
833  for(i=1; i<=s; i++)
834    {
835      K = ideal(1);
836      //...comnpute the intersection K of the other associated prime ideals...
837      for(k=1; k<=s; k++)
838        {
839          if(i!=k)
840            {
841              K = intersect(K,L[k]);
842            }
843        }
844
845      //...and find an element that lies in K but not in P_i...
846      g = randomid(K,1)[1];
847      Z = L[i];
848      Z = std(Z);
849      while(reduce(g,Z)==0)
850        {
851          g = randomid(K,1)[1];
852        }
853      setring base;
854      g = imap(@r,g);
855      //...and compute the corresponding ideal I(i)
856      l[i] = quotient(I,g);
857      setring @r;
858
859    }
860  setring base;
861  return(l);
862}
863
864
865/////////////////////////////////////////////////////////////////////
866static proc equiAssEHV(ideal I)
867"USAGE:   equiAssEHV(I); I equidimensional, radical, and homogeneous ideal
868RETURN:  a list, the associated prime ideals of I"
869{
870  if(printlevel > 2){"Entering equiAssEHV.";}
871  if(attrib(basering,"global")==0)
872    {
873      ERROR("// Not implemented for this ordering, please change to global ordering.");
874    }
875  list L;
876  def base = basering;
877  int n = nvars(basering);
878  int i,j;
879
880  //Compute the normalization of I.
881  if(printlevel > 2){"Entering Normalization.";}
882  list norOut = normal(I, "noDeco");
883  list K = norOut[1];
884  if(printlevel > 2){"Leaving Normalisation.";}
885
886  //The normalization algorithm returns k factors.
887  int k = size(K);
888  if(printlevel > 1){"Normalization algorithm splits ideal in " + string(k) + " factors.";}
889
890  //Examine each factor of the normalization.
891  def P;
892  for(i=1; i<=k; i++)
893    {
894      P = K[i];
895      setring P;
896
897      //Use procedure idempotent to split the i-th factor of
898      //the normalization in a product of integral domains.
899      if(printlevel > 1){"Examining " + string(i) +". factor.";}
900      list l = idempotent(norid);
901      int leng = size(l);
902      if(printlevel > 1){"Idempotent algorithm splits factor " + string(i) + " in " + string(leng) + " factors.";}
903
904      //Intersect the minimal primes corresponding to the
905      //integral domains obtained from idempotent with the groundring,
906      //i.e. compute the preimages w.r.t. the corresponding normalization map
907      ideal J;
908      for(j=1; j<=leng; j++)
909        {
910          J = l[j];
911          setring base;
912          L[size(L)+1] = preimage(P,normap,J);
913          setring P;
914        }
915      kill l, leng, J;
916      setring base;
917    }
918
919  return(L);
920}
921
922
923/////////////////////////////////////////////////////////////////////
924proc AssEHV(ideal I, list #)
925"USAGE:  AssEHV(I [,Strategy]); I Ideal, Strategy list
926RETURN:  a list, the associated prime ideals of I
927NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos.
928         The (optional) second argument determines the strategy used:
929         Strategy[1]        > strategy for the equidimensional part
930                     = 0    : uses equiMaxEHV
931                     = 1    : uses equidimMax
932         Strategy[2]        > strategy for the equidimensional radical
933                     = 0    : uses equiRadEHV
934                     = 1    : uses equiRadical
935         Strategy[3]        > strategy for equiRadEHV
936                     = 0    : combination of strategy 1 and 2
937                     = 1    : computation of the radical just with the
938                              help of regular sequences
939                     = 2    : does not try to find a regular sequence
940         Strategy[4]        > strategy for the computation of ideal quotients
941                     = n    : uses quot(.,.,n) for the ideal quotient computations
942         If no second argument is given, Strategy=(0,0,0,0) is used.
943EXAMPLE: example AssEHV; shows an example"
944{
945  if(attrib(basering,"global")==0)
946    {
947      ERROR("// Not implemented for this ordering, please change to global ordering.");
948    }
949  if(printlevel > 2){"Entering AssEHV";}
950
951  //Specify the strategy to be used.
952  if(size(#)==0)
953    {
954      # = 0,0,0,0;
955    }
956  if(size(#)==1)
957    {
958      # = #[1],0,0,0;
959    }
960  if(size(#)==2)
961    {
962      # = #[1],#[2],0,0;
963    }
964  if(size(#)==3)
965    {
966      # = #[1],#[2],#[3],0;
967    }
968
969  list L;
970  ideal K;
971  def base = basering;
972  int m,j;
973  ideal J = groebner(I);
974  int n = nvars(basering);
975  int d = dim(J);
976
977  //Compute a resolution of I.
978  if(printlevel > 2){"Computing resolution.";}
979  if(homog(I)==1)
980    {
981      list re = sres(J,0);
982      re = minres(re);
983    }
984  else
985    {
986      list re = mres(I,0);
987    }
988  ideal ann;
989  int cod;
990
991  //For 0<=i<= dim(I) compute the intersection of the i-dimensional associated primes of I
992  for(int i=0; i<=d; i++)
993    {
994      if(printlevel > 1){"Are there components of dimension " + string(i) + "?";}
995      ann = AnnExtEHV(n-i,re);
996      cod = n - dim(groebner(ann));
997
998      //If there are associated primes of dimension i...
999      if(cod == n-i)
1000        {
1001          if(printlevel > 1){"Yes. There are components of dimension " + string(i) + ".";}
1002          //...then compute the intersection K of all associated primes of I of dimension i
1003          if(#[2]==0)
1004            {
1005              if(size(#) > 3)
1006                {
1007                  K = equiRadEHV(ann,#[1],#[3],#[4]);
1008                }
1009              if(size(#) > 2)
1010                {
1011                  K = equiRadEHV(ann,#[1],#[3]);
1012                }
1013              else
1014                {
1015                  K = equiRadEHV(ann,#[1]);
1016                }
1017            }
1018          if(#[2]==1)
1019            {
1020              K = equiRadical(ann);
1021            }
1022          attrib(K,"isEquidimensional",1);
1023          attrib(K,"isRadical",1);
1024
1025          //If K is already homogeneous then use equiAssEHV to recover the associated primes of K...
1026           if(homog(K)==1)
1027            {
1028              if(printlevel > 2){"Input already homogeneous.";}
1029              L = L + equiAssEHV(K);
1030            }
1031
1032          //...else...
1033          else
1034            {
1035              //...homogenize K w.r.t. t,...
1036              if(printlevel > 2){"Input not homogeneous; must homogenize.";}
1037              def homoR=changeord(list(list("dp",1:nvars(basering))));
1038              setring homoR;
1039              ideal homoJ = fetch(base,K);
1040              homoJ = groebner(homoJ);
1041              list l2;
1042              for (int ii = 1; ii <= n; ii++)
1043              {
1044                l2[i] = "x("+string(ii)+")";
1045              }
1046              l2[n+1] = "t";
1047              ring newR = create_ring(ringlist(base)[1], l2, "dp", "no_minpoly");
1048              ideal homoK = fetch(homoR,homoJ);
1049              homoK = homog(homoK,t);
1050              attrib(homoK,"isEquidimensional",1);
1051              attrib(homoK,"isRadical",1);
1052
1053              //...compute the associated primes of the homogenization using equiAssEHV,...
1054              list l = equiAssEHV(homoK);
1055
1056              //...and set t=1 in the generators of the associated primes just computed.
1057              ideal Z;
1058              for(j=1; j<=size(l); j++)
1059                {
1060                  Z = subst(l[j],t,1);
1061                  setring base;
1062                  L[size(L)+1] = fetch(newR,Z);
1063                  setring newR;
1064                }
1065              setring base;
1066              kill homoR;
1067              kill newR;
1068            }
1069        }
1070      else
1071        {
1072          if(printlevel > 1){"No. There are no components of dimension " + string(i) + ".";}
1073        }
1074    }
1075  return(L);
1076}
1077
1078example
1079{
1080  "EXAMPLE:";
1081  echo = 2;
1082  ring  r = 0,(x,y,z),dp;
1083  poly  p = z2+1;
1084  poly  q = z3+2;
1085  ideal i = p*q^2,y-z2;
1086  list pr = AssEHV(i);
1087  pr;
1088}
1089
1090/////////////////////////////////////////////////////////////////////
1091proc minAssEHV(ideal I, list #)
1092"USAGE:  minAssEHV(I [,Strategy]); I ideal, Strategy list
1093RETURN:  a list, the minimal associated prime ideals of I
1094NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos.
1095         The (optional) second argument determines the strategy used:
1096         Strategy[1]        > strategy for the equidimensional part
1097                     = 0    : uses equiMaxEHV
1098                     = 1    : uses equidimMax
1099         Strategy[2]        > strategy for the equidimensional radical
1100                     = 0    : uses equiRadEHV, resp. radicalEHV
1101                     = 1    : uses equiRadical, resp. radical
1102         Strategy[3]        > strategy for equiRadEHV
1103                     = 0    : combination of strategy 1 and 2
1104                     = 1    : computation of the radical just with the
1105                              help of regular sequences
1106                     = 2    : does not try to find a regular sequence
1107         Strategy[4]        > strategy for the computation of ideal quotients
1108                     = n    : uses quot(.,.,n) for the ideal quotient computations
1109         If no second argument is given, Strategy=(0,0,0,0) is used.
1110EXAMPLE: example minAssEHV; shows an example"
1111{
1112  if(attrib(basering,"global")==0)
1113    {ERROR("// Not implemented for this ordering, please change to global ordering.");}
1114
1115  //Specify the strategy to be used.
1116  if(size(#)==0)
1117    {
1118      # = 0,0,0,0;
1119    }
1120  if(size(#)==1)
1121    {
1122      # = #[1],0,0,0;
1123    }
1124  if(size(#)==2)
1125    {
1126      # = #[1],#[2],0,0;
1127    }
1128  if(size(#)==3)
1129    {
1130      # = #[1],#[2],#[3],0;
1131    }
1132
1133  //Compute the radical of I.
1134  if(#[2]==0)
1135    {
1136      I = radEHV(I,#);
1137    }
1138  if(#[2]==1)
1139    {
1140      I = radical(I);
1141    }
1142
1143  //Compute the associated primes of the radical.
1144  return(AssEHV(I,#));
1145}
1146
1147example
1148{
1149  "EXAMPLE:";
1150  echo = 2;
1151  ring  r = 0,(x,y,z),dp;
1152  poly  p = z2+1;
1153  poly  q = z3+2;
1154  ideal i = p*q^2,y-z2;
1155  list pr = minAssEHV(i);
1156  pr;
1157}
1158
1159
1160/////////////////////////////////////////////////////////////////////
1161//                                                                 //
1162//        P R I M A R Y   D E C O M P O S I T I O N                //
1163//                                                                 //
1164/////////////////////////////////////////////////////////////////////
1165
1166
1167/////////////////////////////////////////////////////////////////////
1168proc localize(ideal I, ideal P, list l)
1169"USAGE:   localize(I,P,l); I ideal, P an associated prime ideal of I,
1170         l list of all associated primes of I
1171RETURN:  ideal,  the contraction of the ideal generated by I
1172         in the localization w.r.t P
1173EXAMPLE: example localize; shows an example"
1174{
1175  if(attrib(basering,"global")==0)
1176    {
1177      ERROR("// Not implemented for this ordering, please change to global ordering.");
1178    }
1179
1180  ideal Intersection = ideal(1);
1181  int s = size(l);
1182  if(attrib(P,"isSB")!=1)
1183    {
1184      P = groebner(P);
1185    }
1186
1187  //Compute the intersection of all associated primes of I that are not contained in P.
1188  for(int i=1; i<=s; i++)
1189    {
1190      if(isSubset(l[i],P)!=1)
1191        {
1192          Intersection = intersect(Intersection,l[i]);
1193        }
1194    }
1195  Intersection = groebner(Intersection);
1196
1197  //If the intersection is the entire ring...
1198  if(reduce(1,Intersection)==0)
1199    {
1200      //...then return I...
1201      return(I);
1202    }
1203  //...else try to find an element f that lies in the intersection but outside P...
1204  poly f = 0;
1205  while(reduce(f,P) == 0)
1206    {
1207      f = randomid(Intersection,1)[1];
1208    }
1209
1210  //...and saturate I w.r.t. f.
1211  I = sat(I,f)[1];
1212  return(I);
1213}
1214
1215/////////////////////////////////////////////////////////////////////
1216proc componentEHV(ideal I, ideal P, list L, list #)
1217"USAGE:   componentEHV(I,P,L [,Strategy]);
1218         I ideal, P associated prime of I,
1219         L list of all associated primes of I, Strategy list
1220RETURN:  ideal, a P-primary component Q for I
1221NOTE:    The (optional) second argument determines the strategy used:
1222         Strategy[1]        > strategy for equidimensional part
1223                     = 0    : uses equiMaxEHV
1224                     = 1    : uses equidimMax
1225         If no second argument is given then Strategy=0 is used.
1226EXAMPLE: example componentEHV; shows an example"
1227{
1228  if(printlevel > 2){"Entering componentEHV.";}
1229  if(attrib(basering,"global")==0)
1230    {
1231      ERROR("// Not implemented for this ordering, please change to global ordering.");
1232    }
1233
1234  //If no strategy is specified use standard strategy.
1235  if(size(#)==0)
1236    {
1237      # = 0;
1238    }
1239
1240  ideal T = P;
1241  ideal Q;
1242
1243  //Compute the localization of I at P...
1244  ideal IP = groebner(localize(I,P,L));
1245
1246  //...and compute the saturation of the localization w.r.t. P.
1247  ideal IP2 = sat(IP,P)[1];
1248
1249  //As long as we have not found a primary component...
1250  int isPrimaryComponent = 0;
1251  while(isPrimaryComponent!=1)
1252    {
1253      //...compute the equidimensional part Q of I+P^n...
1254      if(#[1]==0)
1255        {
1256          Q = equiMaxEHV(I+T);
1257        }
1258      if(#[1]==1)
1259        {
1260          Q = equidimMax(I+T);
1261        }
1262      //...and check if it is a primary component for P.
1263      if(isSubset(intersect(IP2,Q),IP)==1)
1264        {
1265          isPrimaryComponent = 1;
1266        }
1267      else
1268        {
1269          T = T*P;
1270        }
1271    }
1272  if(printlevel > 2){"Leaving componentEHV.";}
1273  return(Q);
1274}
1275
1276example
1277{
1278  "EXAMPLE:";
1279  echo = 2;
1280  ring  r = 0,(x,y,z),dp;
1281  poly  p = z2+1;
1282  poly  q = z3+2;
1283  ideal i = p*q^2,y-z2;
1284  list pr = AssEHV(i);
1285  componentEHV(i,pr[1],pr);
1286}
1287
1288
1289/////////////////////////////////////////////////////////////////////
1290proc primdecEHV(ideal I, list #)
1291"USAGE:   primdecEHV(I [,Strategy]); I ideal, Strategy list
1292RETURN:  a list pr of primary ideals and their associated primes:
1293                pr[i][1]   the i-th primary component,
1294                pr[i][2]   the i-th prime component.
1295NOTE:    Algorithm of Eisenbud/Huneke/Vasconcelos.
1296         The (optional) second argument determines the strategy used:
1297         Strategy[1]        > strategy for equidimensional part
1298                     = 0    : uses equiMaxEHV
1299                     = 1    : uses equidimMax
1300         Strategy[2]        > strategy for equidimensional radical
1301                     = 0    : uses equiRadEHV, resp. radicalEHV
1302                     = 1    : uses equiRadical, resp. radical
1303         Strategy[3]        > strategy for equiRadEHV
1304                     = 0    : combination of strategy 1 and 2
1305                     = 1    : computation of the radical just with the
1306                              help of regular sequences
1307                     = 2    : does not try to find a regular sequence
1308         Strategy[4]        > strategy for the computation of ideal quotients
1309                     = n    : uses quot(.,.,n) for the ideal quotient computations
1310         If no second argument is given then Strategy=(0,0,0,0) is used.
1311EXAMPLE: example primdecEHV; shows an example"
1312{
1313  if(printlevel > 2){"Entering primdecEHV.";}
1314  if(attrib(basering,"global")==0)
1315    {
1316      ERROR("// Not implemented for this ordering, please change to global ordering.");
1317    }
1318  list L,K;
1319
1320  //Specify the strategy to be used.
1321  if(size(#)==0)
1322    {
1323      # = 0,0,0,0;
1324    }
1325  if(size(#)==1)
1326    {
1327      # = #[1],0,0,0;
1328    }
1329  if(size(#)==2)
1330    {
1331      # = #[1],#[2],0,0;
1332    }
1333  if(size(#)==3)
1334    {
1335      # = #[1],#[2],#[3],0;
1336    }
1337
1338  //Compute the associated primes of I...
1339  L = AssEHV(I,#);
1340  if(printlevel > 0){"We have " + string(size(L)) + " prime components.";}
1341
1342  //...and compute for each associated prime of I a corresponding primary component.
1343  int l = size(L);
1344  for(int i=1; i<=l; i++)
1345    {
1346      K[i] = list();
1347      K[i][2] = L[i];
1348      K[i][1] = componentEHV(I,L[i],L,#[1]);
1349    }
1350  if(printlevel > 2){"Leaving primdecEHV.";}
1351  return(K);
1352
1353}
1354
1355example
1356{
1357  "EXAMPLE:";
1358  echo = 2;
1359  ring  r = 0,(x,y,z),dp;
1360  poly  p = z2+1;
1361  poly  q = z3+2;
1362  ideal i = p*q^2,y-z2;
1363  list pr = primdecEHV(i);
1364  pr;
1365}
1366
1367
1368/////////////////////////////////////////////////////////////////////
1369//                                                                 //
1370//        A L G O R I T H M S   F O R   T E S T I N G              //
1371//                                                                 //
1372/////////////////////////////////////////////////////////////////////
1373proc removeRedundantComponents(list primdecResult)
1374"USAGE:   removeRedundantComponents(L); L a primary decomposition
1375RETURN:  an irredundant primary decomposition.
1376EXAMPLE: example removeRedundantComponents; shows an example"
1377proc removeRedundantComponents(list primdecResult)
1378{
1379   int i,j;
1380
1381   i = 1;
1382   j = 1;
1383   while  ( i<=size(primdecResult) )
1384   {
1385     ASSUME(1, i>0 &&  i<=size(primdecResult) );
1386     ASSUME(1, j>0 &&  j<=size(primdecResult) );
1387
1388      if (i!=j)
1389      {
1390          if ( isSubset( std(primdecResult[j][1]) , std(primdecResult[i][1]) ) )
1391          {
1392              primdecResult[i] =  primdecResult[j];
1393              primdecResult = delete(primdecResult,j);
1394              j = j - 1; /*correct j due to deletion */
1395              if ( j<i )  {    i = i -1;   } /*correct i */
1396          }
1397      }
1398      j = j + 1;
1399      if ( j > size(primdecResult) )
1400      {
1401         i = i + 1;
1402         j = 1;
1403      }
1404   }
1405
1406   return (primdecResult);
1407}
1408example
1409{
1410  "EXAMPLE:";
1411  ring rng = 0,(x,y),dp;
1412  def L1 = list(list(ideal(x-1),ideal(x-1)),list(ideal(x-1),ideal(x-1)),list(ideal(y-2),ideal(y-2)));
1413  def L2 = list(list(ideal(x-1),ideal(x-1)),list(ideal(y-2),ideal(y-2)));
1414  L1;
1415  L2;
1416  ASSUME(0, primDecsAreEquivalent( L1, L2 ) );
1417}
1418
1419/////////////////////////////////////////////////////////////////////
1420proc primDecsAreEquivalent(list L, list K)
1421"USAGE:   primDecsAreEquivalent(L,K); L,K list of ideals
1422RETURN:  integer, 1 if the lists are the same up to ordering and redundant components; 0 otherwise
1423EXAMPLE: example primDecsAreEquivalent; shows an example"
1424{
1425  L = removeRedundantComponents(L);
1426  K = removeRedundantComponents(K);
1427  int s1 = size(L);
1428  int s2 = size(K);
1429  if(s1!=s2)
1430    {
1431      return(0);
1432    }
1433  list L1, K1;
1434  int i,j,t;
1435  list N;
1436  for(i=1; i<=s1; i++)
1437    {
1438      L1[i]=std(L[i][2]);
1439      K1[i]=std(K[i][2]);
1440    }
1441  for(i=1; i<=s1; i++)
1442    {
1443      for(j=1; j<=s1; j++)
1444        {
1445          if(isSubset(L1[i],K1[j]))
1446            {
1447              if(isSubset(K1[j],L1[i]))
1448                {
1449                  for(t=1; t<=size(N); t++)
1450                    {
1451                      if(N[t]==j)
1452                        {
1453                          return(0);
1454                        }
1455                    }
1456                  N[size(N)+1]=j;
1457                }
1458
1459            }
1460        }
1461    }
1462  if ( size(N) != s1 )
1463  {
1464     return(0);
1465  }
1466  return(1);
1467}
1468
1469example
1470{
1471  "EXAMPLE:";
1472  echo = 2;
1473  ring  r = 0,(x,y),dp;
1474  ideal i = x2,xy;
1475  list L1 = primdecGTZ(i);
1476  list L2 = primdecEHV(i);
1477  primDecsAreEquivalent(L1,L2);
1478}
Note: See TracBrowser for help on using the repository browser.