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

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