source: git/Singular/LIB/primdec.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: 67.6 KB
Line 
1// $Id: primdec.lib,v 1.6 1997-11-05 18:16:57 Singular Exp $
2///////////////////////////////////////////////////////
3// primdec.lib
4// algorithms for primary decomposition based on
5// the ideas of Gianni,Trager,Zacharias
6// written by Gerhard Pfister
7//////////////////////////////////////////////////////
8
9LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I)
10
11  minAssPrimes (ideal I)
12  //computes the minimal associated primes of the ideal I
13
14  decomp (ideal I)
15  // Computes a complete primary decomposition via Gianni,Trager,Zacharias
16
17  radical(ideal I)
18  //computes the radical of the ideal I
19
20  equiRadical(ideal I)
21  //computes the radical of the equidimensional part of the ideal I
22
23  prepareAss(ideal I)
24  //computes the radicals of the equidimensional parts of I
25
26LIB "random.lib";
27LIB "elim.lib";
28///////////////////////////////////////////////////////////////////////////////
29
30proc sat1 (ideal id, poly p)
31USAGE:   sat1(id,j);  id ideal, j polynomial
32RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
33NOTE:    result is a std basis in the basering
34EXAMPLE: example sat; shows an example
35{
36   int @k;
37   ideal inew=std(id);
38   ideal iold;
39   option(returnSB);
40   while(specialIdealsEqual(iold,inew)==0 )
41   {
42      iold=inew;
43      inew=quotient(iold,p);
44      @k++;
45   }
46   @k--;
47   option(noreturnSB);
48   list L =inew,p^@k;
49   return (L);
50}
51
52///////////////////////////////////////////////////////////////////////////////
53
54proc sat2 (ideal id, ideal h)
55USAGE:   sat2(id,j);  id ideal, j polynomial
56RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
57NOTE:    result is a std basis in the basering
58EXAMPLE: example sat2; shows an example
59{
60   int @k,@i;
61   def @P= basering;
62   if(ordstr(basering)[1,2]!="dp")
63   {
64      execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),dp;";
65      ideal inew=std(imap(@P,id));
66      ideal  @h=imap(@P,h);
67   }
68   else
69   {
70      ideal @h=h;
71      ideal inew=std(id);
72   }
73   ideal fac;
74
75   for(@i=1;@i<=ncols(@h);@i++)
76   {
77     if(deg(@h[@i])>0)
78     {
79        fac=fac+factorize(@h[@i],1);
80     }
81   }
82   fac=simplify(fac,4);
83   poly @f=1;
84   if(deg(fac[1])>0)
85   {
86      ideal iold;
87
88      for(@i=1;@i<=size(fac);@i++)
89      {
90        @f=@f*fac[@i];
91      }
92      option(returnSB);
93      while(specialIdealsEqual(iold,inew)==0 )
94      {
95         iold=inew;
96         if(deg(iold[size(iold)])!=1)
97         {
98            inew=quotient(iold,@f);
99         }
100         else
101         {
102            inew=iold;
103         }
104         @k++;
105      }
106      option(noreturnSB);
107      @k--;
108   }
109
110   if(ordstr(@P)[1,2]!="dp")
111   {
112      setring @P;
113      ideal inew=std(imap(@Phelp,inew));
114      poly @f=imap(@Phelp,@f);
115   }
116   list L =inew,@f^@k;
117   return (L);
118}
119
120///////////////////////////////////////////////////////////////////////////////
121
122proc minSat(ideal inew, ideal h)
123{
124   int i,k;
125   poly f=1;
126   ideal iold,fac;
127   list quotM,l;
128
129   for(i=1;i<=ncols(h);i++)
130   {
131      if(deg(h[i])>0)
132      {
133         fac=fac+factorize(h[i],1);
134      }
135   }
136   fac=simplify(fac,4);
137   if(size(fac)==0)
138   {
139      l=inew,1;
140      return(l);
141   }
142   fac=sort(fac)[1];
143   for(i=1;i<=size(fac);i++)
144   {
145      f=f*fac[i];
146   }
147   quotM[1]=inew;
148   quotM[2]=fac;
149   quotM[3]=f;
150   f=1;
151   option(returnSB);
152   while(specialIdealsEqual(iold,quotM[1])==0)
153   {
154      if(k>0)
155      {
156         f=f*quotM[3];
157      }
158      iold=quotM[1];
159      quotM=quotMin(quotM);
160      k++;
161   }
162   option(noreturnSB);
163   l=quotM[1],f;
164   return(l);
165}
166
167proc quotMin(list tsil)
168{
169   int i,j,k,action;
170   ideal verg;
171   list l;
172   poly g;
173
174   ideal laedi=tsil[1];
175   ideal fac=tsil[2];
176   poly f=tsil[3];
177
178   ideal star=quotient(laedi,f);
179
180   action=1;
181
182   while(action==1)
183   {
184      if(size(fac)==1)
185      {
186         action=0;
187         break;
188      }
189      for(i=1;i<=size(fac);i++)
190      {
191        g=1;
192         for(j=1;j<=size(fac);j++)
193         {
194            if(i!=j)
195            {
196               g=g*fac[j];
197            }
198         }
199         verg=quotient(laedi,g);
200         if(specialIdealsEqual(verg,star)==1)
201         {
202            f=g;
203            fac[i]=0;
204            fac=simplify(fac,2);
205            break;
206         }
207         if(i==size(fac))
208         {
209            action=0;
210         }
211      }
212   }
213   l=star,fac,f;
214   return(l);
215}
216
217
218////////////////////////////////////////////////////////////////////////////////
219proc testFactor(list act,poly p)
220{
221   int i;
222   poly q=act[1][1]^act[2][1];
223   for(i=2;i<=size(act[1]);i++)
224   {
225      q=q*act[1][i]^act[2][i];
226   }
227   q=1/leadcoef(q)*q;
228   p=1/leadcoef(p)*p;
229   if(p-q!=0)
230   {
231      "ERROR IN FACTOR";
232      act;
233      p;
234      q;
235      pause;
236   }
237}
238////////////////////////////////////////////////////////////////////////////////
239
240proc factor(poly p)
241USAGE:   factor(p) p poly
242RETURN:  list=;
243NOTE:
244EXAMPLE: example factor; shows an example
245{
246
247  ideal @i;
248  list @l;
249  intvec @v,@w;
250  int @j,@k,@n;
251
252  if(deg(p)<=1)
253  {
254     @i=ideal(p);
255     @v=1;
256     @l[1]=@i;
257     @l[2]=@v;
258     return(@l);
259  }
260  if (size(p)==1)
261  {
262     @w=leadexp(p);
263     for(@j=1;@j<=nvars(basering);@j++)
264     {
265        if(@w[@j]!=0)
266        {
267           @k++;
268           @v[@k]=@w[@j];
269           @i=@i+ideal(var(@j));
270        }
271     }
272     @l[1]=@i;
273     @l[2]=@v;
274     return(@l);
275  }
276  @l=factorize(p,2);
277  if(npars(basering)>0)
278  {
279     for(@j=1;@j<=size(@l[1]);@j++)
280     {
281        if(deg(@l[1][@j])==0)
282        {
283           @n++;
284        }
285     }
286     if(@n>0)
287     {
288        if(@n==size(@l[1]))
289        {
290           @l[1]=ideal(1);
291           @v=1;
292           @l[2]=@v;
293        }
294        else
295        {
296           @k=0;
297           int pleh;
298           for(@j=1;@j<=size(@l[1]);@j++)
299           {
300              if(deg(@l[1][@j])!=0)
301              {
302                 @k++;
303                 @i=@i+ideal(@l[1][@j]);
304                 if(size(@i)==pleh)
305                 {
306"factorization error";
307@l;
308                    @k--;
309                    @v[@k]=@v[@k]+@l[2][@j];
310                 }
311                 else
312                 {
313                    pleh++;
314                    @v[@k]=@l[2][@j];
315                 }
316              }
317           }
318           @l[1]=@i;
319           @l[2]=@v;
320        }
321     }
322  }
323  return(@l);
324}
325example
326{ "EXAMPLE:"; echo = 2;
327   ring  r = 0,(x,y,z),lp;
328   poly  p = (x+y)^2*(y-z)^3;
329   list  l = factor(p);
330   l;
331   ring r1 =(0,b,d,f,g),(a,c,e),lp;
332   poly p  =(1*d)*e^2+(1*d*f^2*g);
333   list  l = factor(p);
334   l;
335   ring r2 =(0,b,f,g),(a,c,e,d),lp;
336   poly p  =(1*d)*e^2+(1*d*f^2*g);
337   list  l = factor(p);
338   l;
339
340}
341
342
343
344////////////////////////////////////////////////////////////////////////////////
345
346proc idealsEqual( ideal k, ideal j)
347{
348   return(stdIdealsEqual(std(k),std(j)));
349}
350
351proc specialIdealsEqual( ideal k1, ideal k2)
352{
353   int j;
354
355   if(size(k1)==size(k2))
356   {
357      for(j=1;j<=size(k1);j++)
358      {
359         if(leadexp(k1[j])!=leadexp(k2[j]))
360         {
361            return(0);
362         }
363      }
364      return(1);
365   }
366   return(0);
367}
368
369proc stdIdealsEqual( ideal k1, ideal k2)
370{
371   int j;
372
373   if(size(k1)==size(k2))
374   {
375      for(j=1;j<=size(k1);j++)
376      {
377         if(leadexp(k1[j])!=leadexp(k2[j]))
378         {
379            return(0);
380         }
381      }
382      attrib(k2,"isSB",1);
383      if(size(reduce(k1,k2))==0)
384      {
385         return(1);
386      }
387   }
388   return(0);
389}
390
391
392////////////////////////////////////////////////////////////////////////////////
393
394proc testPrimary(list pr, ideal k)
395USAGE:   testPrimary(pr,k) pr list, k ideal;
396RETURN:  int = 1, if the intersection of the ideals in pr is k, 0 if not
397NOTE:
398EEXAMPLE: example testPrimary ; shows an example
399{
400   int i;
401   ideal j=pr[1];
402   for (i=2;i<=size(pr)/2;i++)
403   {
404       j=intersect(j,pr[2*i-1]);
405   }
406   return(idealsEqual(j,k));
407}
408example
409{ "EXAMPLE:"; echo = 2;
410   ring  s = 0,(x,y,z),lp;
411   ideal i=x3-x2-x+1,xy-y;
412   ideal i1=x-1;
413   ideal i2=x-1;
414   ideal i3=y,x2-2x+1;
415   ideal i4=y,x-1;
416   ideal i5=y,x+1;
417   ideal i6=y,x+1;
418   list pr=i1,i2,i3,i4,i5,i6;
419   testPrimary(pr,i);
420   pr[5]=y+1,x+1;
421   testPrimary(pr,i);
422}
423
424////////////////////////////////////////////////////////////////////////////////
425proc printPrimary( list l, list #)
426USAGE:   printPrimary(l) l list;
427RETURN:  nothing
428NOTE:
429EXAMPLE: example printPrimary; shows an example
430{
431   if(size(#)>0)
432   {
433      "                                            ";
434      "  The primary decomposition of the ideal    ";
435      #[1];
436      "                                            ";
437      "  is:                                       ";
438      "                                            ";
439   }
440   int k;
441   for (k=1;k<=size(l)/2;k=k+1)
442   {
443      "                                            ";
444      "primary ideal:                              ";
445      l[2*k-1];
446      "                                            ";
447      "associated prime ideal                      ";
448      l[2*k];
449      "                                            ";
450   }
451}
452example
453{ "EXAMPLE:"; echo = 2;
454}
455
456////////////////////////////////////////////////////////////////////////////////
457
458
459proc randomLast(int b)
460USAGE:   randomLast
461RETURN:  ideal = maxideal(1) but the last variable exchanged by
462         a sum of it with a linear random combination of the other
463         variables
464NOTE:
465EXAMPLE: example randomLast; shows an example
466{
467
468  ideal i=maxideal(1);
469  int k=size(i);
470  i[k]=0;
471  i=randomid(i,size(i),b);
472  ideal ires=maxideal(1);
473  ires[k]=i[1]+var(k);
474  return(ires);
475}
476example
477{ "EXAMPLE:"; echo = 2;
478   ring  r = 0,(x,y,z),lp;
479   ideal i = randomLast(10);
480   i;
481}
482
483
484////////////////////////////////////////////////////////////////////////////////
485
486
487proc primaryTest (ideal i, poly p)
488{
489   int m=1;
490   int n=nvars(basering);
491   int e;
492   poly t;
493   ideal h;
494
495   ideal prm=p;
496   attrib(prm,"isSB",1);
497
498   while (n>1)
499   {
500      n=n-1;
501      m=m+1;
502
503      //search for i[m] which has a power of var(n) as leading term
504      if (n==1)
505      {
506         m=size(i);
507      }
508      else
509      {
510         while (lead(i[m])/var(n-1)==0)
511        {
512            m=m+1;
513         }
514         m=m-1;
515      }
516      //check whether i[m] =(c*var(n)+h)^e modulo prm for some
517      //h in K[var(n+1),...,var(nvars(basering))], c in K
518      //if not (0) is returned, else var(n)+h is added to prm
519
520         e=deg(lead(i[m]));
521        // t=hilfe1(i,prm,m,n);
522         t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
523
524         i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
525
526         if (reduce(i[m]-t^e,prm) !=0)
527         {
528           return(ideal(0));
529         }
530         h=interred(t);
531         t=h[1];
532
533      prm = prm,t;
534      attrib(prm,"isSB",1);
535   }
536   return(prm);
537}
538
539proc hilfe(ideal i,ideal prm,int m)
540{
541   poly t;
542   int e;
543
544   if(size(i[m])==1)
545   {
546      t=var(n);
547   }
548   else
549   {
550      e=deg(lead(i[m]));
551
552      if(deg(poly(e))>=0)
553      {
554           t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
555           i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
556      {
557      else
558      {
559         i[m]=i[m]/leadcoef(i[m]);
560         t=reduce(coef(i[m],var(n))[2,e+1],prm);
561         t=var(n)+factorize(t,1)[1];
562      }
563   }
564   return(t);
565}
566proc hilfe1(ideal i,ideal prm,int m,int n)
567{
568   poly t;
569   int e;
570   if(size(i[m])==1)
571   {
572      t=var(n);
573   }
574   else
575   {
576      e=deg(lead(i[m]));
577      i[m]=i[m]/leadcoef(i[m]);
578      t=reduce(coeffs(i[m],var(n))[1,1],prm);
579      if(size(t)==0){return(var(n));}
580      t=var(n)+factorize(t,1)[1];
581   }
582   return(t);
583}
584
585///////////////////////////////////////////////////////////////////////////////
586proc splitPrimary(list l,ideal ser,int @wr,list sact)
587{
588   int i,j,k,s,r,w;
589   list keepresult,act,keepprime;
590   poly @f;
591   int sl=size(l);
592
593   for(i=1;i<=sl/2;i++)
594   {
595      if(sact[2][i]>1)
596      {
597         keepprime[i]=l[2*i-1]+ideal(sact[1][i]);
598      }
599      else
600      {
601         keepprime[i]=l[2*i-1];
602      }
603  }
604   i=0;
605   while(i<size(l)/2)
606   {
607      i++;
608      if((size(ser)>0)&&(size(reduce(ser,l[2*i-1]))==0))
609      {
610         l[2*i-1]=ideal(1);
611         l[2*i]=ideal(1);
612         continue;
613      }
614
615
616      if(size(l[2*i])==0)
617      {
618         if(homog(l[2*i-1])==1)
619         {
620            l[2*i]=maxideal(1);
621            continue;
622         }
623         j=0;
624         if(i<=sl/2)
625         {
626            j=1;
627         }
628         while(j<size(l[2*i-1]))
629         {
630            j++;
631            act=factor(l[2*i-1][j]);
632            r=size(act[1]);
633            attrib(l[2*i-1],"isSB",1);
634            if((r==1)&&(vdim(l[2*i-1])==deg(l[2*i-1][j])))
635            {
636              l[2*i]=std(l[2*i-1],act[1][1]);
637              break;
638            }
639            if((r==1)&&(act[2][1]>1))
640            {
641               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
642               if(homog(keepprime[i])==1)
643               {
644                  l[2*i]=maxideal(1);
645                  break;
646               }
647            }
648            if(gcdTest(act[1])==1)
649            {
650               for(k=2;k<=r;k++)
651               {
652                  keepprime[size(l)/2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
653               }
654               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
655               for(k=1;k<=r;k++)
656               {
657                  if(@wr==0)
658                  {
659                     keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]);
660                  }
661                  else
662                  {
663                     keepresult[k]=std(l[2*i-1],act[1][k]);
664                  }
665               }
666               l[2*i-1]=keepresult[1];
667               if(vdim(keepresult[1])==deg(act[1][1]))
668               {
669                  l[2*i]=keepresult[1];
670               }
671               if((homog(keepresult[1])==1)||(homog(keepprime[i])==1))
672               {
673                  l[2*i]=maxideal(1);
674               }
675               s=size(l)-2;
676               for(k=2;k<=r;k++)
677               {
678                  l[s+2*k-1]=keepresult[k];
679                  keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k]));
680                  if(vdim(keepresult[k])==deg(act[1][k]))
681                  {
682                     l[s+2*k]=keepresult[k];
683                  }
684                  else
685                  {
686                     l[s+2*k]=ideal(0);
687                  }
688                  if((homog(keepresult[k])==1)||(homog(keepprime[s/2+k])==1))
689                  {
690                    l[s+2*k]=maxideal(1);
691                  }
692               }
693               i--;
694               break;
695            }
696            if(r>=2)
697            {
698               s=size(l);
699               @f=act[1][1];
700               act=sat1(l[2*i-1],act[1][1]);
701               if(deg(act[1][1])>0)
702               {
703                  l[s+1]=std(l[2*i-1],act[2]);
704                  if(homog(l[s+1])==1)
705                  {
706                     l[s+2]=maxideal(1);
707                  }
708                  else
709                  {
710                     l[s+2]=ideal(0);
711                  }
712                  keepprime[s/2+1]=interred(keepprime[i]+ideal(@f));
713                  if(homog(keepprime[s/2+1])==1)
714                  {
715                     l[s+2]=maxideal(1);
716                  }
717                  keepprime[i]=act[1];
718                  l[2*i-1]=act[1];
719                  attrib(l[2*i-1],"isSB",1);
720                  if(homog(l[2*i-1])==1)
721                  {
722                     l[2*i]=maxideal(1);
723                  }
724
725                  i--;
726                  break;
727               }
728            }
729         }
730      }
731   }
732   if(sl==size(l))
733   {
734      return(l);
735   }
736   for(i=1;i<=size(l)/2;i++)
737   {
738      if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
739      {
740         keepprime[i]=std(keepprime[i]);
741         if(homog(keepprime[i])==1)
742         {
743             l[2*i]=maxideal(1);
744         }
745         else
746         {
747            act=zero_decomp(keepprime[i],ideal(0),@wr,1);
748            if(size(act)==2)
749            {
750               l[2*i]=act[2];
751            }
752         }
753      }
754   }
755   return(l);
756}
757example
758{ "EXAMPLE:"; echo=2;
759   LIB "primdec.lib";
760   ring  r = 32003,(x,y,z),lp;
761   ideal i1=x*(x+1),yz,(z+1)*(z-1);
762   ideal i2=xy,yz,(x-2)*(x+3);
763   list l=i1,ideal(0),i2,ideal(0),i2,ideal(1);
764   list l1=splitPrimary(l,ideal(0),0);
765   l1;
766}
767
768////////////////////////////////////////////////////////////////////////////////
769
770proc zero_decomp (ideal j,ideal ser,int @wr,list #)
771USAGE:   zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
772         (@wr=0 for primary decomposition, @wr=1 for computaion of associated
773         primes)
774RETURN:  list = list of primary ideals and their radicals (at even positions
775         in the list) if the input is zero-dimensional and a standardbases
776         with respect to lex-ordering
777         If ser!=(0) and ser is contained in j or if j is not zero-dimen-
778         sional then ideal(1),ideal(1) is returned
779NOTE:    Algorithm of Gianni, Traeger, Zacharias
780EXAMPLE: example zero_decomp; shows an example
781{
782  def   @P = basering;
783  int nva = nvars(basering);
784  int @k,@s,@n,@k1;
785  list primary,lres,act,@lh,@wh;
786  map phi,psi;
787  ideal jmap,helpprim,@qh,@qht;
788  intvec @vh,@hilb;
789  string @ri;
790  poly @f;
791
792  if (dim(j)>0)
793  {
794     primary[1]=ideal(1);
795     primary[2]=ideal(1);
796     return(primary);
797  }
798  j=interred(j);
799  attrib(j,"isSB",1);
800  if(vdim(j)==deg(j[1]))
801  {
802     if((size(ser)>0)&&(size(reduce(ser,j))==0))
803     {
804          primary[1]=ideal(1);
805          primary[2]=ideal(1);
806          return(primary);
807     }
808     act=factor(j[1]);
809     for(@k=1;@k<=size(act[1]);@k++)
810     {
811       @qh=j;
812       if(@wr==0)
813       {
814          @qh[1]=act[1][@k]^act[2][@k];
815       }
816       else
817       {
818          @qh[1]=act[1][@k];
819       }
820       primary[2*@k-1]=interred(@qh);
821       @qh=j;
822       @qh[1]=act[1][@k];
823       primary[2*@k]=interred(@qh);
824     }
825        return(primary);
826  }
827
828  if(homog(j)==1)
829  {
830     primary[1]=j;
831     if((size(ser)>0)&&(size(reduce(ser,j))==0))
832     {
833          primary[1]=ideal(1);
834          primary[2]=ideal(1);
835          return(primary);
836     }
837     if(dim(j)==-1)
838     {
839        primary[1]=ideal(1);
840        primary[2]=ideal(1);
841     }
842     else
843     {
844        primary[2]=maxideal(1);
845     }
846     return(primary);
847  }
848
849//the first element in the standardbase is factorized
850  if(deg(j[1])>0)
851  {
852    act=factor(j[1]);
853    testFactor(act,j[1]);
854  }
855  else
856  {
857     primary[1]=ideal(1);
858     primary[2]=ideal(1);
859     return(primary);
860  }
861
862//withe the factors new ideals (hopefully the primary decomposition)
863//are created
864
865  if(size(act[1])>1)
866  {
867     if(size(#)>1)
868     {
869        primary[1]=ideal(1);
870        primary[2]=ideal(1);
871        primary[3]=ideal(1);
872        primary[4]=ideal(1);
873        return(primary);
874     }
875     for(@k=1;@k<=size(act[1]);@k++)
876     {
877        if(@wr==0)
878        {
879           primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
880        }
881        else
882        {
883           primary[2*@k-1]=std(j,act[1][@k]);
884        }
885        if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
886        {
887           primary[2*@k]   = primary[2*@k-1];
888        }
889        else
890        {
891           primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
892        }
893     }
894  }
895  else
896  {
897     primary[1]=j;
898     if((size(#)>0)&&(act[2][1]>1))
899     {
900        act[2]=1;
901        primary[1]=std(primary[1],act[1][1]);
902     }
903
904     if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
905     {
906        primary[2]=primary[1];
907     }
908     else
909     {
910        primary[2]=primaryTest(primary[1],act[1][1]);
911     }
912  }
913  if(size(#)==0)
914  {
915     primary=splitPrimary(primary,ser,@wr,act);
916  }
917
918//test whether all ideals in the decomposition are primary and
919//in general position
920//if not after a random coordinate transformation of the last
921//variable the corresponding ideal is decomposed again.
922
923  @k=0;
924  int zz;
925  while(@k<(size(primary)/2))
926  {
927    @k++;
928    if (size(primary[2*@k])==0)
929    {
930       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
931       {
932          if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
933          {
934             primary[2*@k]=primary[2*@k-1];
935          }
936       }
937    }
938  }
939  @k=0;
940  while(@k<(size(primary)/2))
941  {
942    @k++;
943    if (size(primary[2*@k])==0)
944    {
945//      "the univariat polynomials";
946//       int qwe=timer;
947//       system("finduni",primary[2*@k-1]);
948
949       jmap=randomLast(100);
950//       timer-qwe;
951//       primary[2*@k-1];
952//       pause;
953       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
954       {
955          if(deg(lead(primary[2*@k-1][@n]))==1)
956          {
957             jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0);
958          }
959       }
960       phi=@P,jmap;
961       jmap[nva]=-(jmap[nva]-2*var(nva));
962       psi=@P,jmap;
963       @qht=primary[2*@k-1];
964       @qh=phi(@qht);
965       if(npars(@P)>0)
966       {
967          @ri= "ring @Phelp ="
968                  +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;";
969       }
970       else
971       {
972          @ri= "ring @Phelp ="
973                  +string(char(@P))+",("+varstr(@P)+",@t),dp;";
974       }
975       execute(@ri);
976       ideal @qh=homog(imap(@P,@qht),@t);
977
978       ideal @qh1=std(@qh);
979       @hilb=hilb(@qh1,1);
980       @ri= "ring @Phelp1 ="
981                  +string(char(@P))+",("+varstr(@Phelp)+"),lp;";
982       execute(@ri);
983       ideal @qh=homog(imap(@P,@qh),@t);
984       kill @Phelp;
985       @qh=std(@qh,@hilb);
986       @qh=subst(@qh,@t,1);
987       setring @P;
988       @qh=imap(@Phelp1,@qh);
989       kill @Phelp1;
990       @qh=clearSB(@qh);
991       attrib(@qh,"isSB",1);
992
993       @lh=zero_decomp (@qh,psi(ser),@wr);
994
995       kill lres;
996       list lres;
997       if(size(@lh)==2)
998       {
999          helpprim=@lh[2];
1000          lres[1]=primary[2*@k-1];
1001          lres[2]=psi(helpprim);
1002          if(size(reduce(lres[2],lres[1]))==0)
1003          {
1004             primary[2*@k]=primary[2*@k-1];
1005             continue;
1006          }
1007       }
1008       else
1009       {
1010          act=factor(@qh[1]);
1011          if(2*size(act[1])==size(@lh))
1012          {
1013             for(@n=1;@n<=size(act[1]);@n++)
1014             {
1015                @f=act[1][@n]^act[2][@n];
1016                lres[2*@n-1]=interred(primary[2*@k-1]+psi(@f));
1017                helpprim=@lh[2*@n];
1018                lres[2*@n]=psi(helpprim);
1019             }
1020          }
1021          else
1022          {
1023             lres=psi(@lh);
1024          }
1025       }
1026       if(npars(@P)>0)
1027       {
1028          @ri= "ring @Phelp ="
1029                  +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;";
1030       }
1031       else
1032       {
1033          @ri= "ring @Phelp ="
1034                  +string(char(@P))+",("+varstr(@P)+",@t),dp;";
1035       }
1036       execute(@ri);
1037       list @lvec;
1038       list @lr=imap(@P,lres);
1039       ideal @lr1;
1040
1041       if(size(@lr)==2)
1042       {
1043          @lr[2]=homog(@lr[2],@t);
1044          @lr1=std(@lr[2]);
1045          @lvec[2]=hilb(@lr1,1);
1046       }
1047       else
1048       {
1049          for(@n=1;@n<=size(@lr)/2;@n++)
1050          {
1051             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1052             {
1053                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1054                @lr1=std(@lr[2*@n-1]);
1055                @lvec[2*@n-1]=hilb(@lr1,1);
1056                @lvec[2*@n]=@lvec[2*@n-1];
1057             }
1058             else
1059             {
1060                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1061                @lr1=std(@lr[2*@n-1]);
1062                @lvec[2*@n-1]=hilb(@lr1,1);
1063                @lr[2*@n]=homog(@lr[2*@n],@t);
1064                @lr1=std(@lr[2*@n]);
1065                @lvec[2*@n]=hilb(@lr1,1);
1066
1067             }
1068         }
1069       }
1070       @ri= "ring @Phelp1 ="
1071                  +string(char(@P))+",("+varstr(@Phelp)+"),lp;";
1072       execute(@ri);
1073       list @lr=imap(@Phelp,@lr);
1074
1075       kill @Phelp;
1076       if(size(@lr)==2)
1077       {
1078          @lr[2]=std(@lr[2],@lvec[2]);
1079          @lr[2]=subst(@lr[2],@t,1);
1080
1081       }
1082       else
1083       {
1084          for(@n=1;@n<=size(@lr)/2;@n++)
1085          {
1086             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1087             {
1088                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1089                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1090                @lr[2*@n]=@lr[2*@n-1];
1091                attrib(@lr[2*@n],"isSB",1);
1092             }
1093             else
1094             {
1095                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1096                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1097                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
1098                @lr[2*@n]=subst(@lr[2*@n],@t,1);
1099             }
1100          }
1101       }
1102       kill @lvec;
1103       setring @P;
1104       lres=imap(@Phelp1,@lr);
1105       kill @Phelp1;
1106       for(@n=1;@n<=size(lres);@n++)
1107       {
1108          lres[@n]=clearSB(lres[@n]);
1109          attrib(lres[@n],"isSB",1);
1110       }
1111
1112       primary[2*@k-1]=lres[1];
1113       primary[2*@k]=lres[2];
1114       @s=size(primary)/2;
1115       for(@n=1;@n<=size(lres)/2-1;@n++)
1116       {
1117         primary[2*@s+2*@n-1]=lres[2*@n+1];
1118         primary[2*@s+2*@n]=lres[2*@n+2];
1119       }
1120       @k--;
1121     }
1122  }
1123  return(primary);
1124}
1125example
1126{ "EXAMPLE:"; echo = 2;
1127   ring  r = 0,(x,y,z),lp;
1128   poly  p = z2+1;
1129   poly  q = z4+2;
1130   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1131   i=std(i);
1132   list  pr= zero_decomp(i,ideal(0),0);
1133   pr;
1134}
1135
1136////////////////////////////////////////////////////////////////////////////////
1137
1138proc ggt (ideal i)
1139USAGE:   ggt(i); i list of polynomials
1140RETURN:  poly = ggt(i[1],...,i[size(i)])
1141NOTE:
1142EXAMPLE: example ggt; shows an example
1143{
1144  int k;
1145  poly p=i[1];
1146  if(deg(p)==0)
1147  {
1148    return(1);
1149  }
1150
1151
1152  for (k=2;k<=size(i);k++)
1153  {
1154     if(deg(i[k])==0)
1155     {
1156        return(1)
1157     }
1158     p=GCD(p,i[k]);
1159     if(deg(p)==0)
1160     {
1161        return(1);
1162     }
1163  }
1164  return(p);
1165}
1166example
1167{ "EXAMPLE:"; echo = 2;
1168   ring  r = 0,(x,y,z),lp;
1169   poly  p = (x+y)*(y+z);
1170   poly  q = (z4+2)*(y+z);
1171   ideal l=p,q;
1172   poly  pr= ggt(l);
1173   pr;
1174}
1175///////////////////////////////////////////////////////////////////////////////
1176proc gcdTest(ideal act)
1177{
1178  int i,j;
1179  if(size(act)<=1)
1180  {
1181     return(0);
1182  }
1183  for (i=1;i<=size(act)-1;i++)
1184  {
1185     for(j=i+1;j<=size(act);j++)
1186     {
1187        if(deg(std(ideal(act[i],act[j]))[1])>0)
1188        {
1189           return(0);
1190        }
1191     }
1192  }
1193  return(1);
1194}
1195
1196///////////////////////////////////////////////////////////////////////////////
1197proc coeffLcm(ideal h)
1198{
1199   string @pa=parstr(basering);
1200   if(size(@pa)==0)
1201   {
1202      return(lcmP(h));
1203   }
1204   def bsr= basering;
1205   string @id=string(h);
1206   execute "ring @r=0,("+@pa+","+varstr(bsr)+"),dp;";
1207   execute "ideal @i="+@id+";";
1208   poly @p=lcmP(@i);
1209   string @ps=string(@p);
1210   setring bsr;
1211   execute "poly @p="+@ps+";";
1212   return(@p);
1213}
1214example
1215{
1216   "EXAMPLE:"; echo = 2;
1217   ring  r =( 0,a,b),(x,y,z),lp;
1218   poly  p = (a+b)*(y-z);
1219   poly  q = (a+b)*(y+z);
1220   ideal l=p,q;
1221   poly  pr= coeffLcm(l);
1222   pr;
1223}
1224
1225///////////////////////////////////////////////////////////////////////////////
1226
1227proc lcmP(ideal i)
1228USAGE:   lcm(i); i list of polynomials
1229RETURN:  poly = lcm(i[1],...,i[size(i)])
1230NOTE:
1231EXAMPLE: example lcm; shows an example
1232{
1233  int k,j;
1234   poly p,q;
1235  i=simplify(i,10);
1236  for(j=1;j<=size(i);j++)
1237  {
1238    if(deg(i[j])>0)
1239    {
1240      p=i[j];
1241      break;
1242    }
1243  }
1244  if(deg(p)==-1)
1245  {
1246    return(1);
1247  }
1248  for (k=j+1;k<=size(i);k++)
1249  {
1250     if(deg(i[k])!=0)
1251     {
1252        q=GCD(p,i[k]);
1253        if(deg(q)==0)
1254        {
1255           p=p*i[k];
1256        }
1257        else
1258        {
1259           p=p/q;
1260           p=p*i[k];
1261        }
1262     }
1263   }
1264  return(p);
1265}
1266example
1267{ "EXAMPLE:"; echo = 2;
1268   ring  r = 0,(x,y,z),lp;
1269   poly  p = (x+y)*(y+z);
1270   poly  q = (z4+2)*(y+z);
1271   ideal l=p,q;
1272   poly  pr= lcmP(l);
1273   pr;
1274   l=1,-1,p,1,-1,q,1;
1275   pr=lcmP(l);
1276   pr;
1277}
1278
1279///////////////////////////////////////////////////////////////////////////////
1280proc clearSB (ideal i,list #)
1281USAGE:   clearSB(i); i ideal which is SB ordered by monomial ordering
1282RETURN:  ideal = minimal SB
1283NOTE:
1284EXAMPLE: example clearSB; shows an example
1285{
1286  int k,j;
1287  poly m;
1288  int c=size(i);
1289
1290  if(size(#)==0)
1291  {
1292    for(j=1;j<c;j++)
1293    {
1294      if(deg(i[j])==0)
1295      {
1296        i=ideal(1);
1297        return(i);
1298      }
1299      if(deg(i[j])>0)
1300      {
1301        m=lead(i[j]);
1302        for(k=j+1;k<=c;k++)
1303        {
1304           if(size(lead(i[k])/m)>0)
1305           {
1306             i[k]=0;
1307           }
1308        }
1309      }
1310    }
1311  }
1312  else
1313  {
1314    j=0;
1315    while(j<c-1)
1316    {
1317      j++;
1318      if(deg(i[j])==0)
1319      {
1320        i=ideal(1);
1321        return(i);
1322      }
1323      if(deg(i[j])>0)
1324      {
1325        m=lead(i[j]);
1326        for(k=j+1;k<=c;k++)
1327        {
1328           if(size(lead(i[k])/m)>0)
1329           {
1330             if((leadexp(m)!=leadexp(i[k]))||(#[j]<=#[k]))
1331             {
1332                i[k]=0;
1333             }
1334             else
1335             {
1336                i[j]=0;
1337                break;
1338             }
1339           }
1340        }
1341      }
1342    }
1343  }
1344  return(simplify(i,2));
1345}
1346example
1347{ "EXAMPLE:"; echo = 2;
1348   LIB "primdec.lib";
1349   ring  r = (0,a,b),(x,y,z),dp;
1350   ideal i=ax2+y,a2x+y,bx;
1351   list l=1,2,1;
1352   ideal j=clearSB(i,l);
1353   j;
1354}
1355
1356///////////////////////////////////////////////////////////////////////////////
1357
1358proc independSet (ideal j)
1359USAGE:   independentSet(i); i ideal
1360RETURN:  list = new varstring with the independent set at the end,
1361                ordstring with the corresponding block ordering,
1362                the integer where the independent set starts in the varstring
1363NOTE:
1364EXAMPLE: example independentSet; shows an example
1365{
1366   int n,k,di;
1367   list resu,hilf;
1368   string var1,var2;
1369   list v=system("indsetall",j,1);
1370
1371   for(n=1;n<=size(v);n++)
1372   {
1373      di=0;
1374      var1="";
1375      var2="";
1376      for(k=1;k<=size(v[n]);k++)
1377      {
1378         if(v[n][k]!=0)
1379         {
1380            di++;
1381            var2=var2+"var("+string(k)+"),";
1382         }
1383         else
1384         {
1385            var1=var1+"var("+string(k)+"),";
1386         }
1387      }
1388      if(di>0)
1389      {
1390         var1=var1+var2;
1391         var1=var1[1..size(var1)-1];
1392         hilf[1]=var1;
1393         hilf[2]="lp";
1394         //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
1395         hilf[3]=di;
1396         resu[n]=hilf;
1397      }
1398      else
1399      {
1400         resu[n]=varstr(basering),ordstr(basering),0;
1401      }
1402   }
1403   return(resu);
1404}
1405example
1406{ "EXAMPLE:"; echo = 2;
1407   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1408   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1409   i=std(i);
1410   list  l=independSet(i);
1411   l;
1412   i=i,g;
1413   l=independSet(i);
1414   l;
1415
1416   ring s=0,(x,y,z),lp;
1417   ideal i=z,yx;
1418   list l=independSet(i);
1419   l;
1420
1421
1422}
1423///////////////////////////////////////////////////////////////////////////////
1424
1425proc maxIndependSet (ideal j)
1426USAGE:   maxIndependentSet(i); i ideal
1427RETURN:  list = new varstring with the maximal independent set at the end,
1428                ordstring with the corresponding block ordering,
1429                the integer where the independent set starts in the varstring
1430NOTE:
1431EXAMPLE: example maxIndependentSet; shows an example
1432{
1433   int n,k,di;
1434   list resu,hilf;
1435   string var1,var2;
1436   list v=system("indsetall",j,0);
1437
1438   for(n=1;n<=size(v);n++)
1439   {
1440      di=0;
1441      var1="";
1442      var2="";
1443      for(k=1;k<=size(v[n]);k++)
1444      {
1445         if(v[n][k]!=0)
1446         {
1447            di++;
1448            var2=var2+"var("+string(k)+"),";
1449         }
1450         else
1451         {
1452            var1=var1+"var("+string(k)+"),";
1453         }
1454      }
1455      if(di>0)
1456      {
1457         var1=var1+var2;
1458         var1=var1[1..size(var1)-1];
1459         hilf[1]=var1;
1460         hilf[2]="lp";
1461         hilf[3]=di;
1462         resu[n]=hilf;
1463      }
1464      else
1465      {
1466         resu[n]=varstr(basering),ordstr(basering),0;
1467      }
1468   }
1469   return(resu);
1470}
1471example
1472{ "EXAMPLE:"; echo = 2;
1473   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1474   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1475   i=std(i);
1476   list  l=maxIndependSet(i);
1477   l;
1478   i=i,g;
1479   l=maxIndependSet(i);
1480   l;
1481
1482   ring s=0,(x,y,z),lp;
1483   ideal i=z,yx;
1484   list l=maxIndependSet(i);
1485   l;
1486
1487
1488}
1489
1490///////////////////////////////////////////////////////////////////////////////
1491
1492proc prepareQuotientring (int nnp)
1493USAGE:   prepareQuotientring(nnp); nnp int
1494RETURN:  string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
1495NOTE:
1496EXAMPLE: example independentSet; shows an example
1497{
1498  ideal @ih,@jh;
1499  int npar=npars(basering);
1500  int @n;
1501
1502  string quotring= "ring quring = ("+charstr(basering);
1503  for(@n=nnp+1;@n<=nvars(basering);@n++)
1504  {
1505     quotring=quotring+",var("+string(@n)+")";
1506     @ih=@ih+var(@n);
1507  }
1508
1509  quotring=quotring+"),(var(1)";
1510  @jh=@jh+var(1);
1511  for(@n=2;@n<=nnp;@n++)
1512  {
1513    quotring=quotring+",var("+string(@n)+")";
1514    @jh=@jh+var(@n);
1515  }
1516  quotring=quotring+"),lp;";
1517
1518  return(quotring);
1519
1520}
1521example
1522{ "EXAMPLE:"; echo = 2;
1523   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
1524   def @Q=basering;
1525   list l= prepareQuotientring(3);
1526   l;
1527   execute l[1];
1528   execute l[2];
1529   basering;
1530   phi;
1531   setring @Q;
1532
1533}
1534
1535///////////////////////////////////////////////////////////////////////
1536
1537proc projdim(list l)
1538{
1539   int i=size(l)+1;
1540
1541   while(i>2)
1542   {
1543      i--;
1544      if((size(l[i])>0)&&(deg(l[i][1])>0))
1545      {
1546         return(i);
1547      }
1548   }
1549}
1550
1551///////////////////////////////////////////////////////////////////////////////
1552proc cleanPrimary(list l)
1553{
1554   int i,j;
1555   list lh;
1556   for(i=1;i<=size(l)/2;i++)
1557   {
1558      if(deg(l[2*i-1][1])>0)
1559      {
1560         j++;
1561         lh[j]=l[2*i-1];
1562         j++;
1563         lh[j]=l[2*i];
1564      }
1565   }
1566   return(lh);
1567}
1568///////////////////////////////////////////////////////////////////////////////
1569
1570proc minAssPrimes(ideal i, list #)
1571USAGE:   minAssPrimes(i); i ideal
1572         minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
1573RETURN:  list = the minimal associated prime ideals of i
1574NOTE:
1575EXAMPLE: example minAssPrimes; shows an example
1576{
1577   #[1]=1;
1578   def @P=basering;
1579   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
1580             +ordstr(basering)+");";
1581
1582
1583   ideal i=fetch(@P,i);
1584   if(size(#)==0)
1585   {
1586      int @wr;
1587      list tluser,@res;
1588      list primary=decomp(i,2);
1589
1590      @res[1]=primary;
1591
1592      tluser=union(@res);
1593      setring @P;
1594      list @res=imap(gnir,tluser);
1595      return(@res);
1596   }
1597   list @res,empty;
1598   ideal ser;
1599   option(redSB);
1600   list @pr=facstd(i);
1601   if(size(@pr)==1)
1602   {
1603      attrib(@pr[1],"isSB",1);
1604      if((dim(@pr[1])==0)&&(homog(@pr[1])==1))
1605      {
1606         setring @P;
1607         list @res=maxideal(1);
1608         return(@res);
1609      }
1610      if(dim(@pr[1])>1)
1611      {
1612         setring @P;
1613         kill gnir;
1614         execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
1615         ideal i=fetch(@P,i);
1616         list @pr=facstd(i);
1617         ideal ser;
1618      }
1619   }
1620   option(noredSB);
1621   int j,k,odim,ndim,count;
1622   attrib(@pr[1],"isSB",1);
1623   if(#[1]==77)
1624   {
1625     odim=dim(@pr[1]);
1626     count=1;
1627     intvec pos;
1628     pos[size(@pr)]=0;
1629     for(j=2;j<=size(@pr);j++)
1630     {
1631        attrib(@pr[j],"isSB",1);
1632        ndim=dim(@pr[j]);
1633        if(ndim>odim)
1634        {
1635           for(k=count;k<=j-1;k++)
1636           {
1637              pos[k]=1;
1638           }
1639           count=j;
1640           odim=ndim;
1641        }
1642        if(ndim<odim)
1643        {
1644           pos[j]=1;
1645        }
1646     }
1647     for(j=1;j<=size(@pr);j++)
1648     {
1649        if(pos[j]!=1)
1650        {
1651            @res[j]=decomp(@pr[j],2);
1652        }
1653        else
1654        {
1655           @res[j]=empty;
1656        }
1657     }
1658   }
1659   else
1660   {
1661     ser=ideal(1);
1662     for(j=1;j<=size(@pr);j++)
1663     {
1664         @res[j]=decomp(@pr[j],2);
1665 //      @res[j]=decomp(@pr[j],2,@pr[j],ser);
1666 //      for(k=1;k<=size(@res[j]);k++)
1667 //      {
1668 //         ser=intersect1(ser,@res[j][k]);
1669 //      }
1670     }
1671   }
1672
1673   @res=union(@res);
1674   setring @P;
1675   list @res=imap(gnir,@res);
1676   return(@res);
1677}
1678example
1679{ "EXAMPLE:"; echo = 2;
1680   ring  r = 32003,(x,y,z),lp;
1681   poly  p = z2+1;
1682   poly  q = z4+2;
1683   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1684   LIB "primaryDecomposition.lib";
1685   list pr= minAssPrimes(i);
1686   pr;
1687   pr= minAssPrimes(i,1);
1688}
1689
1690///////////////////////////////////////////////////////////////////////////////
1691
1692proc union(list li)
1693{
1694  int i,j,k;
1695
1696  def P=basering;
1697
1698  execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
1699  list l=fetch(P,li);
1700  list @erg;
1701
1702  for(k=1;k<=size(l);k++)
1703  {
1704     for(j=1;j<=size(l[k])/2;j++)
1705     {
1706        if(deg(l[k][2*j][1])!=0)
1707        {
1708           i++;
1709           @erg[i]=l[k][2*j];
1710        }
1711     }
1712  }
1713
1714  list @wos;
1715  i=0;
1716  ideal i1,i2;
1717  while(i<size(@erg)-1)
1718  {
1719     i++;
1720     k=i+1;
1721     i1=lead(@erg[i]);
1722      attrib(i1,"isSB",1);
1723      attrib(@erg[i],"isSB",1);
1724
1725     while(k<=size(@erg))
1726     {
1727        if(deg(@erg[i][1])==0)
1728        {
1729           break;
1730        }
1731        i2=lead(@erg[k]);
1732        attrib(@erg[k],"isSB",1);
1733        attrib(i2,"isSB",1);
1734
1735        if(size(reduce(i1,i2,1))==0)
1736        {
1737           if(size(reduce(@erg[i],@erg[k]))==0)
1738           {
1739              @erg[k]=ideal(1);
1740              i2=ideal(1);
1741           }
1742        }
1743        if(size(reduce(i2,i1,1))==0)
1744        {
1745           if(size(reduce(@erg[k],@erg[i]))==0)
1746           {
1747              break;
1748           }
1749        }
1750        k++;
1751        if(k>size(@erg))
1752        {
1753           @wos[size(@wos)+1]=@erg[i];
1754        }
1755     }
1756  }
1757  if(deg(@erg[size(@erg)][1])!=0)
1758  {
1759     @wos[size(@wos)+1]=@erg[size(@erg)];
1760  }
1761  setring P;
1762  list @ser=fetch(ir,@wos);
1763  return(@ser);
1764}
1765///////////////////////////////////////////////////////////////////////////////
1766proc radicalOld(ideal i)
1767{
1768   list pr=minAssPrimes(i,1);
1769   int j;
1770   ideal k=pr[1];
1771   for(j=2;j<=size(pr);j++)
1772   {
1773      k=intersect(k,pr[j]);
1774   }
1775   return(k);
1776}
1777///////////////////////////////////////////////////////////////////////////////
1778proc equiRadical(ideal i)
1779{
1780   return(radical(i,1));
1781}
1782///////////////////////////////////////////////////////////////////////////////
1783proc decomp (ideal i,list #)
1784USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
1785         decomp(i,1);        (for the minimal associated primes) )
1786RETURN:  list = list of primary ideals and their associated primes
1787         (at even positions in the list)
1788         (resp. a list of the minimal associated primes)
1789NOTE:    Algorithm of Gianni, Traeger, Zacharias
1790EXAMPLE: example decomp; shows an example
1791{
1792  def  @P = basering;
1793  list primary,indep,ltras;
1794  intvec @vh,isat;
1795  int @wr,@k,@n,@m,@n1,@n2,@n3,homo;
1796  ideal peek=i;
1797  ideal ser,tras;
1798
1799  int @aa=timer;
1800
1801  homo=homog(i);
1802  if(size(#)>0)
1803  {
1804     if((#[1]==1)||(#[1]==2))
1805     {
1806        @wr=#[1];
1807        if(size(#)>1)
1808        {
1809          peek=#[2];
1810          ser=#[3];
1811        }
1812      }
1813      else
1814      {
1815         peek=#[1];
1816         ser=#[2];
1817      }
1818  }
1819
1820  if(homo==1)
1821  {
1822     ltras=mstd(i);
1823     attrib(ltras[1],"isSB",1);
1824     tras=ltras[1];
1825     if(dim(tras)==0)
1826     {
1827        primary[1]=ltras[2];
1828        primary[2]=maxideal(1);
1829        if(@wr>0)
1830        {
1831           list l;
1832           l[1]=maxideal(1);
1833           l[2]=maxideal(1);
1834           return(l);
1835        }
1836        return(primary);
1837     }
1838      intvec @hilb=hilb(tras,1);
1839  }
1840
1841  //----------------------------------------------------------------
1842  //i is the zero-ideal
1843  //----------------------------------------------------------------
1844
1845  if(size(i)==0)
1846  {
1847     primary=i,i;
1848     return(primary);
1849  }
1850
1851  //----------------------------------------------------------------
1852  //pass to the lexicographical ordering and compute a standardbasis
1853  //----------------------------------------------------------------
1854
1855  execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
1856
1857  option(redSB);
1858  ideal ser=fetch(@P,ser);
1859  ideal peek=std(fetch(@P,peek));
1860  homo=homog(peek);
1861
1862  if(homo==1)
1863  {
1864     if(ordstr(@P)[1,2]!="lp")
1865     {
1866        ideal @j=std(fetch(@P,i),@hilb);
1867     }
1868     else
1869     {
1870        ideal @j=fetch(@P,tras);
1871        attrib(@j,"isSB",1);
1872     }
1873  }
1874  else
1875  {
1876     ideal @j=std(fetch(@P,i));
1877  }
1878
1879  //----------------------------------------------------------------
1880  //j is the ring
1881  //----------------------------------------------------------------
1882
1883  if (dim(@j)==-1)
1884  {
1885     setring @P;
1886     option(noredSB);
1887     return(ideal(0));
1888  }
1889
1890  //----------------------------------------------------------------
1891  //  the case of one variable
1892  //----------------------------------------------------------------
1893
1894  if(nvars(basering)==1)
1895  {
1896
1897     list fac=factor(@j[1]);
1898     list gprimary;
1899     for(@k=1;@k<=size(fac[1]);@k++)
1900     {
1901        if(@wr==0)
1902        {
1903           gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
1904           gprimary[2*@k]=ideal(fac[1][@k]);
1905        }
1906        else
1907        {
1908           gprimary[2*@k-1]=ideal(fac[1][@k]);
1909           gprimary[2*@k]=ideal(fac[1][@k]);
1910        }
1911     }
1912     setring @P;
1913     option(noredSB);
1914     primary=fetch(gnir,gprimary);
1915
1916     return(primary);
1917  }
1918
1919 //------------------------------------------------------------------
1920 //the zero-dimensional case
1921 //------------------------------------------------------------------
1922
1923  if (dim(@j)==0)
1924  {
1925      list gprimary= zero_decomp(@j,ser,@wr);
1926
1927      setring @P;
1928      option(noredSB);
1929      primary=fetch(gnir,gprimary);
1930      if(size(ser)>0)
1931      {
1932         primary=cleanPrimary(primary);
1933      }
1934      return(primary);
1935   }
1936
1937  poly @gs,@gh,@p;
1938  string @va,quotring;
1939  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
1940  ideal @h;
1941  int jdim=dim(@j);
1942  list fett;
1943  int lauf,di;
1944  //------------------------------------------------------------------
1945  //search for a maximal independent set indep,i.e.
1946  //look for subring such that the intersection with the ideal is zero
1947  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
1948  //indep[1] is the new varstring and indep[2] the string for the block-ordering
1949  //------------------------------------------------------------------
1950
1951  if(@wr!=1)
1952  {
1953     allindep=independSet(@j);
1954     for(@m=1;@m<=size(allindep);@m++)
1955     {
1956        if(allindep[@m][3]==jdim)
1957        {
1958           di++;
1959           indep[di]=allindep[@m];
1960        }
1961        else
1962        {
1963           lauf++;
1964           restindep[lauf]=allindep[@m];
1965        }
1966     }
1967   }
1968   else
1969   {
1970      indep=maxIndependSet(@j);
1971   }
1972
1973  ideal jkeep=@j;
1974
1975  if(ordstr(@P)[1]=="w")
1976  {
1977     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");";
1978  }
1979  else
1980  {
1981     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),dp;";
1982  }
1983  ideal jwork=std(imap(gnir,@j));
1984  poly @p,@q;
1985  ideal @h,fac;
1986  di=dim(jwork);
1987  setring gnir;
1988  for(@m=1;@m<=size(indep);@m++)
1989  {
1990     isat=0;
1991     @n2=0;
1992     if((indep[@m][1]==varstr(basering))&&(@m==1))
1993     //this is the good case, nothing to do, just to have the same notations
1994     //change the ring
1995     {
1996        execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
1997                              +ordstr(basering)+");";
1998        ideal @j=fetch(gnir,@j);
1999        attrib(@j,"isSB",1);
2000     }
2001     else
2002     {
2003        @va=string(maxideal(1));
2004        execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2005                              +indep[@m][2]+");";
2006        execute "map phi=gnir,"+@va+";";
2007        if(homo==1)
2008        {
2009           ideal @j=std(phi(@j),@hilb);
2010        }
2011        else
2012        {
2013          ideal @j=std(phi(@j));
2014        }
2015      }
2016     if((deg(@j[1])==0)||(dim(@j)<jdim))
2017     {
2018       setring gnir;
2019       break;
2020     }
2021     for (lauf=1;lauf<=size(@j);lauf++)
2022     {
2023         fett[lauf]=size(@j[lauf]);
2024     }
2025     //------------------------------------------------------------------------------------
2026     //we have now the following situation:
2027     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2028     //to this quotientring, j is their still a standardbasis, the
2029     //leading coefficients of the polynomials  there (polynomials in
2030     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2031     //we need their ggt, gh, because of the following:
2032     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2033     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2034     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2035
2036     //------------------------------------------------------------------------------------
2037
2038     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2039     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
2040     //-------------------------------------------------------------------------------------
2041
2042     quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2043
2044     //---------------------------------------------------------------------
2045     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2046     //---------------------------------------------------------------------
2047
2048     execute quotring;
2049
2050    // @j considered in the quotientring
2051     ideal @j=imap(gnir1,@j);
2052     ideal ser=imap(gnir,ser);
2053
2054     kill gnir1;
2055
2056     //j is a standardbasis in the quotientring but usually not minimal
2057     //here it becomes minimal
2058
2059     @j=clearSB(@j,fett);
2060     attrib(@j,"isSB",1);
2061
2062     //we need later ggt(h[1],...)=gh for saturation
2063     ideal @h;
2064     if(deg(@j[1])>0)
2065     {
2066        for(@n=1;@n<=size(@j);@n++)
2067        {
2068           @h[@n]=leadcoef(@j[@n]);
2069        }
2070        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2071
2072        list uprimary= zero_decomp(@j,ser,@wr);
2073
2074     }
2075     else
2076     {
2077       list uprimary;
2078       uprimary[1]=ideal(1);
2079       uprimary[2]=ideal(1);
2080     }
2081
2082     //we need the intersection of the ideals in the list quprimary with the
2083     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2084     //but fi polynomials, then the intersection of q with the polynomialring
2085     //is the saturation of the ideal generated by f1,...,fr with respect to
2086     //h which is the lcm of the leading coefficients of the fi considered in the
2087     //quotientring: this is coded in saturn
2088
2089     list saturn;
2090     ideal hpl;
2091
2092     for(@n=1;@n<=size(uprimary);@n++)
2093     {
2094        hpl=0;
2095        for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2096        {
2097           hpl=hpl,leadcoef(uprimary[@n][@n1]);
2098        }
2099        saturn[@n]=hpl;
2100     }
2101     //--------------------------------------------------------------------
2102     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2103     //back to the polynomialring
2104     //---------------------------------------------------------------------
2105     setring gnir;
2106
2107     collectprimary=imap(quring,uprimary);
2108     lsau=imap(quring,saturn);
2109     @h=imap(quring,@h);
2110
2111     kill quring;
2112
2113
2114     @n2=size(quprimary);
2115     @n3=@n2;
2116
2117     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2118     {
2119        if(deg(collectprimary[2*@n1][1])>0)
2120        {
2121           @n2++;
2122           quprimary[@n2]=collectprimary[2*@n1-1];
2123           lnew[@n2]=lsau[2*@n1-1];
2124           @n2++;
2125           lnew[@n2]=lsau[2*@n1];
2126           quprimary[@n2]=collectprimary[2*@n1];
2127        }
2128     }
2129
2130     //here the intersection with the polynomialring
2131     //mentioned above is really computed
2132
2133    for(@n=@n3/2+1;@n<=@n2/2;@n++)
2134     {
2135        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
2136        {
2137           quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2138           quprimary[2*@n]=quprimary[2*@n-1];
2139        }
2140        else
2141        {
2142           if(@wr==0)
2143           {
2144              quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2145           }
2146           quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
2147        }
2148     }
2149     if(size(@h)>0)
2150     {
2151           //---------------------------------------------------------------
2152           //we change to @Phelp to have the ordering dp for saturation
2153           //---------------------------------------------------------------
2154        setring @Phelp;
2155        @h=imap(gnir,@h);
2156        if(@wr!=1)
2157//        if(@wr==0)
2158        {
2159           @q=minSat(jwork,@h)[2];
2160        }
2161        else
2162        {
2163            fac=ideal(0);
2164            for(lauf=1;lauf<=ncols(@h);lauf++)
2165           {
2166              if(deg(@h[lauf])>0)
2167              {
2168                 fac=fac+factorize(@h[lauf],1);
2169              }
2170           }
2171           fac=simplify(fac,4);
2172           @q=1;
2173           for(lauf=1;lauf<=size(fac);lauf++)
2174           {
2175             @q=@q*fac[lauf];
2176           }
2177        }
2178        jwork=std(jwork,@q);
2179        if(dim(jwork)<di)
2180        {
2181           setring gnir;
2182           @j=imap(@Phelp,jwork);
2183           break;
2184        }
2185        if(homo==1)
2186        {
2187              @hilb=hilb(jwork,1);
2188        }
2189
2190        setring gnir;
2191        @j=imap(@Phelp,jwork);
2192     }
2193  }
2194  if((size(quprimary)==0)&&(@wr>0))
2195  {
2196     @j=ideal(1);
2197     quprimary[1]=ideal(1);
2198     quprimary[2]=ideal(1);
2199  }
2200
2201  //---------------------------------------------------------------
2202  //notice that j=sat(j,gh) intersected with (j,gh^n)
2203  //we finished with sat(j,gh) and have to start with (j,gh^n)
2204  //---------------------------------------------------------------
2205  if((deg(@j[1])!=0)&&(@wr!=1))
2206  {
2207     int uq=size(quprimary);
2208     if(uq>0)
2209     {
2210        if(@wr==0)
2211        {
2212           ideal htest=quprimary[1];
2213
2214           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
2215           {
2216              htest=intersect(htest,quprimary[2*@n1-1]);
2217           }
2218        }
2219        else
2220        {
2221           ideal htest=quprimary[2];
2222
2223           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
2224           {
2225              htest=intersect(htest,quprimary[2*@n1]);
2226           }
2227        }
2228        if(size(ser)>0)
2229        {
2230           htest=intersect(htest,ser);
2231        }
2232        ser=std(htest);
2233     }
2234     //we are not ready yet
2235     if (specialIdealsEqual(ser,peek)!=1)
2236     {
2237        for(@m=1;@m<=size(restindep);@m++)
2238        {
2239           isat=0;
2240           @n2=0;
2241           if(restindep[@m][1]==varstr(basering))
2242           //this is the good case, nothing to do, just to have the same notations
2243           //change the ring
2244           {
2245              execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2246                              +ordstr(basering)+");";
2247              ideal @j=fetch(gnir,jkeep);
2248              attrib(@j,"isSB",1);
2249           }
2250           else
2251           {
2252              @va=string(maxideal(1));
2253              execute "ring gnir1 = ("+charstr(basering)+"),("+restindep[@m][1]+"),("
2254                              +restindep[@m][2]+");";
2255              execute "map phi=gnir,"+@va+";";
2256              if(homo==1)
2257              {
2258                 ideal @j=std(phi(jkeep),@hilb);
2259              }
2260              else
2261              {
2262                ideal @j=std(phi(jkeep));
2263              }
2264           }
2265
2266           for (lauf=1;lauf<=size(@j);lauf++)
2267           {
2268              fett[lauf]=size(@j[lauf]);
2269           }
2270           //------------------------------------------------------------------------------------
2271           //we have now the following situation:
2272           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2273           //to this quotientring, j is their still a standardbasis, the
2274           //leading coefficients of the polynomials  there (polynomials in
2275           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2276           //we need their ggt, gh, because of the following:
2277           //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2278           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2279           //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2280
2281           //------------------------------------------------------------------------------------
2282
2283           //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2284           //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
2285           //-------------------------------------------------------------------------------------
2286
2287           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
2288
2289           //---------------------------------------------------------------------
2290           //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2291           //---------------------------------------------------------------------
2292
2293           execute quotring;
2294
2295           // @j considered in the quotientring
2296           ideal @j=imap(gnir1,@j);
2297           ideal ser=imap(gnir,ser);
2298
2299           kill gnir1;
2300
2301           //j is a standardbasis in the quotientring but usually not minimal
2302           //here it becomes minimal
2303           @j=clearSB(@j,fett);
2304           attrib(@j,"isSB",1);
2305
2306           //we need later ggt(h[1],...)=gh for saturation
2307           ideal @h;
2308
2309           for(@n=1;@n<=size(@j);@n++)
2310           {
2311              @h[@n]=leadcoef(@j[@n]);
2312           }
2313           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2314
2315            list uprimary= zero_decomp(@j,ser,@wr);
2316
2317           //we need the intersection of the ideals in the list quprimary with the
2318           //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2319           //but fi polynomials, then the intersection of q with the polynomialring
2320           //is the saturation of the ideal generated by f1,...,fr with respect to
2321           //h which is the lcm of the leading coefficients of the fi considered in the
2322           //quotientring: this is coded in saturn
2323
2324           list saturn;
2325           ideal hpl;
2326
2327           for(@n=1;@n<=size(uprimary);@n++)
2328           {
2329              hpl=0;
2330              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2331              {
2332                 hpl=hpl,leadcoef(uprimary[@n][@n1]);
2333              }
2334               saturn[@n]=hpl;
2335           }
2336           //--------------------------------------------------------------------
2337           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2338           //back to the polynomialring
2339           //---------------------------------------------------------------------
2340           setring gnir;
2341
2342           collectprimary=imap(quring,uprimary);
2343           lsau=imap(quring,saturn);
2344           @h=imap(quring,@h);
2345
2346           kill quring;
2347
2348
2349           @n2=size(quprimary);
2350           @n3=@n2;
2351
2352           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2353           {
2354              if(deg(collectprimary[2*@n1][1])>0)
2355              {
2356                 @n2++;
2357                 quprimary[@n2]=collectprimary[2*@n1-1];
2358                 lnew[@n2]=lsau[2*@n1-1];
2359                 @n2++;
2360                 lnew[@n2]=lsau[2*@n1];
2361                 quprimary[@n2]=collectprimary[2*@n1];
2362              }
2363           }
2364
2365           //here the intersection with the polynomialring
2366           //mentioned above is really computed
2367
2368           for(@n=@n3/2+1;@n<=@n2/2;@n++)
2369           {
2370              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
2371              {
2372                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2373                 quprimary[2*@n]=quprimary[2*@n-1];
2374              }
2375              else
2376              {
2377                 if(@wr==0)
2378                 {
2379                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2380                 }
2381                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
2382             }
2383             if(@wr==0)
2384             {
2385                ser=std(intersect(ser,quprimary[2*@n-1]));
2386             }
2387             else
2388             {
2389                ser=std(intersect(ser,quprimary[2*@n]));
2390             }
2391           }
2392        }
2393        //we are not ready yet
2394        if (specialIdealsEqual(ser,peek)!=1)
2395        {
2396           if(@wr>0)
2397           {
2398              htprimary=decomp(@j,@wr,peek,ser);
2399           }
2400           else
2401           {
2402              htprimary=decomp(@j,peek,ser);
2403           }
2404           // here we collect now both results primary(sat(j,gh))
2405           // and primary(j,gh^n)
2406
2407           @n=size(quprimary);
2408           for (@k=1;@k<=size(htprimary);@k++)
2409           {
2410              quprimary[@n+@k]=htprimary[@k];
2411           }
2412        }
2413     }
2414   }
2415  //------------------------------------------------------------
2416  //back to the ring we started with
2417  //the final result: primary
2418  //------------------------------------------------------------
2419  setring @P;
2420  primary=imap(gnir,quprimary);
2421
2422  option(noredSB);
2423  return(primary);
2424}
2425
2426
2427example
2428{ "EXAMPLE:"; echo = 2;
2429   ring  r = 32003,(x,y,z),lp;
2430   poly  p = z2+1;
2431   poly  q = z4+2;
2432   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
2433   LIB "primdec.lib";
2434   list pr= decomp(i);
2435   pr;
2436   testPrimary( pr, i);
2437}
2438
2439///////////////////////////////////////////////////////////////////////////////
2440proc radical(ideal i,list #)
2441{
2442   def @P=basering;
2443   int j,il;
2444   if(size(#)>0)
2445   {
2446     il=#[1];
2447   }
2448   ideal re=1;
2449   option(redSB);
2450   list pr=facstd(i);
2451
2452   if(size(pr)==1)
2453   {
2454      attrib(pr[1],"isSB",1);
2455      if((dim(pr[1])==0)&&(homog(pr[1])==1))
2456      {
2457         ideal @res=maxideal(1);
2458         return(@res);
2459      }
2460      if(dim(pr[1])>1)
2461      {
2462         execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
2463         ideal i=fetch(@P,i);
2464         list @pr=facstd(i);
2465         setring @P;
2466         pr=fetch(gnir,@pr);
2467      }
2468   }
2469   option(noredSB);
2470   int s=size(pr);
2471   if(s==1)
2472   {
2473      return(radicalEHV(i,ideal(1),il));
2474   }
2475   intvec pos;
2476   pos[s]=0;
2477
2478   if(il==1)
2479   {
2480     int ndim,k;
2481     attrib(pr[1],"isSB",1);
2482     int odim=dim(pr[1]);
2483     int count=1;
2484
2485     for(j=2;j<=s;j++)
2486     {
2487        attrib(pr[j],"isSB",1);
2488        ndim=dim(pr[j]);
2489        if(ndim>odim)
2490        {
2491           for(k=count;k<=j-1;k++)
2492           {
2493              pos[k]=1;
2494           }
2495           count=j;
2496           odim=ndim;
2497        }
2498        if(ndim<odim)
2499        {
2500           pos[j]=1;
2501        }
2502     }
2503   }
2504
2505   for(j=1;j<=s;j++)
2506   {
2507      if(pos[j]==0)
2508      {
2509         re=intersect1(re,radicalEHV(pr[s+1-j],re,il));
2510      }
2511   }
2512   return(re);
2513}
2514
2515proc intersect1(ideal i,ideal j)
2516{
2517   def R=basering;
2518   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+",@t),dp;";
2519   ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
2520   ideal j=eliminate(i,var(nvars(basering)));
2521   setring R;
2522   map phi=gnir,maxideal(1);
2523   return(phi(j));
2524}
2525
2526
2527///////////////////////////////////////////////////////////////////////////////
2528proc radicalKL (list m,ideal ser,list #)
2529USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
2530         decomp(i,1);        (for the minimal associated primes) )
2531RETURN:  list = list of primary ideals and their associated primes
2532         (at even positions in the list)
2533         (resp. a list of the minimal associated primes)
2534NOTE:    Algorithm of Gianni, Traeger, Zacharias
2535EXAMPLE: example decomp; shows an example
2536{
2537   ideal i=m[2];
2538  //----------------------------------------------------------------
2539  //i is the zero-ideal
2540  //----------------------------------------------------------------
2541
2542   if(size(i)==0)
2543   {
2544      return(ideal(0));
2545   }
2546
2547   def  @P = basering;
2548   list indep,allindep,restindep,fett,@mu;
2549   intvec @vh,isat;
2550   int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di;
2551   ideal @j,@j1,fac,@h,collectrad,htrad,lsau;
2552   ideal rad=ideal(1);
2553   ideal te=ser;
2554   poly @p,@q;
2555   string @va,quotring;
2556   int  homo=homog(i);
2557
2558   if(size(#)>0)
2559   {
2560      @wr=#[1];
2561   }
2562   @j=m[1];
2563   @j1=m[2];
2564   int jdim=dim(@j);
2565   if(size(reduce(ser,@j))==0)
2566   {
2567      return(ser);
2568   }
2569   if(homo==1)
2570   {
2571      if(jdim==0)
2572      {
2573         option(noredSB);
2574         return(maxideal(1));
2575      }
2576      intvec @hilb=hilb(@j,1);
2577   }
2578
2579
2580  //----------------------------------------------------------------
2581  //j is the ring
2582  //----------------------------------------------------------------
2583
2584   if (jdim==-1)
2585   {
2586      option(noredSB);
2587      return(ideal(0));
2588   }
2589
2590  //----------------------------------------------------------------
2591  //  the case of one variable
2592  //----------------------------------------------------------------
2593
2594   if(nvars(basering)==1)
2595   {
2596      fac=factorize(@j[1],1);
2597      poly @p=1;
2598      for(@k=1;@k<=size(fac);@k++)
2599      {
2600         @p=@p*fac[@k];
2601      }
2602      option(noredSB);
2603
2604      return(ideal(@p));
2605   }
2606   //------------------------------------------------------------------
2607   //the case of a complete intersection
2608   //------------------------------------------------------------------
2609    if(jdim+size(@j1)==nvars(basering))
2610    {
2611      // ideal jac=minor(jacob(@j1),size(@j1));
2612      // return(quotient(@j1,jac));
2613    }
2614
2615   //------------------------------------------------------------------
2616   //the zero-dimensional case
2617   //------------------------------------------------------------------
2618
2619   if (jdim==0)
2620   {
2621      @j1=system("finduni",@j);
2622      for(@k=1;@k<=size(@j1);@k++)
2623      {
2624         fac=factorize(cleardenom(@j1[@k]),1);
2625         @p=fac[1];
2626         for(@n=2;@n<=size(fac);@n++)
2627         {
2628            @p=@p*fac[@n];
2629         }
2630         @j=@j,@p;
2631      }
2632      @j=std(@j);
2633      option(noredSB);
2634      return(@j);
2635   }
2636
2637   //------------------------------------------------------------------
2638   //search for a maximal independent set indep,i.e.
2639   //look for subring such that the intersection with the ideal is zero
2640   //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
2641   //indep[1] is the new varstring and indep[2] the string for the block-ordering
2642   //------------------------------------------------------------------
2643
2644   indep=maxIndependSet(@j);
2645
2646   di=dim(@j);
2647
2648   for(@m=1;@m<=size(indep);@m++)
2649   {
2650      if((indep[@m][1]==varstr(basering))&&(@m==1))
2651      //this is the good case, nothing to do, just to have the same notations
2652      //change the ring
2653      {
2654         execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2655                              +ordstr(basering)+");";
2656         ideal @j=fetch(@P,@j);
2657         attrib(@j,"isSB",1);
2658      }
2659      else
2660      {
2661         @va=string(maxideal(1));
2662         execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2663                              +indep[@m][2]+");";
2664         execute "map phi=@P,"+@va+";";
2665         if(homo==1)
2666         {
2667            ideal @j=std(phi(@j),@hilb);
2668         }
2669         else
2670         {
2671           ideal @j=std(phi(@j));
2672         }
2673      }
2674      if((deg(@j[1])==0)||(dim(@j)<jdim))
2675      {
2676         setring @P;
2677         break;
2678      }
2679      for (lauf=1;lauf<=size(@j);lauf++)
2680      {
2681         fett[lauf]=size(@j[lauf]);
2682      }
2683     //------------------------------------------------------------------------------------
2684     //we have now the following situation:
2685     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2686     //to this quotientring, j is their still a standardbasis, the
2687     //leading coefficients of the polynomials  there (polynomials in
2688     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2689     //we need their ggt, gh, because of the following:
2690     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2691     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2692     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2693
2694     //------------------------------------------------------------------------------------
2695
2696     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2697     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
2698     //-------------------------------------------------------------------------------------
2699
2700      quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2701
2702     //---------------------------------------------------------------------
2703     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2704     //---------------------------------------------------------------------
2705
2706      execute quotring;
2707
2708    // @j considered in the quotientring
2709      ideal @j=imap(gnir1,@j);
2710
2711      kill gnir1;
2712
2713     //j is a standardbasis in the quotientring but usually not minimal
2714     //here it becomes minimal
2715
2716      @j=clearSB(@j,fett);
2717      attrib(@j,"isSB",1);
2718
2719     //we need later ggt(h[1],...)=gh for saturation
2720      ideal @h;
2721      if(deg(@j[1])>0)
2722      {
2723         for(@n=1;@n<=size(@j);@n++)
2724         {
2725            @h[@n]=leadcoef(@j[@n]);
2726         }
2727        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2728         option(redSB);
2729         @j=interred(@j);
2730         attrib(@j,"isSB",1);
2731         list  @mo=@j,@j;
2732         ideal zero_rad= radicalKL(@mo,ideal(1));
2733      }
2734      else
2735      {
2736         ideal zero_rad=ideal(1);
2737      }
2738
2739     //we need the intersection of the ideals in the list quprimary with the
2740     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2741     //but fi polynomials, then the intersection of q with the polynomialring
2742     //is the saturation of the ideal generated by f1,...,fr with respect to
2743     //h which is the lcm of the leading coefficients of the fi considered in the
2744     //quotientring: this is coded in saturn
2745
2746      ideal hpl;
2747
2748      for(@n=1;@n<=size(zero_rad);@n++)
2749      {
2750         hpl=hpl,leadcoef(zero_rad[@n]);
2751      }
2752
2753     //--------------------------------------------------------------------
2754     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2755     //back to the polynomialring
2756     //---------------------------------------------------------------------
2757      setring @P;
2758
2759      collectrad=imap(quring,zero_rad);
2760      lsau=simplify(imap(quring,hpl),2);
2761      @h=imap(quring,@h);
2762
2763      kill quring;
2764
2765
2766     //here the intersection with the polynomialring
2767     //mentioned above is really computed
2768
2769      collectrad=sat2(collectrad,lsau)[1];
2770
2771      if(deg(@h[1])>=0)
2772      {
2773         fac=ideal(0);
2774         for(lauf=1;lauf<=ncols(@h);lauf++)
2775         {
2776            if(deg(@h[lauf])>0)
2777            {
2778                fac=fac+factorize(@h[lauf],1);
2779            }
2780         }
2781         fac=simplify(fac,4);
2782         @q=1;
2783         for(lauf=1;lauf<=size(fac);lauf++)
2784         {
2785            @q=@q*fac[lauf];
2786         }
2787
2788
2789         @mu=mstd(quotient(@j+ideal(@q),rad));
2790         @j=@mu[1];
2791         attrib(@j,"isSB",1);
2792
2793      }
2794      if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
2795      {
2796int xyz=timer;
2797"bei collecterad";
2798         rad=intersect(rad,collectrad);
2799timer-xyz;
2800      }
2801      else
2802      {
2803         if(deg(collectrad[1])>0)
2804         {
2805            rad=collectrad;
2806         }
2807      }
2808
2809      te=simplify(reduce(te*rad,@j),2);
2810
2811      if((dim(@j)<di)||(size(te)==0))
2812      {
2813         break;
2814      }
2815      if(homo==1)
2816      {
2817         @hilb=hilb(@j,1);
2818      }
2819   }
2820
2821   if(((@wr==1)&&(dim(@j)<di))||(deg(@j[1])==0)||(size(te)==0))
2822   {
2823      return(rad);
2824   }
2825   ideal tes=radicalKL(@mu,rad,@wr);
2826int sml=timer;
2827"bei rad";
2828   rad=intersect(rad,tes);
2829timer-sml;
2830  // rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
2831
2832
2833   option(noredSB);
2834   return(rad);
2835}
2836
2837
2838example
2839{ "EXAMPLE:"; echo = 2;
2840}
2841
2842
2843proc radicalEHV(ideal i,ideal re,list #)
2844{
2845   ideal J,I,I0,radI0,L,radI1,I2,radI2;
2846   int l,il;
2847   if(size(#)>0)
2848   {
2849      il=#[1];
2850   }
2851   option(redSB);
2852   list m=mstd(i);
2853        I=m[2];
2854   option(noredSB);
2855   if(size(reduce(re,m[1]))==0)
2856   {
2857      return(re);
2858   }
2859   int cod=nvars(basering)-dim(m[1]);
2860   if(nvars(basering)<9)
2861   {
2862      if(cod==size(m[2]))
2863      {
2864         J=minor(jacob(I),cod);
2865         return(quotient(I,J));
2866      }
2867
2868      for(l=1;l<=cod;l++)
2869      {
2870         I0[l]=I[l];
2871      }
2872      if(dim(std(I0))+cod==nvars(basering))
2873      {
2874         J=minor(jacob(I0),cod);
2875         radI0=quotient(I0,J);
2876         L=quotient(radI0,I);
2877         radI1=quotient(radI0,L);
2878
2879         if(size(reduce(radI1,m[1]))==0)
2880         {
2881            return(I);
2882         }
2883         if(il==1)
2884         {
2885            return(radI1);
2886         }
2887
2888         I2=sat(I,radI1)[1];
2889
2890         if(deg(I2[1])<=0)
2891         {
2892            return(radI1);
2893         }
2894         return(intersect(radI1,radicalEHV(I2,re,il)));
2895      }
2896   }
2897   return(radicalKL(m,re,il));
2898}
2899
2900proc Ann(module M)
2901{
2902  M=prune(M);  //to obtain a small embedding
2903  return(quotient(M,freemodule(nrows(M))));
2904}
2905
2906//computes the equidimensional part of the ideal i of codimension e
2907proc int_ass_primary_e(ideal i, int e)
2908{
2909  if(homog(i)!=1)
2910  {
2911     i=std(i);
2912  }
2913  list re=sres(i,0);                   //the resolution
2914  re=minres(re);                       //minimized resolution
2915  ideal ann=AnnExt_R(e,re);
2916  if(nvars(basering)-dim(std(ann))!=e)
2917  {
2918    return(ideal(1));
2919  }
2920  return(ann);
2921}
2922
2923//computes all equidimensional parts of the ideal i
2924proc prepareAss(ideal i)
2925{
2926  ideal j=std(i);
2927  int cod=nvars(basering)-dim(j);
2928  int e;
2929  list er;
2930  ideal ann;
2931  if(homog(i)==1)
2932  {
2933     list re=sres(i,0);                   //the resolution
2934     re=minres(re);                       //minimized resolution
2935  }
2936  else
2937  {
2938     list re=mres(i,0);             //fehler in sres
2939  }
2940  for(e=cod;e<=nvars(basering);e++)
2941  {
2942     ann=AnnExt_R(e,re);
2943     if(nvars(basering)-dim(std(ann))==e)
2944     {
2945        er[size(er)+1]=equiRadical(ann);
2946     }
2947  }
2948  return(er);
2949}
2950
2951//computes the annihilator of Ext^n(R/i,R) with given resolution re
2952//n is not necessarily the number of variables
2953proc AnnExt_R(int n,list re)
2954{
2955  if(n<nvars(basering))
2956  {
2957     matrix f=transpose(re[n+1]);      //Hom(_,R)
2958     module k=res(f,2)[2];             //the kernel
2959     matrix g=transpose(re[n]);        //the image of Hom(_,R)
2960     ideal ann=quotient(g,k);           //the anihilator
2961  }
2962  else
2963  {
2964     ideal ann=Ann(transpose(re[n]));
2965  }
2966  return(ann);
2967}
2968
Note: See TracBrowser for help on using the repository browser.