source: git/Singular/LIB/primdec.lib @ 0b59f5

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