source: git/Singular/LIB/primdec.lib @ 76aca2

spielwiese
Last change on this file since 76aca2 was 76aca2, checked in by Hans Schönemann <hannes@…>, 19 years ago
*lossen: Ann, headStand: help, example git-svn-id: file:///usr/local/Singular/svn/trunk@8080 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 136.8 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: primdec.lib,v 1.105 2005-05-06 16:01:50 Singular Exp $";
3category="Commutative Algebra";
4info="
5LIBRARY: primdec.lib   Primary Decomposition and Radical of Ideals
6AUTHORS:  Gerhard Pfister, pfister@mathematik.uni-kl.de (GTZ)
7@*        Wolfram Decker, decker@math.uni-sb.de         (SY)
8@*        Hans Schoenemann, hannes@mathematik.uni-kl.de (SY)
9
10OVERVIEW:
11    Algorithms for primary decomposition based on the ideas of
12    Gianni, Trager and Zacharias (implementation by Gerhard Pfister),
13    respectively based on the ideas of Shimoyama and Yokoyama (implementation
14    by Wolfram Decker and Hans Schoenemann).
15 @* The procedures are implemented to be used in characteristic 0.
16 @* They also work in positive characteristic >> 0.
17 @* In small characteristic and for algebraic extensions, primdecGTZ
18    may not terminate.
19    Algorithms for the computation of the radical based on the ideas of
20    Krick, Logar and Kemper (implementation by Gerhard Pfister).
21
22PROCEDURES:
23 Ann(M);           annihilator of R^n/M, R=basering, M in R^n
24 primdecGTZ(I);    complete primary decomposition via Gianni,Trager,Zacharias
25 primdecSY(I...);  complete primary decomposition via Shimoyama-Yokoyama
26 minAssGTZ(I);     the minimal associated primes via Gianni,Trager,Zacharias
27 minAssChar(I...); the minimal associated primes using characteristic sets
28 testPrimary(L,k); tests the result of the primary decomposition
29 radical(I);       computes the radical of I via  Krick/Logar and Kemper
30 radicalEHV(I);    computes the radical of I via  Eisenbud,Huneke,Vasconcelos
31 equiRadical(I);   the radical of the equidimensional part of the ideal I
32 prepareAss(I);    list of radicals of the equidimensional components of I
33 equidim(I);       weak equidimensional decomposition of I
34 equidimMax(I);    equidimensional locus of I
35 equidimMaxEHV(I); equidimensional locus of I via Eisenbud,Huneke,Vasconcelos
36 zerodec(I);       zerodimensional decomposition via Monico
37";
38
39LIB "general.lib";
40LIB "elim.lib";
41LIB "poly.lib";
42LIB "random.lib";
43LIB "inout.lib";
44LIB "matrix.lib";
45LIB "triang.lib";
46///////////////////////////////////////////////////////////////////////////////
47//
48//                      Gianni/Trager/Zacharias
49//
50///////////////////////////////////////////////////////////////////////////////
51
52static proc sat1 (ideal id, poly p)
53"USAGE:   sat1(id,j);  id ideal, j polynomial
54RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
55NOTE:    result is a std basis in the basering
56"
57{
58   int @k;
59   ideal inew=std(id);
60   ideal iold;
61   intvec op=option(get);
62   option(returnSB);
63   while(specialIdealsEqual(iold,inew)==0 )
64   {
65      iold=inew;
66      inew=quotient(iold,p);
67      @k++;
68   }
69   @k--;
70   option(set,op);
71   list L =inew,p^@k;
72   return (L);
73}
74
75///////////////////////////////////////////////////////////////////////////////
76
77static proc sat2 (ideal id, ideal h)
78"USAGE:   sat2(id,j);  id ideal, j polynomial
79RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
80NOTE:    result is a std basis in the basering
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      for(@i=1;@i<=size(fac);@i++)
111      {
112        @f=@f*fac[@i];
113     }
114     intvec op = option(get);
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(set,op);
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
145
146proc minSat(ideal inew, ideal h)
147{
148   int i,k;
149   poly f=1;
150   ideal iold,fac;
151   list quotM,l;
152
153   for(i=1;i<=ncols(h);i++)
154   {
155      if(deg(h[i])>0)
156      {
157         fac=fac+factorize(h[i],1);
158      }
159   }
160   fac=simplify(fac,4);
161   if(size(fac)==0)
162   {
163      l=inew,1;
164      return(l);
165   }
166   fac=sort(fac)[1];
167   for(i=1;i<=size(fac);i++)
168   {
169      f=f*fac[i];
170   }
171   quotM[1]=inew;
172   quotM[2]=fac;
173   quotM[3]=f;
174   f=1;
175   intvec op = option(get);
176   option(returnSB);
177   while(specialIdealsEqual(iold,quotM[1])==0)
178   {
179      if(k>0)
180      {
181         f=f*quotM[3];
182      }
183      iold=quotM[1];
184      quotM=quotMin(quotM);
185      k++;
186   }
187   option(set,op);
188   l=quotM[1],f;
189   return(l);
190}
191
192static proc quotMin(list tsil)
193{
194   int i,j,k,action;
195   ideal verg;
196   list l;
197   poly g;
198
199   ideal laedi=tsil[1];
200   ideal fac=tsil[2];
201   poly f=tsil[3];
202
203   ideal star=quotient(laedi,f);
204
205   if(specialIdealsEqual(star,laedi))
206   {
207      l=star,fac,f;
208      return(l);
209   }
210
211   action=1;
212
213   while(action==1)
214   {
215      if(size(fac)==1)
216      {
217         action=0;
218         break;
219      }
220      for(i=1;i<=size(fac);i++)
221      {
222        g=1;
223        verg=laedi;
224
225         for(j=1;j<=size(fac);j++)
226         {
227            if(i!=j)
228            {
229               g=g*fac[j];
230            }
231         }
232
233         verg=quotient(laedi,g);
234
235         if(specialIdealsEqual(verg,star)==1)
236         {
237            f=g;
238            fac[i]=0;
239            fac=simplify(fac,2);
240            break;
241         }
242         if(i==size(fac))
243         {
244            action=0;
245         }
246      }
247   }
248   l=star,fac,f;
249   return(l);
250}
251
252///////////////////////////////////////////////////////////////////////////////
253
254static proc testFactor(list act,poly p)
255{
256  poly keep=p;
257
258   int i;
259   poly q=act[1][1]^act[2][1];
260   for(i=2;i<=size(act[1]);i++)
261   {
262      q=q*act[1][i]^act[2][i];
263   }
264   q=1/leadcoef(q)*q;
265   p=1/leadcoef(p)*p;
266   if(p-q!=0)
267   {
268      "ERROR IN FACTOR, please inform the authors";
269   }
270}
271///////////////////////////////////////////////////////////////////////////////
272
273static proc factor(poly p)
274"USAGE:   factor(p) p poly
275RETURN:  list=;
276NOTE:
277EXAMPLE: example factor; shows an example
278"
279{
280
281  ideal @i;
282  list @l;
283  intvec @v,@w;
284  int @j,@k,@n;
285
286  if(deg(p)<=1)
287  {
288     @i=ideal(p);
289     @v=1;
290     @l[1]=@i;
291     @l[2]=@v;
292     return(@l);
293  }
294  if (size(p)==1)
295  {
296     @w=leadexp(p);
297     for(@j=1;@j<=nvars(basering);@j++)
298     {
299        if(@w[@j]!=0)
300        {
301           @k++;
302           @v[@k]=@w[@j];
303           @i=@i+ideal(var(@j));
304        }
305     }
306     @l[1]=@i;
307     @l[2]=@v;
308     return(@l);
309  }
310  //  @l=factorize(p,2);
311  @l=factorize(p);
312  // if(npars(basering)>0)
313  // {
314     for(@j=1;@j<=size(@l[1]);@j++)
315     {
316        if(deg(@l[1][@j])==0)
317        {
318           @n++;
319        }
320     }
321     if(@n>0)
322     {
323        if(@n==size(@l[1]))
324        {
325           @l[1]=ideal(1);
326           @v=1;
327           @l[2]=@v;
328        }
329        else
330        {
331           @k=0;
332           int pleh;
333           for(@j=1;@j<=size(@l[1]);@j++)
334           {
335              if(deg(@l[1][@j])!=0)
336              {
337                 @k++;
338                 @i=@i+ideal(@l[1][@j]);
339                 if(size(@i)==pleh)
340                 {
341                   "//factorization error";
342                   @l;
343                    @k--;
344                    @v[@k]=@v[@k]+@l[2][@j];
345                 }
346                 else
347                 {
348                    pleh++;
349                    @v[@k]=@l[2][@j];
350                 }
351              }
352           }
353           @l[1]=@i;
354           @l[2]=@v;
355        }
356     }
357     // }
358  return(@l);
359}
360example
361{ "EXAMPLE:"; echo = 2;
362   ring  r = 0,(x,y,z),lp;
363   poly  p = (x+y)^2*(y-z)^3;
364   list  l = factor(p);
365   l;
366   ring r1 =(0,b,d,f,g),(a,c,e),lp;
367   poly p  =(1*d)*e^2+(1*d*f^2*g);
368   list  l = factor(p);
369   l;
370   ring r2 =(0,b,f,g),(a,c,e,d),lp;
371   poly p  =(1*d)*e^2+(1*d*f^2*g);
372   list  l = factor(p);
373   l;
374}
375
376///////////////////////////////////////////////////////////////////////////////
377
378proc idealsEqual( ideal k, ideal j)
379{
380   return(stdIdealsEqual(std(k),std(j)));
381}
382
383static proc specialIdealsEqual( ideal k1, ideal k2)
384{
385   int j;
386
387   if(size(k1)==size(k2))
388   {
389      for(j=1;j<=size(k1);j++)
390      {
391         if(leadexp(k1[j])!=leadexp(k2[j]))
392         {
393            return(0);
394         }
395      }
396      return(1);
397   }
398   return(0);
399}
400
401static proc stdIdealsEqual( ideal k1, ideal k2)
402{
403   int j;
404
405   if(size(k1)==size(k2))
406   {
407      for(j=1;j<=size(k1);j++)
408      {
409         if(leadexp(k1[j])!=leadexp(k2[j]))
410         {
411            return(0);
412         }
413      }
414      attrib(k2,"isSB",1);
415      if(size(reduce(k1,k2,1))==0)
416      {
417         return(1);
418      }
419   }
420   return(0);
421}
422///////////////////////////////////////////////////////////////////////////////
423
424proc primaryTest (ideal i, poly p)
425{
426   int m=1;
427   int n=nvars(basering);
428   int e,f;
429   poly t;
430   ideal h;
431   list act;
432
433   ideal prm=p;
434   attrib(prm,"isSB",1);
435
436   while (n>1)
437   {
438      n=n-1;
439      m=m+1;
440
441      //search for i[m] which has a power of var(n) as leading term
442      if (n==1)
443      {
444         m=size(i);
445      }
446      else
447      {
448        while (lead(i[m])/var(n-1)==0)
449        {
450          m=m+1;
451        }
452        m=m-1;
453      }
454      //check whether i[m] =(c*var(n)+h)^e modulo prm for some
455      //h in K[var(n+1),...,var(nvars(basering))], c in K
456      //if not (0) is returned, else var(n)+h is added to prm
457
458      e=deg(lead(i[m]));
459      if(char(basering)!=0)
460      {
461         f=1;
462         if(e mod char(basering)==0)
463         {
464           if ( voice >=15 )
465           {
466              "// WARNING: The characteristic is perhaps too small to use";
467              "// the algorithm of Gianni/Trager/Zacharias.";
468              "// This may result in an infinte loop";
469              "// loop in primaryTest, voice:",voice;"";
470           }
471           while(e mod char(basering)==0)
472           {
473              f=f*char(basering);
474              e=e/char(basering);
475           }
476
477         }
478         t=leadcoef(i[m])*e*var(n)^f+(i[m]-lead(i[m]))/var(n)^((e-1)*f);
479         i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
480         if (reduce(i[m]-t^e,prm,1) !=0)
481         {
482           return(ideal(0));
483         }
484         if(f>1)
485         {
486            act=factorize(t);
487            if(size(act[1])>2)
488            {
489              return(ideal(0));
490            }
491            if(deg(act[1][2])>1)
492            {
493              return(ideal(0));
494            }
495            t=act[1][2];
496         }
497      }
498      else
499      {
500         t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
501         i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
502         if (reduce(i[m]-t^e,prm,1) !=0)
503         {
504           return(ideal(0));
505         }
506      }
507
508      h=interred(t);
509      t=h[1];
510
511      prm = prm,t;
512      attrib(prm,"isSB",1);
513   }
514   return(prm);
515}
516
517///////////////////////////////////////////////////////////////////////////////
518proc gcdTest(ideal act)
519{
520  int i,j;
521  if(size(act)<=1)
522  {
523     return(0);
524  }
525  for (i=1;i<=size(act)-1;i++)
526  {
527     for(j=i+1;j<=size(act);j++)
528     {
529        if(deg(std(ideal(act[i],act[j]))[1])>0)
530        {
531           return(0);
532        }
533     }
534  }
535  return(1);
536}
537
538///////////////////////////////////////////////////////////////////////////////
539static proc splitPrimary(list l,ideal ser,int @wr,list sact)
540{
541   int i,j,k,s,r,w;
542   list keepresult,act,keepprime;
543   poly @f;
544   int sl=size(l);
545   for(i=1;i<=sl/2;i++)
546   {
547      if(sact[2][i]>1)
548      {
549         keepprime[i]=l[2*i-1]+ideal(sact[1][i]);
550      }
551      else
552      {
553         keepprime[i]=l[2*i-1];
554      }
555  }
556   i=0;
557   while(i<size(l)/2)
558   {
559      i++;
560      if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0))
561      {
562         l[2*i-1]=ideal(1);
563         l[2*i]=ideal(1);
564         continue;
565      }
566
567      if(size(l[2*i])==0)
568      {
569         if(homog(l[2*i-1])==1)
570         {
571            l[2*i]=maxideal(1);
572            continue;
573         }
574         j=0;
575         if(i<=sl/2)
576         {
577            j=1;
578         }
579         while(j<size(l[2*i-1]))
580         {
581            j++;
582            act=factor(l[2*i-1][j]);
583            r=size(act[1]);
584            attrib(l[2*i-1],"isSB",1);
585            if((r==1)&&(vdim(l[2*i-1])==deg(l[2*i-1][j])))
586            {
587              l[2*i]=std(l[2*i-1],act[1][1]);
588              break;
589            }
590            if((r==1)&&(act[2][1]>1))
591            {
592               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
593               if(homog(keepprime[i])==1)
594               {
595                  l[2*i]=maxideal(1);
596                  break;
597               }
598            }
599            if(gcdTest(act[1])==1)
600            {
601               for(k=2;k<=r;k++)
602               {
603                  keepprime[size(l)/2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
604               }
605               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
606               for(k=1;k<=r;k++)
607               {
608                  if(@wr==0)
609                  {
610                     keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]);
611                  }
612                  else
613                  {
614                     keepresult[k]=std(l[2*i-1],act[1][k]);
615                  }
616               }
617               l[2*i-1]=keepresult[1];
618               if(vdim(keepresult[1])==deg(act[1][1]))
619               {
620                  l[2*i]=keepresult[1];
621               }
622               if((homog(keepresult[1])==1)||(homog(keepprime[i])==1))
623               {
624                  l[2*i]=maxideal(1);
625               }
626               s=size(l)-2;
627               for(k=2;k<=r;k++)
628               {
629                  l[s+2*k-1]=keepresult[k];
630                  keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k]));
631                  if(vdim(keepresult[k])==deg(act[1][k]))
632                  {
633                     l[s+2*k]=keepresult[k];
634                  }
635                  else
636                  {
637                     l[s+2*k]=ideal(0);
638                  }
639                  if((homog(keepresult[k])==1)||(homog(keepprime[s/2+k])==1))
640                  {
641                    l[s+2*k]=maxideal(1);
642                  }
643               }
644               i--;
645               break;
646            }
647            if(r>=2)
648            {
649               s=size(l);
650               @f=act[1][1];
651               act=sat1(l[2*i-1],act[1][1]);
652               if(deg(act[1][1])>0)
653               {
654                  l[s+1]=std(l[2*i-1],act[2]);
655                  if(homog(l[s+1])==1)
656                  {
657                     l[s+2]=maxideal(1);
658                  }
659                  else
660                  {
661                     l[s+2]=ideal(0);
662                  }
663                  keepprime[s/2+1]=interred(keepprime[i]+ideal(@f));
664                  if(homog(keepprime[s/2+1])==1)
665                  {
666                     l[s+2]=maxideal(1);
667                  }
668                  keepprime[i]=act[1];
669                  l[2*i-1]=act[1];
670                  attrib(l[2*i-1],"isSB",1);
671                  if(homog(l[2*i-1])==1)
672                  {
673                     l[2*i]=maxideal(1);
674                  }
675
676                  i--;
677                  break;
678               }
679            }
680         }
681      }
682   }
683   if(sl==size(l))
684   {
685      return(l);
686   }
687   for(i=1;i<=size(l)/2;i++)
688   {
689     attrib(l[2*i-1],"isSB",1);
690
691     if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0)&&(deg(l[2*i-1][1])>0))
692     {
693       "Achtung in split";
694
695         l[2*i-1]=ideal(1);
696         l[2*i]=ideal(1);
697     }
698     if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
699      {
700         keepprime[i]=std(keepprime[i]);
701         if(homog(keepprime[i])==1)
702         {
703             l[2*i]=maxideal(1);
704         }
705         else
706         {
707            act=zero_decomp(keepprime[i],ideal(0),@wr,1);
708            if(size(act)==2)
709            {
710               l[2*i]=act[2];
711            }
712         }
713      }
714   }
715   return(l);
716}
717example
718{ "EXAMPLE:"; echo=2;
719   ring  r = 32003,(x,y,z),lp;
720   ideal i1=x*(x+1),yz,(z+1)*(z-1);
721   ideal i2=xy,yz,(x-2)*(x+3);
722   list l=i1,ideal(0),i2,ideal(0),i2,ideal(1);
723   list l1=splitPrimary(l,ideal(0),0);
724   l1;
725}
726///////////////////////////////////////////////////////////////////////////////
727static proc splitCharp(list l)
728{
729  if((char(basering)==0)||(npars(basering)>0))
730  {
731     return(l);
732  }
733  def P=basering;
734  int i,j,k,m,q,d,o;
735  int n=nvars(basering);
736  ideal s,t,u,sact;
737  poly ni;
738  string minp,gnir,va;
739  list sa,keep,rp,keep1;
740  for(i=1;i<=size(l)/2;i++)
741  {
742    if(size(l[2*i])==0)
743    {
744       if(deg(l[2*i-1][1])==vdim(l[2*i-1]))
745       {
746          l[2*i]=l[2*i-1];
747       }
748    }
749  }
750  for(i=1;i<=size(l)/2;i++)
751  {
752    if(size(l[2*i])==0)
753    {
754      s=factorize(l[2*i-1][1],1);   //vermeiden!!!
755      t=l[2*i-1];
756      m=size(t);
757      ni=s[1];
758      if(deg(ni)>1)
759      {
760        va=varstr(P);
761        j=size(va);
762        while(va[j]!=","){j--;}
763        va=va[1..j-1];
764        gnir="ring RL=("+string(char(P))+","+string(var(n))+"),("+va+"),lp;";
765        execute(gnir);
766        minpoly=leadcoef(imap(P,ni));
767        ideal act;
768        ideal t=imap(P,t);
769
770        for(k=2;k<=m;k++)
771        {
772           act=factorize(t[k],1);
773           if(size(act)>1){break;}
774        }
775        setring P;
776        sact=imap(RL,act);
777
778        if(size(sact)>1)
779        {
780           sa=sat1(l[2*i-1],sact[1]);
781           keep[size(keep)+1]=std(l[2*i-1],sa[2]);
782           l[2*i-1]=std(sa[1]);
783           l[2*i]=primaryTest(sa[1],sa[1][1]);
784        }
785        if((size(sact)==1)&&(m==2))
786        {
787           l[2*i]=l[2*i-1];
788           attrib(l[2*i],"isSB",1);
789        }
790        if((size(sact)==1)&&(m>2))
791        {
792           setring RL;
793           option(redSB);
794           t=std(t);
795
796           list sp=zero_decomp(t,0,0);
797
798           setring P;
799           rp=imap(RL,sp);
800           for(o=1;o<=size(rp);o++)
801           {
802              rp[o]=interred(simplify(rp[o],1)+ideal(ni));
803           }
804           l[2*i-1]=rp[1];
805           l[2*i]=rp[2];
806           rp=delete(rp,1);
807           rp=delete(rp,1);
808           keep1=keep1+rp;
809           option(noredSB);
810        }
811        kill RL;
812      }
813    }
814  }
815  if(size(keep)>0)
816  {
817    for(i=1;i<=size(keep);i++)
818    {
819       if(deg(keep[i][1])>0)
820       {
821          l[size(l)+1]=keep[i];
822          l[size(l)+1]=primaryTest(keep[i],keep[i][1]);
823       }
824    }
825  }
826  l=l+keep1;
827  return(l);
828}
829
830///////////////////////////////////////////////////////////////////////////////
831
832proc zero_decomp (ideal j,ideal ser,int @wr,list #)
833"USAGE:   zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
834         (@wr=0 for primary decomposition, @wr=1 for computaion of associated
835         primes)
836RETURN:  list = list of primary ideals and their radicals (at even positions
837         in the list) if the input is zero-dimensional and a standardbases
838         with respect to lex-ordering
839         If ser!=(0) and ser is contained in j or if j is not zero-dimen-
840         sional then ideal(1),ideal(1) is returned
841NOTE:    Algorithm of Gianni/Trager/Zacharias
842EXAMPLE: example zero_decomp; shows an example
843"
844{
845  def   @P = basering;
846  int uytrewq;
847  int nva = nvars(basering);
848  int @k,@s,@n,@k1,zz;
849  list primary,lres0,lres1,act,@lh,@wh;
850  map phi,psi,phi1,psi1;
851  ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1;
852  intvec @vh,@hilb;
853  string @ri;
854  poly @f;
855  if (dim(j)>0)
856  {
857     primary[1]=ideal(1);
858     primary[2]=ideal(1);
859     return(primary);
860  }
861  j=interred(j);
862
863  attrib(j,"isSB",1);
864
865  if(vdim(j)==deg(j[1]))
866  {
867     act=factor(j[1]);
868     for(@k=1;@k<=size(act[1]);@k++)
869     {
870       @qh=j;
871       if(@wr==0)
872       {
873          @qh[1]=act[1][@k]^act[2][@k];
874       }
875       else
876       {
877          @qh[1]=act[1][@k];
878       }
879       primary[2*@k-1]=interred(@qh);
880       @qh=j;
881       @qh[1]=act[1][@k];
882       primary[2*@k]=interred(@qh);
883       attrib( primary[2*@k-1],"isSB",1);
884
885       if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
886       {
887          primary[2*@k-1]=ideal(1);
888          primary[2*@k]=ideal(1);
889       }
890     }
891     return(primary);
892  }
893
894  if(homog(j)==1)
895  {
896     primary[1]=j;
897     if((size(ser)>0)&&(size(reduce(ser,j,1))==0))
898     {
899          primary[1]=ideal(1);
900          primary[2]=ideal(1);
901          return(primary);
902     }
903     if(dim(j)==-1)
904     {
905        primary[1]=ideal(1);
906        primary[2]=ideal(1);
907     }
908     else
909     {
910        primary[2]=maxideal(1);
911     }
912     return(primary);
913  }
914
915//the first element in the standardbase is factorized
916  if(deg(j[1])>0)
917  {
918    act=factor(j[1]);
919    testFactor(act,j[1]);
920  }
921  else
922  {
923     primary[1]=ideal(1);
924     primary[2]=ideal(1);
925     return(primary);
926  }
927
928//with the factors new ideals (hopefully the primary decomposition)
929//are created
930  if(size(act[1])>1)
931  {
932     if(size(#)>1)
933     {
934        primary[1]=ideal(1);
935        primary[2]=ideal(1);
936        primary[3]=ideal(1);
937        primary[4]=ideal(1);
938        return(primary);
939     }
940     for(@k=1;@k<=size(act[1]);@k++)
941     {
942        if(@wr==0)
943        {
944           primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
945
946        }
947        else
948        {
949           primary[2*@k-1]=std(j,act[1][@k]);
950        }
951        if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
952        {
953           primary[2*@k]   = primary[2*@k-1];
954        }
955        else
956        {
957
958           primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
959
960        }
961     }
962  }
963  else
964  {
965     primary[1]=j;
966     if((size(#)>0)&&(act[2][1]>1))
967     {
968        act[2]=1;
969        primary[1]=std(primary[1],act[1][1]);
970     }
971
972     if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
973     {
974        primary[2]=primary[1];
975     }
976     else
977     {
978        primary[2]=primaryTest(primary[1],act[1][1]);
979     }
980  }
981
982  if(size(#)==0)
983  {
984     primary=splitPrimary(primary,ser,@wr,act);
985  }
986
987  if((voice>=6)&&(char(basering)<=181))
988  {
989     primary=splitCharp(primary);
990  }
991
992  if((@wr==2)&&(npars(basering)>0)&&(voice>=6)&&(char(basering)>0))
993  {
994  //the prime decomposition of Yokoyama in characteristic p
995     list ke,ek;
996     @k=0;
997     while(@k<size(primary)/2)
998     {
999        @k++;
1000        if(size(primary[2*@k])==0)
1001        {
1002           ek=insepDecomp(primary[2*@k-1]);
1003           primary=delete(primary,2*@k);
1004           primary=delete(primary,2*@k-1);
1005           @k--;
1006        }
1007        ke=ke+ek;
1008     }
1009     for(@k=1;@k<=size(ke);@k++)
1010     {
1011        primary[size(primary)+1]=ke[@k];
1012        primary[size(primary)+1]=ke[@k];
1013     }
1014  }
1015
1016  if(voice>=8){primary=extF(primary)};
1017
1018//test whether all ideals in the decomposition are primary and
1019//in general position
1020//if not after a random coordinate transformation of the last
1021//variable the corresponding ideal is decomposed again.
1022  if((npars(basering)>0)&&(voice>=8))
1023  {
1024     poly randp;
1025     for(zz=1;zz<nvars(basering);zz++)
1026     {
1027        randp=randp
1028              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(zz);
1029     }
1030     randp=randp+var(nvars(basering));
1031  }
1032  @k=0;
1033  while(@k<(size(primary)/2))
1034  {
1035    @k++;
1036    if (size(primary[2*@k])==0)
1037    {
1038       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
1039       {
1040          attrib(primary[2*@k-1],"isSB",1);
1041          if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
1042          {
1043             primary[2*@k]=primary[2*@k-1];
1044          }
1045       }
1046    }
1047  }
1048
1049  @k=0;
1050  ideal keep;
1051  while(@k<(size(primary)/2))
1052  {
1053    @k++;
1054    if (size(primary[2*@k])==0)
1055    {
1056
1057       jmap=randomLast(100);
1058       jmap1=maxideal(1);
1059       jmap2=maxideal(1);
1060       @qht=primary[2*@k-1];
1061       if((npars(basering)>0)&&(voice>=10))
1062       {
1063          jmap[size(jmap)]=randp;
1064       }
1065
1066       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
1067       {
1068          if(deg(lead(primary[2*@k-1][@n]))==1)
1069          {
1070             for(zz=1;zz<=nva;zz++)
1071             {
1072                if(lead(primary[2*@k-1][@n])/var(zz)!=0)
1073                {
1074                 jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
1075                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
1076                    jmap2[zz]=primary[2*@k-1][@n];
1077                    @qht[@n]=var(zz);
1078
1079                }
1080             }
1081             jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0);
1082          }
1083       }
1084       if(size(subst(jmap[nva],var(1),0)-var(nva))!=0)
1085       {
1086         // jmap[nva]=subst(jmap[nva],var(1),0);
1087         //hier geaendert +untersuchen!!!!!!!!!!!!!!
1088       }
1089       phi1=@P,jmap1;
1090       phi=@P,jmap;
1091       for(@n=1;@n<=nva;@n++)
1092       {
1093          jmap[@n]=-(jmap[@n]-2*var(@n));
1094       }
1095       psi=@P,jmap;
1096       psi1=@P,jmap2;
1097       @qh=phi(@qht);
1098
1099//=================== the new part ============================
1100
1101       @qh=groebner(@qh);
1102
1103//=============================================================
1104//       if(npars(@P)>0)
1105//       {
1106//          @ri= "ring @Phelp ="
1107//                  +string(char(@P))+",
1108//                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
1109//       }
1110//       else
1111//       {
1112//          @ri= "ring @Phelp ="
1113//                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
1114//       }
1115//       execute(@ri);
1116//       ideal @qh=homog(imap(@P,@qht),@t);
1117//
1118//       ideal @qh1=std(@qh);
1119//       @hilb=hilb(@qh1,1);
1120//       @ri= "ring @Phelp1 ="
1121//                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
1122//       execute(@ri);
1123//       ideal @qh=homog(imap(@P,@qh),@t);
1124//       kill @Phelp;
1125//       @qh=std(@qh,@hilb);
1126//       @qh=subst(@qh,@t,1);
1127//       setring @P;
1128//       @qh=imap(@Phelp1,@qh);
1129//       kill @Phelp1;
1130//       @qh=clearSB(@qh);
1131//       attrib(@qh,"isSB",1);
1132//=============================================================
1133
1134       ser1=phi1(ser);
1135
1136       @lh=zero_decomp (@qh,phi(ser1),@wr);
1137
1138
1139       kill lres0;
1140       list lres0;
1141       if(size(@lh)==2)
1142       {
1143          helpprim=@lh[2];
1144          lres0[1]=primary[2*@k-1];
1145          ser1=psi(helpprim);
1146          lres0[2]=psi1(ser1);
1147          if(size(reduce(lres0[2],lres0[1],1))==0)
1148          {
1149             primary[2*@k]=primary[2*@k-1];
1150             continue;
1151          }
1152       }
1153       else
1154       {
1155
1156          lres1=psi(@lh);
1157          lres0=psi1(lres1);
1158       }
1159
1160//=================== the new part ============================
1161
1162       primary=delete(primary,2*@k-1);
1163       primary=delete(primary,2*@k-1);
1164       @k--;
1165       if(size(lres0)==2)
1166       {
1167          lres0[2]=groebner(lres0[2]);
1168       }
1169       else
1170       {
1171          for(@n=1;@n<=size(lres0)/2;@n++)
1172          {
1173             if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1)
1174             {
1175                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
1176                lres0[2*@n]=lres0[2*@n-1];
1177                attrib(lres0[2*@n],"isSB",1);
1178             }
1179             else
1180             {
1181                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
1182                lres0[2*@n]=groebner(lres0[2*@n]);
1183             }
1184          }
1185       }
1186       primary=primary+lres0;
1187
1188//=============================================================
1189//       if(npars(@P)>0)
1190//       {
1191//          @ri= "ring @Phelp ="
1192//                  +string(char(@P))+",
1193//                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
1194//       }
1195//       else
1196//       {
1197//          @ri= "ring @Phelp ="
1198//                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
1199//       }
1200//       execute(@ri);
1201//       list @lvec;
1202//       list @lr=imap(@P,lres0);
1203//       ideal @lr1;
1204//
1205//       if(size(@lr)==2)
1206//       {
1207//          @lr[2]=homog(@lr[2],@t);
1208//          @lr1=std(@lr[2]);
1209//          @lvec[2]=hilb(@lr1,1);
1210//       }
1211//       else
1212//       {
1213//          for(@n=1;@n<=size(@lr)/2;@n++)
1214//          {
1215//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1216//             {
1217//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1218//                @lr1=std(@lr[2*@n-1]);
1219//                @lvec[2*@n-1]=hilb(@lr1,1);
1220//                @lvec[2*@n]=@lvec[2*@n-1];
1221//             }
1222//             else
1223//             {
1224//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1225//                @lr1=std(@lr[2*@n-1]);
1226//                @lvec[2*@n-1]=hilb(@lr1,1);
1227//                @lr[2*@n]=homog(@lr[2*@n],@t);
1228//                @lr1=std(@lr[2*@n]);
1229//                @lvec[2*@n]=hilb(@lr1,1);
1230//
1231//             }
1232//         }
1233//       }
1234//       @ri= "ring @Phelp1 ="
1235//                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
1236//       execute(@ri);
1237//       list @lr=imap(@Phelp,@lr);
1238//
1239//       kill @Phelp;
1240//       if(size(@lr)==2)
1241//      {
1242//          @lr[2]=std(@lr[2],@lvec[2]);
1243//          @lr[2]=subst(@lr[2],@t,1);
1244//
1245//       }
1246//       else
1247//       {
1248//          for(@n=1;@n<=size(@lr)/2;@n++)
1249//          {
1250//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1251//             {
1252//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1253//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1254//                @lr[2*@n]=@lr[2*@n-1];
1255//                attrib(@lr[2*@n],"isSB",1);
1256//             }
1257//             else
1258//             {
1259//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1260//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1261//                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
1262//                @lr[2*@n]=subst(@lr[2*@n],@t,1);
1263//             }
1264//          }
1265//       }
1266//       kill @lvec;
1267//       setring @P;
1268//       lres0=imap(@Phelp1,@lr);
1269//       kill @Phelp1;
1270//       for(@n=1;@n<=size(lres0);@n++)
1271//       {
1272//          lres0[@n]=clearSB(lres0[@n]);
1273//          attrib(lres0[@n],"isSB",1);
1274//       }
1275//
1276//       primary[2*@k-1]=lres0[1];
1277//       primary[2*@k]=lres0[2];
1278//       @s=size(primary)/2;
1279//       for(@n=1;@n<=size(lres0)/2-1;@n++)
1280//       {
1281//         primary[2*@s+2*@n-1]=lres0[2*@n+1];
1282//         primary[2*@s+2*@n]=lres0[2*@n+2];
1283//       }
1284//       @k--;
1285//=============================================================
1286     }
1287  }
1288  return(primary);
1289}
1290example
1291{ "EXAMPLE:"; echo = 2;
1292   ring  r = 0,(x,y,z),lp;
1293   poly  p = z2+1;
1294   poly  q = z4+2;
1295   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1296   i=std(i);
1297   list  pr= zero_decomp(i,ideal(0),0);
1298   pr;
1299}
1300///////////////////////////////////////////////////////////////////////////////
1301proc extF(list l,list #)
1302{
1303//zero_dimensional primary decomposition after finite field extension
1304   def R=basering;
1305   int p=char(R);
1306
1307   if((p==0)||(p>13)||(npars(R)>0)){return(l);}
1308
1309   int ex=3;
1310   if(size(#)>0){ex=#[1];}
1311
1312   list peek,peek1;
1313   while(size(l)>0)
1314   {
1315      if(size(l[2])==0)
1316      {
1317         peek[size(peek)+1]=l[1];
1318      }
1319      else
1320      {
1321         peek1[size(peek1)+1]=l[1];
1322         peek1[size(peek1)+1]=l[2];
1323      }
1324      l=delete(l,1);
1325      l=delete(l,1);
1326   }
1327   if(size(peek)==0){return(peek1);}
1328
1329   string gnir="ring RH=("+string(p)+"^"+string(ex)+",a),("+varstr(R)+"),lp;";
1330   execute(gnir);
1331   string mp="minpoly="+string(minpoly)+";";
1332   gnir="ring RL=("+string(p)+",a),("+varstr(R)+"),lp;";
1333   execute(gnir);
1334   execute(mp);
1335   list L=imap(R,peek);
1336   list pr, keep;
1337   int i;
1338   for(i=1;i<=size(L);i++)
1339   {
1340      attrib(L[i],"isSB",1);
1341      pr=zero_decomp(L[i],0,0);
1342      keep=keep+pr;
1343   }
1344   for(i=1;i<=size(keep);i++)
1345   {
1346      keep[i]=simplify(keep[i],1);
1347   }
1348   mp="poly pp="+string(minpoly)+";";
1349
1350   string gnir1="ring RS="+string(p)+",("+varstr(R)+",a),lp;";
1351   execute(gnir1);
1352   execute(mp);
1353   list L=imap(RL,keep);
1354
1355   for(i=1;i<=size(L);i++)
1356   {
1357      L[i]=eliminate(L[i]+ideal(pp),a);
1358   }
1359   i=0;
1360   int j;
1361   while(i<size(L)/2-1)
1362   {
1363      i++;
1364      j=i;
1365      while(j<size(L)/2)
1366      {
1367         j++;
1368         if(idealsEqual(L[2*i-1],L[2*j-1]))
1369         {
1370            L=delete(L,2*j-1);
1371            L=delete(L,2*j-1);
1372            j--;
1373         }
1374      }
1375   }
1376   setring R;
1377   list re=imap(RS,L);
1378   re=re+peek1;
1379
1380   return(extF(re,ex+1));
1381}
1382
1383///////////////////////////////////////////////////////////////////////////////
1384proc zeroSp(ideal i)
1385{
1386//preparation for the separable closure
1387//decomposition into ideals of special type
1388//i.e. the minimal polynomials of every variable mod i are irreducible
1389//returns a list of 2 lists: rr=pe,qe
1390//the ideals in pe[l] are special, their special elements are in qe[l]
1391//pe[l] is a dp-Groebnerbasis
1392//the radical of the intersection of the pe[l] is equal to the radical of i
1393
1394   def R=basering;
1395
1396   //i has to be a reduced groebner basis
1397   ideal F=finduni(i);
1398
1399   int j,k,l,ready;
1400   list fa;
1401   fa[1]=factorize(F[1],1);
1402   poly te,ti;
1403   ideal tj;
1404   //avoid factorization of the same polynomial
1405   for(j=2;j<=size(F);j++)
1406   {
1407      for(k=1;k<=j-1;k++)
1408      {
1409         ti=F[k];
1410         te=subst(ti,var(k),var(j));
1411         if(te==F[j])
1412         {
1413            tj=fa[k];
1414            fa[j]=subst(tj,var(k),var(j));
1415            ready=1;
1416            break;
1417         }
1418      }
1419      if(!ready)
1420      {
1421         fa[j]=factorize(F[j],1);
1422      }
1423      ready=0;
1424   }
1425   execute( "ring P=("+charstr(R)+"),("+varstr(R)+"),(C,dp);");
1426   ideal i=imap(R,i);
1427   if(npars(basering)==0)
1428   {
1429      ideal J=fglm(R,i);
1430   }
1431   else
1432   {
1433      ideal J=groebner(i);
1434   }
1435   list fa=imap(R,fa);
1436   list qe=J;          //collects a dp-Groebnerbasis of the special ideals
1437   list keep=ideal(0); //collects the special elements
1438
1439   list re,em,ke;
1440   ideal K,L;
1441
1442   for(j=1;j<=nvars(basering);j++)
1443   {
1444      for(l=1;l<=size(qe);l++)
1445      {
1446         for(k=1;k<=size(fa[j]);k++)
1447         {
1448            L=std(qe[l],fa[j][k]);
1449            K=keep[l],fa[j][k];
1450            if(deg(L[1])>0)
1451            {
1452               re[size(re)+1]=L;
1453               ke[size(ke)+1]=K;
1454            }
1455         }
1456      }
1457      qe=re;
1458      re=em;
1459      keep=ke;
1460      ke=em;
1461   }
1462
1463   setring R;
1464   list qe=imap(P,keep);
1465   list pe=imap(P,qe);
1466   for(l=1;l<=size(qe);l++)
1467   {
1468      qe[l]=simplify(qe[l],2);
1469    }
1470    list rr=pe,qe;
1471    return(rr);
1472}
1473///////////////////////////////////////////////////////////////////////////////
1474
1475proc zeroSepClos(ideal I,ideal F)
1476{
1477//computes the separable closure of the special ideal I
1478//F is the set of special elements of I
1479//returns the separable closure sc(I) of I and an intvec v
1480//such that sc(I)=preimage(frobenius definde by v)
1481//i.e. var(i)----->var(i)^(p^v[i])
1482
1483   if(homog(I)==1){return(maxideal(1));}
1484
1485   //assume F[i] irreducible in I and depending only on var(i)
1486
1487   def R=basering;
1488   int n=nvars(R);
1489   int p=char(R);
1490   intvec v;
1491   v[n]=0;
1492   int i,k;
1493   list l;
1494
1495   for(i=1;i<=n;i++)
1496   {
1497      l[i]=sep(F[i],i);
1498      F[i]=l[i][1];
1499      if(l[i][2]>k){k=l[i][2];}
1500   }
1501
1502   if(k==0){return(list(I,v));}        //the separable case
1503   ideal m;
1504
1505   for(i=1;i<=n;i++)
1506   {
1507      m[i]=var(i)^(p^l[i][2]);
1508      v[i]=l[i][2];
1509   }
1510   map phi=R,m;
1511   ideal J=preimage(R,phi,I);
1512   return(list(J,v));
1513}
1514///////////////////////////////////////////////////////////////////////////////
1515
1516proc insepDecomp(ideal i)
1517{
1518//decomposes i into special ideals
1519//computes the prime decomposition of the special ideals
1520//and transforms it back to a decomposition of i
1521
1522   def R=basering;
1523   list pr=zeroSp(i);
1524   int l,k;
1525   list re,wo,qr;
1526   ideal m=maxideal(1);
1527   ideal K;
1528   map phi=R,m;
1529   int p=char(R);
1530   intvec op=option(get);
1531
1532   for(l=1;l<=size(pr[1]);l++)
1533   {
1534      wo=zeroSepClos(pr[1][l],pr[2][l]);
1535      for(k=1;k<=nvars(basering);k++)
1536      {
1537         m[k]=var(k)^(p^wo[2][k]);
1538      }
1539      phi=R,m;
1540      qr=decomp(wo[1],2);
1541
1542      option(redSB);
1543      for(k=1;k<=size(qr)/2;k++)
1544      {
1545         K=qr[2*k];
1546         K=phi(K);
1547         K=groebner(K);
1548         re[size(re)+1]=zeroRad(K);
1549      }
1550      option(noredSB);
1551   }
1552   option(set,op);
1553   return(re);
1554}
1555
1556
1557///////////////////////////////////////////////////////////////////////////////
1558
1559static proc clearSB (ideal i,list #)
1560"USAGE:   clearSB(i); i ideal which is SB ordered by monomial ordering
1561RETURN:  ideal = minimal SB
1562NOTE:
1563EXAMPLE: example clearSB; shows an example
1564"
1565{
1566  int k,j;
1567  poly m;
1568  int c=size(i);
1569
1570  if(size(#)==0)
1571  {
1572    for(j=1;j<c;j++)
1573    {
1574      if(deg(i[j])==0)
1575      {
1576        i=ideal(1);
1577        return(i);
1578      }
1579      if(deg(i[j])>0)
1580      {
1581        m=lead(i[j]);
1582        for(k=j+1;k<=c;k++)
1583        {
1584           if(size(lead(i[k])/m)>0)
1585           {
1586             i[k]=0;
1587           }
1588        }
1589      }
1590    }
1591  }
1592  else
1593  {
1594    j=0;
1595    while(j<c-1)
1596    {
1597      j++;
1598      if(deg(i[j])==0)
1599      {
1600        i=ideal(1);
1601        return(i);
1602      }
1603      if(deg(i[j])>0)
1604      {
1605        m=lead(i[j]);
1606        for(k=j+1;k<=c;k++)
1607        {
1608           if(size(lead(i[k])/m)>0)
1609           {
1610             if((leadexp(m)!=leadexp(i[k]))||(#[j]<=#[k]))
1611             {
1612                i[k]=0;
1613             }
1614             else
1615             {
1616                i[j]=0;
1617                break;
1618             }
1619           }
1620        }
1621      }
1622    }
1623  }
1624  return(simplify(i,2));
1625}
1626example
1627{ "EXAMPLE:"; echo = 2;
1628   ring  r = (0,a,b),(x,y,z),dp;
1629   ideal i=ax2+y,a2x+y,bx;
1630   list l=1,2,1;
1631   ideal j=clearSB(i,l);
1632   j;
1633}
1634
1635///////////////////////////////////////////////////////////////////////////////
1636
1637static proc independSet (ideal j)
1638"USAGE:   independentSet(i); i ideal
1639RETURN:  list = new varstring with the independent set at the end,
1640                ordstring with the corresponding block ordering,
1641                the integer where the independent set starts in the varstring
1642NOTE:
1643EXAMPLE: example independentSet; shows an example
1644"
1645{
1646   int n,k,di;
1647   list resu,hilf;
1648   string var1,var2;
1649   list v=indepSet(j,1);
1650
1651   for(n=1;n<=size(v);n++)
1652   {
1653      di=0;
1654      var1="";
1655      var2="";
1656      for(k=1;k<=size(v[n]);k++)
1657      {
1658         if(v[n][k]!=0)
1659         {
1660            di++;
1661            var2=var2+"var("+string(k)+"),";
1662         }
1663         else
1664         {
1665            var1=var1+"var("+string(k)+"),";
1666         }
1667      }
1668      if(di>0)
1669      {
1670         var1=var1+var2;
1671         var1=var1[1..size(var1)-1];
1672         hilf[1]=var1;
1673         hilf[2]="lp";
1674         //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
1675         hilf[3]=di;
1676         resu[n]=hilf;
1677      }
1678      else
1679      {
1680         resu[n]=varstr(basering),ordstr(basering),0;
1681      }
1682   }
1683   return(resu);
1684}
1685example
1686{ "EXAMPLE:"; echo = 2;
1687   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1688   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1689   i=std(i);
1690   list  l=independSet(i);
1691   l;
1692   i=i,g;
1693   l=independSet(i);
1694   l;
1695
1696   ring s=0,(x,y,z),lp;
1697   ideal i=z,yx;
1698   list l=independSet(i);
1699   l;
1700
1701
1702}
1703///////////////////////////////////////////////////////////////////////////////
1704
1705static proc maxIndependSet (ideal j)
1706"USAGE:   maxIndependentSet(i); i ideal
1707RETURN:  list = new varstring with the maximal independent set at the end,
1708                ordstring with the corresponding block ordering,
1709                the integer where the independent set starts in the varstring
1710NOTE:
1711EXAMPLE: example maxIndependentSet; shows an example
1712"
1713{
1714   int n,k,di;
1715   list resu,hilf;
1716   string var1,var2;
1717   list v=indepSet(j,0);
1718
1719   for(n=1;n<=size(v);n++)
1720   {
1721      di=0;
1722      var1="";
1723      var2="";
1724      for(k=1;k<=size(v[n]);k++)
1725      {
1726         if(v[n][k]!=0)
1727         {
1728            di++;
1729            var2=var2+"var("+string(k)+"),";
1730         }
1731         else
1732         {
1733            var1=var1+"var("+string(k)+"),";
1734         }
1735      }
1736      if(di>0)
1737      {
1738         var1=var1+var2;
1739         var1=var1[1..size(var1)-1];
1740         hilf[1]=var1;
1741         hilf[2]="lp";
1742         hilf[3]=di;
1743         resu[n]=hilf;
1744      }
1745      else
1746      {
1747         resu[n]=varstr(basering),ordstr(basering),0;
1748      }
1749   }
1750   return(resu);
1751}
1752example
1753{ "EXAMPLE:"; echo = 2;
1754   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1755   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1756   i=std(i);
1757   list  l=maxIndependSet(i);
1758   l;
1759   i=i,g;
1760   l=maxIndependSet(i);
1761   l;
1762
1763   ring s=0,(x,y,z),lp;
1764   ideal i=z,yx;
1765   list l=maxIndependSet(i);
1766   l;
1767
1768
1769}
1770
1771///////////////////////////////////////////////////////////////////////////////
1772
1773static proc prepareQuotientring (int nnp)
1774"USAGE:   prepareQuotientring(nnp); nnp int
1775RETURN:  string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
1776NOTE:
1777EXAMPLE: example independentSet; shows an example
1778"
1779{
1780  ideal @ih,@jh;
1781  int npar=npars(basering);
1782  int @n;
1783
1784  string quotring= "ring quring = ("+charstr(basering);
1785  for(@n=nnp+1;@n<=nvars(basering);@n++)
1786  {
1787     quotring=quotring+",var("+string(@n)+")";
1788     @ih=@ih+var(@n);
1789  }
1790
1791  quotring=quotring+"),(var(1)";
1792  @jh=@jh+var(1);
1793  for(@n=2;@n<=nnp;@n++)
1794  {
1795    quotring=quotring+",var("+string(@n)+")";
1796    @jh=@jh+var(@n);
1797  }
1798  quotring=quotring+"),(C,lp);";
1799
1800  return(quotring);
1801
1802}
1803example
1804{ "EXAMPLE:"; echo = 2;
1805   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
1806   def @Q=basering;
1807   list l= prepareQuotientring(3);
1808   l;
1809   execute(l[1]);
1810   execute(l[2]);
1811   basering;
1812   phi;
1813   setring @Q;
1814
1815}
1816
1817///////////////////////////////////////////////////////////////////////////////
1818static proc cleanPrimary(list l)
1819{
1820   int i,j;
1821   list lh;
1822   for(i=1;i<=size(l)/2;i++)
1823   {
1824      if(deg(l[2*i-1][1])>0)
1825      {
1826         j++;
1827         lh[j]=l[2*i-1];
1828         j++;
1829         lh[j]=l[2*i];
1830      }
1831   }
1832   return(lh);
1833}
1834///////////////////////////////////////////////////////////////////////////////
1835
1836
1837proc minAssPrimesold(ideal i, list #)
1838"USAGE:   minAssPrimes(i); i ideal
1839         minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
1840RETURN:  list = the minimal associated prime ideals of i
1841EXAMPLE: example minAssPrimes; shows an example
1842"
1843{
1844   def @P=basering;
1845   if(size(i)==0){return(list(ideal(0)));}
1846   list qr=simplifyIdeal(i);
1847   map phi=@P,qr[2];
1848   i=qr[1];
1849
1850   execute ("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
1851             +ordstr(basering)+");");
1852
1853
1854   ideal i=fetch(@P,i);
1855   if(size(#)==0)
1856   {
1857      int @wr;
1858      list tluser,@res;
1859      list primary=decomp(i,2);
1860
1861      @res[1]=primary;
1862
1863      tluser=union(@res);
1864      setring @P;
1865      list @res=imap(gnir,tluser);
1866      return(phi(@res));
1867   }
1868   list @res,empty;
1869   ideal ser;
1870   option(redSB);
1871   list @pr=facstd(i);
1872   if(size(@pr)==1)
1873   {
1874      attrib(@pr[1],"isSB",1);
1875      if((dim(@pr[1])==0)&&(homog(@pr[1])==1))
1876      {
1877         setring @P;
1878         list @res=maxideal(1);
1879         return(phi(@res));
1880      }
1881      if(dim(@pr[1])>1)
1882      {
1883         setring @P;
1884        // kill gnir;
1885         execute ("ring gnir1 = ("+charstr(basering)+"),
1886                              ("+varstr(basering)+"),(C,lp);");
1887         ideal i=fetch(@P,i);
1888         list @pr=facstd(i);
1889        // ideal ser;
1890         setring gnir;
1891         @pr=fetch(gnir1,@pr);
1892         kill gnir1;
1893      }
1894   }
1895    option(noredSB);
1896   int j,k,odim,ndim,count;
1897   attrib(@pr[1],"isSB",1);
1898   if(#[1]==77)
1899   {
1900     odim=dim(@pr[1]);
1901     count=1;
1902     intvec pos;
1903     pos[size(@pr)]=0;
1904     for(j=2;j<=size(@pr);j++)
1905     {
1906        attrib(@pr[j],"isSB",1);
1907        ndim=dim(@pr[j]);
1908        if(ndim>odim)
1909        {
1910           for(k=count;k<=j-1;k++)
1911           {
1912              pos[k]=1;
1913           }
1914           count=j;
1915           odim=ndim;
1916        }
1917        if(ndim<odim)
1918        {
1919           pos[j]=1;
1920        }
1921     }
1922     for(j=1;j<=size(@pr);j++)
1923     {
1924        if(pos[j]!=1)
1925        {
1926            @res[j]=decomp(@pr[j],2);
1927        }
1928        else
1929        {
1930           @res[j]=empty;
1931        }
1932     }
1933   }
1934   else
1935   {
1936     ser=ideal(1);
1937     for(j=1;j<=size(@pr);j++)
1938     {
1939//@pr[j];
1940//pause();
1941        @res[j]=decomp(@pr[j],2);
1942//       @res[j]=decomp(@pr[j],2,@pr[j],ser);
1943//       for(k=1;k<=size(@res[j]);k++)
1944//       {
1945//          ser=intersect(ser,@res[j][k]);
1946//       }
1947     }
1948   }
1949
1950   @res=union(@res);
1951   setring @P;
1952   list @res=imap(gnir,@res);
1953   return(phi(@res));
1954}
1955example
1956{ "EXAMPLE:"; echo = 2;
1957   ring  r = 32003,(x,y,z),lp;
1958   poly  p = z2+1;
1959   poly  q = z4+2;
1960   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1961   list pr= minAssPrimes(i);  pr;
1962
1963   minAssPrimes(i,1);
1964}
1965
1966static proc primT(ideal i)
1967{
1968   //assumes that all generators of i are irreducible
1969   //i is standard basis
1970
1971   attrib(i,"isSB",1);
1972   int j=size(i);
1973   int k;
1974   while(j>0)
1975   {
1976     if(deg(i[j])>1){break;}
1977     j--;
1978   }
1979   if(j==0){return(1);}
1980   if(deg(i[j])==vdim(i)){return(1);}
1981   return(0);
1982}
1983
1984
1985static proc minAssPrimes(ideal i, list #)
1986"USAGE:   minAssPrimes(i); i ideal
1987         minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
1988RETURN:  list = the minimal associated prime ideals of i
1989EXAMPLE: example minAssPrimes; shows an example
1990"
1991{
1992   def P=basering;
1993   if(size(i)==0){return(list(ideal(0)));}
1994   list q=simplifyIdeal(i);
1995   list re=maxideal(1);
1996   int j,a,k;
1997   intvec op=option(get);
1998   map phi=P,q[2];
1999
2000   if(npars(P)==0){option(redSB);}
2001
2002   i=std(q[1]);
2003   if(dim(i)==-1){re=ideal(1);return(re);}
2004   if((dim(i)==0)&&(npars(P)==0))
2005   {
2006      int di=vdim(i);
2007      execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
2008      ideal J=interred(imap(P,i));
2009      attrib(J,"isSB",1);
2010      if(vdim(J)!=di)
2011      {
2012         J=fglm(P,i);
2013      }
2014      list pr=triangMH(J,2);
2015      list qr,re;
2016
2017      for(k=1;k<=size(pr);k++)
2018      {
2019         if(primT(pr[k]))
2020         {
2021            re[size(re)+1]=pr[k];
2022         }
2023         else
2024         {
2025            attrib(pr[k],"isSB",1);
2026            qr=decomp(pr[k],2);
2027            for(j=1;j<=size(qr)/2;j++)
2028            {
2029               re[size(re)+1]=qr[2*j];
2030            }
2031         }
2032      }
2033      setring P;
2034      re=imap(gnir,re);
2035      option(set,op);
2036      return(phi(re));
2037   }
2038
2039   if((size(#)==0)||(dim(i)==0))
2040   {
2041      re[1]=decomp(i,2);
2042      re=union(re);
2043      option(set,op);
2044      return(phi(re));
2045   }
2046
2047   q=facstd(i);
2048
2049   if((size(q)==1)&&(dim(i)>1))
2050   {
2051      execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
2052
2053      list p=facstd(fetch(P,i));
2054      if(size(p)>1)
2055      {
2056         a=1;
2057         setring P;
2058         q=fetch(gnir,p);
2059      }
2060      else
2061      {
2062         setring P;
2063      }
2064      kill gnir;
2065   }
2066
2067   option(set,op);
2068   for(j=1;j<=size(q);j++)
2069   {
2070      if(a==0){attrib(q[j],"isSB",1);}
2071      re[j]=decomp(q[j],2);
2072   }
2073   re=union(re);
2074   return(phi(re));
2075}
2076example
2077{ "EXAMPLE:"; echo = 2;
2078   ring  r = 32003,(x,y,z),lp;
2079   poly  p = z2+1;
2080   poly  q = z4+2;
2081   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
2082   list pr= minAssPrimes(i);  pr;
2083
2084   minAssPrimes(i,1);
2085}
2086
2087static proc union(list li)
2088{
2089  int i,j,k;
2090
2091  def P=basering;
2092
2093  execute("ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);");
2094  list l=fetch(P,li);
2095  list @erg;
2096
2097  for(k=1;k<=size(l);k++)
2098  {
2099     for(j=1;j<=size(l[k])/2;j++)
2100     {
2101        if(deg(l[k][2*j][1])!=0)
2102        {
2103           i++;
2104           @erg[i]=l[k][2*j];
2105        }
2106     }
2107  }
2108
2109  list @wos;
2110  i=0;
2111  ideal i1,i2;
2112  while(i<size(@erg)-1)
2113  {
2114     i++;
2115     k=i+1;
2116     i1=lead(@erg[i]);
2117      attrib(i1,"isSB",1);
2118      attrib(@erg[i],"isSB",1);
2119
2120     while(k<=size(@erg))
2121     {
2122        if(deg(@erg[i][1])==0)
2123        {
2124           break;
2125        }
2126        i2=lead(@erg[k]);
2127        attrib(@erg[k],"isSB",1);
2128        attrib(i2,"isSB",1);
2129
2130        if(size(reduce(i1,i2,1))==0)
2131        {
2132           if(size(reduce(@erg[i],@erg[k],1))==0)
2133           {
2134              @erg[k]=ideal(1);
2135              i2=ideal(1);
2136           }
2137        }
2138        if(size(reduce(i2,i1,1))==0)
2139        {
2140           if(size(reduce(@erg[k],@erg[i],1))==0)
2141           {
2142              break;
2143           }
2144        }
2145        k++;
2146        if(k>size(@erg))
2147        {
2148           @wos[size(@wos)+1]=@erg[i];
2149        }
2150     }
2151  }
2152  if(deg(@erg[size(@erg)][1])!=0)
2153  {
2154     @wos[size(@wos)+1]=@erg[size(@erg)];
2155  }
2156  setring P;
2157  list @ser=fetch(ir,@wos);
2158  return(@ser);
2159}
2160///////////////////////////////////////////////////////////////////////////////
2161proc equidim(ideal i,list #)
2162"USAGE:  equidim(i) or equidim(i,1) ; i ideal
2163RETURN: list of equidimensional ideals a[1],...,a[s] with:
2164        - a[s] the equidimensional locus of i, i.e. the intersection
2165          of the primary ideals of dimension of i
2166        - a[1],...,a[s-1] the lower dimensional equidimensional loci.
2167NOTE:    An embedded component q (primary ideal) of i can be replaced in the
2168         decomposition by a primary ideal q1 with the same radical as q. @*
2169         @code{equidim(i,1)} uses the algorithm of Eisenbud/Huneke/Vasconcelos.
2170
2171EXAMPLE:example equidim; shows an example
2172"
2173{
2174  if(ord_test(basering)!=1)
2175  {
2176      ERROR(
2177      "// Not implemented for this ordering, please change to global ordering."
2178      );
2179  }
2180  intvec op ;
2181  def  P = basering;
2182  list eq;
2183  intvec w;
2184  int n,m;
2185  int g=size(i);
2186  int a=attrib(i,"isSB");
2187  int homo=homog(i);
2188  if(size(#)!=0)
2189  {
2190     m=1;
2191  }
2192
2193  if(((homo==1)||(a==1))&&(find(ordstr(basering),"l")==0)
2194                                &&(find(ordstr(basering),"s")==0))
2195  {
2196     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
2197                              +ordstr(basering)+");");
2198     ideal i=imap(P,i);
2199     ideal j=i;
2200     if(a==1)
2201     {
2202       attrib(j,"isSB",1);
2203     }
2204     else
2205     {
2206       j=groebner(i);
2207     }
2208  }
2209  else
2210  {
2211     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;");
2212     ideal i=imap(P,i);
2213     ideal j=groebner(i);
2214  }
2215  if(homo==1)
2216  {
2217     for(n=1;n<=nvars(basering);n++)
2218     {
2219        w[n]=ord(var(n));
2220     }
2221     intvec hil=hilb(j,1,w);
2222  }
2223
2224  if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1)
2225                  ||(dim(j)==0)||(dim(j)+g==nvars(basering)))
2226  {
2227    setring P;
2228    eq[1]=i;
2229    return(eq);
2230  }
2231
2232  if(m==0)
2233  {
2234     ideal k=equidimMax(j);
2235  }
2236  else
2237  {
2238     ideal k=equidimMaxEHV(j);
2239  }
2240  if(size(reduce(k,j,1))==0)
2241  {
2242    setring P;
2243    eq[1]=i;
2244    kill gnir;
2245    return(eq);
2246  }
2247  op=option(get);
2248  option(returnSB);
2249  j=quotient(j,k);
2250  option(set,op);
2251
2252  list equi=equidim(j);
2253  if(deg(equi[size(equi)][1])<=0)
2254  {
2255      equi[size(equi)]=k;
2256  }
2257  else
2258  {
2259    equi[size(equi)+1]=k;
2260  }
2261  setring P;
2262  eq=imap(gnir,equi);
2263  kill gnir;
2264  return(eq);
2265}
2266example
2267{ "EXAMPLE:"; echo = 2;
2268   ring  r = 32003,(x,y,z),dp;
2269   ideal i = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
2270   equidim(i);
2271}
2272
2273///////////////////////////////////////////////////////////////////////////////
2274proc equidimMax(ideal i)
2275"USAGE:  equidimMax(i); i ideal
2276RETURN:  ideal of equidimensional locus (of maximal dimension) of i.
2277EXAMPLE: example equidimMax; shows an example
2278"
2279{
2280  if(ord_test(basering)!=1)
2281  {
2282      ERROR(
2283      "// Not implemented for this ordering, please change to global ordering."
2284      );
2285  }
2286  def  P = basering;
2287  ideal eq;
2288  intvec w;
2289  int n;
2290  int g=size(i);
2291  int a=attrib(i,"isSB");
2292  int homo=homog(i);
2293
2294  if(((homo==1)||(a==1))&&(find(ordstr(basering),"l")==0)
2295                                &&(find(ordstr(basering),"s")==0))
2296  {
2297     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
2298                              +ordstr(basering)+");");
2299     ideal i=imap(P,i);
2300     ideal j=i;
2301     if(a==1)
2302     {
2303       attrib(j,"isSB",1);
2304     }
2305     else
2306     {
2307       j=groebner(i);
2308     }
2309  }
2310  else
2311  {
2312     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;");
2313     ideal i=imap(P,i);
2314     ideal j=groebner(i);
2315  }
2316  list indep;
2317  ideal equ,equi;
2318  if(homo==1)
2319  {
2320     for(n=1;n<=nvars(basering);n++)
2321     {
2322        w[n]=ord(var(n));
2323     }
2324     intvec hil=hilb(j,1,w);
2325  }
2326  if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1)
2327                  ||(dim(j)==0)||(dim(j)+g==nvars(basering)))
2328  {
2329    setring P;
2330    return(i);
2331  }
2332
2333  indep=maxIndependSet(j);
2334
2335  execute("ring gnir1 = ("+charstr(basering)+"),("+indep[1][1]+"),("
2336                              +indep[1][2]+");");
2337  if(homo==1)
2338  {
2339     ideal j=std(imap(gnir,j),hil,w);
2340  }
2341  else
2342  {
2343     ideal j=groebner(imap(gnir,j));
2344  }
2345  string quotring=prepareQuotientring(nvars(basering)-indep[1][3]);
2346  execute(quotring);
2347  ideal j=imap(gnir1,j);
2348  kill gnir1;
2349  j=clearSB(j);
2350  ideal h;
2351  for(n=1;n<=size(j);n++)
2352  {
2353     h[n]=leadcoef(j[n]);
2354  }
2355  setring gnir;
2356  ideal h=imap(quring,h);
2357  kill quring;
2358
2359  list l=minSat(j,h);
2360
2361  if(deg(l[2])>0)
2362  {
2363    equ=l[1];
2364    attrib(equ,"isSB",1);
2365    j=std(j,l[2]);
2366
2367    if(dim(equ)==dim(j))
2368    {
2369      equi=equidimMax(j);
2370      equ=interred(intersect(equ,equi));
2371    }
2372  }
2373  else
2374  {
2375    equ=i;
2376  }
2377
2378  setring P;
2379  eq=imap(gnir,equ);
2380  kill gnir;
2381  return(eq);
2382}
2383example
2384{ "EXAMPLE:"; echo = 2;
2385   ring  r = 32003,(x,y,z),dp;
2386   ideal i = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
2387   equidimMax(i);
2388}
2389///////////////////////////////////////////////////////////////////////////////
2390static proc islp()
2391{
2392   string s=ordstr(basering);
2393   int n=find(s,"lp");
2394   if(!n){return(0);}
2395   int k=find(s,",");
2396   string t=s[k+1..size(s)];
2397   int l=find(t,",");
2398   t=s[1..k-1];
2399   int m=find(t,",");
2400   if(l+m){return(0);}
2401   return(1);
2402}
2403///////////////////////////////////////////////////////////////////////////////
2404
2405proc algeDeco(ideal i, int w)
2406{
2407//reduces primery decomposition over algebraic extensions to
2408//the other cases
2409   def R=basering;
2410   int n=nvars(R);
2411
2412//---Anfang Provisorium
2413   if(size(i)==2)
2414   {
2415      option(redSB);
2416      ideal J=std(i);
2417      option(noredSB);
2418      if((size(J)==2)&&(deg(J[1])==1))
2419      {
2420         ideal keep;
2421         poly f;
2422         int j;
2423         for(j=1;j<=nvars(basering);j++)
2424         {
2425           f=J[2];
2426           while((f/var(j))*var(j)-f==0)
2427           {
2428             f=f/var(j);
2429             keep=keep,var(j);
2430           }
2431           J[2]=f;
2432         }
2433         ideal K=factorize(J[2],1);
2434         if(deg(K[1])==0){K=0;}
2435         K=K+std(keep);
2436         ideal L;
2437         list resu;
2438         for(j=1;j<=size(K);j++)
2439         {
2440            L=J[1],K[j];
2441            resu[j]=L;
2442         }
2443         return(resu);
2444      }
2445   }
2446//---Ende Provisorium
2447   string mp="poly p="+string(minpoly)+";";
2448   string gnir="ring RH="+string(char(R))+",("+varstr(R)+","+string(par(1))
2449                +"),dp;";
2450   execute(gnir);
2451   execute(mp);
2452   ideal i=imap(R,i);
2453   ideal I=subst(i,var(nvars(basering)),0);
2454   int j;
2455   for(j=1;j<=ncols(i);j++)
2456   {
2457     if(i[j]!=I[j]){break;}
2458   }
2459   if((j>ncols(i))&&(deg(p)==1))
2460   {
2461     setring R;
2462     kill RH;
2463     kill gnir;
2464     string gnir="ring RH="+string(char(R))+",("+varstr(R)+"),dp;";
2465     execute(gnir);
2466     ideal i=imap(R,i);
2467     ideal J;
2468   }
2469   else
2470   {
2471      i=i,p;
2472   }
2473   list pr;
2474
2475   if(w==0)
2476   {
2477      pr=decomp(i);
2478   }
2479   if(w==1)
2480   {
2481      pr=prim_dec(i,1);
2482      pr=reconvList(pr);
2483   }
2484   if(w==2)
2485   {
2486      pr=minAssPrimes(i);
2487   }
2488   if(n<nvars(basering))
2489   {
2490      gnir="ring RS="+string(char(R))+",("+varstr(RH)
2491                +"),(dp("+string(n)+"),lp);";
2492      execute(gnir);
2493      list pr=imap(RH,pr);
2494      ideal K;
2495      for(j=1;j<=size(pr);j++)
2496      {
2497         K=groebner(pr[j]);
2498         K=K[2..size(K)];
2499         pr[j]=K;
2500      }
2501      setring R;
2502      list pr=imap(RS,pr);
2503   }
2504   else
2505   {
2506      setring R;
2507      list pr=imap(RH,pr);
2508   }
2509   list re;
2510   if(w==2)
2511   {
2512      re=pr;
2513   }
2514   else
2515   {
2516      re=convList(pr);
2517   }
2518   return(re);
2519}
2520
2521///////////////////////////////////////////////////////////////////////////////
2522static proc decomp(ideal i,list #)
2523"USAGE:  decomp(i); i ideal  (for primary decomposition)   (resp.
2524         decomp(i,1);        (for the associated primes of dimension of i) )
2525         decomp(i,2);        (for the minimal associated primes) )
2526RETURN:  list = list of primary ideals and their associated primes
2527         (at even positions in the list)
2528         (resp. a list of the minimal associated primes)
2529NOTE:    Algorithm of Gianni/Trager/Zacharias
2530EXAMPLE: example decomp; shows an example
2531"
2532{
2533  intvec op;
2534  def  @P = basering;
2535  list primary,indep,ltras;
2536  intvec @vh,isat,@w;
2537  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi;
2538  ideal peek=i;
2539  ideal ser,tras;
2540  int isS=(attrib(i,"isSB")==1);
2541
2542  if(size(#)>0)
2543  {
2544     if((#[1]==1)||(#[1]==2))
2545     {
2546        @wr=#[1];
2547        if(size(#)>1)
2548        {
2549          seri=1;
2550          peek=#[2];
2551          ser=#[3];
2552        }
2553      }
2554      else
2555      {
2556        seri=1;
2557        peek=#[1];
2558        ser=#[2];
2559      }
2560  }
2561
2562  homo=homog(i);
2563  if(homo==1)
2564  {
2565    if(attrib(i,"isSB")!=1)
2566    {
2567      ltras=mstd(i);
2568      attrib(ltras[1],"isSB",1);
2569    }
2570    else
2571    {
2572      ltras=i,i;
2573      attrib(ltras[1],"isSB",1);
2574    }
2575    tras=ltras[1];
2576    attrib(tras,"isSB",1);
2577    if(dim(tras)==0)
2578    {
2579        primary[1]=ltras[2];
2580        primary[2]=maxideal(1);
2581        if(@wr>0)
2582        {
2583           list l;
2584           l[1]=maxideal(1);
2585           l[2]=maxideal(1);
2586           return(l);
2587        }
2588        return(primary);
2589     }
2590     for(@n=1;@n<=nvars(basering);@n++)
2591     {
2592        @w[@n]=ord(var(@n));
2593     }
2594     intvec @hilb=hilb(tras,1,@w);
2595     intvec keephilb=@hilb;
2596  }
2597
2598  //----------------------------------------------------------------
2599  //i is the zero-ideal
2600  //----------------------------------------------------------------
2601
2602  if(size(i)==0)
2603  {
2604     primary=i,i;
2605     return(primary);
2606  }
2607
2608  //----------------------------------------------------------------
2609  //pass to the lexicographical ordering and compute a standardbasis
2610  //----------------------------------------------------------------
2611
2612  int lp=islp();
2613
2614  execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);");
2615  op=option(get);
2616  option(redSB);
2617
2618  ideal ser=fetch(@P,ser);
2619
2620  if(homo==1)
2621  {
2622     if(!lp)
2623     {
2624       ideal @j=std(fetch(@P,i),@hilb,@w);
2625     }
2626     else
2627     {
2628        ideal @j=fetch(@P,tras);
2629        attrib(@j,"isSB",1);
2630     }
2631  }
2632  else
2633  {
2634      if(lp&&isS)
2635      {
2636        ideal @j=fetch(@P,i);
2637        attrib(@j,"isSB",1);
2638      }
2639      else
2640      {
2641        ideal @j=groebner(fetch(@P,i));
2642      }
2643  }
2644  option(set,op);
2645  if(seri==1)
2646  {
2647    ideal peek=fetch(@P,peek);
2648    attrib(peek,"isSB",1);
2649  }
2650  else
2651  {
2652    ideal peek=@j;
2653  }
2654  if(size(ser)==0)
2655  {
2656    ideal fried;
2657    @n=size(@j);
2658    for(@k=1;@k<=@n;@k++)
2659    {
2660      if(deg(lead(@j[@k]))==1)
2661      {
2662        fried[size(fried)+1]=@j[@k];
2663        @j[@k]=0;
2664      }
2665    }
2666    if(size(fried)==nvars(basering))
2667    {
2668       setring @P;
2669       primary[1]=i;
2670       primary[2]=i;
2671       return(primary);
2672    }
2673    if(size(fried)>0)
2674    {
2675       string newva;
2676       string newma;
2677       for(@k=1;@k<=nvars(basering);@k++)
2678       {
2679          @n1=0;
2680          for(@n=1;@n<=size(fried);@n++)
2681          {
2682             if(leadmonom(fried[@n])==var(@k))
2683             {
2684                @n1=1;
2685                break;
2686             }
2687          }
2688          if(@n1==0)
2689          {
2690            newva=newva+string(var(@k))+",";
2691            newma=newma+string(var(@k))+",";
2692          }
2693          else
2694          {
2695            newma=newma+string(0)+",";
2696          }
2697       }
2698       newva[size(newva)]=")";
2699       newma[size(newma)]=";";
2700       execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;");
2701       execute("map @kappa=gnir,"+newma);
2702       ideal @j= @kappa(@j);
2703       @j=simplify(@j,2);
2704       attrib(@j,"isSB",1);
2705       list pr=decomp(@j);
2706       setring gnir;
2707       list pr=imap(@deirf,pr);
2708       for(@k=1;@k<=size(pr);@k++)
2709       {
2710         @j=pr[@k]+fried;
2711         pr[@k]=@j;
2712       }
2713       setring @P;
2714       return(imap(gnir,pr));
2715    }
2716  }
2717  //----------------------------------------------------------------
2718  //j is the ring
2719  //----------------------------------------------------------------
2720
2721  if (dim(@j)==-1)
2722  {
2723    setring @P;
2724    primary=ideal(1),ideal(1);
2725    return(primary);
2726  }
2727
2728  //----------------------------------------------------------------
2729  //  the case of one variable
2730  //----------------------------------------------------------------
2731
2732  if(nvars(basering)==1)
2733  {
2734
2735     list fac=factor(@j[1]);
2736     list gprimary;
2737     for(@k=1;@k<=size(fac[1]);@k++)
2738     {
2739        if(@wr==0)
2740        {
2741           gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
2742           gprimary[2*@k]=ideal(fac[1][@k]);
2743        }
2744        else
2745        {
2746           gprimary[2*@k-1]=ideal(fac[1][@k]);
2747           gprimary[2*@k]=ideal(fac[1][@k]);
2748        }
2749     }
2750     setring @P;
2751     primary=fetch(gnir,gprimary);
2752
2753     return(primary);
2754  }
2755
2756 //------------------------------------------------------------------
2757 //the zero-dimensional case
2758 //------------------------------------------------------------------
2759  if (dim(@j)==0)
2760  {
2761    op=option(get);
2762    option(redSB);
2763
2764    list gprimary= zero_decomp(@j,ser,@wr);
2765    option(set,op);
2766    setring @P;
2767    primary=fetch(gnir,gprimary);
2768    if(size(ser)>0)
2769    {
2770      primary=cleanPrimary(primary);
2771    }
2772    return(primary);
2773  }
2774
2775  poly @gs,@gh,@p;
2776  string @va,quotring;
2777  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
2778  ideal @h;
2779  int jdim=dim(@j);
2780  list fett;
2781  int lauf,di,newtest;
2782  //------------------------------------------------------------------
2783  //search for a maximal independent set indep,i.e.
2784  //look for subring such that the intersection with the ideal is zero
2785  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
2786  //indep[1] is the new varstring and indep[2] the string for block-ordering
2787  //------------------------------------------------------------------
2788  if(@wr!=1)
2789  {
2790     allindep=independSet(@j);
2791     for(@m=1;@m<=size(allindep);@m++)
2792     {
2793        if(allindep[@m][3]==jdim)
2794        {
2795           di++;
2796           indep[di]=allindep[@m];
2797        }
2798        else
2799        {
2800           lauf++;
2801           restindep[lauf]=allindep[@m];
2802        }
2803     }
2804   }
2805   else
2806   {
2807      indep=maxIndependSet(@j);
2808   }
2809
2810  ideal jkeep=@j;
2811  if(ordstr(@P)[1]=="w")
2812  {
2813     execute("ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");");
2814  }
2815  else
2816  {
2817     execute( "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);");
2818  }
2819
2820  if(homo==1)
2821  {
2822    if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
2823       ||(ordstr(@P)[3]=="w"))
2824    {
2825      ideal jwork=imap(@P,tras);
2826      attrib(jwork,"isSB",1);
2827    }
2828    else
2829    {
2830      ideal jwork=std(imap(gnir,@j),@hilb,@w);
2831    }
2832
2833  }
2834  else
2835  {
2836    ideal jwork=groebner(imap(gnir,@j));
2837  }
2838  list hquprimary;
2839  poly @p,@q;
2840  ideal @h,fac,ser;
2841  di=dim(jwork);
2842  keepdi=di;
2843
2844  setring gnir;
2845  for(@m=1;@m<=size(indep);@m++)
2846  {
2847     isat=0;
2848     @n2=0;
2849     if((indep[@m][1]==varstr(basering))&&(@m==1))
2850     //this is the good case, nothing to do, just to have the same notations
2851     //change the ring
2852     {
2853        execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2854                              +ordstr(basering)+");");
2855        ideal @j=fetch(gnir,@j);
2856        attrib(@j,"isSB",1);
2857        ideal ser=fetch(gnir,ser);
2858
2859     }
2860     else
2861     {
2862        @va=string(maxideal(1));
2863        execute("ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
2864                              +indep[@m][2]+");");
2865        execute("map phi=gnir,"+@va+";");
2866        op=option(get);
2867        option(redSB);
2868        if(homo==1)
2869        {
2870           ideal @j=std(phi(@j),@hilb,@w);
2871        }
2872        else
2873        {
2874          ideal @j=groebner(phi(@j));
2875        }
2876        ideal ser=phi(ser);
2877
2878        option(set,op);
2879     }
2880     if((deg(@j[1])==0)||(dim(@j)<jdim))
2881     {
2882       setring gnir;
2883       break;
2884     }
2885     for (lauf=1;lauf<=size(@j);lauf++)
2886     {
2887         fett[lauf]=size(@j[lauf]);
2888     }
2889     //------------------------------------------------------------------------
2890     //we have now the following situation:
2891     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2892     //to this quotientring, j is their still a standardbasis, the
2893     //leading coefficients of the polynomials  there (polynomials in
2894     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2895     //we need their ggt, gh, because of the following: let
2896     //(j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2897     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2898     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2899
2900     //------------------------------------------------------------------------
2901
2902     //arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and
2903     //map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..rest..]
2904     //------------------------------------------------------------------------
2905
2906     quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
2907
2908     //---------------------------------------------------------------------
2909     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2910     //---------------------------------------------------------------------
2911
2912     execute(quotring);
2913
2914    // @j considered in the quotientring
2915     ideal @j=imap(gnir1,@j);
2916     ideal ser=imap(gnir1,ser);
2917
2918     kill gnir1;
2919
2920     //j is a standardbasis in the quotientring but usually not minimal
2921     //here it becomes minimal
2922
2923     @j=clearSB(@j,fett);
2924     attrib(@j,"isSB",1);
2925
2926     //we need later ggt(h[1],...)=gh for saturation
2927     ideal @h;
2928     if(deg(@j[1])>0)
2929     {
2930        for(@n=1;@n<=size(@j);@n++)
2931        {
2932           @h[@n]=leadcoef(@j[@n]);
2933        }
2934        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2935        op=option(get);
2936        option(redSB);
2937
2938        list uprimary= zero_decomp(@j,ser,@wr);
2939        option(set,op);
2940     }
2941     else
2942     {
2943       list uprimary;
2944       uprimary[1]=ideal(1);
2945       uprimary[2]=ideal(1);
2946     }
2947     //we need the intersection of the ideals in the list quprimary with the
2948     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2949     //but fi polynomials, then the intersection of q with the polynomialring
2950     //is the saturation of the ideal generated by f1,...,fr with respect to
2951     //h which is the lcm of the leading coefficients of the fi considered in
2952     //in the quotientring: this is coded in saturn
2953
2954     list saturn;
2955     ideal hpl;
2956
2957     for(@n=1;@n<=size(uprimary);@n++)
2958     {
2959        uprimary[@n]=interred(uprimary[@n]); // temporary fix
2960        hpl=0;
2961        for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2962        {
2963           hpl=hpl,leadcoef(uprimary[@n][@n1]);
2964        }
2965        saturn[@n]=hpl;
2966     }
2967
2968     //--------------------------------------------------------------------
2969     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2970     //back to the polynomialring
2971     //---------------------------------------------------------------------
2972     setring gnir;
2973
2974     collectprimary=imap(quring,uprimary);
2975     lsau=imap(quring,saturn);
2976     @h=imap(quring,@h);
2977
2978     kill quring;
2979
2980     @n2=size(quprimary);
2981     @n3=@n2;
2982
2983     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2984     {
2985        if(deg(collectprimary[2*@n1][1])>0)
2986        {
2987           @n2++;
2988           quprimary[@n2]=collectprimary[2*@n1-1];
2989           lnew[@n2]=lsau[2*@n1-1];
2990           @n2++;
2991           lnew[@n2]=lsau[2*@n1];
2992           quprimary[@n2]=collectprimary[2*@n1];
2993        }
2994     }
2995
2996     //here the intersection with the polynomialring
2997     //mentioned above is really computed
2998     for(@n=@n3/2+1;@n<=@n2/2;@n++)
2999     {
3000        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
3001        {
3002           quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
3003           quprimary[2*@n]=quprimary[2*@n-1];
3004        }
3005        else
3006        {
3007           if(@wr==0)
3008           {
3009              quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
3010           }
3011           quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
3012        }
3013     }
3014
3015     if(size(@h)>0)
3016     {
3017           //---------------------------------------------------------------
3018           //we change to @Phelp to have the ordering dp for saturation
3019           //---------------------------------------------------------------
3020        setring @Phelp;
3021        @h=imap(gnir,@h);
3022        if(@wr!=1)
3023        {
3024          @q=minSat(jwork,@h)[2];
3025        }
3026        else
3027        {
3028            fac=ideal(0);
3029            for(lauf=1;lauf<=ncols(@h);lauf++)
3030           {
3031              if(deg(@h[lauf])>0)
3032              {
3033                 fac=fac+factorize(@h[lauf],1);
3034              }
3035           }
3036           fac=simplify(fac,4);
3037           @q=1;
3038           for(lauf=1;lauf<=size(fac);lauf++)
3039           {
3040             @q=@q*fac[lauf];
3041           }
3042        }
3043        jwork=std(jwork,@q);
3044        keepdi=dim(jwork);
3045        if(keepdi<di)
3046        {
3047           setring gnir;
3048           @j=imap(@Phelp,jwork);
3049           break;
3050        }
3051        if(homo==1)
3052        {
3053              @hilb=hilb(jwork,1,@w);
3054        }
3055
3056        setring gnir;
3057        @j=imap(@Phelp,jwork);
3058     }
3059  }
3060
3061  if((size(quprimary)==0)&&(@wr==1))
3062  {
3063     @j=ideal(1);
3064     quprimary[1]=ideal(1);
3065     quprimary[2]=ideal(1);
3066  }
3067  if((size(quprimary)==0))
3068  {
3069    keepdi=di-1;
3070  }
3071  //---------------------------------------------------------------
3072  //notice that j=sat(j,gh) intersected with (j,gh^n)
3073  //we finished with sat(j,gh) and have to start with (j,gh^n)
3074  //---------------------------------------------------------------
3075  if((deg(@j[1])!=0)&&(@wr!=1))
3076  {
3077     if(size(quprimary)>0)
3078     {
3079        setring @Phelp;
3080        ser=imap(gnir,ser);
3081        hquprimary=imap(gnir,quprimary);
3082        if(@wr==0)
3083        {
3084           ideal htest=hquprimary[1];
3085           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
3086           {
3087              htest=intersect(htest,hquprimary[2*@n1-1]);
3088           }
3089        }
3090        else
3091        {
3092           ideal htest=hquprimary[2];
3093
3094           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
3095           {
3096              htest=intersect(htest,hquprimary[2*@n1]);
3097           }
3098        }
3099
3100        if(size(ser)>0)
3101        {
3102           ser=intersect(htest,ser);
3103        }
3104        else
3105        {
3106          ser=htest;
3107        }
3108        setring gnir;
3109        ser=imap(@Phelp,ser);
3110     }
3111     if(size(reduce(ser,peek,1))!=0)
3112     {
3113        for(@m=1;@m<=size(restindep);@m++)
3114        {
3115         // if(restindep[@m][3]>=keepdi)
3116         // {
3117           isat=0;
3118           @n2=0;
3119
3120           if(restindep[@m][1]==varstr(basering))
3121           //the good case, nothing to do, just to have the same notations
3122           //change the ring
3123           {
3124              execute("ring gnir1 = ("+charstr(basering)+"),("+
3125                varstr(basering)+"),("+ordstr(basering)+");");
3126              ideal @j=fetch(gnir,jkeep);
3127              attrib(@j,"isSB",1);
3128           }
3129           else
3130           {
3131              @va=string(maxideal(1));
3132              execute("ring gnir1 = ("+charstr(basering)+"),("+
3133                      restindep[@m][1]+"),(" +restindep[@m][2]+");");
3134              execute("map phi=gnir,"+@va+";");
3135              op=option(get);
3136              option(redSB);
3137              if(homo==1)
3138              {
3139                 ideal @j=std(phi(jkeep),keephilb,@w);
3140              }
3141              else
3142              {
3143                ideal @j=groebner(phi(jkeep));
3144              }
3145              ideal ser=phi(ser);
3146              option(set,op);
3147           }
3148
3149           for (lauf=1;lauf<=size(@j);lauf++)
3150           {
3151              fett[lauf]=size(@j[lauf]);
3152           }
3153           //------------------------------------------------------------------
3154           //we have now the following situation:
3155           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may
3156           //pass to this quotientring, j is their still a standardbasis, the
3157           //leading coefficients of the polynomials  there (polynomials in
3158           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
3159           //we need their ggt, gh, because of the following:
3160           //let (j:gh^n)=(j:gh^infinity) then
3161           //j*K(var(nnp+1),..,var(nva))[..the rest..]
3162           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
3163           //on the other hand j=(j,gh^n) intersected with (j:gh^n)
3164
3165           //------------------------------------------------------------------
3166
3167           //the arrangement for the quotientring
3168           // K(var(nnp+1),..,var(nva))[..the rest..]
3169           //and the map phi:K[var(1),...,var(nva)] ---->
3170           //--->K(var(nnpr+1),..,var(nva))[..the rest..]
3171           //------------------------------------------------------------------
3172
3173           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
3174
3175           //------------------------------------------------------------------
3176           //we pass to the quotientring  K(var(nnp+1),..,var(nva))[..rest..]
3177           //------------------------------------------------------------------
3178
3179           execute(quotring);
3180
3181           // @j considered in the quotientring
3182           ideal @j=imap(gnir1,@j);
3183           ideal ser=imap(gnir1,ser);
3184
3185           kill gnir1;
3186
3187           //j is a standardbasis in the quotientring but usually not minimal
3188           //here it becomes minimal
3189           @j=clearSB(@j,fett);
3190           attrib(@j,"isSB",1);
3191
3192           //we need later ggt(h[1],...)=gh for saturation
3193           ideal @h;
3194
3195           for(@n=1;@n<=size(@j);@n++)
3196           {
3197              @h[@n]=leadcoef(@j[@n]);
3198           }
3199           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
3200
3201           op=option(get);
3202           option(redSB);
3203           list uprimary= zero_decomp(@j,ser,@wr);
3204           option(set,op);
3205
3206           //we need the intersection of the ideals in the list quprimary with
3207           //the polynomialring, i.e. let q=(f1,...,fr) in the quotientring
3208           //such an ideal but fi polynomials, then the intersection of q with
3209           //the polynomialring is the saturation of the ideal generated by
3210           //f1,...,fr with respect toh which is the lcm of the leading
3211           //coefficients of the fi considered in the quotientring:
3212           //this is coded in saturn
3213
3214           list saturn;
3215           ideal hpl;
3216
3217           for(@n=1;@n<=size(uprimary);@n++)
3218           {
3219              hpl=0;
3220              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
3221              {
3222                 hpl=hpl,leadcoef(uprimary[@n][@n1]);
3223              }
3224               saturn[@n]=hpl;
3225           }
3226           //------------------------------------------------------------------
3227           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..rest..]
3228           //back to the polynomialring
3229           //------------------------------------------------------------------
3230           setring gnir;
3231
3232           collectprimary=imap(quring,uprimary);
3233           lsau=imap(quring,saturn);
3234           @h=imap(quring,@h);
3235
3236           kill quring;
3237
3238
3239           @n2=size(quprimary);
3240           @n3=@n2;
3241
3242           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
3243           {
3244              if(deg(collectprimary[2*@n1][1])>0)
3245              {
3246                 @n2++;
3247                 quprimary[@n2]=collectprimary[2*@n1-1];
3248                 lnew[@n2]=lsau[2*@n1-1];
3249                 @n2++;
3250                 lnew[@n2]=lsau[2*@n1];
3251                 quprimary[@n2]=collectprimary[2*@n1];
3252              }
3253           }
3254
3255
3256           //here the intersection with the polynomialring
3257           //mentioned above is really computed
3258
3259           for(@n=@n3/2+1;@n<=@n2/2;@n++)
3260           {
3261              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
3262              {
3263                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
3264                 quprimary[2*@n]=quprimary[2*@n-1];
3265              }
3266              else
3267              {
3268                 if(@wr==0)
3269                 {
3270                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
3271                 }
3272                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
3273              }
3274           }
3275           if(@n2>=@n3+2)
3276           {
3277              setring @Phelp;
3278              ser=imap(gnir,ser);
3279              hquprimary=imap(gnir,quprimary);
3280              for(@n=@n3/2+1;@n<=@n2/2;@n++)
3281              {
3282                if(@wr==0)
3283                {
3284                   ser=intersect(ser,hquprimary[2*@n-1]);
3285                }
3286                else
3287                {
3288                   ser=intersect(ser,hquprimary[2*@n]);
3289                }
3290              }
3291              setring gnir;
3292              ser=imap(@Phelp,ser);
3293           }
3294
3295         // }
3296        }
3297        if(size(reduce(ser,peek,1))!=0)
3298        {
3299           if(@wr>0)
3300           {
3301              htprimary=decomp(@j,@wr,peek,ser);
3302           }
3303           else
3304           {
3305              htprimary=decomp(@j,peek,ser);
3306           }
3307           // here we collect now both results primary(sat(j,gh))
3308           // and primary(j,gh^n)
3309           @n=size(quprimary);
3310           for (@k=1;@k<=size(htprimary);@k++)
3311           {
3312              quprimary[@n+@k]=htprimary[@k];
3313           }
3314        }
3315     }
3316
3317   }
3318  //---------------------------------------------------------------------------
3319  //back to the ring we started with
3320  //the final result: primary
3321  //---------------------------------------------------------------------------
3322
3323  setring @P;
3324  primary=imap(gnir,quprimary);
3325  return(primary);
3326}
3327
3328
3329example
3330{ "EXAMPLE:"; echo = 2;
3331   ring  r = 32003,(x,y,z),lp;
3332   poly  p = z2+1;
3333   poly  q = z4+2;
3334   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
3335   list pr= decomp(i);
3336   pr;
3337   testPrimary( pr, i);
3338}
3339
3340///////////////////////////////////////////////////////////////////////////////
3341static proc powerCoeffs(poly f,int e)
3342//computes a polynomial with the same monomials as f but coefficients
3343//the p^e th power of the coefficients of f
3344{
3345   int i;
3346   poly g;
3347   int ex=char(basering)^e;
3348   for(i=1;i<=size(f);i++)
3349   {
3350      g=g+leadcoef(f[i])^ex*leadmonom(f[i]);
3351   }
3352   return(g);
3353}
3354///////////////////////////////////////////////////////////////////////////////
3355
3356proc sep(poly f,int i, list #)
3357"USAGE:  input: a polynomial f depending on the i-th variable and optional
3358         an integer k considering the polynomial f defined over Fp(t1,...,tm)
3359         as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
3360 RETURN: the separabel part of f as polynomial in Fp(t1,...,tm)
3361        and an integer k to indicate that f should be considerd
3362        as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
3363 EXAMPLE: example sep; shows an example
3364{
3365   def R=basering;
3366   int k;
3367   if(size(#)>0){k=#[1];}
3368
3369
3370   poly h=gcd(f,diff(f,var(i)));
3371   if((reduce(f,std(h))!=0)||(reduce(diff(f,var(i)),std(h))!=0))
3372   {
3373      ERROR("FEHLER IN GCD");
3374   }
3375   poly g1=lift(h,f)[1][1];    //  f/h
3376   poly h1;
3377
3378   while(h!=h1)
3379   {
3380      h1=h;
3381      h=gcd(h,diff(h,var(i)));
3382   }
3383
3384   if(deg(h1)==0){return(list(g1,k));} //in characteristic 0 we return here
3385
3386   k++;
3387
3388   ideal ma=maxideal(1);
3389   ma[i]=var(i)^char(R);
3390   map phi=R,ma;
3391   ideal hh=h;    //this is technical because preimage works only for ideals
3392
3393   poly u=preimage(R,phi,hh)[1]; //h=u(x(i)^p)
3394
3395   list g2=sep(u,i,k);           //we consider u(t(1)^(p^-1),...,t(m)^(p^-1))
3396   g1=powerCoeffs(g1,g2[2]-k+1); //to have g1 over the same field as g2[1]
3397
3398   list g3=sep(g1*g2[1],i,g2[2]);
3399   return(g3);
3400}
3401example
3402{ "EXAMPLE:"; echo = 2;
3403   ring R=(5,t,s),(x,y,z),dp;
3404   poly f=(x^25-t*x^5+t)*(x^3+s);
3405   sep(f,1);
3406}
3407
3408///////////////////////////////////////////////////////////////////////////////
3409 proc zeroRad(ideal I,list #)
3410"USAGE:  zeroRad(I) , I a zero-dimensional ideal
3411 RETURN: the radical of I
3412 NOTE:  Algorithm of Kemper
3413 EXAMPLE: example zeroRad; shows an example
3414{
3415   if(homog(I)==1){return(maxideal(1));}
3416   //I needs to be a reduced standard basis
3417   def R=basering;
3418   int m=npars(R);
3419   int n=nvars(R);
3420   int p=char(R);
3421   int d=vdim(I);
3422   int i,k;
3423   list l;
3424   if(((p==0)||(p>d))&&(d==deg(I[1])))
3425   {
3426     intvec e=leadexp(I[1]);
3427     for(i=1;i<=nvars(basering);i++)
3428     {
3429       if(e[i]!=0) break;
3430     }
3431     I[1]=sep(I[1],i)[1];
3432     return(interred(I));
3433   }
3434   intvec op=option(get);
3435
3436   option(redSB);
3437   ideal F=finduni(I);//F[i] generates I intersected with K[var(i)]
3438
3439   option(set,op);
3440   if(size(#)>0){I=#[1];}
3441
3442   for(i=1;i<=n;i++)
3443   {
3444      l[i]=sep(F[i],i);
3445      F[i]=l[i][1];
3446      if(l[i][2]>k){k=l[i][2];}  //computation of the maximal k
3447   }
3448
3449   if((k==0)||(m==0)){return(interred(I+F));}        //the separable case
3450
3451   for(i=1;i<=n;i++)             //consider all polynomials over
3452   {                             //Fp(t(1)^(p^-k),...,t(m)^(p^-k))
3453      F[i]=powerCoeffs(F[i],k-l[i][2]);
3454   }
3455
3456   string cR="ring @R="+string(p)+",("+parstr(R)+","+varstr(R)+"),dp;";
3457   execute(cR);
3458   ideal F=imap(R,F);
3459
3460   string nR="ring @S="+string(p)+",(y(1..m),"+varstr(R)+","+parstr(R)+"),dp;";
3461   execute(nR);
3462
3463   ideal G=fetch(@R,F);    //G[i](t(1)^(p^-k),...,t(m)^(p^-k),x(i))=sep(F[i])
3464
3465   ideal I=imap(R,I);
3466   ideal J=I+G;
3467   poly el=1;
3468   k=p^k;
3469   for(i=1;i<=m;i++)
3470   {
3471      J=J,var(i)^k-var(m+n+i);
3472      el=el*y(i);
3473   }
3474
3475   J=eliminate(J,el);
3476   setring R;
3477   ideal J=imap(@S,J);
3478   return(J);
3479}
3480example
3481{ "EXAMPLE:"; echo = 2;
3482   ring R=(5,t),(x,y),dp;
3483   ideal I=x^5-t,y^5-t;
3484   zeroRad(I);
3485}
3486
3487///////////////////////////////////////////////////////////////////////////////
3488static proc radicalKL (ideal i,ideal ser,list #)
3489{
3490   attrib(i,"isSB",1);   // i needs to be a reduced standard basis
3491   list indep,fett;
3492   intvec @w,@hilb,op;
3493   int @wr,@n,@m,lauf,di;
3494   ideal fac,@h,collectrad,lsau;
3495   poly @q;
3496   string @va,quotring;
3497
3498   def  @P = basering;
3499   int jdim=dim(i);
3500   int  homo=homog(i);
3501   ideal rad=ideal(1);
3502   ideal te=ser;
3503   if(size(#)>0)
3504   {
3505      @wr=#[1];
3506   }
3507   if(homo==1)
3508   {
3509      for(@n=1;@n<=nvars(basering);@n++)
3510      {
3511         @w[@n]=ord(var(@n));
3512      }
3513      @hilb=hilb(i,1,@w);
3514   }
3515
3516
3517  //---------------------------------------------------------------------------
3518  //j is the ring
3519  //---------------------------------------------------------------------------
3520
3521   if (jdim==-1)
3522   {
3523
3524      return(ideal(1));
3525   }
3526
3527  //---------------------------------------------------------------------------
3528  //the zero-dimensional case
3529  //---------------------------------------------------------------------------
3530
3531   if (jdim==0)
3532   {
3533      return(zeroRad(i));
3534   }
3535   //-------------------------------------------------------------------------
3536   //search for a maximal independent set indep,i.e.
3537   //look for subring such that the intersection with the ideal is zero
3538   //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
3539   //indep[1] is the new varstring, indep[2] the string for the block-ordering
3540   //-------------------------------------------------------------------------
3541
3542   indep=maxIndependSet(i);
3543
3544   for(@m=1;@m<=size(indep);@m++)
3545   {
3546      if((indep[@m][1]==varstr(basering))&&(@m==1))
3547      //this is the good case, nothing to do, just to have the same notations
3548      //change the ring
3549      {
3550        execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
3551                              +ordstr(basering)+");");
3552         ideal @j=fetch(@P,i);
3553         attrib(@j,"isSB",1);
3554      }
3555      else
3556      {
3557         @va=string(maxideal(1));
3558         execute("ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
3559                              +indep[@m][2]+");");
3560         execute("map phi=@P,"+@va+";");
3561         if(homo==1)
3562         {
3563            ideal @j=std(phi(i),@hilb,@w);
3564         }
3565         else
3566         {
3567           ideal @j=groebner(phi(i));
3568         }
3569      }
3570      if((deg(@j[1])==0)||(dim(@j)<jdim))
3571      {
3572         setring @P;
3573         break;
3574      }
3575      for (lauf=1;lauf<=size(@j);lauf++)
3576      {
3577         fett[lauf]=size(@j[lauf]);
3578      }
3579     //------------------------------------------------------------------------
3580     //we have now the following situation:
3581     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
3582     //to this quotientring, j is their still a standardbasis, the
3583     //leading coefficients of the polynomials  there (polynomials in
3584     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
3585     //we need their ggt, gh, because of the following:
3586     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
3587     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
3588     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
3589
3590     //------------------------------------------------------------------------
3591
3592     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
3593     //and the map phi:K[var(1),...,var(nva)] ----->
3594     //K(var(nnpr+1),..,var(nva))[..the rest..]
3595     //------------------------------------------------------------------------
3596      quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
3597
3598     //------------------------------------------------------------------------
3599     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
3600     //------------------------------------------------------------------------
3601
3602      execute(quotring);
3603
3604    // @j considered in the quotientring
3605      ideal @j=imap(gnir1,@j);
3606
3607      kill gnir1;
3608
3609     //j is a standardbasis in the quotientring but usually not minimal
3610     //here it becomes minimal
3611
3612      @j=clearSB(@j,fett);
3613
3614     //we need later ggt(h[1],...)=gh for saturation
3615      ideal @h;
3616      if(deg(@j[1])>0)
3617      {
3618         for(@n=1;@n<=size(@j);@n++)
3619         {
3620            @h[@n]=leadcoef(@j[@n]);
3621         }
3622         op=option(get);
3623         option(redSB);
3624         @j=interred(@j);  //to obtain a reduced standardbasis
3625         attrib(@j,"isSB",1);
3626         option(set,op);
3627
3628         ideal zero_rad= zeroRad(@j);
3629      }
3630      else
3631      {
3632         ideal zero_rad=ideal(1);
3633      }
3634
3635     //we need the intersection of the ideals in the list quprimary with the
3636     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
3637     //but fi polynomials, then the intersection of q with the polynomialring
3638     //is the saturation of the ideal generated by f1,...,fr with respect to
3639     //h which is the lcm of the leading coefficients of the fi considered in
3640     //the quotientring: this is coded in saturn
3641
3642      zero_rad=std(zero_rad);
3643
3644      ideal hpl;
3645
3646      for(@n=1;@n<=size(zero_rad);@n++)
3647      {
3648         hpl=hpl,leadcoef(zero_rad[@n]);
3649      }
3650
3651     //------------------------------------------------------------------------
3652     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
3653     //back to the polynomialring
3654     //------------------------------------------------------------------------
3655      setring @P;
3656
3657      collectrad=imap(quring,zero_rad);
3658      lsau=simplify(imap(quring,hpl),2);
3659      @h=imap(quring,@h);
3660
3661      kill quring;
3662
3663
3664     //here the intersection with the polynomialring
3665     //mentioned above is really computed
3666
3667      collectrad=sat2(collectrad,lsau)[1];
3668      if(deg(@h[1])>=0)
3669      {
3670         fac=ideal(0);
3671         for(lauf=1;lauf<=ncols(@h);lauf++)
3672         {
3673            if(deg(@h[lauf])>0)
3674            {
3675                fac=fac+factorize(@h[lauf],1);
3676            }
3677         }
3678         fac=simplify(fac,4);
3679         @q=1;
3680         for(lauf=1;lauf<=size(fac);lauf++)
3681         {
3682            @q=@q*fac[lauf];
3683         }
3684         op=option(get);
3685         option(returnSB);
3686         option(redSB);
3687         i=quotient(i+ideal(@q),rad);
3688         attrib(i,"isSB",1);
3689         option(set,op);
3690
3691      }
3692      if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
3693      {
3694         rad=intersect(rad,collectrad);
3695         te=intersect(te,collectrad);
3696         te=simplify(reduce(te,i,1),2);
3697      }
3698      else
3699      {
3700         if(deg(collectrad[1])>0)
3701         {
3702            rad=collectrad;
3703            te=intersect(te,collectrad);
3704            te=simplify(reduce(te,i,1),2);
3705         }
3706      }
3707
3708      if((dim(i)<jdim)||(size(te)==0))
3709      {
3710         break;
3711      }
3712      if(homo==1)
3713      {
3714         @hilb=hilb(i,1,@w);
3715      }
3716   }
3717   if(((@wr==1)&&(dim(i)<jdim))||(deg(i[1])==0)||(size(te)==0))
3718   {
3719      return(rad);
3720   }
3721   rad=intersect(rad,radicalKL(i,ideal(1),@wr));
3722   return(rad);
3723}
3724
3725///////////////////////////////////////////////////////////////////////////////
3726
3727proc radicalEHV(ideal i)
3728"USAGE:   radicalEHV(i); i ideal.
3729RETURN:  ideal, the radical of i.
3730NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos, which
3731         reduces the computation to the complete intersection case,
3732         by taking, in the general case, a generic linear combination
3733         of the input.
3734         Works only in characteristic 0 or p large.
3735EXAMPLE: example radicalEHV; shows an example
3736"
3737{
3738   if(ord_test(basering)!=1)
3739   {
3740      ERROR(
3741      "// Not implemented for this ordering, please change to global ordering."
3742      );
3743   }
3744   if((char(basering)<100)&&(char(basering)!=0))
3745   {
3746      "WARNING: The characteristic is too small, the result may be wrong";
3747   }
3748   ideal J,I,I0,radI0,L,radI1,I2,radI2;
3749   int l,n;
3750   intvec op=option(get);
3751   matrix M;
3752
3753   option(redSB);
3754   list m=mstd(i);
3755        I=m[2];
3756   option(set,op);
3757
3758   int cod=nvars(basering)-dim(m[1]);
3759   //-------------------complete intersection case:----------------------
3760   if(cod==size(m[2]))
3761   {
3762     J=minor(jacob(I),cod);
3763     return(quotient(I,J));
3764   }
3765   //-----first codim elements of I are a complete intersection:---------
3766   for(l=1;l<=cod;l++)
3767   {
3768      I0[l]=I[l];
3769   }
3770   n=dim(std(I0))+cod-nvars(basering);
3771   //-----last codim elements of I are a complete intersection:----------
3772   if(n!=0)
3773   {
3774      for(l=1;l<=cod;l++)
3775      {
3776         I0[l]=I[size(I)-l+1];
3777      }
3778      n=dim(std(I0))+cod-nvars(basering);
3779   }
3780   //-----taking a generic linear combination of the input:--------------
3781   if(n!=0)
3782   {
3783      M=transpose(sparsetriag(size(m[2]),cod,95,1));
3784      I0=ideal(M*transpose(I));
3785      n=dim(std(I0))+cod-nvars(basering);
3786   }
3787   //-----taking a more generic linear combination of the input:---------
3788   if(n!=0)
3789   {
3790      M=transpose(sparsetriag(size(m[2]),cod,0,100));
3791      I0=ideal(M*transpose(I));
3792      n=dim(std(I0))+cod-nvars(basering);
3793   }
3794   if(n==0)
3795   {
3796      J=minor(jacob(I0),cod);
3797      radI0=quotient(I0,J);
3798      L=quotient(radI0,I);
3799      radI1=quotient(radI0,L);
3800
3801      if(size(reduce(radI1,m[1],1))==0)
3802      {
3803         return(I);
3804      }
3805
3806      I2=sat(I,radI1)[1];
3807
3808      if(deg(I2[1])<=0)
3809      {
3810         return(radI1);
3811      }
3812      return(intersect(radI1,radicalEHV(I2)));
3813   }
3814   //---------------------general case-------------------------------------
3815   return(radical(I));
3816}
3817example
3818{ "EXAMPLE:";  echo = 2;
3819   ring  r = 0,(x,y,z),dp;
3820   poly  p = z2+1;
3821   poly  q = z3+2;
3822   ideal i = p*q^2,y-z2;
3823   ideal pr= radicalEHV(i);
3824   pr;
3825}
3826
3827///////////////////////////////////////////////////////////////////////////////
3828
3829proc Ann(module M)
3830"USAGE:   Ann(M);  M module
3831RETURN:  ideal, the annihilator of coker(M)
3832NOTE:    The output is the ideal of all elements a of the basering R such that
3833         a * R^m is contained in M  (m=number of rows of M).
3834EXAMPLE: example Ann; shows an example
3835"
3836{
3837  M=prune(M);  //to obtain a small embedding
3838  ideal ann=quotient1(M,freemodule(nrows(M)));
3839  return(ann);
3840}
3841example
3842{ "EXAMPLE:"; echo = 2;
3843   ring  r = 0,(x,y,z),lp;
3844   module M = x2-y2,z3;
3845   Ann(M);
3846   M = [1,x2],[y,x];
3847   Ann(M);
3848   qring Q=std(xy-1);
3849   module M=imap(r,M);
3850   Ann(M);
3851}
3852
3853///////////////////////////////////////////////////////////////////////////////
3854
3855//computes the equidimensional part of the ideal i of codimension e
3856static proc int_ass_primary_e(ideal i, int e)
3857{
3858  if(homog(i)!=1)
3859  {
3860     i=std(i);
3861  }
3862  list re=sres(i,0);                   //the resolution
3863  re=minres(re);                       //minimized resolution
3864  ideal ann=AnnExt_R(e,re);
3865  if(nvars(basering)-dim(std(ann))!=e)
3866  {
3867    return(ideal(1));
3868  }
3869  return(ann);
3870}
3871
3872///////////////////////////////////////////////////////////////////////////////
3873
3874//computes the annihilator of Ext^n(R/i,R) with given resolution re
3875//n is not necessarily the number of variables
3876static proc AnnExt_R(int n,list re)
3877{
3878  if(n<nvars(basering))
3879  {
3880     matrix f=transpose(re[n+1]);      //Hom(_,R)
3881     module k=nres(f,2)[2];            //the kernel
3882     matrix g=transpose(re[n]);        //the image of Hom(_,R)
3883
3884     ideal ann=quotient1(g,k);           //the anihilator
3885  }
3886  else
3887  {
3888     ideal ann=Ann(transpose(re[n]));
3889  }
3890  return(ann);
3891}
3892///////////////////////////////////////////////////////////////////////////////
3893
3894static proc analyze(list pr)
3895{
3896   int ii,jj;
3897   for(ii=1;ii<=size(pr)/2;ii++)
3898   {
3899      dim(std(pr[2*ii]));
3900      idealsEqual(pr[2*ii-1],pr[2*ii]);
3901      "===========================";
3902   }
3903
3904   for(ii=size(pr)/2;ii>1;ii--)
3905   {
3906      for(jj=1;jj<ii;jj++)
3907      {
3908         if(size(reduce(pr[2*jj],std(pr[2*ii],1)))==0)
3909         {
3910            "eingebette Komponente";
3911            jj;
3912            ii;
3913         }
3914      }
3915   }
3916}
3917
3918///////////////////////////////////////////////////////////////////////////////
3919//
3920//                  Shimoyama-Yokoyama
3921//
3922///////////////////////////////////////////////////////////////////////////////
3923
3924static proc simplifyIdeal(ideal i)
3925{
3926  def r=basering;
3927
3928  int j,k;
3929  map phi;
3930  poly p;
3931
3932  ideal iwork=i;
3933  ideal imap1=maxideal(1);
3934  ideal imap2=maxideal(1);
3935
3936
3937  for(j=1;j<=nvars(basering);j++)
3938  {
3939    for(k=1;k<=size(i);k++)
3940    {
3941      if(deg(iwork[k]/var(j))==0)
3942      {
3943        p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
3944        imap1[j]=p+2*var(j);
3945        phi=r,imap1;
3946        iwork=phi(iwork);
3947        iwork=subst(iwork,var(j),0);
3948        iwork[k]=var(j);
3949        imap1=maxideal(1);
3950        imap2[j]=-p;
3951        break;
3952      }
3953    }
3954  }
3955  return(iwork,imap2);
3956}
3957
3958
3959///////////////////////////////////////////////////////
3960// ini_mod
3961// input: a polynomial p
3962// output: the initial term of p as needed
3963// in the context of characteristic sets
3964//////////////////////////////////////////////////////
3965
3966static proc ini_mod(poly p)
3967{
3968  if (p==0)
3969  {
3970    return(0);
3971  }
3972  int n; matrix m;
3973  for( n=nvars(basering); n>0; n=n-1)
3974  {
3975    m=coef(p,var(n));
3976    if(m[1,1]!=1)
3977    {
3978      p=m[2,1];
3979      break;
3980    }
3981  }
3982  if(deg(p)==0)
3983  {
3984    p=0;
3985  }
3986  return(p);
3987}
3988///////////////////////////////////////////////////////
3989// min_ass_prim_charsets
3990// input: generators of an ideal PS and an integer cho
3991// If cho=0, the given ordering of the variables is used.
3992// Otherwise, the system tries to find an "optimal ordering",
3993// which in some cases may considerably speed up the algorithm
3994// output: the minimal associated primes of PS
3995// algorithm: via characteriostic sets
3996//////////////////////////////////////////////////////
3997
3998
3999static proc min_ass_prim_charsets (ideal PS, int cho)
4000{
4001  if((cho<0) and (cho>1))
4002    {
4003      "ERROR: <int> must be 0 or 1"
4004      return();
4005    }
4006  if(system("version")>933)
4007  {
4008    option(notWarnSB);
4009  }
4010  if(cho==0)
4011  {
4012    return(min_ass_prim_charsets0(PS));
4013  }
4014  else
4015  {
4016     return(min_ass_prim_charsets1(PS));
4017  }
4018}
4019///////////////////////////////////////////////////////
4020// min_ass_prim_charsets0
4021// input: generators of an ideal PS
4022// output: the minimal associated primes of PS
4023// algorithm: via characteristic sets
4024// the given ordering of the variables is used
4025//////////////////////////////////////////////////////
4026
4027
4028static proc min_ass_prim_charsets0 (ideal PS)
4029{
4030  intvec op;
4031  matrix m=char_series(PS);  // We compute an irreducible
4032                             // characteristic series
4033  int i,j,k;
4034  list PSI;
4035  list PHI;  // the ideals given by the characteristic series
4036  for(i=nrows(m);i>=1; i--)
4037  {
4038     PHI[i]=ideal(m[i,1..ncols(m)]);
4039  }
4040  // We compute the radical of each ideal in PHI
4041  ideal I,JS,II;
4042  int sizeJS, sizeII;
4043  for(i=size(PHI);i>=1; i--)
4044  {
4045     I=0;
4046     for(j=size(PHI[i]);j>0;j=j-1)
4047     {
4048       I=I+ini_mod(PHI[i][j]);
4049     }
4050     JS=std(PHI[i]);
4051     sizeJS=size(JS);
4052     for(j=size(I);j>0;j=j-1)
4053     {
4054       II=0;
4055       sizeII=0;
4056       k=0;
4057       while(k<=sizeII)                  // successive saturation
4058       {
4059         op=option(get);
4060         option(returnSB);
4061         II=quotient(JS,I[j]);
4062         option(set,op);
4063         sizeII=size(II);
4064         if(sizeII==sizeJS)
4065         {
4066            for(k=1;k<=sizeII;k++)
4067            {
4068              if(leadexp(II[k])!=leadexp(JS[k])) break;
4069            }
4070         }
4071         JS=II;
4072         sizeJS=sizeII;
4073       }
4074    }
4075    PSI=insert(PSI,JS);
4076  }
4077  int sizePSI=size(PSI);
4078  // We eliminate redundant ideals
4079  for(i=1;i<sizePSI;i++)
4080  {
4081    for(j=i+1;j<=sizePSI;j++)
4082    {
4083      if(size(PSI[i])!=0)
4084      {
4085        if(size(PSI[j])!=0)
4086        {
4087          if(size(NF(PSI[i],PSI[j],1))==0)
4088          {
4089            PSI[j]=ideal(0);
4090          }
4091          else
4092          {
4093            if(size(NF(PSI[j],PSI[i],1))==0)
4094            {
4095              PSI[i]=ideal(0);
4096            }
4097          }
4098        }
4099      }
4100    }
4101  }
4102  for(i=sizePSI;i>=1;i--)
4103  {
4104    if(size(PSI[i])==0)
4105    {
4106      PSI=delete(PSI,i);
4107    }
4108  }
4109  return (PSI);
4110}
4111
4112///////////////////////////////////////////////////////
4113// min_ass_prim_charsets1
4114// input: generators of an ideal PS
4115// output: the minimal associated primes of PS
4116// algorithm: via characteristic sets
4117// input: generators of an ideal PS and an integer i
4118// The system tries to find an "optimal ordering" of
4119// the variables
4120//////////////////////////////////////////////////////
4121
4122
4123static proc min_ass_prim_charsets1 (ideal PS)
4124{
4125  intvec op;
4126  def oldring=basering;
4127  string n=system("neworder",PS);
4128  execute("ring r=("+charstr(oldring)+"),("+n+"),dp;");
4129  ideal PS=imap(oldring,PS);
4130  matrix m=char_series(PS);  // We compute an irreducible
4131                             // characteristic series
4132  int i,j,k;
4133  ideal I;
4134  list PSI;
4135  list PHI;    // the ideals given by the characteristic series
4136  list ITPHI;  // their initial terms
4137  for(i=nrows(m);i>=1; i--)
4138  {
4139     PHI[i]=ideal(m[i,1..ncols(m)]);
4140     I=0;
4141     for(j=size(PHI[i]);j>0;j=j-1)
4142     {
4143       I=I,ini_mod(PHI[i][j]);
4144     }
4145      I=I[2..ncols(I)];
4146      ITPHI[i]=I;
4147  }
4148  setring oldring;
4149  matrix m=imap(r,m);
4150  list PHI=imap(r,PHI);
4151  list ITPHI=imap(r,ITPHI);
4152  // We compute the radical of each ideal in PHI
4153  ideal I,JS,II;
4154  int sizeJS, sizeII;
4155  for(i=size(PHI);i>=1; i--)
4156  {
4157     I=0;
4158     for(j=size(PHI[i]);j>0;j=j-1)
4159     {
4160       I=I+ITPHI[i][j];
4161     }
4162     JS=std(PHI[i]);
4163     sizeJS=size(JS);
4164     for(j=size(I);j>0;j=j-1)
4165     {
4166       II=0;
4167       sizeII=0;
4168       k=0;
4169       while(k<=sizeII)                  // successive iteration
4170       {
4171         op=option(get);
4172         option(returnSB);
4173         II=quotient(JS,I[j]);
4174         option(set,op);
4175//std
4176//         II=std(II);
4177         sizeII=size(II);
4178         if(sizeII==sizeJS)
4179         {
4180            for(k=1;k<=sizeII;k++)
4181            {
4182              if(leadexp(II[k])!=leadexp(JS[k])) break;
4183            }
4184         }
4185         JS=II;
4186         sizeJS=sizeII;
4187       }
4188    }
4189    PSI=insert(PSI,JS);
4190  }
4191  int sizePSI=size(PSI);
4192  // We eliminate redundant ideals
4193  for(i=1;i<sizePSI;i++)
4194  {
4195    for(j=i+1;j<=sizePSI;j++)
4196    {
4197      if(size(PSI[i])!=0)
4198      {
4199        if(size(PSI[j])!=0)
4200        {
4201          if(size(NF(PSI[i],PSI[j],1))==0)
4202          {
4203            PSI[j]=ideal(0);
4204          }
4205          else
4206          {
4207            if(size(NF(PSI[j],PSI[i],1))==0)
4208            {
4209              PSI[i]=ideal(0);
4210            }
4211          }
4212        }
4213      }
4214    }
4215  }
4216  for(i=sizePSI;i>=1;i--)
4217  {
4218    if(size(PSI[i])==0)
4219    {
4220      PSI=delete(PSI,i);
4221    }
4222  }
4223  return (PSI);
4224}
4225
4226
4227/////////////////////////////////////////////////////
4228// proc prim_dec
4229// input:  generators of an ideal I and an integer choose
4230// If choose=0, min_ass_prim_charsets with the given
4231// ordering of the variables is used.
4232// If choose=1, min_ass_prim_charsets with the "optimized"
4233// ordering of the variables is used.
4234// If choose=2, minAssPrimes from primdec.lib is used
4235// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
4236// output: a primary decomposition of I, i.e., a list
4237// of pairs consisting of a standard basis of a primary component
4238// of I and a standard basis of the corresponding associated prime.
4239// To compute the minimal associated primes of a given ideal
4240// min_ass_prim_l is called, i.e., the minimal associated primes
4241// are computed via characteristic sets.
4242// In the homogeneous case, the performance of the procedure
4243// will be improved if I is already given by a minimal set of
4244// generators. Apply minbase if necessary.
4245//////////////////////////////////////////////////////////
4246
4247
4248static proc prim_dec(ideal I, int choose)
4249{
4250   if((choose<0) or (choose>3))
4251   {
4252     "ERROR: <int> must be 0 or 1 or 2 or 3";
4253      return();
4254   }
4255   if(system("version")>933)
4256   {
4257      option(notWarnSB);
4258   }
4259  ideal H=1; // The intersection of the primary components
4260  list U;    // the leaves of the decomposition tree, i.e.,
4261             // pairs consisting of a primary component of I
4262             // and the corresponding associated prime
4263  list W;    // the non-leaf vertices in the decomposition tree.
4264             // every entry has 6 components:
4265                // 1- the vertex itself , i.e., a standard bais of the
4266                //    given ideal I (type 1), or a standard basis of a
4267                //    pseudo-primary component arising from
4268                //    pseudo-primary decomposition (type 2), or a
4269                //    standard basis of a remaining component arising from
4270                //    pseudo-primary decomposition or extraction (type 3)
4271                // 2- the type of the vertex as indicated above
4272                // 3- the weighted_tree_depth of the vertex
4273                // 4- the tester of the vertex
4274                // 5- a standard basis of the associated prime
4275                //    of a vertex of type 2, or 0 otherwise
4276                // 6- a list of pairs consisting of a standard
4277                //    basis of a minimal associated prime ideal
4278                //    of the father of the vertex and the
4279                //    irreducible factors of the "minimal
4280                //    divisor" of the seperator or extractor
4281                //    corresponding to the prime ideal
4282                //    as computed by the procedure minsat,
4283                //    if the vertex is of type 3, or
4284                //    the empty list otherwise
4285  ideal SI=std(I);
4286  if(SI[1]==1)  // primdecSY(ideal(1))
4287  {
4288    return(list());
4289  }
4290  int ncolsSI=ncols(SI);
4291  int ncolsH=1;
4292  W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
4293  int weighted_tree_depth;
4294  int i,j;
4295  int check;
4296  list V;  // current vertex
4297  list VV; // new vertex
4298  list QQ;
4299  list WI;
4300  ideal Qi,SQ,SRest,fac;
4301  poly tester;
4302
4303  while(1)
4304  {
4305    i=1;
4306    while(1)
4307    {
4308      while(i<=size(W)) // find vertex V of smallest weighted tree-depth
4309      {
4310        if (W[i][3]<=weighted_tree_depth) break;
4311        i++;
4312      }
4313      if (i<=size(W)) break;
4314      i=1;
4315      weighted_tree_depth++;
4316    }
4317    V=W[i];
4318    W=delete(W,i); // delete V from W
4319
4320    // now proceed by type of vertex V
4321
4322    if (V[2]==2)  // extraction needed
4323    {
4324      SQ,SRest,fac=extraction(V[1],V[5]);
4325                        // standard basis of primary component,
4326                        // standard basis of remaining component,
4327                        // irreducible factors of
4328                        // the "minimal divisor" of the extractor
4329                        // as computed by the procedure minsat,
4330      check=0;
4331      for(j=1;j<=ncolsH;j++)
4332      {
4333        if (NF(H[j],SQ,1)!=0) // Q is not redundant
4334        {
4335          check=1;
4336          break;
4337        }
4338      }
4339      if(check==1)             // Q is not redundant
4340      {
4341        QQ=list();
4342        QQ[1]=list(SQ,V[5]);  // primary component, associated prime,
4343                              // i.e., standard bases thereof
4344        U=U+QQ;
4345        H=intersect(H,SQ);
4346        H=std(H);
4347        ncolsH=ncols(H);
4348        check=0;
4349        if(ncolsH==ncolsSI)
4350        {
4351          for(j=1;j<=ncolsSI;j++)
4352          {
4353            if(leadexp(H[j])!=leadexp(SI[j]))
4354            {
4355              check=1;
4356              break;
4357            }
4358          }
4359        }
4360        else
4361        {
4362          check=1;
4363        }
4364        if(check==0) // H==I => U is a primary decomposition
4365        {
4366          return(U);
4367        }
4368      }
4369      if (SRest[1]!=1)        // the remaining component is not
4370                              // the whole ring
4371      {
4372        if (rad_con(V[4],SRest)==0) // the new vertex is not the
4373                                    // root of a redundant subtree
4374        {
4375          VV[1]=SRest;     // remaining component
4376          VV[2]=3;         // pseudoprimdec_special
4377          VV[3]=V[3]+1;    // weighted depth
4378          VV[4]=V[4];      // the tester did not change
4379          VV[5]=ideal(0);
4380          VV[6]=list(list(V[5],fac));
4381          W=insert(W,VV,size(W));
4382        }
4383      }
4384    }
4385    else
4386    {
4387      if (V[2]==3) // pseudo_prim_dec_special is needed
4388      {
4389        QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
4390                         // QQ = quadruples:
4391                         // standard basis of pseudo-primary component,
4392                         // standard basis of corresponding prime,
4393                         // seperator, irreducible factors of
4394                         // the "minimal divisor" of the seperator
4395                         // as computed by the procedure minsat,
4396                         // SRest=standard basis of remaining component
4397      }
4398      else     // V is the root, pseudo_prim_dec is needed
4399      {
4400        QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
4401                         // QQ = quadruples:
4402                         // standard basis of pseudo-primary component,
4403                         // standard basis of corresponding prime,
4404                         // seperator, irreducible factors of
4405                         // the "minimal divisor" of the seperator
4406                         // as computed by the procedure minsat,
4407                         // SRest=standard basis of remaining component
4408
4409      }
4410      //check
4411      for(i=size(QQ);i>=1;i--)
4412      //for(i=1;i<=size(QQ);i++)
4413      {
4414        tester=QQ[i][3]*V[4];
4415        Qi=QQ[i][2];
4416        if(NF(tester,Qi,1)!=0)  // the new vertex is not the
4417                                // root of a redundant subtree
4418        {
4419          VV[1]=QQ[i][1];
4420          VV[2]=2;
4421          VV[3]=V[3]+1;
4422          VV[4]=tester;      // the new tester as computed above
4423          VV[5]=Qi;          // QQ[i][2];
4424          VV[6]=list();
4425          W=insert(W,VV,size(W));
4426        }
4427      }
4428      if (SRest[1]!=1)        // the remaining component is not
4429                              // the whole ring
4430      {
4431        if (rad_con(V[4],SRest)==0) // the vertex is not the root
4432                                    // of a redundant subtree
4433        {
4434          VV[1]=SRest;
4435          VV[2]=3;
4436          VV[3]=V[3]+2;
4437          VV[4]=V[4];      // the tester did not change
4438          VV[5]=ideal(0);
4439          WI=list();
4440          for(i=1;i<=size(QQ);i++)
4441          {
4442            WI=insert(WI,list(QQ[i][2],QQ[i][4]));
4443          }
4444          VV[6]=WI;
4445          W=insert(W,VV,size(W));
4446        }
4447      }
4448    }
4449  }
4450}
4451
4452//////////////////////////////////////////////////////////////////////////
4453// proc pseudo_prim_dec_charsets
4454// input: Generators of an arbitrary ideal I, a standard basis SI of I,
4455// and an integer choo
4456// If choo=0, min_ass_prim_charsets with the given
4457// ordering of the variables is used.
4458// If choo=1, min_ass_prim_charsets with the "optimized"
4459// ordering of the variables is used.
4460// If choo=2, minAssPrimes from primdec.lib is used
4461// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
4462// output: a pseudo primary decomposition of I, i.e., a list
4463// of pseudo primary components together with a standard basis of the
4464// remaining component. Each pseudo primary component is
4465// represented by a quadrupel: A standard basis of the component,
4466// a standard basis of the corresponding associated prime, the
4467// seperator of the component, and the irreducible factors of the
4468// "minimal divisor" of the seperator as computed by the procedure minsat,
4469// calls  proc pseudo_prim_dec_i
4470//////////////////////////////////////////////////////////////////////////
4471
4472
4473static proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
4474{
4475  list L;          // The list of minimal associated primes,
4476                   // each one given by a standard basis
4477  if((choo==0) or (choo==1))
4478    {
4479       L=min_ass_prim_charsets(I,choo);
4480    }
4481  else
4482    {
4483      if(choo==2)
4484      {
4485        L=minAssPrimes(I);
4486      }
4487      else
4488      {
4489        L=minAssPrimes(I,1);
4490      }
4491      for(int i=size(L);i>=1;i=i-1)
4492        {
4493          L[i]=std(L[i]);
4494        }
4495    }
4496  return (pseudo_prim_dec_i(SI,L));
4497}
4498
4499////////////////////////////////////////////////////////////////
4500// proc pseudo_prim_dec_special_charsets
4501// input: a standard basis of an ideal I whose radical is the
4502// intersection of the radicals of ideals generated by one prime ideal
4503// P_i together with one polynomial f_i, the list V6 must be the list of
4504// pairs (standard basis of P_i, irreducible factors of f_i),
4505// and an integer choo
4506// If choo=0, min_ass_prim_charsets with the given
4507// ordering of the variables is used.
4508// If choo=1, min_ass_prim_charsets with the "optimized"
4509// ordering of the variables is used.
4510// If choo=2, minAssPrimes from primdec.lib is used
4511// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
4512// output: a pseudo primary decomposition of I, i.e., a list
4513// of pseudo primary components together with a standard basis of the
4514// remaining component. Each pseudo primary component is
4515// represented by a quadrupel: A standard basis of the component,
4516// a standard basis of the corresponding associated prime, the
4517// seperator of the component, and the irreducible factors of the
4518// "minimal divisor" of the seperator as computed by the procedure minsat,
4519// calls  proc pseudo_prim_dec_i
4520////////////////////////////////////////////////////////////////
4521
4522
4523static proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
4524{
4525  int i,j,l;
4526  list m;
4527  list L;
4528  int sizeL;
4529  ideal P,SP; ideal fac;
4530  int dimSP;
4531  for(l=size(V6);l>=1;l--)   // creates a list of associated primes
4532                             // of I, possibly redundant
4533  {
4534    P=V6[l][1];
4535    fac=V6[l][2];
4536    for(i=ncols(fac);i>=1;i--)
4537    {
4538      SP=P+fac[i];
4539      SP=std(SP);
4540      if(SP[1]!=1)
4541      {
4542        if((choo==0) or (choo==1))
4543        {
4544          m=min_ass_prim_charsets(SP,choo);  // a list of SB
4545        }
4546        else
4547        {
4548          if(choo==2)
4549          {
4550            m=minAssPrimes(SP);
4551          }
4552          else
4553          {
4554            m=minAssPrimes(SP,1);
4555          }
4556          for(j=size(m);j>=1;j=j-1)
4557            {
4558              m[j]=std(m[j]);
4559            }
4560        }
4561        dimSP=dim(SP);
4562        for(j=size(m);j>=1; j--)
4563        {
4564          if(dim(m[j])==dimSP)
4565          {
4566            L=insert(L,m[j],size(L));
4567          }
4568        }
4569      }
4570    }
4571  }
4572  sizeL=size(L);
4573  for(i=1;i<sizeL;i++)     // get rid of redundant primes
4574  {
4575    for(j=i+1;j<=sizeL;j++)
4576    {
4577      if(size(L[i])!=0)
4578      {
4579        if(size(L[j])!=0)
4580        {
4581          if(size(NF(L[i],L[j],1))==0)
4582          {
4583            L[j]=ideal(0);
4584          }
4585          else
4586          {
4587            if(size(NF(L[j],L[i],1))==0)
4588            {
4589              L[i]=ideal(0);
4590            }
4591          }
4592        }
4593      }
4594    }
4595  }
4596  for(i=sizeL;i>=1;i--)
4597  {
4598    if(size(L[i])==0)
4599    {
4600      L=delete(L,i);
4601    }
4602  }
4603  return (pseudo_prim_dec_i(SI,L));
4604}
4605
4606
4607////////////////////////////////////////////////////////////////
4608// proc pseudo_prim_dec_i
4609// input: A standard basis of an arbitrary ideal I, and standard bases
4610// of the minimal associated primes of I
4611// output: a pseudo primary decomposition of I, i.e., a list
4612// of pseudo primary components together with a standard basis of the
4613// remaining component. Each pseudo primary component is
4614// represented by a quadrupel: A standard basis of the component Q_i,
4615// a standard basis of the corresponding associated prime P_i, the
4616// seperator of the component, and the irreducible factors of the
4617// "minimal divisor" of the seperator as computed by the procedure minsat,
4618////////////////////////////////////////////////////////////////
4619
4620
4621static proc pseudo_prim_dec_i (ideal SI, list L)
4622{
4623  list Q;
4624  if (size(L)==1)               // one minimal associated prime only
4625                                // the ideal is already pseudo primary
4626  {
4627    Q=SI,L[1],1;
4628    list QQ;
4629    QQ[1]=Q;
4630    return (QQ,ideal(1));
4631  }
4632
4633  poly f0,f,g;
4634  ideal fac;
4635  int i,j,k,l;
4636  ideal SQi;
4637  ideal I'=SI;
4638  list QP;
4639  int sizeL=size(L);
4640  for(i=1;i<=sizeL;i++)
4641  {
4642    fac=0;
4643    for(j=1;j<=sizeL;j++)           // compute the seperator sep_i
4644                                    // of the i-th component
4645    {
4646      if (i!=j)                       // search g not in L[i], but L[j]
4647      {
4648        for(k=1;k<=ncols(L[j]);k++)
4649        {
4650          if(NF(L[j][k],L[i],1)!=0)
4651          {
4652            break;
4653          }
4654        }
4655        fac=fac+L[j][k];
4656      }
4657    }
4658    // delete superfluous polynomials
4659    fac=simplify(fac,8);
4660    // saturation
4661    SQi,f0,f,fac=minsat_ppd(SI,fac);
4662    I'=I',f;
4663    QP=SQi,L[i],f0,fac;
4664             // the quadrupel:
4665             // a standard basis of Q_i,
4666             // a standard basis of P_i,
4667             // sep_i,
4668             // irreducible factors of
4669             // the "minimal divisor" of the seperator
4670             //  as computed by the procedure minsat,
4671    Q[i]=QP;
4672  }
4673  I'=std(I');
4674  return (Q, I');
4675                   // I' = remaining component
4676}
4677
4678
4679////////////////////////////////////////////////////////////////
4680// proc extraction
4681// input: A standard basis of a pseudo primary ideal I, and a standard
4682// basis of the unique minimal associated prime P of I
4683// output: an extraction of I, i.e., a standard basis of the primary
4684// component Q of I with associated prime P, a standard basis of the
4685// remaining component, and the irreducible factors of the
4686// "minimal divisor" of the extractor as computed by the procedure minsat
4687////////////////////////////////////////////////////////////////
4688
4689
4690static proc extraction (ideal SI, ideal SP)
4691{
4692  list indsets=indepSet(SP,0);
4693  poly f;
4694  if(size(indsets)!=0)      //check, whether dim P != 0
4695  {
4696    intvec v;               // a maximal independent set of variables
4697                            // modulo P
4698    string U;               // the independent variables
4699    string A;               // the dependent variables
4700    int j,k;
4701    int a;                  //  the size of A
4702    int degf;
4703    ideal g;
4704    list polys;
4705    int sizepolys;
4706    list newpoly;
4707    def R=basering;
4708    //intvec hv=hilb(SI,1);
4709    for (k=1;k<=size(indsets);k++)
4710    {
4711      v=indsets[k];
4712      for (j=1;j<=nvars(R);j++)
4713      {
4714        if (v[j]==1)
4715        {
4716          U=U+varstr(j)+",";
4717        }
4718        else
4719        {
4720          A=A+varstr(j)+",";
4721          a++;
4722        }
4723      }
4724
4725      U[size(U)]=")";           // we compute the extractor of I (w.r.t. U)
4726      execute("ring RAU=("+charstr(basering)+"),("+A+U+",(dp("+string(a)+"),dp);");
4727      ideal I=imap(R,SI);
4728      //I=std(I,hv);            // the standard basis in (R[U])[A]
4729      I=std(I);            // the standard basis in (R[U])[A]
4730      A[size(A)]=")";
4731      execute("ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;");
4732      ideal I=imap(RAU,I);
4733      //"std in lokalisierung:"+newline,I;
4734      ideal h;
4735      for(j=ncols(I);j>=1;j--)
4736      {
4737        h[j]=leadcoef(I[j]);  // consider I in (R(U))[A]
4738      }
4739      setring R;
4740      g=imap(Rloc,h);
4741      kill RAU,Rloc;
4742      U="";
4743      A="";
4744      a=0;
4745      f=lcm(g);
4746      newpoly[1]=f;
4747      polys=polys+newpoly;
4748      newpoly=list();
4749    }
4750    f=polys[1];
4751    degf=deg(f);
4752    sizepolys=size(polys);
4753    for (k=2;k<=sizepolys;k++)
4754    {
4755      if (deg(polys[k])<degf)
4756      {
4757        f=polys[k];
4758        degf=deg(f);
4759      }
4760    }
4761  }
4762  else
4763  {
4764    f=1;
4765  }
4766  poly f0,h0; ideal SQ; ideal fac;
4767  if(f!=1)
4768  {
4769    SQ,f0,h0,fac=minsat(SI,f);
4770    return(SQ,std(SI+h0),fac);
4771             // the tripel
4772             // a standard basis of Q,
4773             // a standard basis of remaining component,
4774             // irreducible factors of
4775             // the "minimal divisor" of the extractor
4776             // as computed by the procedure minsat
4777  }
4778  else
4779  {
4780    return(SI,ideal(1),ideal(1));
4781  }
4782}
4783
4784/////////////////////////////////////////////////////
4785// proc minsat
4786// input:  a standard basis of an ideal I and a polynomial p
4787// output: a standard basis IS of the saturation of I w.r. to p,
4788// the maximal squarefree factor f0 of p,
4789// the "minimal divisor" f of f0 such that the saturation of
4790// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
4791// the irreducible factors of f
4792//////////////////////////////////////////////////////////
4793
4794
4795static proc minsat(ideal SI, poly p)
4796{
4797  ideal fac=factorize(p,1);       //the irreducible factors of p
4798  fac=sort(fac)[1];
4799  int i,k;
4800  poly f0=1;
4801  for(i=ncols(fac);i>=1;i--)
4802  {
4803    f0=f0*fac[i];
4804  }
4805  poly f=1;
4806  ideal iold;
4807  list quotM;
4808  quotM[1]=SI;
4809  quotM[2]=fac;
4810  quotM[3]=f0;
4811  // we deal seperately with the first quotient;
4812  // factors, which do not contribute to this one,
4813  // are omitted
4814  iold=quotM[1];
4815  quotM=minquot(quotM);
4816  fac=quotM[2];
4817  if(quotM[3]==1)
4818    {
4819      return(quotM[1],f0,f,fac);
4820    }
4821  while(special_ideals_equal(iold,quotM[1])==0)
4822    {
4823      f=f*quotM[3];
4824      iold=quotM[1];
4825      quotM=minquot(quotM);
4826    }
4827  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
4828}
4829
4830/////////////////////////////////////////////////////
4831// proc minsat_ppd
4832// input:  a standard basis of an ideal I and a polynomial p
4833// output: a standard basis IS of the saturation of I w.r. to p,
4834// the maximal squarefree factor f0 of p,
4835// the "minimal divisor" f of f0 such that the saturation of
4836// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
4837// the irreducible factors of f
4838//////////////////////////////////////////////////////////
4839
4840
4841static proc minsat_ppd(ideal SI, ideal fac)
4842{
4843  fac=sort(fac)[1];
4844  int i,k;
4845  poly f0=1;
4846  for(i=ncols(fac);i>=1;i--)
4847  {
4848    f0=f0*fac[i];
4849  }
4850  poly f=1;
4851  ideal iold;
4852  list quotM;
4853  quotM[1]=SI;
4854  quotM[2]=fac;
4855  quotM[3]=f0;
4856  // we deal seperately with the first quotient;
4857  // factors, which do not contribute to this one,
4858  // are omitted
4859  iold=quotM[1];
4860  quotM=minquot(quotM);
4861  fac=quotM[2];
4862  if(quotM[3]==1)
4863    {
4864      return(quotM[1],f0,f,fac);
4865    }
4866  while(special_ideals_equal(iold,quotM[1])==0)
4867  {
4868    f=f*quotM[3];
4869    iold=quotM[1];
4870    quotM=minquot(quotM);
4871    k++;
4872  }
4873  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
4874}
4875/////////////////////////////////////////////////////////////////
4876// proc minquot
4877// input: a list with 3 components: a standard basis
4878// of an ideal I, a set of irreducible polynomials, and
4879// there product f0
4880// output: a standard basis of the ideal (I:f0), the irreducible
4881// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
4882// the "minimal divisor" f
4883/////////////////////////////////////////////////////////////////
4884
4885static proc minquot(list tsil)
4886{
4887   intvec op;
4888   int i,j,k,action;
4889   ideal verg;
4890   list l;
4891   poly g;
4892   ideal laedi=tsil[1];
4893   ideal fac=tsil[2];
4894   poly f=tsil[3];
4895
4896//std
4897//   ideal star=quotient(laedi,f);
4898//   star=std(star);
4899   op=option(get);
4900   option(returnSB);
4901   ideal star=quotient(laedi,f);
4902   option(set,op);
4903   if(special_ideals_equal(laedi,star)==1)
4904     {
4905       return(laedi,ideal(1),1);
4906     }
4907   action=1;
4908   while(action==1)
4909   {
4910      if(size(fac)==1)
4911      {
4912         action=0;
4913         break;
4914      }
4915      for(i=1;i<=size(fac);i++)
4916      {
4917        g=1;
4918         for(j=1;j<=size(fac);j++)
4919         {
4920            if(i!=j)
4921            {
4922               g=g*fac[j];
4923            }
4924         }
4925//std
4926//         verg=quotient(laedi,g);
4927//         verg=std(verg);
4928         op=option(get);
4929         option(returnSB);
4930         verg=quotient(laedi,g);
4931         option(set,op);
4932         if(special_ideals_equal(verg,star)==1)
4933         {
4934            f=g;
4935            fac[i]=0;
4936            fac=simplify(fac,2);
4937            break;
4938         }
4939         if(i==size(fac))
4940         {
4941            action=0;
4942         }
4943      }
4944   }
4945   l=star,fac,f;
4946   return(l);
4947}
4948/////////////////////////////////////////////////
4949// proc special_ideals_equal
4950// input: standard bases of ideal k1 and k2 such that
4951// k1 is contained in k2, or k2 is contained ink1
4952// output: 1, if k1 equals k2, 0 otherwise
4953//////////////////////////////////////////////////
4954
4955static proc special_ideals_equal( ideal k1, ideal k2)
4956{
4957   int j;
4958   if(size(k1)==size(k2))
4959   {
4960      for(j=1;j<=size(k1);j++)
4961      {
4962         if(leadexp(k1[j])!=leadexp(k2[j]))
4963         {
4964            return(0);
4965         }
4966      }
4967      return(1);
4968   }
4969   return(0);
4970}
4971
4972
4973///////////////////////////////////////////////////////////////////////////////
4974
4975static proc convList(list l)
4976{
4977   int i;
4978   list re,he;
4979   for(i=1;i<=size(l)/2;i++)
4980   {
4981      he=l[2*i-1],l[2*i];
4982      re[i]=he;
4983   }
4984   return(re);
4985}
4986///////////////////////////////////////////////////////////////////////////////
4987
4988static proc reconvList(list l)
4989{
4990   int i;
4991   list re;
4992   for(i=1;i<=size(l);i++)
4993   {
4994      re[2*i-1]=l[i][1];
4995      re[2*i]=l[i][2];
4996   }
4997   return(re);
4998}
4999
5000///////////////////////////////////////////////////////////////////////////////
5001//
5002//     The main procedures
5003//
5004///////////////////////////////////////////////////////////////////////////////
5005
5006proc primdecGTZ(ideal i)
5007"USAGE:   primdecGTZ(i); i ideal
5008RETURN:  a list pr of primary ideals and their associated primes:
5009@format
5010   pr[i][1]   the i-th primary component,
5011   pr[i][2]   the i-th prime component.
5012@end format
5013NOTE:    Algorithm of Gianni/Trager/Zacharias.
5014         Designed for characteristic 0, works also in char k > 0, if it
5015         terminates (may result in an infinite loop in small characteristic!)
5016EXAMPLE: example primdecGTZ; shows an example
5017"
5018{
5019   if(ord_test(basering)!=1)
5020   {
5021      ERROR(
5022      "// Not implemented for this ordering, please change to global ordering."
5023      );
5024   }
5025   if(minpoly!=0)
5026   {
5027      return(algeDeco(i,0));
5028      ERROR(
5029      "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
5030      );
5031   }
5032  return(convList(decomp(i)));
5033}
5034example
5035{ "EXAMPLE:";  echo = 2;
5036   ring  r = 0,(x,y,z),lp;
5037   poly  p = z2+1;
5038   poly  q = z3+2;
5039   ideal i = p*q^2,y-z2;
5040   list pr = primdecGTZ(i);
5041   pr;
5042}
5043///////////////////////////////////////////////////////////////////////////////
5044
5045proc primdecSY(ideal i, list #)
5046"USAGE:   primdecSY(i); i ideal, c int
5047RETURN:  a list pr of primary ideals and their associated primes:
5048@format
5049   pr[i][1]   the i-th primary component,
5050   pr[i][2]   the i-th prime component.
5051@end format
5052NOTE:    Algorithm of Shimoyama/Yokoyama.
5053@format
5054   if c=0,  the given ordering of the variables is used,
5055   if c=1,  minAssChar tries to use an optimal ordering,
5056   if c=2,  minAssGTZ is used,
5057   if c=3,  minAssGTZ and facstd are used.
5058@end format
5059EXAMPLE: example primdecSY; shows an example
5060"
5061{
5062   if(ord_test(basering)!=1)
5063   {
5064      ERROR(
5065      "// Not implemented for this ordering, please change to global ordering."
5066      );
5067   }
5068   i=simplify(i,2);
5069   if ((i[1]==0)||(i[1]==1))
5070   {
5071     list L=list(ideal(i[1]),ideal(i[1]));
5072     return(list(L));
5073   }
5074   if(minpoly!=0)
5075   {
5076      return(algeDeco(i,1));
5077   }
5078   if (size(#)==1)
5079   { return(prim_dec(i,#[1])); }
5080   else
5081   { return(prim_dec(i,1)); }
5082}
5083example
5084{ "EXAMPLE:";  echo = 2;
5085   ring  r = 0,(x,y,z),lp;
5086   poly  p = z2+1;
5087   poly  q = z3+2;
5088   ideal i = p*q^2,y-z2;
5089   list pr = primdecSY(i);
5090   pr;
5091}
5092///////////////////////////////////////////////////////////////////////////////
5093proc minAssGTZ(ideal i,list #)
5094"USAGE:   minAssGTZ(i); i ideal
5095          minAssGTZ(i,1); i ideal  does not use the factorizing Groebner
5096RETURN:  a list, the minimal associated prime ideals of i.
5097NOTE:    Designed for characteristic 0, works also in char k > 0 based
5098         on an algorithm of Yokoyama
5099EXAMPLE: example minAssGTZ; shows an example
5100"
5101{
5102   if(ord_test(basering)!=1)
5103   {
5104      ERROR(
5105      "// Not implemented for this ordering, please change to global ordering."
5106      );
5107   }
5108   if(minpoly!=0)
5109   {
5110      return(algeDeco(i,2));
5111      ERROR(
5112      "// Not implemented for algebraic extensions. Simulate the ring extension by adding the minpoly to the ideal"
5113      );
5114   }
5115  if(size(#)==0)
5116   {
5117      return(minAssPrimes(i,1));
5118   }
5119   return(minAssPrimes(i));
5120}
5121example
5122{ "EXAMPLE:";  echo = 2;
5123   ring  r = 0,(x,y,z),dp;
5124   poly  p = z2+1;
5125   poly  q = z3+2;
5126   ideal i = p*q^2,y-z2;
5127   list pr = minAssGTZ(i);
5128   pr;
5129}
5130
5131///////////////////////////////////////////////////////////////////////////////
5132proc minAssChar(ideal i, list #)
5133"USAGE:   minAssChar(i[,c]); i ideal, c int.
5134RETURN:  list, the minimal associated prime ideals of i.
5135NOTE:    If c=0, the given ordering of the variables is used. @*
5136         Otherwise, the system tries to find an optimal ordering,
5137         which in some cases may considerably speed up the algorithm. @*
5138         Due to a bug in the factorization, the result may be not completely
5139         decomposed in small characteristic.
5140EXAMPLE: example minAssChar; shows an example
5141"
5142{
5143   if(ord_test(basering)!=1)
5144   {
5145      ERROR(
5146      "// Not implemented for this ordering, please change to global ordering."
5147      );
5148   }
5149   if (size(#)==1)
5150   { return(min_ass_prim_charsets(i,#[1])); }
5151   else
5152   { return(min_ass_prim_charsets(i,1)); }
5153}
5154example
5155{ "EXAMPLE:";  echo = 2;
5156   ring  r = 0,(x,y,z),dp;
5157   poly  p = z2+1;
5158   poly  q = z3+2;
5159   ideal i = p*q^2,y-z2;
5160   list pr = minAssChar(i);
5161   pr;
5162}
5163///////////////////////////////////////////////////////////////////////////////
5164proc equiRadical(ideal i)
5165"USAGE:   equiRadical(i); i ideal
5166RETURN:  ideal, intersection of associated primes of i of maximal dimension.
5167NOTE:    A combination of the algorithms of Krick/Logar and Kemper is used.
5168         Works also in positive characteristic (Kempers algorithm).
5169EXAMPLE: example equiRadical; shows an example
5170"
5171{
5172   if(ord_test(basering)!=1)
5173   {
5174      ERROR(
5175      "// Not implemented for this ordering, please change to global ordering."
5176      );
5177   }
5178   return(radical(i,1));
5179}
5180example
5181{ "EXAMPLE:";  echo = 2;
5182   ring  r = 0,(x,y,z),dp;
5183   poly  p = z2+1;
5184   poly  q = z3+2;
5185   ideal i = p*q^2,y-z2;
5186   ideal pr= equiRadical(i);
5187   pr;
5188}
5189
5190///////////////////////////////////////////////////////////////////////////////
5191proc radical(ideal i,list #)
5192"USAGE:   radical(i); i ideal.
5193RETURN:  ideal, the radical of i.
5194NOTE:    A combination of the algorithms of Krick/Logar and Kemper is used.
5195         Works also in positive characteristic (Kempers algorithm).
5196EXAMPLE: example radical; shows an example
5197"
5198{
5199   if(ord_test(basering)!=1)
5200   {
5201      ERROR(
5202      "// Not implemented for this ordering, please change to global ordering."
5203      );
5204   }
5205   def @P=basering;
5206   int j,il;
5207   if(size(#)>0){il=#[1];}
5208   if(size(i)==0){return(ideal(0));}
5209   ideal re=1;
5210   intvec op = option(get);
5211   list qr=simplifyIdeal(i);
5212   ideal isave=i;
5213   map phi=@P,qr[2];
5214
5215   option(redSB);
5216   i=groebner(qr[1]);
5217   option(set,op);
5218   int di=dim(i);
5219
5220   if(di==0)
5221   {
5222     i=zeroRad(i,qr[1]);
5223     return(interred(phi(i)));
5224   }
5225
5226   option(redSB);
5227   list pr=i;
5228   if (!homog(i))
5229   {
5230     pr=facstd(i);
5231   }
5232   option(set,op);
5233   int s=size(pr);
5234
5235   for(j=1;j<=s;j++)
5236   {
5237      attrib(pr[s+1-j],"isSB",1);
5238      if((size(reduce(re,pr[s+1-j],1))!=0)&&((dim(pr[s+1-j])==di)||!il))
5239      {
5240         re=intersect(re,radicalKL(pr[s+1-j],re,il));
5241      }
5242   }
5243   return(interred(phi(re)));
5244}
5245example
5246{ "EXAMPLE:";  echo = 2;
5247   ring  r = 0,(x,y,z),dp;
5248   poly  p = z2+1;
5249   poly  q = z3+2;
5250   ideal i = p*q^2,y-z2;
5251   ideal pr= radical(i);
5252   pr;
5253}
5254///////////////////////////////////////////////////////////////////////////////
5255proc prepareAss(ideal i)
5256"USAGE:   prepareAss(i); i ideal
5257RETURN:  list, the radicals of the maximal dimensional components of i.
5258NOTE:    Uses algorithm of Eisenbud/Huneke/Vasconcelos.
5259EXAMPLE: example prepareAss; shows an example
5260"
5261{
5262  if(ord_test(basering)!=1)
5263  {
5264      ERROR(
5265      "// Not implemented for this ordering, please change to global ordering."
5266      );
5267  }
5268  ideal j=std(i);
5269  int cod=nvars(basering)-dim(j);
5270  int e;
5271  list er;
5272  ideal ann;
5273  if(homog(i)==1)
5274  {
5275     list re=sres(j,0);                   //the resolution
5276     re=minres(re);                       //minimized resolution
5277  }
5278  else
5279  {
5280    list re=mres(i,0);
5281  }
5282  for(e=cod;e<=nvars(basering);e++)
5283  {
5284     ann=AnnExt_R(e,re);
5285
5286     if(nvars(basering)-dim(std(ann))==e)
5287     {
5288        er[size(er)+1]=equiRadical(ann);
5289     }
5290  }
5291  return(er);
5292}
5293example
5294{ "EXAMPLE:";  echo = 2;
5295   ring  r = 0,(x,y,z),dp;
5296   poly  p = z2+1;
5297   poly  q = z3+2;
5298   ideal i = p*q^2,y-z2;
5299   list pr = prepareAss(i);
5300   pr;
5301}
5302///////////////////////////////////////////////////////////////////////////////
5303proc equidimMaxEHV(ideal i)
5304"USAGE:  equidimMaxEHV(i); i ideal
5305RETURN:  ideal, the equidimensional component (of maximal dimension) of i.
5306NOTE:    Uses algorithm of Eisenbud, Huneke and Vasconcelos.
5307EXAMPLE: example equidimMaxEHV; shows an example
5308"
5309{
5310  if(ord_test(basering)!=1)
5311  {
5312      ERROR(
5313      "// Not implemented for this ordering, please change to global ordering."
5314      );
5315  }
5316  ideal j=groebner(i);
5317  int cod=nvars(basering)-dim(j);
5318  int e;
5319  ideal ann;
5320  if(homog(i)==1)
5321  {
5322     list re=sres(j,0);                   //the resolution
5323     re=minres(re);                       //minimized resolution
5324  }
5325  else
5326  {
5327    list re=mres(i,0);
5328  }
5329  ann=AnnExt_R(cod,re);
5330  return(ann);
5331}
5332example
5333{ "EXAMPLE:";  echo = 2;
5334   ring  r = 0,(x,y,z),dp;
5335   ideal i=intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
5336   equidimMaxEHV(i);
5337}
5338
5339proc testPrimary(list pr, ideal k)
5340"USAGE:   testPrimary(pr,k); pr a list, k an ideal.
5341ASSUME:  pr is the result of primdecGTZ(k) or primdecSY(k).
5342RETURN:  int, 1 if the intersection of the ideals in pr is k, 0 if not
5343EXAMPLE: example testPrimary; shows an example
5344"
5345{
5346   int i;
5347   pr=reconvList(pr);
5348   ideal j=pr[1];
5349   for (i=2;i<=size(pr)/2;i++)
5350   {
5351       j=intersect(j,pr[2*i-1]);
5352   }
5353   return(idealsEqual(j,k));
5354}
5355example
5356{ "EXAMPLE:";  echo = 2;
5357   ring  r = 32003,(x,y,z),dp;
5358   poly  p = z2+1;
5359   poly  q = z4+2;
5360   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
5361   list pr = primdecGTZ(i);
5362   testPrimary(pr,i);
5363}
5364
5365///////////////////////////////////////////////////////////////////////////////
5366proc zerodec(ideal I)
5367"USAGE:   zerodec(I); I ideal
5368ASSUME:  I is zero-dimensional, the characteristic of the ground field is 0
5369RETURN:  list of primary ideals, the zero-dimensional decomposition of I
5370NOTE:    The algorithm (of Monico), works well only for a small total number
5371         of solutions (@code{vdim(std(I))} should be < 100) and without
5372         parameters. In practice, it works also in large characteristic p>0
5373         but may fail for small p.
5374@*       If printlevel > 0 (default = 0) additional information is displayed.
5375EXAMPLE: example zerodec; shows an example
5376"
5377{
5378   if(ord_test(basering)!=1)
5379   {
5380      ERROR(
5381      "// Not implemented for this ordering, please change to global ordering."
5382      );
5383   }
5384   def R=basering;
5385   poly q;
5386   int j,time;
5387   matrix m;
5388   list re;
5389   poly va=var(1);
5390   ideal J=groebner(I);
5391   ideal ba=kbase(J);
5392   int d=vdim(J);
5393   dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d));
5394//------ compute matrix of multiplication on R/I with generic element p -----
5395   int e=nvars(basering);
5396   poly p=randomLast(100)[e]+random(-50,50);     //the generic element
5397   matrix n[d][d];
5398   time = timer;
5399   for(j=2;j<=e;j++)
5400   {
5401      va=va*var(j);
5402   }
5403   for(j=1;j<=d;j++)
5404   {
5405      q=reduce(p*ba[j],J);
5406      m=coeffs(q,ba,va);
5407      n[j,1..d]=m[1..d,1];
5408   }
5409   dbprint(printlevel-voice+2,
5410      "// time for computing multiplication matrix (with generic element) : "+
5411      string(timer-time));
5412//---------------- compute characteristic polynomial of matrix --------------
5413   execute("ring P1=("+charstr(R)+"),T,dp;");
5414   matrix n=imap(R,n);
5415   time = timer;
5416   poly charpol=det(n-T*freemodule(d));
5417   dbprint(printlevel-voice+2,"// time for computing char poly: "+
5418           string(timer-time));
5419//------------------- factorize characteristic polynomial -------------------
5420//check first if constant term of charpoly is != 0 (which is true for
5421//sufficiently generic element)
5422   if(charpol[size(charpol)]!=0)
5423   {
5424     time = timer;
5425     list fac=factor(charpol);
5426     testFactor(fac,charpol);
5427     dbprint(printlevel-voice+2,"// time for factorizing char poly: "+
5428             string(timer-time));
5429     int f=size(fac[1]);
5430//--------------------------- the irreducible case --------------------------
5431     if(f==1)
5432     {
5433       setring R;
5434       re=I;
5435       return(re);
5436     }
5437//---------------------------- the reducible case ---------------------------
5438//if f_i are the irreducible factors of charpoly, mult=ri, then <I,g_i^ri>
5439//are the primary components where g_i = f_i(p). However, substituting p in
5440//f_i may result in a huge object although the final result may be small.
5441//Hence it is better to simultaneously reduce with I. For this we need a new
5442//ring.
5443     execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);");
5444     list rfac=imap(P1,fac);
5445     intvec ov=option(get);;
5446     option(redSB);
5447     list re1;
5448     ideal new = T-imap(R,p),imap(R,J);
5449     attrib(new, "isSB",1);    //we know that new is a standard basis
5450     for(j=1;j<=f;j++)
5451     {
5452        re1[j]=reduce(rfac[1][j]^rfac[2][j],new);
5453     }
5454     setring R;
5455     re = imap(P,re1);
5456     for(j=1;j<=f;j++)
5457     {
5458        J=I,re[j];
5459        re[j]=interred(J);
5460     }
5461     option(set,ov);
5462     return(re);
5463  }
5464  else
5465//------------------- choice of generic element failed -------------------
5466  {
5467     dbprint(printlevel-voice+2,"// try new generic element!");
5468     setring R;
5469     return(zerodec(I));
5470  }
5471}
5472example
5473{ "EXAMPLE:";  echo = 2;
5474   ring r  = 0,(x,y),dp;
5475   ideal i = x2-2,y2-2;
5476   list pr = zerodec(i);
5477   pr;
5478}
5479////////////////////////////////////////////////////////////////////////////
5480/*
5481//Beispiele Wenk-Dipl (in ~/Texfiles/Diplom/Wenk/Examples/)
5482//Zeiten: Singular/Singular/Singular -r123456789 -v :wilde13 (PentiumPro200)
5483//Singular for HPUX-9 version 1-3-8  (2000060214)  Jun  2 2000 15:31:26
5484//(wilde13)
5485
5486//1. vdim=20, 3  Komponenten
5487//zerodec-time:2(1)  (matrix:1 charpoly:0 factor:1)
5488//primdecGTZ-time: 1(0)
5489ring rs= 0,(a,b,c),dp;
5490poly f1= a^2*b*c + a*b^2*c + a*b*c^2 + a*b*c + a*b + a*c + b*c;
5491poly f2= a^2*b^2*c + a*b^2*c^2 + a^2*b*c + a*b*c + b*c + a + c;
5492poly f3= a^2*b^2*c^2 + a^2*b^2*c + a*b^2*c + a*b*c + a*c + c + 1;
5493ideal gls=f1,f2,f3;
5494int time=timer;
5495printlevel =1;
5496time=timer; list pr1=zerodec(gls); timer-time;size(pr1);
5497time=timer; list pr =primdecGTZ(gls); timer-time;size(pr);
5498time=timer; ideal ra =radical(gls); timer-time;size(pr);
5499
5500//2.cyclic5  vdim=70, 20 Komponenten
5501//zerodec-time:36(28)  (matrix:1(0) charpoly:18(19) factor:17(9)
5502//primdecGTZ-time: 28(5)
5503//radical : 0
5504ring rs= 0,(a,b,c,d,e),dp;
5505poly f0= a + b + c + d + e + 1;
5506poly f1= a + b + c + d + e;
5507poly f2= a*b + b*c + c*d + a*e + d*e;
5508poly f3= a*b*c + b*c*d + a*b*e + a*d*e + c*d*e;
5509poly f4= a*b*c*d + a*b*c*e + a*b*d*e + a*c*d*e + b*c*d*e;
5510poly f5= a*b*c*d*e - 1;
5511ideal gls= f1,f2,f3,f4,f5;
5512
5513//3. random vdim=40, 1 Komponente
5514//zerodec-time:126(304)  (matrix:1 charpoly:115(298) factor:10(5))
5515//primdecGTZ-time:17 (11)
5516ring rs=0,(x,y,z),dp;
5517poly f1=2*x^2 + 4*x + 3*y^2 + 7*x*z + 9*y*z + 5*z^2;
5518poly f2=7*x^3 + 8*x*y + 12*y^2 + 18*x*z + 3*y^4*z + 10*z^3 + 12;
5519poly f3=3*x^4 + 1*x*y*z + 6*y^3 + 3*x*z^2 + 2*y*z^2 + 4*z^2 + 5;
5520ideal gls=f1,f2,f3;
5521
5522//4. introduction into resultants, sturmfels, vdim=28, 1 Komponente
5523//zerodec-time:4  (matrix:0 charpoly:0 factor:4)
5524//primdecGTZ-time:1
5525ring rs=0,(x,y),dp;
5526poly f1= x4+y4-1;
5527poly f2= x5y2-4x3y3+x2y5-1;
5528ideal gls=f1,f2;
5529
5530//5. 3 quadratic equations with random coeffs, vdim=8, 1 Komponente
5531//zerodec-time:0(0)  (matrix:0 charpoly:0 factor:0)
5532//primdecGTZ-time:1(0)
5533ring rs=0,(x,y,z),dp;
5534poly f1=2*x^2 + 4*x*y + 3*y^2 + 7*x*z + 9*y*z + 5*z^2 + 2;
5535poly f2=7*x^2 + 8*x*y + 12*y^2 + 18*x*z + 3*y*z + 10*z^2 + 12;
5536poly f3=3*x^2 + 1*x*y + 6*y^2 + 3*x*z + 2*y*z + 4*z^2 + 5;
5537ideal gls=f1,f2,f3;
5538
5539//6. 3 polys    vdim=24, 1 Komponente
5540// run("ex14",2);
5541//zerodec-time:5(4)  (matrix:0 charpoly:3(3) factor:2(1))
5542//primdecGTZ-time:4 (2)
5543ring rs=0,(x1,x2,x3,x4),dp;
5544poly f1=16*x1^2 + 3*x2^2 + 5*x3^4 - 1 - 4*x4 + x4^3;
5545poly f2=5*x1^3 + 3*x2^2 + 4*x3^2*x4 + 2*x1*x4 - 1 + x4 + 4*x1 + x2 + x3 + x4;
5546poly f3=-4*x1^2 + x2^2 + x3^2 - 3 + x4^2 + 4*x1 + x2 + x3 + x4;
5547poly f4=-4*x1 + x2 + x3 + x4;
5548ideal gls=f1,f2,f3,f4;
5549
5550//7. ex43, PoSSo, caprasse   vdim=56, 16 Komponenten
5551//zerodec-time:23(15)  (matrix:0 charpoly:16(13) factor:3(2))
5552//primdecGTZ-time:3 (2)
5553ring rs= 0,(y,z,x,t),dp;
5554ideal gls=y^2*z+2*y*x*t-z-2*x,
55554*y^2*z*x-z*x^3+2*y^3*t+4*y*x^2*t-10*y^2+4*z*x+4*x^2-10*y*t+2,
55562*y*z*t+x*t^2-2*z-x,
5557-z^3*x+4*y*z^2*t+4*z*x*t^2+2*y*t^3+4*z^2+4*z*x-10*y*t-10*t^2+2;
5558
5559//8. Arnborg-System, n=6 (II),    vdim=156, 90 Komponenten
5560//zerodec-time (char32003):127(45)(matrix:2(0) charpoly:106(37) factor:16(7))
5561//primdecGTZ-time(char32003) :81 (18)
5562//ring rs= 0,(a,b,c,d,x,f),dp;
5563ring rs= 32003,(a,b,c,d,x,f),dp;
5564ideal gls=a+b+c+d+x+f, ab+bc+cd+dx+xf+af, abc+bcd+cdx+d*xf+axf+abf,
5565abcd+bcdx+cd*xf+ad*xf+abxf+abcf, abcdx+bcd*xf+acd*xf+abd*xf+abcxf+abcdf,
5566abcd*xf-1;
5567
5568//9. ex42, PoSSo, Methan6_1, vdim=27, 2 Komponenten
5569//zerodec-time:610  (matrix:10 charpoly:557 factor:26)
5570//primdecGTZ-time: 118
5571//zerodec-time(char32003):2
5572//primdecGTZ-time(char32003):4
5573//ring rs= 0,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp;
5574ring rs= 32003,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp;
5575ideal gls=64*x2*x7-10*x1*x8+10*x7*x9+11*x7*x10-320000*x1,
5576-32*x2*x7-5*x2*x8-5*x2*x10+160000*x1-5000*x2,
5577-x3*x8+x6*x8+x9*x10+210*x6+1300000,
5578-x4*x8+700000,
5579x10^2-2*x5,
5580-x6*x8+x7*x9-210*x6,
5581-64*x2*x7-10*x7*x9-11*x7*x10+320000*x1-16*x7+7000000,
5582-10*x1*x8-10*x2*x8-10*x3*x8-10*x4*x8-10*x6*x8+10*x2*x10+11*x7*x10
5583    +20000*x2+14*x5,
5584x4*x8-x7*x9-x9*x10-410*x9,
558510*x2*x8+10*x3*x8+10*x6*x8+10*x7*x9-10*x2*x10-11*x7*x10-10*x9*x10
5586    -10*x10^2+1400*x6-4200*x10;
5587
5588//10. ex33, walk-s7, Diplomarbeit von Tim, vdim=114
5589//zerfaellt in unterschiedlich viele Komponenten in versch. Charkteristiken:
5590//char32003:30, char0:3(2xdeg1,1xdeg112!), char181:4(2xdeg1,1xdeg28,1xdeg84)
5591//char 0: zerodec-time:10075 (ca 3h) (matrix:3 charpoly:9367, factor:680
5592//        + 24 sec fuer Normalform (anstatt einsetzen), total [29623k])
5593//        primdecGTZ-time: 214
5594//char 32003:zerodec-time:197(68) (matrix:2(1) charpoly:173(60) factor:15(6))
5595//        primdecGTZ-time:14 (5)
5596//char 181:zerodec-time:(87) (matrix:(1) charpoly:(58) factor:(25))
5597//        primdecGTZ-time:(2)
5598//in char181 stimmen Ergebnisse von zerodec und primdecGTZ ueberein (gecheckt)
5599
5600//ring rs= 0,(a,b,c,d,e,f,g),dp;
5601ring rs= 32003,(a,b,c,d,e,f,g),dp;
5602poly f1= 2gb + 2fc + 2ed + a2 + a;
5603poly f2= 2gc + 2fd + e2 + 2ba + b;
5604poly f3= 2gd + 2fe + 2ca + c + b2;
5605poly f4= 2ge + f2 + 2da + d + 2cb;
5606poly f5= 2fg + 2ea + e + 2db + c2;
5607poly f6= g2 + 2fa + f + 2eb + 2dc;
5608poly f7= 2ga + g + 2fb + 2ec + d2;
5609ideal gls= f1,f2,f3,f4,f5,f6,f7;
5610
5611~/Singular/Singular/Singular -r123456789 -v
5612LIB"./primdec.lib";
5613timer=1;
5614int time=timer;
5615printlevel =1;
5616option(prot,mem);
5617time=timer; list pr1=zerodec(gls); timer-time;
5618
5619time=timer; list pr =primdecGTZ(gls); timer-time;
5620time=timer; list pr =primdecSY(gls); timer-time;
5621time=timer; ideal ra =radical(gls); timer-time;size(pr);
5622LIB"all.lib";
5623
5624ring R=0,(a,b,c,d,e,f),dp;
5625ideal I=cyclic(6);
5626minAssGTZ(I);
5627
5628
5629ring S=(2,a,b),(x,y),lp;
5630ideal I=x8-b,y4+a;
5631minAssGTZ(I);
5632
5633ring S1=2,(x,y,a,b),lp;
5634ideal I=x8-b,y4+a;
5635minAssGTZ(I);
5636
5637
5638ring S2=(2,z),(x,y),dp;
5639minpoly=z2+z+1;
5640ideal I=y3+y+1,x4+x+1;
5641primdecGTZ(I);
5642minAssGTZ(I);
5643
5644ring S3=2,(x,y,z),dp;
5645ideal I=y3+y+1,x4+x+1,z2+z+1;
5646primdecGTZ(I);
5647minAssGTZ(I);
5648
5649
5650ring R1=2,(x,y,z),lp;
5651ideal I=y6+y5+y3+y2+1,x4+x+1,z2+z+1;
5652primdecGTZ(I);
5653minAssGTZ(I);
5654
5655
5656ring R2=(2,z),(x,y),lp;
5657minpoly=z3+z+1;
5658ideal I=y2+y+(z2+z+1),x4+x+1;
5659primdecGTZ(I);
5660minAssGTZ(I);
5661
5662*/
Note: See TracBrowser for help on using the repository browser.