source: git/Singular/LIB/primdec.lib @ 8677912

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