# source:git/Singular/LIB/primdec.lib@82716e

spielwiese
Last change on this file since 82716e was 82716e, checked in by Hans Schönemann <hannes@…>, 25 years ago
*hannes: typos in the info-help-string git-svn-id: file:///usr/local/Singular/svn/trunk@1773 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 100.6 KB
Line
1// \$Id: primdec.lib,v 1.18 1998-05-14 18:45:13 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// algorithms for primary decomposition based on
9// the ideas of Shimoyama/Yokoyama
10// written by Wolfram Decker and Hans Schoenemann
11//////////////////////////////////////////////////////
12
13version="\$Id: primdec.lib,v 1.18 1998-05-14 18:45:13 Singular Exp \$";
14info="
15LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I)
16
17  minAssPrimes (ideal I)
18  //computes the minimal associated primes of the ideal I
19
20  primdecGTZ (ideal I)
21  // Computes a complete primary decomposition via Gianni,Trager,Zacharias
22
24  //computes the radical of the ideal I
25
27  //computes the radical of the equidimensional part of the ideal I
28
29  prepareAss(ideal I)
30  //computes the radicals of the equidimensional parts of I
31
32  min_ass_prim_charsets (ideal I, int choose)
33  // minimal associated primes via the  characteristic set
34  // package written by Michael Messollen.
35  // The integer choose must be either 0 or 1.
36  // If choose=0, the given ordering of the variables is used.
37  // If choose=1, the system tries to find an \"optimal ordering\",
38  // which in some cases may considerably speed up the algorithm
39
40  // You may also may want to try one of the algorithms for
41  // minimal associated primes in the the library
42  // primdec.lib, written by Gerhard Pfister.
43  // These algorithms are variants of the algorithm
44  // of Gianni-Trager-Zacharias
45
46  primdecSY (ideal I, int choose)
47
48  // Computes a complete primary decomposition via
49  // a variant of the pseudoprimary approach of
50  // Shimoyama-Yokoyama.
51  // The integer choose must be either 0, 1, 2 or 3.
52  // If choose=0, min_ass_prim_charsets with the given
53  // ordering of the variables is used.
54  // If choose=1, min_ass_prim_charsets with the \"optimized\"
55  // ordering of the variables is used.
56  // If choose=2, minAssPrimes from primdec.lib is used
57  // If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
58";
59
60LIB "general.lib";
61LIB "elim.lib";
62LIB "poly.lib";
63LIB "random.lib";
64///////////////////////////////////////////////////////////////////////////////
65
66proc sat1 (ideal id, poly p)
67"USAGE:   sat1(id,j);  id ideal, j polynomial
68RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
69NOTE:    result is a std basis in the basering
70EXAMPLE: example sat; shows an example
71"
72{
73   int @k;
74   ideal inew=std(id);
75   ideal iold;
76   option(returnSB);
77   while(specialIdealsEqual(iold,inew)==0 )
78   {
79      iold=inew;
80      inew=quotient(iold,p);
81      @k++;
82   }
83   @k--;
84   option(noreturnSB);
85   list L =inew,p^@k;
86   return (L);
87}
88
89///////////////////////////////////////////////////////////////////////////////
90
91proc sat2 (ideal id, ideal h)
92"USAGE:   sat2(id,j);  id ideal, j polynomial
93RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
94NOTE:    result is a std basis in the basering
95EXAMPLE: example sat2; shows an example
96"
97{
98   int @k,@i;
99   def @P= basering;
100   if(ordstr(basering)[1,2]!="dp")
101   {
102      execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),(C,dp);";
103      ideal inew=std(imap(@P,id));
104      ideal  @h=imap(@P,h);
105   }
106   else
107   {
108      ideal @h=h;
109      ideal inew=std(id);
110   }
111   ideal fac;
112
113   for(@i=1;@i<=ncols(@h);@i++)
114   {
115     if(deg(@h[@i])>0)
116     {
117        fac=fac+factorize(@h[@i],1);
118     }
119   }
120   fac=simplify(fac,4);
121   poly @f=1;
122   if(deg(fac[1])>0)
123   {
124      ideal iold;
125
126      for(@i=1;@i<=size(fac);@i++)
127      {
128        @f=@f*fac[@i];
129      }
130      option(returnSB);
131      while(specialIdealsEqual(iold,inew)==0 )
132      {
133         iold=inew;
134         if(deg(iold[size(iold)])!=1)
135         {
136            inew=quotient(iold,@f);
137         }
138         else
139         {
140            inew=iold;
141         }
142         @k++;
143      }
144      option(noreturnSB);
145      @k--;
146   }
147
148   if(ordstr(@P)[1,2]!="dp")
149   {
150      setring @P;
151      ideal inew=std(imap(@Phelp,inew));
152      poly @f=imap(@Phelp,@f);
153   }
154   list L =inew,@f^@k;
155   return (L);
156}
157
158///////////////////////////////////////////////////////////////////////////////
159
160proc minSat(ideal inew, ideal h)
161{
162   int i,k;
163   poly f=1;
164   ideal iold,fac;
165   list quotM,l;
166
167   for(i=1;i<=ncols(h);i++)
168   {
169      if(deg(h[i])>0)
170      {
171         fac=fac+factorize(h[i],1);
172      }
173   }
174   fac=simplify(fac,4);
175   if(size(fac)==0)
176   {
177      l=inew,1;
178      return(l);
179   }
180   fac=sort(fac)[1];
181   for(i=1;i<=size(fac);i++)
182   {
183      f=f*fac[i];
184   }
185   quotM[1]=inew;
186   quotM[2]=fac;
187   quotM[3]=f;
188   f=1;
189   option(returnSB);
190   while(specialIdealsEqual(iold,quotM[1])==0)
191   {
192      if(k>0)
193      {
194         f=f*quotM[3];
195      }
196      iold=quotM[1];
197      quotM=quotMin(quotM);
198      k++;
199   }
200   option(noreturnSB);
201   l=quotM[1],f;
202   return(l);
203}
204
205proc quotMin(list tsil)
206{
207   int i,j,k,action;
208   ideal verg;
209   list l;
210   poly g;
211
212   ideal laedi=tsil[1];
213   ideal fac=tsil[2];
214   poly f=tsil[3];
215
216   ideal star=quotient(laedi,f);
217   action=1;
218
219   while(action==1)
220   {
221      if(size(fac)==1)
222      {
223         action=0;
224         break;
225      }
226      for(i=1;i<=size(fac);i++)
227      {
228        g=1;
229        verg=laedi;
230
231         for(j=1;j<=size(fac);j++)
232         {
233            if(i!=j)
234            {
235               g=g*fac[j];
236            }
237         }
238         verg=quotient(laedi,g);
239
240         if(specialIdealsEqual(verg,star)==1)
241         {
242            f=g;
243            fac[i]=0;
244            fac=simplify(fac,2);
245            break;
246         }
247         if(i==size(fac))
248         {
249            action=0;
250         }
251      }
252   }
253   l=star,fac,f;
254   return(l);
255}
256
257////////////////////////////////////////////////////////////////////////////////
258proc testFactor(list act,poly p)
259{
260  poly keep=p;
261
262   int i;
263   poly q=act[1][1]^act[2][1];
264   for(i=2;i<=size(act[1]);i++)
265   {
266      q=q*act[1][i]^act[2][i];
267   }
270   if(p-q!=0)
271   {
272      "ERROR IN FACTOR";
273      basering;
274
275      act;
276      keep;
277      pause;
278
279      p;
280      q;
281      pause;
282   }
283}
284////////////////////////////////////////////////////////////////////////////////
285
286proc factor(poly p)
287"USAGE:   factor(p) p poly
288RETURN:  list=;
289NOTE:
290EXAMPLE: example factor; shows an example
291"
292{
293
294  ideal @i;
295  list @l;
296  intvec @v,@w;
297  int @j,@k,@n;
298
299  if(deg(p)<=1)
300  {
301     @i=ideal(p);
302     @v=1;
303     @l[1]=@i;
304     @l[2]=@v;
305     return(@l);
306  }
307  if (size(p)==1)
308  {
310     for(@j=1;@j<=nvars(basering);@j++)
311     {
312        if(@w[@j]!=0)
313        {
314           @k++;
315           @v[@k]=@w[@j];
316           @i=@i+ideal(var(@j));
317        }
318     }
319     @l[1]=@i;
320     @l[2]=@v;
321     return(@l);
322  }
323//  @l=factorize(p,2);
324  @l=factorize(p);
325  // if(npars(basering)>0)
326  // {
327     for(@j=1;@j<=size(@l[1]);@j++)
328     {
329        if(deg(@l[1][@j])==0)
330        {
331           @n++;
332        }
333     }
334     if(@n>0)
335     {
336        if(@n==size(@l[1]))
337        {
338           @l[1]=ideal(1);
339           @v=1;
340           @l[2]=@v;
341        }
342        else
343        {
344           @k=0;
345           int pleh;
346           for(@j=1;@j<=size(@l[1]);@j++)
347           {
348              if(deg(@l[1][@j])!=0)
349              {
350                 @k++;
351                 @i=@i+ideal(@l[1][@j]);
352                 if(size(@i)==pleh)
353                 {
354"factorization error";
355@l;
356                    @k--;
357                    @v[@k]=@v[@k]+@l[2][@j];
358                 }
359                 else
360                 {
361                    pleh++;
362                    @v[@k]=@l[2][@j];
363                 }
364              }
365           }
366           @l[1]=@i;
367           @l[2]=@v;
368        }
369     }
370     // }
371  return(@l);
372}
373example
374{ "EXAMPLE:"; echo = 2;
375   ring  r = 0,(x,y,z),lp;
376   poly  p = (x+y)^2*(y-z)^3;
377   list  l = factor(p);
378   l;
379   ring r1 =(0,b,d,f,g),(a,c,e),lp;
380   poly p  =(1*d)*e^2+(1*d*f^2*g);
381   list  l = factor(p);
382   l;
383   ring r2 =(0,b,f,g),(a,c,e,d),lp;
384   poly p  =(1*d)*e^2+(1*d*f^2*g);
385   list  l = factor(p);
386   l;
387
388}
389
390
391
392////////////////////////////////////////////////////////////////////////////////
393
394proc idealsEqual( ideal k, ideal j)
395{
396   return(stdIdealsEqual(std(k),std(j)));
397}
398
399proc specialIdealsEqual( ideal k1, ideal k2)
400{
401   int j;
402
403   if(size(k1)==size(k2))
404   {
405      for(j=1;j<=size(k1);j++)
406      {
408         {
409            return(0);
410         }
411      }
412      return(1);
413   }
414   return(0);
415}
416
417proc stdIdealsEqual( ideal k1, ideal k2)
418{
419   int j;
420
421   if(size(k1)==size(k2))
422   {
423      for(j=1;j<=size(k1);j++)
424      {
426         {
427            return(0);
428         }
429      }
430      attrib(k2,"isSB",1);
431      if(size(reduce(k1,k2,1))==0)
432      {
433         return(1);
434      }
435   }
436   return(0);
437}
438
439
440////////////////////////////////////////////////////////////////////////////////
441
442proc testPrimary(list pr, ideal k)
443"USAGE:   testPrimary(pr,k) pr list, k ideal;
444RETURN:  int = 1, if the intersection of the ideals in pr is k, 0 if not
445NOTE:
446EXAMPLE: example testPrimary ; shows an example
447"
448{
449   int i;
450   ideal j=pr[1];
451   for (i=2;i<=size(pr)/2;i++)
452   {
453       j=intersect(j,pr[2*i-1]);
454   }
455   return(idealsEqual(j,k));
456}
457example
458{ "EXAMPLE:"; echo = 2;
459   ring  s = 0,(x,y,z),lp;
460   ideal i=x3-x2-x+1,xy-y;
461   ideal i1=x-1;
462   ideal i2=x-1;
463   ideal i3=y,x2-2x+1;
464   ideal i4=y,x-1;
465   ideal i5=y,x+1;
466   ideal i6=y,x+1;
467   list pr=i1,i2,i3,i4,i5,i6;
468   testPrimary(pr,i);
469   pr[5]=y+1,x+1;
470   testPrimary(pr,i);
471}
472
473////////////////////////////////////////////////////////////////////////////////
474proc printPrimary( list l, list #)
475"USAGE:   printPrimary(l) l list;
476RETURN:  nothing
477NOTE:
478EXAMPLE: example printPrimary; shows an example
479"
480{
481   if(size(#)>0)
482   {
483      "                                            ";
484      "  The primary decomposition of the ideal    ";
485      #[1];
486      "                                            ";
487      "  is:                                       ";
488      "                                            ";
489   }
490   int k;
491   for (k=1;k<=size(l)/2;k=k+1)
492   {
493      "                                            ";
494      "primary ideal:                              ";
495      l[2*k-1];
496      "                                            ";
497      "associated prime ideal                      ";
498      l[2*k];
499      "                                            ";
500   }
501}
502example
503{ "EXAMPLE:"; echo = 2;
504}
505
506////////////////////////////////////////////////////////////////////////////////
507
508
509proc randomLast(int b)
510"USAGE:   randomLast
511RETURN:  ideal = maxideal(1) but the last variable exchanged by
512         a sum of it with a linear random combination of the other
513         variables
514NOTE:
515EXAMPLE: example randomLast; shows an example
516"
517{
518
519  ideal i=maxideal(1);
520  int k=size(i);
521  i[k]=0;
522  i=randomid(i,size(i),b);
523  ideal ires=maxideal(1);
524  ires[k]=i[1]+var(k);
525  return(ires);
526}
527example
528{ "EXAMPLE:"; echo = 2;
529   ring  r = 0,(x,y,z),lp;
530   ideal i = randomLast(10);
531   i;
532}
533
534
535////////////////////////////////////////////////////////////////////////////////
536
537
538proc primaryTest (ideal i, poly p)
539{
540   int m=1;
541   int n=nvars(basering);
542   int e;
543   poly t;
544   ideal h;
545
546   ideal prm=p;
547   attrib(prm,"isSB",1);
548
549   while (n>1)
550   {
551      n=n-1;
552      m=m+1;
553
554      //search for i[m] which has a power of var(n) as leading term
555      if (n==1)
556      {
557         m=size(i);
558      }
559      else
560      {
562        {
563            m=m+1;
564         }
565         m=m-1;
566      }
567      //check whether i[m] =(c*var(n)+h)^e modulo prm for some
568      //h in K[var(n+1),...,var(nvars(basering))], c in K
569      //if not (0) is returned, else var(n)+h is added to prm
570
572        // t=hilfe1(i,prm,m,n);
574
576
577         if (reduce(i[m]-t^e,prm,1) !=0)
578         {
579           return(ideal(0));
580         }
581         h=interred(t);
582         t=h[1];
583
584      prm = prm,t;
585      attrib(prm,"isSB",1);
586   }
587   return(prm);
588}
589
590proc hilfe(ideal i,ideal prm,int m)
591{
592   poly t;
593   int e;
594
595   if(size(i[m])==1)
596   {
597      t=var(n);
598   }
599   else
600   {
602
603      if(deg(poly(e))>=0)
604      {
607      }
608      else
609      {
611         t=reduce(coef(i[m],var(n))[2,e+1],prm);
612         t=var(n)+factorize(t,1)[1];
613      }
614   }
615   return(t);
616}
617proc hilfe1(ideal i,ideal prm,int m,int n)
618{
619   poly t;
620   int e;
621   if(size(i[m])==1)
622   {
623      t=var(n);
624   }
625   else
626   {
629      t=reduce(coeffs(i[m],var(n))[1,1],prm);
630      if(size(t)==0){return(var(n));}
631      t=var(n)+factorize(t,1)[1];
632   }
633   return(t);
634}
635
636///////////////////////////////////////////////////////////////////////////////
637proc splitPrimary(list l,ideal ser,int @wr,list sact)
638{
639   int i,j,k,s,r,w;
640   list keepresult,act,keepprime;
641   poly @f;
642   int sl=size(l);
643
644   for(i=1;i<=sl/2;i++)
645   {
646      if(sact[2][i]>1)
647      {
648         keepprime[i]=l[2*i-1]+ideal(sact[1][i]);
649      }
650      else
651      {
652         keepprime[i]=l[2*i-1];
653      }
654  }
655   i=0;
656   while(i<size(l)/2)
657   {
658      i++;
659      if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0))
660      {
661         l[2*i-1]=ideal(1);
662         l[2*i]=ideal(1);
663         continue;
664      }
665
666
667      if(size(l[2*i])==0)
668      {
669         if(homog(l[2*i-1])==1)
670         {
671            l[2*i]=maxideal(1);
672            continue;
673         }
674         j=0;
675         if(i<=sl/2)
676         {
677            j=1;
678         }
679         while(j<size(l[2*i-1]))
680         {
681            j++;
682            act=factor(l[2*i-1][j]);
683            r=size(act[1]);
684            attrib(l[2*i-1],"isSB",1);
685            if((r==1)&&(vdim(l[2*i-1])==deg(l[2*i-1][j])))
686            {
687              l[2*i]=std(l[2*i-1],act[1][1]);
688              break;
689            }
690            if((r==1)&&(act[2][1]>1))
691            {
692               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
693               if(homog(keepprime[i])==1)
694               {
695                  l[2*i]=maxideal(1);
696                  break;
697               }
698            }
699            if(gcdTest(act[1])==1)
700            {
701               for(k=2;k<=r;k++)
702               {
703                  keepprime[size(l)/2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
704               }
705               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
706               for(k=1;k<=r;k++)
707               {
708                  if(@wr==0)
709                  {
710                     keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]);
711                  }
712                  else
713                  {
714                     keepresult[k]=std(l[2*i-1],act[1][k]);
715                  }
716               }
717               l[2*i-1]=keepresult[1];
718               if(vdim(keepresult[1])==deg(act[1][1]))
719               {
720                  l[2*i]=keepresult[1];
721               }
722               if((homog(keepresult[1])==1)||(homog(keepprime[i])==1))
723               {
724                  l[2*i]=maxideal(1);
725               }
726               s=size(l)-2;
727               for(k=2;k<=r;k++)
728               {
729                  l[s+2*k-1]=keepresult[k];
730                  keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k]));
731                  if(vdim(keepresult[k])==deg(act[1][k]))
732                  {
733                     l[s+2*k]=keepresult[k];
734                  }
735                  else
736                  {
737                     l[s+2*k]=ideal(0);
738                  }
739                  if((homog(keepresult[k])==1)||(homog(keepprime[s/2+k])==1))
740                  {
741                    l[s+2*k]=maxideal(1);
742                  }
743               }
744               i--;
745               break;
746            }
747            if(r>=2)
748            {
749               s=size(l);
750               @f=act[1][1];
751               act=sat1(l[2*i-1],act[1][1]);
752               if(deg(act[1][1])>0)
753               {
754                  l[s+1]=std(l[2*i-1],act[2]);
755                  if(homog(l[s+1])==1)
756                  {
757                     l[s+2]=maxideal(1);
758                  }
759                  else
760                  {
761                     l[s+2]=ideal(0);
762                  }
763                  keepprime[s/2+1]=interred(keepprime[i]+ideal(@f));
764                  if(homog(keepprime[s/2+1])==1)
765                  {
766                     l[s+2]=maxideal(1);
767                  }
768                  keepprime[i]=act[1];
769                  l[2*i-1]=act[1];
770                  attrib(l[2*i-1],"isSB",1);
771                  if(homog(l[2*i-1])==1)
772                  {
773                     l[2*i]=maxideal(1);
774                  }
775
776                  i--;
777                  break;
778               }
779            }
780         }
781      }
782   }
783   if(sl==size(l))
784   {
785      return(l);
786   }
787   for(i=1;i<=size(l)/2;i++)
788   {
789     attrib(l[2*i-1],"isSB",1);
790
791     if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0)&&(deg(l[2*i-1][1])>0))
792     {
793       "Achtung in split";
794
795         l[2*i-1]=ideal(1);
796         l[2*i]=ideal(1);
797     }
798     if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
799      {
800         keepprime[i]=std(keepprime[i]);
801         if(homog(keepprime[i])==1)
802         {
803             l[2*i]=maxideal(1);
804         }
805         else
806         {
807            act=zero_decomp(keepprime[i],ideal(0),@wr,1);
808            if(size(act)==2)
809            {
810               l[2*i]=act[2];
811            }
812         }
813      }
814   }
815   return(l);
816}
817example
818{ "EXAMPLE:"; echo=2;
819   LIB "primdec.lib";
820   ring  r = 32003,(x,y,z),lp;
821   ideal i1=x*(x+1),yz,(z+1)*(z-1);
822   ideal i2=xy,yz,(x-2)*(x+3);
823   list l=i1,ideal(0),i2,ideal(0),i2,ideal(1);
824   list l1=splitPrimary(l,ideal(0),0);
825   l1;
826}
827
828////////////////////////////////////////////////////////////////////////////////
829
830proc zero_decomp (ideal j,ideal ser,int @wr,list #)
831"USAGE:   zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
832         (@wr=0 for primary decomposition, @wr=1 for computaion of associated
833         primes)
834RETURN:  list = list of primary ideals and their radicals (at even positions
835         in the list) if the input is zero-dimensional and a standardbases
836         with respect to lex-ordering
837         If ser!=(0) and ser is contained in j or if j is not zero-dimen-
838         sional then ideal(1),ideal(1) is returned
839NOTE:    Algorithm of Gianni, Traeger, Zacharias
840EXAMPLE: example zero_decomp; shows an example
841"
842{
843  def   @P = basering;
844  int nva = nvars(basering);
845  int @k,@s,@n,@k1,zz;
846  list primary,lres0,lres1,act,@lh,@wh;
847  map phi,psi,phi1,psi1;
848  ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1;
849  intvec @vh,@hilb;
850  string @ri;
851  poly @f;
852
853  if (dim(j)>0)
854  {
855     primary[1]=ideal(1);
856     primary[2]=ideal(1);
857     return(primary);
858  }
859
860  j=interred(j);
861  attrib(j,"isSB",1);
862  if(vdim(j)==deg(j[1]))
863  {
864     act=factor(j[1]);
865     for(@k=1;@k<=size(act[1]);@k++)
866     {
867       @qh=j;
868       if(@wr==0)
869       {
870          @qh[1]=act[1][@k]^act[2][@k];
871       }
872       else
873       {
874          @qh[1]=act[1][@k];
875       }
876       primary[2*@k-1]=interred(@qh);
877       @qh=j;
878       @qh[1]=act[1][@k];
879       primary[2*@k]=interred(@qh);
880       attrib( primary[2*@k-1],"isSB",1);
881
882       if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
883       {
884          primary[2*@k-1]=ideal(1);
885          primary[2*@k]=ideal(1);
886       }
887     }
888     return(primary);
889  }
890
891  if(homog(j)==1)
892  {
893     primary[1]=j;
894     if((size(ser)>0)&&(size(reduce(ser,j,1))==0))
895     {
896          primary[1]=ideal(1);
897          primary[2]=ideal(1);
898          return(primary);
899     }
900     if(dim(j)==-1)
901     {
902        primary[1]=ideal(1);
903        primary[2]=ideal(1);
904     }
905     else
906     {
907        primary[2]=maxideal(1);
908     }
909     return(primary);
910  }
911
912//the first element in the standardbase is factorized
913  if(deg(j[1])>0)
914  {
915    act=factor(j[1]);
916    testFactor(act,j[1]);
917  }
918  else
919  {
920     primary[1]=ideal(1);
921     primary[2]=ideal(1);
922     return(primary);
923  }
924
925//withe the factors new ideals (hopefully the primary decomposition)
926//are created
927
928  if(size(act[1])>1)
929  {
930     if(size(#)>1)
931     {
932        primary[1]=ideal(1);
933        primary[2]=ideal(1);
934        primary[3]=ideal(1);
935        primary[4]=ideal(1);
936        return(primary);
937     }
938     for(@k=1;@k<=size(act[1]);@k++)
939     {
940        if(@wr==0)
941        {
942           primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
943        }
944        else
945        {
946           primary[2*@k-1]=std(j,act[1][@k]);
947        }
948        if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
949        {
950           primary[2*@k]   = primary[2*@k-1];
951        }
952        else
953        {
954           primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
955        }
956     }
957  }
958  else
959  {
960     primary[1]=j;
961     if((size(#)>0)&&(act[2][1]>1))
962     {
963        act[2]=1;
964        primary[1]=std(primary[1],act[1][1]);
965     }
966
967     if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
968     {
969        primary[2]=primary[1];
970     }
971     else
972     {
973        primary[2]=primaryTest(primary[1],act[1][1]);
974     }
975  }
976  if(size(#)==0)
977  {
978
979     primary=splitPrimary(primary,ser,@wr,act);
980
981  }
982
983//test whether all ideals in the decomposition are primary and
984//in general position
985//if not after a random coordinate transformation of the last
986//variable the corresponding ideal is decomposed again.
987
988  @k=0;
989  while(@k<(size(primary)/2))
990  {
991    @k++;
992    if (size(primary[2*@k])==0)
993    {
994       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
995       {
996          if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
997          {
998             primary[2*@k]=primary[2*@k-1];
999          }
1000       }
1001    }
1002  }
1003
1004  @k=0;
1005  ideal keep;
1006  while(@k<(size(primary)/2))
1007  {
1008    @k++;
1009    if (size(primary[2*@k])==0)
1010    {
1011
1012       jmap=randomLast(100);
1013       jmap1=maxideal(1);
1014       jmap2=maxideal(1);
1015       @qht=primary[2*@k-1];
1016
1017       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
1018       {
1020          {
1021             for(zz=1;zz<=nva;zz++)
1022             {
1024                {
1027                    jmap2[zz]=primary[2*@k-1][@n];
1028                    @qht[@n]=var(zz);
1029
1030                }
1031             }
1033          }
1034       }
1035       if(size(subst(jmap[nva],var(1),0)-var(nva))!=0)
1036       {
1037          jmap[nva]=subst(jmap[nva],var(1),0);
1038       }
1039       phi1=@P,jmap1;
1040       phi=@P,jmap;
1041
1042       for(@n=1;@n<=nva;@n++)
1043       {
1044          jmap[@n]=-(jmap[@n]-2*var(@n));
1045       }
1046       psi=@P,jmap;
1047       psi1=@P,jmap2;
1048
1049       @qh=phi(@qht);
1050
1051       if(npars(@P)>0)
1052       {
1053          @ri= "ring @Phelp ="
1054                  +string(char(@P))+",
1055                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
1056       }
1057       else
1058       {
1059          @ri= "ring @Phelp ="
1060                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
1061       }
1062       execute(@ri);
1063       ideal @qh=homog(imap(@P,@qht),@t);
1064
1065       ideal @qh1=std(@qh);
1066       @hilb=hilb(@qh1,1);
1067       @ri= "ring @Phelp1 ="
1068                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
1069       execute(@ri);
1070       ideal @qh=homog(imap(@P,@qh),@t);
1071       kill @Phelp;
1072       @qh=std(@qh,@hilb);
1073       @qh=subst(@qh,@t,1);
1074       setring @P;
1075       @qh=imap(@Phelp1,@qh);
1076       kill @Phelp1;
1077       @qh=clearSB(@qh);
1078       attrib(@qh,"isSB",1);
1079       ser1=phi1(ser);
1080       @lh=zero_decomp (@qh,phi(ser1),@wr);
1081//       @lh=zero_decomp (@qh,psi(ser),@wr);
1082
1083
1084       kill lres0;
1085       list lres0;
1086       if(size(@lh)==2)
1087       {
1088          helpprim=@lh[2];
1089          lres0[1]=primary[2*@k-1];
1090          ser1=psi(helpprim);
1091          lres0[2]=psi1(ser1);
1092          if(size(reduce(lres0[2],lres0[1],1))==0)
1093          {
1094             primary[2*@k]=primary[2*@k-1];
1095             continue;
1096          }
1097       }
1098       else
1099       {
1100          act=factor(@qh[1]);
1101          if(2*size(act[1])==size(@lh))
1102          {
1103             for(@n=1;@n<=size(act[1]);@n++)
1104             {
1105                @f=act[1][@n]^act[2][@n];
1106                ser1=psi(@f);
1107                lres0[2*@n-1]=interred(primary[2*@k-1]+psi1(ser1));
1108                helpprim=@lh[2*@n];
1109                ser1=psi(helpprim);
1110                lres0[2*@n]=psi1(ser1);
1111             }
1112          }
1113          else
1114          {
1115             lres1=psi(@lh);
1116             lres0=psi1(lres1);
1117          }
1118       }
1119       if(npars(@P)>0)
1120       {
1121          @ri= "ring @Phelp ="
1122                  +string(char(@P))+",
1123                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
1124       }
1125       else
1126       {
1127          @ri= "ring @Phelp ="
1128                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
1129       }
1130       execute(@ri);
1131       list @lvec;
1132       list @lr=imap(@P,lres0);
1133       ideal @lr1;
1134
1135       if(size(@lr)==2)
1136       {
1137          @lr[2]=homog(@lr[2],@t);
1138          @lr1=std(@lr[2]);
1139          @lvec[2]=hilb(@lr1,1);
1140       }
1141       else
1142       {
1143          for(@n=1;@n<=size(@lr)/2;@n++)
1144          {
1145             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1146             {
1147                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1148                @lr1=std(@lr[2*@n-1]);
1149                @lvec[2*@n-1]=hilb(@lr1,1);
1150                @lvec[2*@n]=@lvec[2*@n-1];
1151             }
1152             else
1153             {
1154                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1155                @lr1=std(@lr[2*@n-1]);
1156                @lvec[2*@n-1]=hilb(@lr1,1);
1157                @lr[2*@n]=homog(@lr[2*@n],@t);
1158                @lr1=std(@lr[2*@n]);
1159                @lvec[2*@n]=hilb(@lr1,1);
1160
1161             }
1162         }
1163       }
1164       @ri= "ring @Phelp1 ="
1165                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
1166       execute(@ri);
1167       list @lr=imap(@Phelp,@lr);
1168
1169       kill @Phelp;
1170       if(size(@lr)==2)
1171       {
1172          @lr[2]=std(@lr[2],@lvec[2]);
1173          @lr[2]=subst(@lr[2],@t,1);
1174
1175       }
1176       else
1177       {
1178          for(@n=1;@n<=size(@lr)/2;@n++)
1179          {
1180             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1181             {
1182                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1183                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1184                @lr[2*@n]=@lr[2*@n-1];
1185                attrib(@lr[2*@n],"isSB",1);
1186             }
1187             else
1188             {
1189                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1190                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1191                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
1192                @lr[2*@n]=subst(@lr[2*@n],@t,1);
1193             }
1194          }
1195       }
1196       kill @lvec;
1197       setring @P;
1198       lres0=imap(@Phelp1,@lr);
1199       kill @Phelp1;
1200       for(@n=1;@n<=size(lres0);@n++)
1201       {
1202          lres0[@n]=clearSB(lres0[@n]);
1203          attrib(lres0[@n],"isSB",1);
1204       }
1205
1206       primary[2*@k-1]=lres0[1];
1207       primary[2*@k]=lres0[2];
1208       @s=size(primary)/2;
1209       for(@n=1;@n<=size(lres0)/2-1;@n++)
1210       {
1211         primary[2*@s+2*@n-1]=lres0[2*@n+1];
1212         primary[2*@s+2*@n]=lres0[2*@n+2];
1213       }
1214       @k--;
1215     }
1216  }
1217  return(primary);
1218}
1219example
1220{ "EXAMPLE:"; echo = 2;
1221   ring  r = 0,(x,y,z),lp;
1222   poly  p = z2+1;
1223   poly  q = z4+2;
1224   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1225   i=std(i);
1226   list  pr= zero_decomp(i,ideal(0),0);
1227   pr;
1228}
1229
1230////////////////////////////////////////////////////////////////////////////////
1231
1232proc ggt (ideal i)
1233"USAGE:   ggt(i); i list of polynomials
1234RETURN:  poly = ggt(i[1],...,i[size(i)])
1235NOTE:
1236EXAMPLE: example ggt; shows an example
1237"
1238{
1239  int k;
1240  poly p=i[1];
1241  if(deg(p)==0)
1242  {
1243    return(1);
1244  }
1245
1246
1247  for (k=2;k<=size(i);k++)
1248  {
1249     if(deg(i[k])==0)
1250     {
1251        return(1)
1252     }
1253     p=GCD(p,i[k]);
1254     if(deg(p)==0)
1255     {
1256        return(1);
1257     }
1258  }
1259  return(p);
1260}
1261example
1262{ "EXAMPLE:"; echo = 2;
1263   ring  r = 0,(x,y,z),lp;
1264   poly  p = (x+y)*(y+z);
1265   poly  q = (z4+2)*(y+z);
1266   ideal l=p,q;
1267   poly  pr= ggt(l);
1268   pr;
1269}
1270///////////////////////////////////////////////////////////////////////////////
1271proc gcdTest(ideal act)
1272{
1273  int i,j;
1274  if(size(act)<=1)
1275  {
1276     return(0);
1277  }
1278  for (i=1;i<=size(act)-1;i++)
1279  {
1280     for(j=i+1;j<=size(act);j++)
1281     {
1282        if(deg(std(ideal(act[i],act[j]))[1])>0)
1283        {
1284           return(0);
1285        }
1286     }
1287  }
1288  return(1);
1289}
1290
1291///////////////////////////////////////////////////////////////////////////////
1292proc coeffLcm(ideal h)
1293{
1294   string @pa=parstr(basering);
1295   if(size(@pa)==0)
1296   {
1297      return(lcmP(h));
1298   }
1299   def bsr= basering;
1300   string @id=string(h);
1301   execute "ring @r=0,("+@pa+","+varstr(bsr)+"),(C,dp);";
1302   execute "ideal @i="+@id+";";
1303   poly @p=lcmP(@i);
1304   string @ps=string(@p);
1305   setring bsr;
1306   execute "poly @p="+@ps+";";
1307   return(@p);
1308}
1309example
1310{
1311   "EXAMPLE:"; echo = 2;
1312   ring  r =( 0,a,b),(x,y,z),lp;
1313   poly  p = (a+b)*(y-z);
1314   poly  q = (a+b)*(y+z);
1315   ideal l=p,q;
1316   poly  pr= coeffLcm(l);
1317   pr;
1318}
1319
1320///////////////////////////////////////////////////////////////////////////////
1321
1322proc lcmP(ideal i)
1323"USAGE:   lcm(i); i list of polynomials
1324RETURN:  poly = lcm(i[1],...,i[size(i)])
1325NOTE:
1326EXAMPLE: example lcm; shows an example
1327"
1328{
1329  int k,j;
1330   poly p,q;
1331  i=simplify(i,10);
1332  for(j=1;j<=size(i);j++)
1333  {
1334    if(deg(i[j])>0)
1335    {
1336      p=i[j];
1337      break;
1338    }
1339  }
1340  if(deg(p)==-1)
1341  {
1342    return(1);
1343  }
1344  for (k=j+1;k<=size(i);k++)
1345  {
1346     if(deg(i[k])!=0)
1347     {
1348        q=GCD(p,i[k]);
1349        if(deg(q)==0)
1350        {
1351           p=p*i[k];
1352        }
1353        else
1354        {
1355           p=p/q;
1356           p=p*i[k];
1357        }
1358     }
1359   }
1360  return(p);
1361}
1362example
1363{ "EXAMPLE:"; echo = 2;
1364   ring  r = 0,(x,y,z),lp;
1365   poly  p = (x+y)*(y+z);
1366   poly  q = (z4+2)*(y+z);
1367   ideal l=p,q;
1368   poly  pr= lcmP(l);
1369   pr;
1370   l=1,-1,p,1,-1,q,1;
1371   pr=lcmP(l);
1372   pr;
1373}
1374
1375///////////////////////////////////////////////////////////////////////////////
1376proc clearSB (ideal i,list #)
1377"USAGE:   clearSB(i); i ideal which is SB ordered by monomial ordering
1378RETURN:  ideal = minimal SB
1379NOTE:
1380EXAMPLE: example clearSB; shows an example
1381"
1382{
1383  int k,j;
1384  poly m;
1385  int c=size(i);
1386
1387  if(size(#)==0)
1388  {
1389    for(j=1;j<c;j++)
1390    {
1391      if(deg(i[j])==0)
1392      {
1393        i=ideal(1);
1394        return(i);
1395      }
1396      if(deg(i[j])>0)
1397      {
1399        for(k=j+1;k<=c;k++)
1400        {
1402           {
1403             i[k]=0;
1404           }
1405        }
1406      }
1407    }
1408  }
1409  else
1410  {
1411    j=0;
1412    while(j<c-1)
1413    {
1414      j++;
1415      if(deg(i[j])==0)
1416      {
1417        i=ideal(1);
1418        return(i);
1419      }
1420      if(deg(i[j])>0)
1421      {
1423        for(k=j+1;k<=c;k++)
1424        {
1426           {
1428             {
1429                i[k]=0;
1430             }
1431             else
1432             {
1433                i[j]=0;
1434                break;
1435             }
1436           }
1437        }
1438      }
1439    }
1440  }
1441  return(simplify(i,2));
1442}
1443example
1444{ "EXAMPLE:"; echo = 2;
1445   LIB "primdec.lib";
1446   ring  r = (0,a,b),(x,y,z),dp;
1447   ideal i=ax2+y,a2x+y,bx;
1448   list l=1,2,1;
1449   ideal j=clearSB(i,l);
1450   j;
1451}
1452
1453///////////////////////////////////////////////////////////////////////////////
1454
1455proc independSet (ideal j)
1456"USAGE:   independentSet(i); i ideal
1457RETURN:  list = new varstring with the independent set at the end,
1458                ordstring with the corresponding block ordering,
1459                the integer where the independent set starts in the varstring
1460NOTE:
1461EXAMPLE: example independentSet; shows an example
1462"
1463{
1464   int n,k,di;
1465   list resu,hilf;
1466   string var1,var2;
1467   list v=indepSet(j,1);
1468
1469   for(n=1;n<=size(v);n++)
1470   {
1471      di=0;
1472      var1="";
1473      var2="";
1474      for(k=1;k<=size(v[n]);k++)
1475      {
1476         if(v[n][k]!=0)
1477         {
1478            di++;
1479            var2=var2+"var("+string(k)+"),";
1480         }
1481         else
1482         {
1483            var1=var1+"var("+string(k)+"),";
1484         }
1485      }
1486      if(di>0)
1487      {
1488         var1=var1+var2;
1489         var1=var1[1..size(var1)-1];
1490         hilf[1]=var1;
1491         hilf[2]="lp";
1492         //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
1493         hilf[3]=di;
1494         resu[n]=hilf;
1495      }
1496      else
1497      {
1498         resu[n]=varstr(basering),ordstr(basering),0;
1499      }
1500   }
1501   return(resu);
1502}
1503example
1504{ "EXAMPLE:"; echo = 2;
1505   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1506   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1507   i=std(i);
1508   list  l=independSet(i);
1509   l;
1510   i=i,g;
1511   l=independSet(i);
1512   l;
1513
1514   ring s=0,(x,y,z),lp;
1515   ideal i=z,yx;
1516   list l=independSet(i);
1517   l;
1518
1519
1520}
1521///////////////////////////////////////////////////////////////////////////////
1522
1523proc maxIndependSet (ideal j)
1524"USAGE:   maxIndependentSet(i); i ideal
1525RETURN:  list = new varstring with the maximal independent set at the end,
1526                ordstring with the corresponding block ordering,
1527                the integer where the independent set starts in the varstring
1528NOTE:
1529EXAMPLE: example maxIndependentSet; shows an example
1530"
1531{
1532   int n,k,di;
1533   list resu,hilf;
1534   string var1,var2;
1535   list v=indepSet(j,0);
1536
1537   for(n=1;n<=size(v);n++)
1538   {
1539      di=0;
1540      var1="";
1541      var2="";
1542      for(k=1;k<=size(v[n]);k++)
1543      {
1544         if(v[n][k]!=0)
1545         {
1546            di++;
1547            var2=var2+"var("+string(k)+"),";
1548         }
1549         else
1550         {
1551            var1=var1+"var("+string(k)+"),";
1552         }
1553      }
1554      if(di>0)
1555      {
1556         var1=var1+var2;
1557         var1=var1[1..size(var1)-1];
1558         hilf[1]=var1;
1559         hilf[2]="lp";
1560         hilf[3]=di;
1561         resu[n]=hilf;
1562      }
1563      else
1564      {
1565         resu[n]=varstr(basering),ordstr(basering),0;
1566      }
1567   }
1568   return(resu);
1569}
1570example
1571{ "EXAMPLE:"; echo = 2;
1572   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1573   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1574   i=std(i);
1575   list  l=maxIndependSet(i);
1576   l;
1577   i=i,g;
1578   l=maxIndependSet(i);
1579   l;
1580
1581   ring s=0,(x,y,z),lp;
1582   ideal i=z,yx;
1583   list l=maxIndependSet(i);
1584   l;
1585
1586
1587}
1588
1589///////////////////////////////////////////////////////////////////////////////
1590
1591proc prepareQuotientring (int nnp)
1592"USAGE:   prepareQuotientring(nnp); nnp int
1593RETURN:  string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
1594NOTE:
1595EXAMPLE: example independentSet; shows an example
1596"
1597{
1598  ideal @ih,@jh;
1599  int npar=npars(basering);
1600  int @n;
1601
1602  string quotring= "ring quring = ("+charstr(basering);
1603  for(@n=nnp+1;@n<=nvars(basering);@n++)
1604  {
1605     quotring=quotring+",var("+string(@n)+")";
1606     @ih=@ih+var(@n);
1607  }
1608
1609  quotring=quotring+"),(var(1)";
1610  @jh=@jh+var(1);
1611  for(@n=2;@n<=nnp;@n++)
1612  {
1613    quotring=quotring+",var("+string(@n)+")";
1614    @jh=@jh+var(@n);
1615  }
1616  quotring=quotring+"),(C,lp);";
1617
1618  return(quotring);
1619
1620}
1621example
1622{ "EXAMPLE:"; echo = 2;
1623   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
1624   def @Q=basering;
1625   list l= prepareQuotientring(3);
1626   l;
1627   execute l[1];
1628   execute l[2];
1629   basering;
1630   phi;
1631   setring @Q;
1632
1633}
1634
1635///////////////////////////////////////////////////////////////////////
1636
1637proc projdim(list l)
1638{
1639   int i=size(l)+1;
1640
1641   while(i>2)
1642   {
1643      i--;
1644      if((size(l[i])>0)&&(deg(l[i][1])>0))
1645      {
1646         return(i);
1647      }
1648   }
1649}
1650
1651///////////////////////////////////////////////////////////////////////////////
1652proc cleanPrimary(list l)
1653{
1654   int i,j;
1655   list lh;
1656   for(i=1;i<=size(l)/2;i++)
1657   {
1658      if(deg(l[2*i-1][1])>0)
1659      {
1660         j++;
1661         lh[j]=l[2*i-1];
1662         j++;
1663         lh[j]=l[2*i];
1664      }
1665   }
1666   return(lh);
1667}
1668///////////////////////////////////////////////////////////////////////////////
1669
1670proc minAssPrimes(ideal i, list #)
1671"USAGE:   minAssPrimes(i); i ideal
1672         minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
1673RETURN:  list = the minimal associated prime ideals of i
1674NOTE:
1675EXAMPLE: example minAssPrimes; shows an example
1676"
1677{
1678   #[1]=1;
1679   def @P=basering;
1680   list qr=simplifyIdeal(i);
1681   map phi=@P,qr[2];
1682   i=qr[1];
1683
1684   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
1685             +ordstr(basering)+");";
1686
1687
1688   ideal i=fetch(@P,i);
1689   if(size(#)==0)
1690   {
1691      int @wr;
1692      list tluser,@res;
1693      list primary=decomp(i,2);
1694
1695      @res[1]=primary;
1696
1697      tluser=union(@res);
1698      setring @P;
1699      list @res=imap(gnir,tluser);
1700      return(phi(@res));
1701   }
1702   list @res,empty;
1703   ideal ser;
1704   option(redSB);
1705   list @pr=facstd(i);
1706   if(size(@pr)==1)
1707   {
1708      attrib(@pr[1],"isSB",1);
1709      if((dim(@pr[1])==0)&&(homog(@pr[1])==1))
1710      {
1711         setring @P;
1712         list @res=maxideal(1);
1713         return(phi(@res));
1714      }
1715      if(dim(@pr[1])>1)
1716      {
1717         setring @P;
1718        // kill gnir;
1719         execute "ring gnir1 = ("+charstr(basering)+"),
1720                              ("+varstr(basering)+"),(C,lp);";
1721         ideal i=fetch(@P,i);
1722         list @pr=facstd(i);
1723        // ideal ser;
1724         setring gnir;
1725         @pr=fetch(gnir1,@pr);
1726         kill gnir1;
1727      }
1728   }
1729    option(noredSB);
1730   int j,k,odim,ndim,count;
1731   attrib(@pr[1],"isSB",1);
1732   if(#[1]==77)
1733   {
1734     odim=dim(@pr[1]);
1735     count=1;
1736     intvec pos;
1737     pos[size(@pr)]=0;
1738     for(j=2;j<=size(@pr);j++)
1739     {
1740        attrib(@pr[j],"isSB",1);
1741        ndim=dim(@pr[j]);
1742        if(ndim>odim)
1743        {
1744           for(k=count;k<=j-1;k++)
1745           {
1746              pos[k]=1;
1747           }
1748           count=j;
1749           odim=ndim;
1750        }
1751        if(ndim<odim)
1752        {
1753           pos[j]=1;
1754        }
1755     }
1756     for(j=1;j<=size(@pr);j++)
1757     {
1758        if(pos[j]!=1)
1759        {
1760            @res[j]=decomp(@pr[j],2);
1761        }
1762        else
1763        {
1764           @res[j]=empty;
1765        }
1766     }
1767   }
1768   else
1769   {
1770     ser=ideal(1);
1771     for(j=1;j<=size(@pr);j++)
1772     {
1773//@pr[j];
1774//pause;
1775        @res[j]=decomp(@pr[j],2);
1776//       @res[j]=decomp(@pr[j],2,@pr[j],ser);
1777//       for(k=1;k<=size(@res[j]);k++)
1778//       {
1779//          ser=intersect(ser,@res[j][k]);
1780//       }
1781     }
1782   }
1783
1784   @res=union(@res);
1785   setring @P;
1786   list @res=imap(gnir,@res);
1787   return(phi(@res));
1788}
1789example
1790{ "EXAMPLE:"; echo = 2;
1791   ring  r = 32003,(x,y,z),lp;
1792   poly  p = z2+1;
1793   poly  q = z4+2;
1794   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1795   LIB "primaryDecomposition.lib";
1796   list pr= minAssPrimes(i);
1797   pr;
1798   pr= minAssPrimes(i,1);
1799}
1800
1801///////////////////////////////////////////////////////////////////////////////
1802
1803proc union(list li)
1804{
1805  int i,j,k;
1806
1807  def P=basering;
1808
1809  execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
1810  list l=fetch(P,li);
1811  list @erg;
1812
1813  for(k=1;k<=size(l);k++)
1814  {
1815     for(j=1;j<=size(l[k])/2;j++)
1816     {
1817        if(deg(l[k][2*j][1])!=0)
1818        {
1819           i++;
1820           @erg[i]=l[k][2*j];
1821        }
1822     }
1823  }
1824
1825  list @wos;
1826  i=0;
1827  ideal i1,i2;
1828  while(i<size(@erg)-1)
1829  {
1830     i++;
1831     k=i+1;
1833      attrib(i1,"isSB",1);
1834      attrib(@erg[i],"isSB",1);
1835
1836     while(k<=size(@erg))
1837     {
1838        if(deg(@erg[i][1])==0)
1839        {
1840           break;
1841        }
1843        attrib(@erg[k],"isSB",1);
1844        attrib(i2,"isSB",1);
1845
1846        if(size(reduce(i1,i2,1))==0)
1847        {
1848           if(size(reduce(@erg[i],@erg[k],1))==0)
1849           {
1850              @erg[k]=ideal(1);
1851              i2=ideal(1);
1852           }
1853        }
1854        if(size(reduce(i2,i1,1))==0)
1855        {
1856           if(size(reduce(@erg[k],@erg[i],1))==0)
1857           {
1858              break;
1859           }
1860        }
1861        k++;
1862        if(k>size(@erg))
1863        {
1864           @wos[size(@wos)+1]=@erg[i];
1865        }
1866     }
1867  }
1868  if(deg(@erg[size(@erg)][1])!=0)
1869  {
1870     @wos[size(@wos)+1]=@erg[size(@erg)];
1871  }
1872  setring P;
1873  list @ser=fetch(ir,@wos);
1874  return(@ser);
1875}
1876///////////////////////////////////////////////////////////////////////////////
1878{
1879   list pr=minAssPrimes(i,1);
1880   int j;
1881   ideal k=pr[1];
1882   for(j=2;j<=size(pr);j++)
1883   {
1884      k=intersect(k,pr[j]);
1885   }
1886   return(k);
1887}
1888///////////////////////////////////////////////////////////////////////////////
1890{
1892}
1893
1894///////////////////////////////////////////////////////////////////////////////
1895proc decomp(ideal i,list #)
1896"USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
1897         decomp(i,1);        (for the minimal associated primes) )
1898RETURN:  list = list of primary ideals and their associated primes
1899         (at even positions in the list)
1900         (resp. a list of the minimal associated primes)
1901NOTE:    Algorithm of Gianni, Traeger, Zacharias
1902EXAMPLE: example decomp; shows an example
1903"
1904{
1905  def  @P = basering;
1906  list primary,indep,ltras;
1907  intvec @vh,isat;
1908  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi;
1909  ideal peek=i;
1910  ideal ser,tras;
1911
1912  if(size(#)>0)
1913  {
1914     if((#[1]==1)||(#[1]==2))
1915     {
1916        @wr=#[1];
1917        if(size(#)>1)
1918        {
1919          seri=1;
1920          peek=#[2];
1921          ser=#[3];
1922        }
1923      }
1924      else
1925      {
1926        seri=1;
1927        peek=#[1];
1928        ser=#[2];
1929      }
1930  }
1931
1932  homo=homog(i);
1933
1934  if(homo==1)
1935  {
1936    if(attrib(i,"isSB")!=1)
1937    {
1938      ltras=mstd(i);
1939      attrib(ltras[1],"isSB",1);
1940    }
1941    else
1942    {
1943      ltras=i,i;
1944    }
1945    tras=ltras[1];
1946    if(dim(tras)==0)
1947    {
1948        primary[1]=ltras[2];
1949        primary[2]=maxideal(1);
1950        if(@wr>0)
1951        {
1952           list l;
1953           l[1]=maxideal(1);
1954           l[2]=maxideal(1);
1955           return(l);
1956        }
1957        return(primary);
1958     }
1959     intvec @hilb=hilb(tras,1);
1960     intvec keephilb=@hilb;
1961  }
1962
1963  //----------------------------------------------------------------
1964  //i is the zero-ideal
1965  //----------------------------------------------------------------
1966
1967  if(size(i)==0)
1968  {
1969     primary=i,i;
1970     return(primary);
1971  }
1972
1973  //----------------------------------------------------------------
1974  //pass to the lexicographical ordering and compute a standardbasis
1975  //----------------------------------------------------------------
1976
1977  execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
1978  option(redSB);
1979
1980  ideal ser=fetch(@P,ser);
1981
1982  if(homo==1)
1983  {
1984     if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)"))
1985     {
1986       ideal @j=std(fetch(@P,i),@hilb);
1987     }
1988     else
1989     {
1990        ideal @j=fetch(@P,tras);
1991        attrib(@j,"isSB",1);
1992     }
1993  }
1994  else
1995  {
1996     ideal @j=std(fetch(@P,i));
1997  }
1998  option(noredSB);
1999  if(seri==1)
2000  {
2001    ideal peek=fetch(@P,peek);
2002    attrib(peek,"isSB",1);
2003  }
2004  else
2005  {
2006    ideal peek=@j;
2007  }
2008  if(size(ser)==0)
2009  {
2010    ideal fried;
2011    @n=size(@j);
2012    for(@k=1;@k<=@n;@k++)
2013    {
2015      {
2016        fried[size(fried)+1]=@j[@k];
2017        @j[@k]=0;
2018      }
2019    }
2020    if(size(fried)>0)
2021    {
2022       @j=simplify(@j,2);
2023       attrib(@j,"isSB",1);
2024       list pr=decomp(@j);
2025       for(@k=1;@k<=size(pr);@k++)
2026       {
2027         @j=pr[@k]+fried;
2028         pr[@k]=@j;
2029       }
2030       setring @P;
2031       return(fetch(gnir,pr));
2032    }
2033  }
2034
2035  //----------------------------------------------------------------
2036  //j is the ring
2037  //----------------------------------------------------------------
2038
2039  if (dim(@j)==-1)
2040  {
2041    setring @P;
2042    return(ideal(0));
2043  }
2044
2045  //----------------------------------------------------------------
2046  //  the case of one variable
2047  //----------------------------------------------------------------
2048
2049  if(nvars(basering)==1)
2050  {
2051
2052     list fac=factor(@j[1]);
2053     list gprimary;
2054     for(@k=1;@k<=size(fac[1]);@k++)
2055     {
2056        if(@wr==0)
2057        {
2058           gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
2059           gprimary[2*@k]=ideal(fac[1][@k]);
2060        }
2061        else
2062        {
2063           gprimary[2*@k-1]=ideal(fac[1][@k]);
2064           gprimary[2*@k]=ideal(fac[1][@k]);
2065        }
2066     }
2067     setring @P;
2068     primary=fetch(gnir,gprimary);
2069
2070     return(primary);
2071  }
2072
2073 //------------------------------------------------------------------
2074 //the zero-dimensional case
2075 //------------------------------------------------------------------
2076
2077  if (dim(@j)==0)
2078  {
2079    option(redSB);
2080    list gprimary= zero_decomp(@j,ser,@wr);
2081    option(noredSB);
2082    setring @P;
2083    primary=fetch(gnir,gprimary);
2084    if(size(ser)>0)
2085    {
2086      primary=cleanPrimary(primary);
2087    }
2088    return(primary);
2089  }
2090
2091  poly @gs,@gh,@p;
2092  string @va,quotring;
2093  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
2094  ideal @h;
2095  int jdim=dim(@j);
2096  list fett;
2097  int lauf,di,newtest;
2098  //------------------------------------------------------------------
2099  //search for a maximal independent set indep,i.e.
2100  //look for subring such that the intersection with the ideal is zero
2101  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
2102  //indep[1] is the new varstring and indep[2] the string for the block-ordering
2103  //------------------------------------------------------------------
2104
2105  if(@wr!=1)
2106  {
2107     allindep=independSet(@j);
2108     for(@m=1;@m<=size(allindep);@m++)
2109     {
2110        if(allindep[@m][3]==jdim)
2111        {
2112           di++;
2113           indep[di]=allindep[@m];
2114        }
2115        else
2116        {
2117           lauf++;
2118           restindep[lauf]=allindep[@m];
2119        }
2120     }
2121   }
2122   else
2123   {
2124      indep=maxIndependSet(@j);
2125   }
2126
2127  ideal jkeep=@j;
2128
2129  if(ordstr(@P)[1]=="w")
2130  {
2131     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");";
2132  }
2133  else
2134  {
2135     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);";
2136  }
2137
2138  if(homo==1)
2139  {
2140    if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
2141       ||(ordstr(@P)[3]=="w"))
2142    {
2143      ideal jwork=imap(@P,tras);
2144      attrib(jwork,"isSB",1);
2145    }
2146    else
2147    {
2148      ideal jwork=std(imap(gnir,@j),@hilb);
2149    }
2150
2151  }
2152  else
2153  {
2154    ideal jwork=std(imap(gnir,@j));
2155  }
2156  list hquprimary;
2157  poly @p,@q;
2158  ideal @h,fac,ser;
2159  di=dim(jwork);
2160  keepdi=di;
2161
2162  setring gnir;
2163  for(@m=1;@m<=size(indep);@m++)
2164  {
2165     isat=0;
2166     @n2=0;
2167     option(redSB);
2168     if((indep[@m][1]==varstr(basering))&&(@m==1))
2169     //this is the good case, nothing to do, just to have the same notations
2170     //change the ring
2171     {
2172        execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2173                              +ordstr(basering)+");";
2174        ideal @j=fetch(gnir,@j);
2175        attrib(@j,"isSB",1);
2176        ideal ser=fetch(gnir,ser);
2177
2178     }
2179     else
2180     {
2181        @va=string(maxideal(1));
2182        execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2183                              +indep[@m][2]+");";
2184        execute "map phi=gnir,"+@va+";";
2185        if(homo==1)
2186        {
2187           ideal @j=std(phi(@j),@hilb);
2188        }
2189        else
2190        {
2191          ideal @j=std(phi(@j));
2192        }
2193        ideal ser=phi(ser);
2194
2195     }
2196     option(noredSB);
2197     if((deg(@j[1])==0)||(dim(@j)<jdim))
2198     {
2199       setring gnir;
2200       break;
2201     }
2202     for (lauf=1;lauf<=size(@j);lauf++)
2203     {
2204         fett[lauf]=size(@j[lauf]);
2205     }
2206     //------------------------------------------------------------------------------------
2207     //we have now the following situation:
2208     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2209     //to this quotientring, j is their still a standardbasis, the
2210     //leading coefficients of the polynomials  there (polynomials in
2211     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2212     //we need their ggt, gh, because of the following:
2213     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2214     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2215     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2216
2217     //------------------------------------------------------------------------------------
2218
2219     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2220     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
2221     //-------------------------------------------------------------------------------------
2222
2223     quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2224
2225     //---------------------------------------------------------------------
2226     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2227     //---------------------------------------------------------------------
2228
2229     execute quotring;
2230
2231    // @j considered in the quotientring
2232     ideal @j=imap(gnir1,@j);
2233     ideal ser=imap(gnir1,ser);
2234
2235     kill gnir1;
2236
2237     //j is a standardbasis in the quotientring but usually not minimal
2238     //here it becomes minimal
2239
2240     @j=clearSB(@j,fett);
2241     attrib(@j,"isSB",1);
2242
2243     //we need later ggt(h[1],...)=gh for saturation
2244     ideal @h;
2245     if(deg(@j[1])>0)
2246     {
2247        for(@n=1;@n<=size(@j);@n++)
2248        {
2250        }
2251        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2252        option(redSB);
2253        list uprimary= zero_decomp(@j,ser,@wr);
2254        option(noredSB);
2255     }
2256     else
2257     {
2258       list uprimary;
2259       uprimary[1]=ideal(1);
2260       uprimary[2]=ideal(1);
2261     }
2262
2263     //we need the intersection of the ideals in the list quprimary with the
2264     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2265     //but fi polynomials, then the intersection of q with the polynomialring
2266     //is the saturation of the ideal generated by f1,...,fr with respect to
2267     //h which is the lcm of the leading coefficients of the fi considered in the
2268     //quotientring: this is coded in saturn
2269
2270     list saturn;
2271     ideal hpl;
2272
2273     for(@n=1;@n<=size(uprimary);@n++)
2274     {
2275        hpl=0;
2276        for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2277        {
2279        }
2280        saturn[@n]=hpl;
2281     }
2282     //--------------------------------------------------------------------
2283     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2284     //back to the polynomialring
2285     //---------------------------------------------------------------------
2286     setring gnir;
2287
2288     collectprimary=imap(quring,uprimary);
2289     lsau=imap(quring,saturn);
2290     @h=imap(quring,@h);
2291
2292     kill quring;
2293
2294
2295     @n2=size(quprimary);
2296     @n3=@n2;
2297
2298     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2299     {
2300        if(deg(collectprimary[2*@n1][1])>0)
2301        {
2302           @n2++;
2303           quprimary[@n2]=collectprimary[2*@n1-1];
2304           lnew[@n2]=lsau[2*@n1-1];
2305           @n2++;
2306           lnew[@n2]=lsau[2*@n1];
2307           quprimary[@n2]=collectprimary[2*@n1];
2308        }
2309     }
2310
2311     //here the intersection with the polynomialring
2312     //mentioned above is really computed
2313     for(@n=@n3/2+1;@n<=@n2/2;@n++)
2314     {
2315        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
2316        {
2317           quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2318           quprimary[2*@n]=quprimary[2*@n-1];
2319        }
2320        else
2321        {
2322           if(@wr==0)
2323           {
2324              quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2325           }
2326           quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
2327        }
2328     }
2329
2330     if(size(@h)>0)
2331     {
2332           //---------------------------------------------------------------
2333           //we change to @Phelp to have the ordering dp for saturation
2334           //---------------------------------------------------------------
2335        setring @Phelp;
2336        @h=imap(gnir,@h);
2337        if(@wr!=1)
2338        {
2339          @q=minSat(jwork,@h)[2];
2340        }
2341        else
2342        {
2343            fac=ideal(0);
2344            for(lauf=1;lauf<=ncols(@h);lauf++)
2345           {
2346              if(deg(@h[lauf])>0)
2347              {
2348                 fac=fac+factorize(@h[lauf],1);
2349              }
2350           }
2351           fac=simplify(fac,4);
2352           @q=1;
2353           for(lauf=1;lauf<=size(fac);lauf++)
2354           {
2355             @q=@q*fac[lauf];
2356           }
2357        }
2358        jwork=std(jwork,@q);
2359        keepdi=dim(jwork);
2360        if(keepdi<di)
2361        {
2362           setring gnir;
2363           @j=imap(@Phelp,jwork);
2364           break;
2365        }
2366        if(homo==1)
2367        {
2368              @hilb=hilb(jwork,1);
2369        }
2370
2371        setring gnir;
2372        @j=imap(@Phelp,jwork);
2373     }
2374  }
2375  if((size(quprimary)==0)&&(@wr>0))
2376  {
2377     @j=ideal(1);
2378     quprimary[1]=ideal(1);
2379     quprimary[2]=ideal(1);
2380  }
2381  if((size(quprimary)==0))
2382  {
2383    keepdi=di-1;
2384  }
2385  //---------------------------------------------------------------
2386  //notice that j=sat(j,gh) intersected with (j,gh^n)
2388  //---------------------------------------------------------------
2389  if((deg(@j[1])!=0)&&(@wr!=1))
2390  {
2391     if(size(quprimary)>0)
2392     {
2393        setring @Phelp;
2394        ser=imap(gnir,ser);
2395        hquprimary=imap(gnir,quprimary);
2396        if(@wr==0)
2397        {
2398           ideal htest=hquprimary[1];
2399           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
2400           {
2401              htest=intersect(htest,hquprimary[2*@n1-1]);
2402           }
2403        }
2404        else
2405        {
2406           ideal htest=hquprimary[2];
2407
2408           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
2409           {
2410              htest=intersect(htest,hquprimary[2*@n1]);
2411           }
2412        }
2413
2414        if(size(ser)>0)
2415        {
2416           ser=intersect(htest,ser);
2417        }
2418        else
2419        {
2420          ser=htest;
2421        }
2422        setring gnir;
2423        ser=imap(@Phelp,ser);
2424     }
2425     if(size(reduce(ser,peek,1))!=0)
2426     {
2427        for(@m=1;@m<=size(restindep);@m++)
2428        {
2429         // if(restindep[@m][3]>=keepdi)
2430         // {
2431           isat=0;
2432           @n2=0;
2433           option(redSB);
2434
2435           if(restindep[@m][1]==varstr(basering))
2436           //this is the good case, nothing to do, just to have the same notations
2437           //change the ring
2438           {
2439              execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2440                              +ordstr(basering)+");";
2441              ideal @j=fetch(gnir,jkeep);
2442              attrib(@j,"isSB",1);
2443           }
2444           else
2445           {
2446              @va=string(maxideal(1));
2447              execute "ring gnir1 = ("+charstr(basering)+"),("+restindep[@m][1]+"),("
2448                              +restindep[@m][2]+");";
2449              execute "map phi=gnir,"+@va+";";
2450              if(homo==1)
2451              {
2452                 ideal @j=std(phi(jkeep),keephilb);
2453              }
2454              else
2455              {
2456                ideal @j=std(phi(jkeep));
2457              }
2458              ideal ser=phi(ser);
2459           }
2460           option(noredSB);
2461
2462           for (lauf=1;lauf<=size(@j);lauf++)
2463           {
2464              fett[lauf]=size(@j[lauf]);
2465           }
2466           //------------------------------------------------------------------------------------
2467           //we have now the following situation:
2468           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2469           //to this quotientring, j is their still a standardbasis, the
2470           //leading coefficients of the polynomials  there (polynomials in
2471           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2472           //we need their ggt, gh, because of the following:
2473           //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2474           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2475           //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2476
2477           //------------------------------------------------------------------------------------
2478
2479           //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2480           //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
2481           //-------------------------------------------------------------------------------------
2482
2483           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
2484
2485           //---------------------------------------------------------------------
2486           //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2487           //---------------------------------------------------------------------
2488
2489           execute quotring;
2490
2491           // @j considered in the quotientring
2492           ideal @j=imap(gnir1,@j);
2493           ideal ser=imap(gnir1,ser);
2494
2495           kill gnir1;
2496
2497           //j is a standardbasis in the quotientring but usually not minimal
2498           //here it becomes minimal
2499           @j=clearSB(@j,fett);
2500           attrib(@j,"isSB",1);
2501
2502           //we need later ggt(h[1],...)=gh for saturation
2503           ideal @h;
2504
2505           for(@n=1;@n<=size(@j);@n++)
2506           {
2508           }
2509           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2510
2511           option(redSB);
2512           list uprimary= zero_decomp(@j,ser,@wr);
2513           option(noredSB);
2514
2515
2516           //we need the intersection of the ideals in the list quprimary with the
2517           //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2518           //but fi polynomials, then the intersection of q with the polynomialring
2519           //is the saturation of the ideal generated by f1,...,fr with respect to
2520           //h which is the lcm of the leading coefficients of the fi considered in the
2521           //quotientring: this is coded in saturn
2522
2523           list saturn;
2524           ideal hpl;
2525
2526           for(@n=1;@n<=size(uprimary);@n++)
2527           {
2528              hpl=0;
2529              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2530              {
2532              }
2533               saturn[@n]=hpl;
2534           }
2535           //--------------------------------------------------------------------
2536           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2537           //back to the polynomialring
2538           //---------------------------------------------------------------------
2539           setring gnir;
2540
2541           collectprimary=imap(quring,uprimary);
2542           lsau=imap(quring,saturn);
2543           @h=imap(quring,@h);
2544
2545           kill quring;
2546
2547
2548           @n2=size(quprimary);
2549           @n3=@n2;
2550
2551           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2552           {
2553              if(deg(collectprimary[2*@n1][1])>0)
2554              {
2555                 @n2++;
2556                 quprimary[@n2]=collectprimary[2*@n1-1];
2557                 lnew[@n2]=lsau[2*@n1-1];
2558                 @n2++;
2559                 lnew[@n2]=lsau[2*@n1];
2560                 quprimary[@n2]=collectprimary[2*@n1];
2561              }
2562           }
2563
2564
2565           //here the intersection with the polynomialring
2566           //mentioned above is really computed
2567
2568           for(@n=@n3/2+1;@n<=@n2/2;@n++)
2569           {
2570              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
2571              {
2572                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2573                 quprimary[2*@n]=quprimary[2*@n-1];
2574              }
2575              else
2576              {
2577                 if(@wr==0)
2578                 {
2579                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2580                 }
2581                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
2582              }
2583           }
2584           if(@n2>=@n3+2)
2585           {
2586              setring @Phelp;
2587              ser=imap(gnir,ser);
2588              hquprimary=imap(gnir,quprimary);
2589              for(@n=@n3/2+1;@n<=@n2/2;@n++)
2590              {
2591                if(@wr==0)
2592                {
2593                   ser=intersect(ser,hquprimary[2*@n-1]);
2594                }
2595                else
2596                {
2597                   ser=intersect(ser,hquprimary[2*@n]);
2598                }
2599              }
2600              setring gnir;
2601              ser=imap(@Phelp,ser);
2602           }
2603
2604         // }
2605        }
2606        if(size(reduce(ser,peek,1))!=0)
2607        {
2608           if(@wr>0)
2609           {
2610              htprimary=decomp(@j,@wr,peek,ser);
2611           }
2612           else
2613           {
2614              htprimary=decomp(@j,peek,ser);
2615           }
2616           // here we collect now both results primary(sat(j,gh))
2617           // and primary(j,gh^n)
2618           @n=size(quprimary);
2619           for (@k=1;@k<=size(htprimary);@k++)
2620           {
2621              quprimary[@n+@k]=htprimary[@k];
2622           }
2623        }
2624     }
2625
2626   }
2627  //------------------------------------------------------------
2628  //back to the ring we started with
2629  //the final result: primary
2630  //------------------------------------------------------------
2631
2632  setring @P;
2633  primary=imap(gnir,quprimary);
2634  return(primary);
2635}
2636
2637
2638example
2639{ "EXAMPLE:"; echo = 2;
2640   ring  r = 32003,(x,y,z),lp;
2641   poly  p = z2+1;
2642   poly  q = z4+2;
2643   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
2644   LIB "primdec.lib";
2645   list pr= decomp(i);
2646   pr;
2647   testPrimary( pr, i);
2648}
2649
2650///////////////////////////////////////////////////////////////////////////////
2651proc primdecGTZ(ideal i)
2652{
2653   return(decomp(i));
2654}
2655///////////////////////////////////////////////////////////////////////////////
2656proc primdecSY(ideal i,int j)
2657{
2658   return(prim_dec(i,j));
2659}
2660///////////////////////////////////////////////////////////////////////////////
2661
2663{
2664   def @P=basering;
2665   int j,il;
2666   if(size(#)>0)
2667   {
2668     il=#[1];
2669   }
2670   ideal re=1;
2671   option(redSB);
2672   list qr=simplifyIdeal(i);
2673   map phi=@P,qr[2];
2674   i=qr[1];
2675
2676   list pr=facstd(i);
2677
2678   if(size(pr)==1)
2679   {
2680      attrib(pr[1],"isSB",1);
2681      if((dim(pr[1])==0)&&(homog(pr[1])==1))
2682      {
2683         ideal @res=maxideal(1);
2684         return(phi(@res));
2685      }
2686      if(dim(pr[1])>1)
2687      {
2688         execute "ring gnir = ("+charstr(basering)+"),
2689                              ("+varstr(basering)+"),(C,lp);";
2690         ideal i=fetch(@P,i);
2691         list @pr=facstd(i);
2692         setring @P;
2693         pr=fetch(gnir,@pr);
2694      }
2695   }
2696   option(noredSB);
2697   int s=size(pr);
2698
2699   if(s==1)
2700   {
2702     return(phi(i));
2703   }
2704   intvec pos;
2705   pos[s]=0;
2706   if(il==1)
2707   {
2708     int ndim,k;
2709     attrib(pr[1],"isSB",1);
2710     int odim=dim(pr[1]);
2711     int count=1;
2712
2713     for(j=2;j<=s;j++)
2714     {
2715        attrib(pr[j],"isSB",1);
2716        ndim=dim(pr[j]);
2717        if(ndim>odim)
2718        {
2719           for(k=count;k<=j-1;k++)
2720           {
2721              pos[k]=1;
2722           }
2723           count=j;
2724           odim=ndim;
2725        }
2726        if(ndim<odim)
2727        {
2728           pos[j]=1;
2729        }
2730     }
2731   }
2732   for(j=1;j<=s;j++)
2733   {
2734      if(pos[j]==0)
2735      {
2737      }
2738   }
2739   return(phi(re));
2740}
2741
2742proc intersect1(ideal i,ideal j)
2743{
2744   def R=basering;
2745   execute "ring gnir = ("+charstr(basering)+"),
2746                        ("+varstr(basering)+",@t),(C,dp);";
2747   ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
2748   ideal j=eliminate(i,var(nvars(basering)));
2749   setring R;
2750   map phi=gnir,maxideal(1);
2751   return(phi(j));
2752}
2753
2754
2755///////////////////////////////////////////////////////////////////////////////
2756proc radicalKL (list m,ideal ser,list #)
2757"USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
2758         decomp(i,1);        (for the minimal associated primes) )
2759RETURN:  list = list of primary ideals and their associated primes
2760         (at even positions in the list)
2761         (resp. a list of the minimal associated primes)
2762NOTE:    Algorithm of Gianni, Traeger, Zacharias
2763EXAMPLE: example decomp; shows an example
2764"
2765{
2766   ideal i=m[2];
2767  //----------------------------------------------------------------
2768  //i is the zero-ideal
2769  //----------------------------------------------------------------
2770
2771   if(size(i)==0)
2772   {
2773      return(ideal(0));
2774   }
2775
2776   def  @P = basering;
2777   list indep,allindep,restindep,fett,@mu;
2778   intvec @vh,isat;
2779   int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di;
2782   ideal te=ser;
2783
2784   poly @p,@q;
2785   string @va,quotring;
2786   int  homo=homog(i);
2787
2788   if(size(#)>0)
2789   {
2790      @wr=#[1];
2791   }
2792   @j=m[1];
2793   @j1=m[2];
2794   int jdim=dim(@j);
2795   if(size(reduce(ser,@j,1))==0)
2796   {
2797      return(ser);
2798   }
2799   if(homo==1)
2800   {
2801      if(jdim==0)
2802      {
2803         option(noredSB);
2804         return(maxideal(1));
2805      }
2806      intvec @hilb=hilb(@j,1);
2807   }
2808
2809
2810  //----------------------------------------------------------------
2811  //j is the ring
2812  //----------------------------------------------------------------
2813
2814   if (jdim==-1)
2815   {
2816      option(noredSB);
2817      return(ideal(0));
2818   }
2819
2820  //----------------------------------------------------------------
2821  //  the case of one variable
2822  //----------------------------------------------------------------
2823
2824   if(nvars(basering)==1)
2825   {
2826      fac=factorize(@j[1],1);
2827      @p=1;
2828      for(@k=1;@k<=size(fac);@k++)
2829      {
2830         @p=@p*fac[@k];
2831      }
2832      option(noredSB);
2833
2834      return(ideal(@p));
2835   }
2836   //------------------------------------------------------------------
2837   //the case of a complete intersection
2838   //------------------------------------------------------------------
2839    if(jdim+size(@j1)==nvars(basering))
2840    {
2841      // ideal jac=minor(jacob(@j1),size(@j1));
2842      // return(quotient(@j1,jac));
2843    }
2844
2845   //------------------------------------------------------------------
2846   //the zero-dimensional case
2847   //------------------------------------------------------------------
2848
2849   if (jdim==0)
2850   {
2851      @j1=finduni(@j);
2852      for(@k=1;@k<=size(@j1);@k++)
2853      {
2854         fac=factorize(cleardenom(@j1[@k]),1);
2855         @p=fac[1];
2856         for(@n=2;@n<=size(fac);@n++)
2857         {
2858            @p=@p*fac[@n];
2859         }
2860         @j=@j,@p;
2861      }
2862      @j=std(@j);
2863      option(noredSB);
2864      return(@j);
2865   }
2866
2867   //------------------------------------------------------------------
2868   //search for a maximal independent set indep,i.e.
2869   //look for subring such that the intersection with the ideal is zero
2870   //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
2871   //indep[1] is the new varstring and indep[2] the string for the block-ordering
2872   //------------------------------------------------------------------
2873
2874   indep=maxIndependSet(@j);
2875
2876   di=dim(@j);
2877
2878   for(@m=1;@m<=size(indep);@m++)
2879   {
2880      if((indep[@m][1]==varstr(basering))&&(@m==1))
2881      //this is the good case, nothing to do, just to have the same notations
2882      //change the ring
2883      {
2884         execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2885                              +ordstr(basering)+");";
2886         ideal @j=fetch(@P,@j);
2887         attrib(@j,"isSB",1);
2888      }
2889      else
2890      {
2891         @va=string(maxideal(1));
2892         execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2893                              +indep[@m][2]+");";
2894         execute "map phi=@P,"+@va+";";
2895         if(homo==1)
2896         {
2897            ideal @j=std(phi(@j),@hilb);
2898         }
2899         else
2900         {
2901           ideal @j=std(phi(@j));
2902         }
2903      }
2904      if((deg(@j[1])==0)||(dim(@j)<jdim))
2905      {
2906         setring @P;
2907         break;
2908      }
2909      for (lauf=1;lauf<=size(@j);lauf++)
2910      {
2911         fett[lauf]=size(@j[lauf]);
2912      }
2913     //------------------------------------------------------------------------------------
2914     //we have now the following situation:
2915     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2916     //to this quotientring, j is their still a standardbasis, the
2917     //leading coefficients of the polynomials  there (polynomials in
2918     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2919     //we need their ggt, gh, because of the following:
2920     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2921     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2922     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2923
2924     //------------------------------------------------------------------------------------
2925
2926     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2927     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
2928     //-------------------------------------------------------------------------------------
2929
2930      quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2931
2932     //---------------------------------------------------------------------
2933     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2934     //---------------------------------------------------------------------
2935
2936      execute quotring;
2937
2938    // @j considered in the quotientring
2939      ideal @j=imap(gnir1,@j);
2940
2941      kill gnir1;
2942
2943     //j is a standardbasis in the quotientring but usually not minimal
2944     //here it becomes minimal
2945
2946      @j=clearSB(@j,fett);
2947      attrib(@j,"isSB",1);
2948
2949     //we need later ggt(h[1],...)=gh for saturation
2950      ideal @h;
2951      if(deg(@j[1])>0)
2952      {
2953         for(@n=1;@n<=size(@j);@n++)
2954         {
2956         }
2957        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2958         option(redSB);
2959         @j=interred(@j);
2960         attrib(@j,"isSB",1);
2961         list  @mo=@j,@j;
2963      }
2964      else
2965      {
2967      }
2968
2969     //we need the intersection of the ideals in the list quprimary with the
2970     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2971     //but fi polynomials, then the intersection of q with the polynomialring
2972     //is the saturation of the ideal generated by f1,...,fr with respect to
2973     //h which is the lcm of the leading coefficients of the fi considered in the
2974     //quotientring: this is coded in saturn
2975
2976      ideal hpl;
2977
2979      {
2981      }
2982
2983     //--------------------------------------------------------------------
2984     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2985     //back to the polynomialring
2986     //---------------------------------------------------------------------
2987      setring @P;
2988
2990      lsau=simplify(imap(quring,hpl),2);
2991      @h=imap(quring,@h);
2992
2993      kill quring;
2994
2995
2996     //here the intersection with the polynomialring
2997     //mentioned above is really computed
2998
3000
3001      if(deg(@h[1])>=0)
3002      {
3003         fac=ideal(0);
3004         for(lauf=1;lauf<=ncols(@h);lauf++)
3005         {
3006            if(deg(@h[lauf])>0)
3007            {
3008                fac=fac+factorize(@h[lauf],1);
3009            }
3010         }
3011         fac=simplify(fac,4);
3012         @q=1;
3013         for(lauf=1;lauf<=size(fac);lauf++)
3014         {
3015            @q=@q*fac[lauf];
3016         }
3017
3018
3020         @j=@mu[1];
3021         attrib(@j,"isSB",1);
3022
3023      }
3025      {
3027      }
3028      else
3029      {
3031         {
3033         }
3034      }
3035
3037
3038      if((dim(@j)<di)||(size(te)==0))
3039      {
3040         break;
3041      }
3042      if(homo==1)
3043      {
3044         @hilb=hilb(@j,1);
3045      }
3046   }
3047
3048   if(((@wr==1)&&(dim(@j)<di))||(deg(@j[1])==0)||(size(te)==0))
3049   {
3051   }
3054
3055
3056   option(noredSB);
3058}
3059
3060
3061example
3062{ "EXAMPLE:"; echo = 2;
3063}
3064
3065
3067{
3069   int l,il;
3070   if(size(#)>0)
3071   {
3072      il=#[1];
3073   }
3074
3075   option(redSB);
3076   list m=mstd(i);
3077        I=m[2];
3078   option(noredSB);
3079   if(size(reduce(re,m[1],1))==0)
3080   {
3081      return(re);
3082   }
3083   int cod=nvars(basering)-dim(m[1]);
3084   if((nvars(basering)<=5)&&(size(m[2])<=5))
3085   {
3086      if(cod==size(m[2]))
3087      {
3088        J=minor(jacob(I),cod);
3089        return(quotient(I,J));
3090      }
3091
3092      for(l=1;l<=cod;l++)
3093      {
3094         I0[l]=I[l];
3095      }
3096      if(dim(std(I0))+cod==nvars(basering))
3097      {
3098         J=minor(jacob(I0),cod);
3102
3104         {
3105            return(I);
3106         }
3107         if(il==1)
3108         {
3109
3111         }
3112
3114
3115         if(deg(I2[1])<=0)
3116         {
3118         }
3120      }
3121   }
3123}
3124
3125proc Ann(module M)
3126{
3127  M=prune(M);  //to obtain a small embedding
3128  ideal ann=quotient1(M,freemodule(nrows(M)));
3129  return(ann);
3130}
3131
3132//computes the equidimensional part of the ideal i of codimension e
3133proc int_ass_primary_e(ideal i, int e)
3134{
3135  if(homog(i)!=1)
3136  {
3137     i=std(i);
3138  }
3139  list re=sres(i,0);                   //the resolution
3140  re=minres(re);                       //minimized resolution
3141  ideal ann=AnnExt_R(e,re);
3142  if(nvars(basering)-dim(std(ann))!=e)
3143  {
3144    return(ideal(1));
3145  }
3146  return(ann);
3147}
3148
3149//computes all equidimensional parts of the ideal i
3150proc prepareAss(ideal i)
3151{
3152  ideal j=std(i);
3153  int cod=nvars(basering)-dim(j);
3154  int e;
3155  list er;
3156  ideal ann;
3157  if(homog(i)==1)
3158  {
3159     list re=sres(i,0);                   //the resolution
3160     re=minres(re);                       //minimized resolution
3161  }
3162  else
3163  {
3164     list re=mres(i,0);             //fehler in sres
3165  }
3166  for(e=cod;e<=nvars(basering);e++)
3167  {
3168     ann=AnnExt_R(e,re);
3169
3170     if(nvars(basering)-dim(std(ann))==e)
3171     {
3173     }
3174  }
3175  return(er);
3176}
3177
3178//computes the annihilator of Ext^n(R/i,R) with given resolution re
3179//n is not necessarily the number of variables
3180proc AnnExt_R(int n,list re)
3181{
3182  if(n<nvars(basering))
3183  {
3184     matrix f=transpose(re[n+1]);      //Hom(_,R)
3185     module k=res(f,2)[2];             //the kernel
3186     matrix g=transpose(re[n]);        //the image of Hom(_,R)
3187
3188     ideal ann=quotient1(g,k);           //the anihilator
3189  }
3190  else
3191  {
3192     ideal ann=Ann(transpose(re[n]));
3193  }
3194  return(ann);
3195}
3196
3197proc quotient1(module a,module b)
3198{
3199   int i;
3200   a=std(a);     //<== neue Zeile
3201   ideal re=quotient(a,module(b[1]));
3202   for(i=2;i<=size(b);i++)
3203   {
3204      re=intersect1(re,quotient(a,module(b[i])));
3205   }
3206   return(re);
3207}
3208
3209
3210
3211proc analyze(list pr)
3212{
3213   int ii,jj;
3214   for(ii=1;ii<=size(pr)/2;ii++)
3215   {
3216      dim(std(pr[2*ii]));
3217      idealsEqual(pr[2*ii-1],pr[2*ii]);
3218      "===========================";
3219   }
3220
3221   for(ii=size(pr)/2;ii>1;ii--)
3222   {
3223      for(jj=1;jj<ii;jj++)
3224      {
3225         if(size(reduce(pr[2*jj],std(pr[2*ii],1)))==0)
3226         {
3227            "eingebette Komponente";
3228            jj;
3229            ii;
3230         }
3231      }
3232   }
3233}
3234
3235
3236proc simplifyIdeal(ideal i)
3237{
3238  def r=basering;
3239
3240  int j,k;
3241  map phi;
3242  poly p;
3243
3244  ideal iwork=i;
3245  ideal imap1=maxideal(1);
3246  ideal imap2=maxideal(1);
3247
3248
3249  for(j=1;j<=nvars(basering);j++)
3250  {
3251    for(k=1;k<=size(i);k++)
3252    {
3253      if(deg(iwork[k]/var(j))==0)
3254      {
3256        imap1[j]=p+2*var(j);
3257        phi=r,imap1;
3258        iwork=phi(iwork);
3259        iwork=subst(iwork,var(j),0);
3260        iwork[k]=var(j);
3261        imap1=maxideal(1);
3262        imap2[j]=-p;
3263        break;
3264      }
3265    }
3266  }
3267  return(iwork,imap2);
3268}
3269
3270
3271///////////////////////////////////////////////////////
3272// ini_mod
3273// input: a polynomial p
3274// output: the initial term of p as needed
3275// in the context of characteristic sets
3276//////////////////////////////////////////////////////
3277
3278proc ini_mod(poly p)
3279{
3280  if (p==0)
3281  {
3282    return(0);
3283  }
3284  int n; matrix m;
3285  for( n=nvars(basering); n>0; n=n-1)
3286  {
3287    m=coef(p,var(n));
3288    if(m[1,1]!=1)
3289    {
3290      p=m[2,1];
3291      break;
3292    }
3293  }
3294  if(deg(p)==0)
3295  {
3296    p=0;
3297  }
3298  return(p);
3299}
3300///////////////////////////////////////////////////////
3301// min_ass_prim_charsets
3302// input: generators of an ideal PS and an integer cho
3303// If cho=0, the given ordering of the variables is used.
3304// Otherwise, the system tries to find an "optimal ordering",
3305// which in some cases may considerably speed up the algorithm
3306// output: the minimal associated primes of PS
3307// algorithm: via characteriostic sets
3308//////////////////////////////////////////////////////
3309
3310
3311proc min_ass_prim_charsets (ideal PS, int cho)
3312{
3313  if((cho<0) and (cho>1))
3314    {
3315      "ERROR: <int> must be 0 or 1"
3316      return();
3317    }
3318  if(system("version")>933)
3319  {
3320    option(notWarnSB);
3321  }
3322  if(cho==0)
3323  {
3324    return(min_ass_prim_charsets0(PS));
3325  }
3326  else
3327  {
3328     return(min_ass_prim_charsets1(PS));
3329  }
3330}
3331///////////////////////////////////////////////////////
3332// min_ass_prim_charsets0
3333// input: generators of an ideal PS
3334// output: the minimal associated primes of PS
3335// algorithm: via characteristic sets
3336// the given ordering of the variables is used
3337//////////////////////////////////////////////////////
3338
3339
3340proc min_ass_prim_charsets0 (ideal PS)
3341{
3342
3343  matrix m=char_series(PS);  // We compute an irreducible
3344                             // characteristic series
3345  int i,j,k;
3346  list PSI;
3347  list PHI;  // the ideals given by the characteristic series
3348  for(i=nrows(m);i>=1; i--)
3349  {
3350     PHI[i]=ideal(m[i,1..ncols(m)]);
3351  }
3352  // We compute the radical of each ideal in PHI
3353  ideal I,JS,II;
3354  int sizeJS, sizeII;
3355  for(i=size(PHI);i>=1; i--)
3356  {
3357     I=0;
3358     for(j=size(PHI[i]);j>0;j=j-1)
3359     {
3360       I=I+ini_mod(PHI[i][j]);
3361     }
3362     JS=std(PHI[i]);
3363     sizeJS=size(JS);
3364     for(j=size(I);j>0;j=j-1)
3365     {
3366       II=0;
3367       sizeII=0;
3368       k=0;
3369       while(k<=sizeII)                  // successive saturation
3370       {
3371         option(returnSB);
3372         II=quotient(JS,I[j]);
3373         option(noreturnSB);
3374//std
3375//         II=std(II);
3376         sizeII=size(II);
3377         if(sizeII==sizeJS)
3378         {
3379            for(k=1;k<=sizeII;k++)
3380            {
3382            }
3383         }
3384         JS=II;
3385         sizeJS=sizeII;
3386       }
3387    }
3388    PSI=insert(PSI,JS);
3389  }
3390  int sizePSI=size(PSI);
3391  // We eliminate redundant ideals
3392  for(i=1;i<sizePSI;i++)
3393  {
3394    for(j=i+1;j<=sizePSI;j++)
3395    {
3396      if(size(PSI[i])!=0)
3397      {
3398        if(size(PSI[j])!=0)
3399        {
3400          if(size(NF(PSI[i],PSI[j],1))==0)
3401          {
3402            PSI[j]=ideal(0);
3403          }
3404          else
3405          {
3406            if(size(NF(PSI[j],PSI[i],1))==0)
3407            {
3408              PSI[i]=ideal(0);
3409            }
3410          }
3411        }
3412      }
3413    }
3414  }
3415  for(i=sizePSI;i>=1;i--)
3416  {
3417    if(size(PSI[i])==0)
3418    {
3419      PSI=delete(PSI,i);
3420    }
3421  }
3422  return (PSI);
3423}
3424
3425///////////////////////////////////////////////////////
3426// min_ass_prim_charsets1
3427// input: generators of an ideal PS
3428// output: the minimal associated primes of PS
3429// algorithm: via characteristic sets
3430// input: generators of an ideal PS and an integer i
3431// The system tries to find an "optimal ordering" of
3432// the variables
3433//////////////////////////////////////////////////////
3434
3435
3436proc min_ass_prim_charsets1 (ideal PS)
3437{
3438  def oldring=basering;
3439  string n=system("neworder",PS);
3440  execute "ring r="+charstr(oldring)+",("+n+"),dp;";
3441  ideal PS=imap(oldring,PS);
3442  matrix m=char_series(PS);  // We compute an irreducible
3443                             // characteristic series
3444  int i,j,k;
3445  ideal I;
3446  list PSI;
3447  list PHI;    // the ideals given by the characteristic series
3448  list ITPHI;  // their initial terms
3449  for(i=nrows(m);i>=1; i--)
3450  {
3451     PHI[i]=ideal(m[i,1..ncols(m)]);
3452     I=0;
3453     for(j=size(PHI[i]);j>0;j=j-1)
3454     {
3455       I=I,ini_mod(PHI[i][j]);
3456     }
3457      I=I[2..ncols(I)];
3458      ITPHI[i]=I;
3459  }
3460  setring oldring;
3461  matrix m=imap(r,m);
3462  list PHI=imap(r,PHI);
3463  list ITPHI=imap(r,ITPHI);
3464  // We compute the radical of each ideal in PHI
3465  ideal I,JS,II;
3466  int sizeJS, sizeII;
3467  for(i=size(PHI);i>=1; i--)
3468  {
3469     I=0;
3470     for(j=size(PHI[i]);j>0;j=j-1)
3471     {
3472       I=I+ITPHI[i][j];
3473     }
3474     JS=std(PHI[i]);
3475     sizeJS=size(JS);
3476     for(j=size(I);j>0;j=j-1)
3477     {
3478       II=0;
3479       sizeII=0;
3480       k=0;
3481       while(k<=sizeII)                  // successive iteration
3482       {
3483         option(returnSB);
3484         II=quotient(JS,I[j]);
3485         option(noreturnSB);
3486//std
3487//         II=std(II);
3488         sizeII=size(II);
3489         if(sizeII==sizeJS)
3490         {
3491            for(k=1;k<=sizeII;k++)
3492            {
3494            }
3495         }
3496         JS=II;
3497         sizeJS=sizeII;
3498       }
3499    }
3500    PSI=insert(PSI,JS);
3501  }
3502  int sizePSI=size(PSI);
3503  // We eliminate redundant ideals
3504  for(i=1;i<sizePSI;i++)
3505  {
3506    for(j=i+1;j<=sizePSI;j++)
3507    {
3508      if(size(PSI[i])!=0)
3509      {
3510        if(size(PSI[j])!=0)
3511        {
3512          if(size(NF(PSI[i],PSI[j],1))==0)
3513          {
3514            PSI[j]=ideal(0);
3515          }
3516          else
3517          {
3518            if(size(NF(PSI[j],PSI[i],1))==0)
3519            {
3520              PSI[i]=ideal(0);
3521            }
3522          }
3523        }
3524      }
3525    }
3526  }
3527  for(i=sizePSI;i>=1;i--)
3528  {
3529    if(size(PSI[i])==0)
3530    {
3531      PSI=delete(PSI,i);
3532    }
3533  }
3534  return (PSI);
3535}
3536
3537
3538/////////////////////////////////////////////////////
3539// proc prim_dec
3540// input:  generators of an ideal I and an integer choose
3541// If choose=0, min_ass_prim_charsets with the given
3542// ordering of the variables is used.
3543// If choose=1, min_ass_prim_charsets with the "optimized"
3544// ordering of the variables is used.
3545// If choose=2, minAssPrimes from primdec.lib is used
3546// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3547// output: a primary decomposition of I, i.e., a list
3548// of pairs consisting of a standard basis of a primary component
3549// of I and a standard basis of the corresponding associated prime.
3550// To compute the minimal associated primes of a given ideal
3551// min_ass_prim_l is called, i.e., the minimal associated primes
3552// are computed via characteristic sets.
3553// In the homogeneous case, the performance of the procedure
3554// will be improved if I is already given by a minimal set of
3555// generators. Apply minbase if necessary.
3556//////////////////////////////////////////////////////////
3557
3558
3559proc prim_dec(ideal I, int choose)
3560{
3561   if((choose<0) or (choose>3))
3562   {
3563     "ERROR: <int> must be 0 or 1 or 2 or 3";
3564      return();
3565   }
3566   if(system("version")>933)
3567   {
3568      option(notWarnSB);
3569    }
3570  ideal H=1; // The intersection of the primary components
3571  list U;    // the leaves of the decomposition tree, i.e.,
3572             // pairs consisting of a primary component of I
3573             // and the corresponding associated prime
3574  list W;    // the non-leaf vertices in the decomposition tree.
3575             // every entry has 6 components:
3576                // 1- the vertex itself , i.e., a standard bais of the
3577                //    given ideal I (type 1), or a standard basis of a
3578                //    pseudo-primary component arising from
3579                //    pseudo-primary decomposition (type 2), or a
3580                //    standard basis of a remaining component arising from
3581                //    pseudo-primary decomposition or extraction (type 3)
3582                // 2- the type of the vertex as indicated above
3583                // 3- the weighted_tree_depth of the vertex
3584                // 4- the tester of the vertex
3585                // 5- a standard basis of the associated prime
3586                //    of a vertex of type 2, or 0 otherwise
3587                // 6- a list of pairs consisting of a standard
3588                //    basis of a minimal associated prime ideal
3589                //    of the father of the vertex and the
3590                //    irreducible factors of the "minimal
3591                //    divisor" of the seperator or extractor
3592                //    corresponding to the prime ideal
3593                //    as computed by the procedure minsat,
3594                //    if the vertex is of type 3, or
3595                //    the empty list otherwise
3596  ideal SI=std(I);
3597  int ncolsSI=ncols(SI);
3598  int ncolsH=1;
3599  W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
3600  int weighted_tree_depth;
3601  int i,j;
3602  int check;
3603  list V;  // current vertex
3604  list VV; // new vertex
3605  list QQ;
3606  list WI;
3607  ideal Qi,SQ,SRest,fac;
3608  poly tester;
3609
3610  while(1)
3611  {
3612    i=1;
3613    while(1)
3614    {
3615      while(i<=size(W)) // find vertex V of smallest weighted tree-depth
3616      {
3617        if (W[i][3]<=weighted_tree_depth) break;
3618        i++;
3619      }
3620      if (i<=size(W)) break;
3621      i=1;
3622      weighted_tree_depth++;
3623    }
3624    V=W[i];
3625    W=delete(W,i); // delete V from W
3626
3627    // now proceed by type of vertex V
3628
3629    if (V[2]==2)  // extraction needed
3630    {
3631      SQ,SRest,fac=extraction(V[1],V[5]);
3632                        // standard basis of primary component,
3633                        // standard basis of remaining component,
3634                        // irreducible factors of
3635                        // the "minimal divisor" of the extractor
3636                        // as computed by the procedure minsat,
3637      check=0;
3638      for(j=1;j<=ncolsH;j++)
3639      {
3640        if (NF(H[j],SQ,1)!=0) // Q is not redundant
3641        {
3642          check=1;
3643          break;
3644        }
3645      }
3646      if(check==1)             // Q is not redundant
3647      {
3648        QQ=list();
3649        QQ[1]=list(SQ,V[5]);  // primary component, associated prime,
3650                              // i.e., standard bases thereof
3651        U=U+QQ;
3652        H=intersect(H,SQ);
3653        H=std(H);
3654        ncolsH=ncols(H);
3655        check=0;
3656        if(ncolsH==ncolsSI)
3657        {
3658          for(j=1;j<=ncolsSI;j++)
3659          {
3661            {
3662              check=1;
3663              break;
3664            }
3665          }
3666        }
3667        else
3668        {
3669          check=1;
3670        }
3671        if(check==0) // H==I => U is a primary decomposition
3672        {
3673          return(U);
3674        }
3675      }
3676      if (SRest[1]!=1)        // the remaining component is not
3677                              // the whole ring
3678      {
3679        if (rad_con(V[4],SRest)==0) // the new vertex is not the
3680                                    // root of a redundant subtree
3681        {
3682          VV[1]=SRest;     // remaining component
3683          VV[2]=3;         // pseudoprimdec_special
3684          VV[3]=V[3]+1;    // weighted depth
3685          VV[4]=V[4];      // the tester did not change
3686          VV[5]=ideal(0);
3687          VV[6]=list(list(V[5],fac));
3688          W=insert(W,VV,size(W));
3689        }
3690      }
3691    }
3692    else
3693    {
3694      if (V[2]==3) // pseudo_prim_dec_special is needed
3695      {
3696        QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
3698                         // standard basis of pseudo-primary component,
3699                         // standard basis of corresponding prime,
3700                         // seperator, irreducible factors of
3701                         // the "minimal divisor" of the seperator
3702                         // as computed by the procedure minsat,
3703                         // SRest=standard basis of remaining component
3704      }
3705      else     // V is the root, pseudo_prim_dec is needed
3706      {
3707        QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
3709                         // standard basis of pseudo-primary component,
3710                         // standard basis of corresponding prime,
3711                         // seperator, irreducible factors of
3712                         // the "minimal divisor" of the seperator
3713                         // as computed by the procedure minsat,
3714                         // SRest=standard basis of remaining component
3715
3716      }
3717//check
3718      for(i=size(QQ);i>=1;i--)
3719      //for(i=1;i<=size(QQ);i++)
3720      {
3721        tester=QQ[i][3]*V[4];
3722        Qi=QQ[i][2];
3723        if(NF(tester,Qi,1)!=0)  // the new vertex is not the
3724                                // root of a redundant subtree
3725        {
3726          VV[1]=QQ[i][1];
3727          VV[2]=2;
3728          VV[3]=V[3]+1;
3729          VV[4]=tester;      // the new tester as computed above
3730          VV[5]=Qi;          // QQ[i][2];
3731          VV[6]=list();
3732          W=insert(W,VV,size(W));
3733        }
3734      }
3735      if (SRest[1]!=1)        // the remaining component is not
3736                              // the whole ring
3737      {
3738        if (rad_con(V[4],SRest)==0) // the vertex is not the root
3739                                    // of a redundant subtree
3740        {
3741          VV[1]=SRest;
3742          VV[2]=3;
3743          VV[3]=V[3]+2;
3744          VV[4]=V[4];      // the tester did not change
3745          VV[5]=ideal(0);
3746          WI=list();
3747          for(i=1;i<=size(QQ);i++)
3748          {
3749            WI=insert(WI,list(QQ[i][2],QQ[i][4]));
3750          }
3751          VV[6]=WI;
3752          W=insert(W,VV,size(W));
3753        }
3754      }
3755    }
3756  }
3757}
3758
3759//////////////////////////////////////////////////////////////////////////
3760// proc pseudo_prim_dec_charsets
3761// input: Generators of an arbitrary ideal I, a standard basis SI of I,
3762// and an integer choo
3763// If choo=0, min_ass_prim_charsets with the given
3764// ordering of the variables is used.
3765// If choo=1, min_ass_prim_charsets with the "optimized"
3766// ordering of the variables is used.
3767// If choo=2, minAssPrimes from primdec.lib is used
3768// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3769// output: a pseudo primary decomposition of I, i.e., a list
3770// of pseudo primary components together with a standard basis of the
3771// remaining component. Each pseudo primary component is
3772// represented by a quadrupel: A standard basis of the component,
3773// a standard basis of the corresponding associated prime, the
3774// seperator of the component, and the irreducible factors of the
3775// "minimal divisor" of the seperator as computed by the procedure minsat,
3776// calls  proc pseudo_prim_dec_i
3777//////////////////////////////////////////////////////////////////////////
3778
3779
3780proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
3781{
3782  list L;          // The list of minimal associated primes,
3783                   // each one given by a standard basis
3784  if((choo==0) or (choo==1))
3785    {
3786       L=min_ass_prim_charsets(I,choo);
3787    }
3788  else
3789    {
3790      if(choo==2)
3791      {
3792        L=minAssPrimes(I);
3793      }
3794      else
3795      {
3796        L=minAssPrimes(I,1);
3797      }
3798      for(int i=size(L);i>=1;i=i-1)
3799        {
3800          L[i]=std(L[i]);
3801        }
3802    }
3803  return (pseudo_prim_dec_i(SI,L));
3804}
3805
3806////////////////////////////////////////////////////////////////
3807// proc pseudo_prim_dec_special_charsets
3808// input: a standard basis of an ideal I whose radical is the
3809// intersection of the radicals of ideals generated by one prime ideal
3810// P_i together with one polynomial f_i, the list V6 must be the list of
3811// pairs (standard basis of P_i, irreducible factors of f_i),
3812// and an integer choo
3813// If choo=0, min_ass_prim_charsets with the given
3814// ordering of the variables is used.
3815// If choo=1, min_ass_prim_charsets with the "optimized"
3816// ordering of the variables is used.
3817// If choo=2, minAssPrimes from primdec.lib is used
3818// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
3819// output: a pseudo primary decomposition of I, i.e., a list
3820// of pseudo primary components together with a standard basis of the
3821// remaining component. Each pseudo primary component is
3822// represented by a quadrupel: A standard basis of the component,
3823// a standard basis of the corresponding associated prime, the
3824// seperator of the component, and the irreducible factors of the
3825// "minimal divisor" of the seperator as computed by the procedure minsat,
3826// calls  proc pseudo_prim_dec_i
3827////////////////////////////////////////////////////////////////
3828
3829
3830proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
3831{
3832  int i,j,l;
3833  list m;
3834  list L;
3835  int sizeL;
3836  ideal P,SP; ideal fac;
3837  int dimSP;
3838  for(l=size(V6);l>=1;l--)   // creates a list of associated primes
3839                             // of I, possibly redundant
3840  {
3841    P=V6[l][1];
3842    fac=V6[l][2];
3843    for(i=ncols(fac);i>=1;i--)
3844    {
3845      SP=P+fac[i];
3846      SP=std(SP);
3847      if(SP[1]!=1)
3848      {
3849        if((choo==0) or (choo==1))
3850        {
3851          m=min_ass_prim_charsets(SP,choo);  // a list of SB
3852        }
3853        else
3854        {
3855          if(choo==2)
3856          {
3857            m=minAssPrimes(SP);
3858          }
3859          else
3860          {
3861            m=minAssPrimes(SP,1);
3862          }
3863          for(j=size(m);j>=1;j=j-1)
3864            {
3865              m[j]=std(m[j]);
3866            }
3867        }
3868        dimSP=dim(SP);
3869        for(j=size(m);j>=1; j--)
3870        {
3871          if(dim(m[j])==dimSP)
3872          {
3873            L=insert(L,m[j],size(L));
3874          }
3875        }
3876      }
3877    }
3878  }
3879  sizeL=size(L);
3880  for(i=1;i<sizeL;i++)     // get rid of redundant primes
3881  {
3882    for(j=i+1;j<=sizeL;j++)
3883    {
3884      if(size(L[i])!=0)
3885      {
3886        if(size(L[j])!=0)
3887        {
3888          if(size(NF(L[i],L[j],1))==0)
3889          {
3890            L[j]=ideal(0);
3891          }
3892          else
3893          {
3894            if(size(NF(L[j],L[i],1))==0)
3895            {
3896              L[i]=ideal(0);
3897            }
3898          }
3899        }
3900      }
3901    }
3902  }
3903  for(i=sizeL;i>=1;i--)
3904  {
3905    if(size(L[i])==0)
3906    {
3907      L=delete(L,i);
3908    }
3909  }
3910  return (pseudo_prim_dec_i(SI,L));
3911}
3912
3913
3914////////////////////////////////////////////////////////////////
3915// proc pseudo_prim_dec_i
3916// input: A standard basis of an arbitrary ideal I, and standard bases
3917// of the minimal associated primes of I
3918// output: a pseudo primary decomposition of I, i.e., a list
3919// of pseudo primary components together with a standard basis of the
3920// remaining component. Each pseudo primary component is
3921// represented by a quadrupel: A standard basis of the component Q_i,
3922// a standard basis of the corresponding associated prime P_i, the
3923// seperator of the component, and the irreducible factors of the
3924// "minimal divisor" of the seperator as computed by the procedure minsat,
3925////////////////////////////////////////////////////////////////
3926
3927
3928proc pseudo_prim_dec_i (ideal SI, list L)
3929{
3930  list Q;
3931  if (size(L)==1)               // one minimal associated prime only
3932                                // the ideal is already pseudo primary
3933  {
3934    Q=SI,L[1],1;
3935    list QQ;
3936    QQ[1]=Q;
3937    return (QQ,ideal(1));
3938  }
3939
3940  poly f0,f,g;
3941  ideal fac;
3942  int i,j,k,l;
3943  ideal SQi;
3944  ideal I'=SI;
3945  list QP;
3946  int sizeL=size(L);
3947  for(i=1;i<=sizeL;i++)
3948  {
3949    fac=0;
3950    for(j=1;j<=sizeL;j++)           // compute the seperator sep_i
3951                                    // of the i-th component
3952    {
3953      if (i!=j)                       // search g not in L[i], but L[j]
3954      {
3955        for(k=1;k<=ncols(L[j]);k++)
3956        {
3957          if(NF(L[j][k],L[i],1)!=0)
3958          {
3959            break;
3960          }
3961        }
3962        fac=fac+L[j][k];
3963      }
3964    }
3965    // delete superfluous polynomials
3966    fac=simplify(fac,8);
3967    // saturation
3968    SQi,f0,f,fac=minsat_ppd(SI,fac);
3969    I'=I',f;
3970    QP=SQi,L[i],f0,fac;
3972             // a standard basis of Q_i,
3973             // a standard basis of P_i,
3974             // sep_i,
3975             // irreducible factors of
3976             // the "minimal divisor" of the seperator
3977             //  as computed by the procedure minsat,
3978    Q[i]=QP;
3979  }
3980  I'=std(I');
3981  return (Q, I');
3982                   // I' = remaining component
3983}
3984
3985
3986////////////////////////////////////////////////////////////////
3987// proc extraction
3988// input: A standard basis of a pseudo primary ideal I, and a standard
3989// basis of the unique minimal associated prime P of I
3990// output: an extraction of I, i.e., a standard basis of the primary
3991// component Q of I with associated prime P, a standard basis of the
3992// remaining component, and the irreducible factors of the
3993// "minimal divisor" of the extractor as computed by the procedure minsat
3994////////////////////////////////////////////////////////////////
3995
3996
3997proc extraction (ideal SI, ideal SP)
3998{
3999  list indsets=indepSet(SP,0);
4000  poly f;
4001  if(size(indsets)!=0)      //check, whether dim P != 0
4002  {
4003    intvec v;               // a maximal independent set of variables
4004                            // modulo P
4005    string U;               // the independent variables
4006    string A;               // the dependent variables
4007    int j,k;
4008    int a;                  //  the size of A
4009    int degf;
4010    ideal g;
4011    list polys;
4012    int sizepolys;
4013    list newpoly;
4014    def R=basering;
4015    //intvec hv=hilb(SI,1);
4016    for (k=1;k<=size(indsets);k++)
4017    {
4018      v=indsets[k];
4019      for (j=1;j<=nvars(R);j++)
4020      {
4021        if (v[j]==1)
4022        {
4023          U=U+varstr(j)+",";
4024        }
4025        else
4026        {
4027          A=A+varstr(j)+",";
4028          a++;
4029        }
4030      }
4031
4032      U[size(U)]=")";           // we compute the extractor of I (w.r.t. U)
4033      execute "ring RAU="+charstr(basering)+",("+A+U+",(dp("+string(a)+"),dp);";
4034      ideal I=imap(R,SI);
4035      //I=std(I,hv);            // the standard basis in (R[U])[A]
4036      I=std(I);            // the standard basis in (R[U])[A]
4037      A[size(A)]=")";
4038      execute "ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;";
4039      ideal I=imap(RAU,I);
4040      //"std in lokalisierung:"+newline,I;
4041      ideal h;
4042      for(j=ncols(I);j>=1;j--)
4043      {
4044        h[j]=leadcoef(I[j]);  // consider I in (R(U))[A]
4045      }
4046      setring R;
4047      g=imap(Rloc,h);
4048      kill RAU,Rloc;
4049      U="";
4050      A="";
4051      a=0;
4052      f=lcm(g);
4053      newpoly[1]=f;
4054      polys=polys+newpoly;
4055      newpoly=list();
4056    }
4057    f=polys[1];
4058    degf=deg(f);
4059    sizepolys=size(polys);
4060    for (k=2;k<=sizepolys;k++)
4061    {
4062      if (deg(polys[k])<degf)
4063      {
4064        f=polys[k];
4065        degf=deg(f);
4066      }
4067    }
4068  }
4069  else
4070  {
4071    f=1;
4072  }
4073  poly f0,h0; ideal SQ; ideal fac;
4074  if(f!=1)
4075  {
4076    SQ,f0,h0,fac=minsat(SI,f);
4077    return(SQ,std(SI+h0),fac);
4078             // the tripel
4079             // a standard basis of Q,
4080             // a standard basis of remaining component,
4081             // irreducible factors of
4082             // the "minimal divisor" of the extractor
4083             // as computed by the procedure minsat
4084  }
4085  else
4086  {
4087    return(SI,ideal(1),ideal(1));
4088  }
4089}
4090
4091/////////////////////////////////////////////////////
4092// proc minsat
4093// input:  a standard basis of an ideal I and a polynomial p
4094// output: a standard basis IS of the saturation of I w.r. to p,
4095// the maximal squarefree factor f0 of p,
4096// the "minimal divisor" f of f0 such that the saturation of
4097// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
4098// the irreducible factors of f
4099//////////////////////////////////////////////////////////
4100
4101
4102proc minsat(ideal SI, poly p)
4103{
4104  ideal fac=factorize(p,1);       //the irreducible factors of p
4105  fac=sort(fac)[1];
4106  int i,k;
4107  poly f0=1;
4108  for(i=ncols(fac);i>=1;i--)
4109  {
4110    f0=f0*fac[i];
4111  }
4112  poly f=1;
4113  ideal iold;
4114  list quotM;
4115  quotM[1]=SI;
4116  quotM[2]=fac;
4117  quotM[3]=f0;
4118  // we deal seperately with the first quotient;
4119  // factors, which do not contribute to this one,
4120  // are omitted
4121  iold=quotM[1];
4122  quotM=minquot(quotM);
4123  fac=quotM[2];
4124  if(quotM[3]==1)
4125    {
4126      return(quotM[1],f0,f,fac);
4127    }
4128  while(special_ideals_equal(iold,quotM[1])==0)
4129    {
4130      f=f*quotM[3];
4131      iold=quotM[1];
4132      quotM=minquot(quotM);
4133    }
4134  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
4135}
4136
4137/////////////////////////////////////////////////////
4138// proc minsat_ppd
4139// input:  a standard basis of an ideal I and a polynomial p
4140// output: a standard basis IS of the saturation of I w.r. to p,
4141// the maximal squarefree factor f0 of p,
4142// the "minimal divisor" f of f0 such that the saturation of
4143// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
4144// the irreducible factors of f
4145//////////////////////////////////////////////////////////
4146
4147
4148proc minsat_ppd(ideal SI, ideal fac)
4149{
4150  fac=sort(fac)[1];
4151  int i,k;
4152  poly f0=1;
4153  for(i=ncols(fac);i>=1;i--)
4154  {
4155    f0=f0*fac[i];
4156  }
4157  poly f=1;
4158  ideal iold;
4159  list quotM;
4160  quotM[1]=SI;
4161  quotM[2]=fac;
4162  quotM[3]=f0;
4163  // we deal seperately with the first quotient;
4164  // factors, which do not contribute to this one,
4165  // are omitted
4166  iold=quotM[1];
4167  quotM=minquot(quotM);
4168  fac=quotM[2];
4169  if(quotM[3]==1)
4170    {
4171      return(quotM[1],f0,f,fac);
4172    }
4173  while(special_ideals_equal(iold,quotM[1])==0)
4174  {
4175    f=f*quotM[3];
4176    iold=quotM[1];
4177    quotM=minquot(quotM);
4178    k++;
4179  }
4180  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
4181}
4182/////////////////////////////////////////////////////////////////
4183// proc minquot
4184// input: a list with 3 components: a standard basis
4185// of an ideal I, a set of irreducible polynomials, and
4186// there product f0
4187// output: a standard basis of the ideal (I:f0), the irreducible
4188// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
4189// the "minimal divisor" f
4190/////////////////////////////////////////////////////////////////
4191
4192proc minquot(list tsil)
4193{
4194   int i,j,k,action;
4195   ideal verg;
4196   list l;
4197   poly g;
4198   ideal laedi=tsil[1];
4199   ideal fac=tsil[2];
4200   poly f=tsil[3];
4201
4202//std
4203//   ideal star=quotient(laedi,f);
4204//   star=std(star);
4205   option(returnSB);
4206   ideal star=quotient(laedi,f);
4207   option(noreturnSB);
4208   if(special_ideals_equal(laedi,star)==1)
4209     {
4210       return(laedi,ideal(1),1);
4211     }
4212   action=1;
4213   while(action==1)
4214   {
4215      if(size(fac)==1)
4216      {
4217         action=0;
4218         break;
4219      }
4220      for(i=1;i<=size(fac);i++)
4221      {
4222        g=1;
4223         for(j=1;j<=size(fac);j++)
4224         {
4225            if(i!=j)
4226            {
4227               g=g*fac[j];
4228            }
4229         }
4230//std
4231//         verg=quotient(laedi,g);
4232//         verg=std(verg);
4233         option(returnSB);
4234         verg=quotient(laedi,g);
4235         option(noreturnSB);
4236         if(special_ideals_equal(verg,star)==1)
4237         {
4238            f=g;
4239            fac[i]=0;
4240            fac=simplify(fac,2);
4241            break;
4242         }
4243         if(i==size(fac))
4244         {
4245            action=0;
4246         }
4247      }
4248   }
4249   l=star,fac,f;
4250   return(l);
4251}
4252/////////////////////////////////////////////////
4253// proc special_ideals_equal
4254// input: standard bases of ideal k1 and k2 such that
4255// k1 is contained in k2, or k2 is contained ink1
4256// output: 1, if k1 equals k2, 0 otherwise
4257//////////////////////////////////////////////////
4258
4259proc special_ideals_equal( ideal k1, ideal k2)
4260{
4261   int j;
4262   if(size(k1)==size(k2))
4263   {
4264      for(j=1;j<=size(k1);j++)
4265      {
4267         {
4268            return(0);
4269         }
4270      }
4271      return(1);
4272   }
4273   return(0);
4274}
4275
4276
4277 proc quotient2(module a,module b)
4278{
4279  string s="ring @newr=("+charstr(basering)+"),(@t,"+varstr(basering)+"),dp;";
4280  def @newP=basering;
4281  execute s;
4282  module b=imap(@newP,b);
4283  module a=imap(@newP,a);
4284 int i;
4285 vector @b=b[1];
4286 for(i=2;i<=size(b);i++)
4287 {
4288    @b=@t^(i-1)*b[i];
4289 }
4290 ideal re=quotient(a,module(@b));
4291  setring @newP;
4292 ideal re=imap(@newr,re);
4293 return(re);
4294}
4295//Im homogenen Fall system("LaScala",i) verwenden
Note: See TracBrowser for help on using the repository browser.