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

spielwiese
Last change on this file since d6db1f2 was d6db1f2, checked in by Hans Schönemann <hannes@…>, 27 years ago
* hannes: changed libs, removed old libs git-svn-id: file:///usr/local/Singular/svn/trunk@240 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 53.7 KB
Line 
1// $Id: primdec.lib,v 1.1 1997-05-05 12:03:01 Singular Exp $
2///////////////////////////////////////////////////////
3// primdec.lib
4// algorithms for primary decomposition based on
5// the ideas of Gianni,Trager,Zacharias
6// written by Gerhard Pfister
7//////////////////////////////////////////////////////
8
9LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I)
10
11  minAssPrimes (ideal I, list choose)
12  // minimal associated primes
13  // The list choose must be either emty (minAssPrimes(I)) or 1
14  // (minAssPrimes(I,1))
15  // In the second case the factorizing Buchberger Algorithm is used
16  // which in most cases may considerably speed up the algorithm
17   
18  primdec (ideal I)
19 
20  // Computes a complete primary decomposition via
21
22  radical(ideal I)
23  //computes the radical of the ideal I
24
25LIB "random.lib";
26///////////////////////////////////////////////////////////////////////////////
27
28proc sat1 (ideal id, poly p)
29USAGE:   sat1(id,j);  id ideal, j polynomial
30RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
31NOTE:    result is a std basis in the basering
32EXAMPLE: example sat; shows an example
33{
34   int @k;
35   ideal inew=std(id);
36   ideal iold;
37   option(returnSB);
38   while(specialIdealsEqual(iold,inew)==0 )
39   {
40      iold=inew;
41      inew=quotient(iold,p);
42      @k++;
43   }
44   @k--;
45   option(noreturnSB);
46   list L =inew,p^@k;
47   return (L);
48}
49
50///////////////////////////////////////////////////////////////////////////////
51
52proc sat2 (ideal id, ideal h)
53USAGE:   sat2(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
56EXAMPLE: example sat2; shows an example
57{
58   int @k,@i;
59   def @P= basering;
60   if(ordstr(basering)[1,2]!="dp")
61   {
62      execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),dp;";
63      ideal inew=std(imap(@P,id));
64      ideal  @h=imap(@P,h);
65   }
66   else
67   {
68      ideal @h=h;
69      ideal inew=std(id);
70   }
71   ideal fac;
72
73   for(@i=1;@i<=ncols(@h);@i++)
74   {
75     if(deg(@h[@i])>0)
76     {
77        fac=fac+factorize(@h[@i],1);
78     }
79   }
80   fac=simplify(fac,4);
81   poly @f=1;
82   if(deg(fac[1])>0)
83   {
84      ideal iold;     
85
86      for(@i=1;@i<=size(fac);@i++)
87      {
88        @f=@f*fac[@i];
89      }
90      option(returnSB);
91      while(specialIdealsEqual(iold,inew)==0 )
92      {
93         iold=inew;
94         if(deg(iold[size(iold)])!=1)
95         {
96            inew=quotient(iold,@f);
97         }
98         else
99         {
100            inew=iold;
101         }
102         @k++;
103      }
104      option(noreturnSB);
105      @k--;
106   }
107
108   if(ordstr(@P)[1,2]!="dp")
109   {
110      setring @P;
111      ideal inew=std(imap(@Phelp,inew));
112      poly @f=imap(@Phelp,@f);
113   }
114   list L =inew,@f^@k;
115   return (L);
116}
117
118///////////////////////////////////////////////////////////////////////////////
119
120proc minSat(ideal inew, ideal h)
121{
122   int i,k;
123   poly f=1;
124   ideal iold,fac;
125   list quotM,l;
126
127   for(i=1;i<=ncols(h);i++)
128   {
129      if(deg(h[i])>0)
130      {
131         fac=fac+factorize(h[i],1);
132      }
133   }
134   fac=simplify(fac,4);
135   if(size(fac)==0)
136   {
137      l=inew,1;
138      return(l); 
139   }
140   fac=sort(fac)[1];
141   for(i=1;i<=size(fac);i++)
142   {
143      f=f*fac[i];
144   }
145   quotM[1]=inew;   
146   quotM[2]=fac;
147   quotM[3]=f;
148   f=1;
149   option(returnSB);   
150   while(specialIdealsEqual(iold,quotM[1])==0)
151   {
152      if(k>0)
153      {
154         f=f*quotM[3];
155      }
156      iold=quotM[1];
157      quotM=quotMin(quotM);
158      k++;
159   }
160   option(noreturnSB);
161   l=quotM[1],f;
162   return(l);
163}
164
165proc quotMin(list tsil)
166{
167   int i,j,k,action;
168   ideal verg;
169   list l;
170   poly g;
171
172   ideal laedi=tsil[1];
173   ideal fac=tsil[2];
174   poly f=tsil[3];
175
176   ideal star=quotient(laedi,f);
177
178   action=1;
179 
180   while(action==1)
181   {
182      if(size(fac)==1)
183      {
184         action=0;
185         break;
186      }
187      for(i=1;i<=size(fac);i++)
188      {
189        g=1;
190         for(j=1;j<=size(fac);j++)
191         {
192            if(i!=j)
193            {
194               g=g*fac[j];
195            }   
196         }
197         verg=quotient(laedi,g);
198         if(specialIdealsEqual(verg,star)==1)
199         {
200            f=g;
201            fac[i]=0;
202            fac=simplify(fac,2);
203            break;
204         }
205         if(i==size(fac))
206         {
207            action=0;
208         }
209      }
210   }
211   l=star,fac,f;
212   return(l); 
213}
214
215
216////////////////////////////////////////////////////////////////////////////////
217proc testFactor(list act,poly p)
218{
219   int i;
220   poly q=act[1][1]^act[2][1];
221   for(i=2;i<=size(act[1]);i++)
222   {
223      q=q*act[1][i]^act[2][i];
224   }
225   q=1/leadcoef(q)*q;
226   p=1/leadcoef(p)*p;
227   if(p-q!=0)
228   {
229      "ERROR IN FACTOR";
230      act;
231      p;
232      q;
233      pause;
234   }
235}
236////////////////////////////////////////////////////////////////////////////////
237
238proc factor(poly p)
239USAGE:   factor(p) p poly
240RETURN:  list=;
241NOTE:   
242EXAMPLE: example factor; shows an example
243{
244 
245  ideal @i;
246  list @l;
247  intvec @v,@w;
248  int @j,@k,@n;
249
250  if(deg(p)<=1)
251  {
252     @i=ideal(p);
253     @v=1;
254     @l[1]=@i;
255     @l[2]=@v;
256     return(@l);
257  }
258  if (size(p)==1)
259  {
260     @w=leadexp(p);
261     for(@j=1;@j<=nvars(basering);@j++)
262     {
263        if(@w[@j]!=0)
264        {
265           @k++;
266           @v[@k]=@w[@j];
267           @i=@i+ideal(var(@j));
268        }
269     }
270     @l[1]=@i;
271     @l[2]=@v;
272     return(@l);
273  }
274  @l=factorize(p,2);
275  if(npars(basering)>0)
276  {
277     for(@j=1;@j<=size(@l[1]);@j++)
278     {
279        if(deg(@l[1][@j])==0)
280        {
281           @n++;
282        }
283     }
284     if(@n>0)
285     {
286        if(@n==size(@l[1]))
287        {
288           @l[1]=ideal(1);
289           @v=1;
290           @l[2]=@v;
291        }
292        else
293        {
294           @k=0;
295           int pleh;
296           for(@j=1;@j<=size(@l[1]);@j++)
297           {
298              if(deg(@l[1][@j])!=0)
299              {
300                 @k++;
301                 @i=@i+ideal(@l[1][@j]);
302                 if(size(@i)==pleh)
303                 {
304"factorization error";
305@l;
306                    @k--;
307                    @v[@k]=@v[@k]+@l[2][@j];
308                 }
309                 else
310                 {
311                    pleh++;
312                    @v[@k]=@l[2][@j];
313                 }
314              }
315           }
316           @l[1]=@i;
317           @l[2]=@v;
318        }
319     }
320  }
321  return(@l);
322}
323example
324{ "EXAMPLE:"; echo = 2;
325   ring  r = 0,(x,y,z),lp;
326   poly  p = (x+y)^2*(y-z)^3;
327   list  l = factor(p);
328   l;
329   ring r1 =(0,b,d,f,g),(a,c,e),lp;
330   poly p  =(1*d)*e^2+(1*d*f^2*g);
331   list  l = factor(p);
332   l;
333   ring r2 =(0,b,f,g),(a,c,e,d),lp;
334   poly p  =(1*d)*e^2+(1*d*f^2*g);
335   list  l = factor(p);
336   l;
337
338}
339
340
341
342////////////////////////////////////////////////////////////////////////////////
343
344proc idealsEqual( ideal k, ideal j)
345{
346   return(stdIdealsEqual(std(k),std(j)));
347}
348
349proc specialIdealsEqual( ideal k1, ideal k2)
350{
351   int j;
352
353   if(size(k1)==size(k2))
354   {
355      for(j=1;j<=size(k1);j++)
356      {
357         if(leadexp(k1[j])!=leadexp(k2[j]))
358         {
359            return(0);
360         }
361      }
362      return(1);
363   }
364   return(0);
365}
366
367proc stdIdealsEqual( ideal k1, ideal k2)
368{
369   int j;
370
371   if(size(k1)==size(k2))
372   {
373      for(j=1;j<=size(k1);j++)
374      {
375         if(leadexp(k1[j])!=leadexp(k2[j]))
376         {
377            return(0);
378         }
379      }
380      attrib(k2,"isSB",1);
381      if(size(reduce(k1,k2))==0)
382      {
383         return(1);
384      }
385   }
386   return(0);
387}
388
389
390////////////////////////////////////////////////////////////////////////////////
391
392proc testPrimary(list pr, ideal k)
393USAGE:   testPrimary(pr,k) pr list, k ideal;
394RETURN:  int = 1, if the intersection of the ideals in pr is k, 0 if not
395NOTE:     
396EEXAMPLE: example testPrimary ; shows an example
397{
398   int i;
399   ideal j=pr[1];
400   for (i=2;i<=size(pr)/2;i++)
401   {
402       j=intersect(j,pr[2*i-1]);
403   }
404   return(idealsEqual(j,k));
405}
406example
407{ "EXAMPLE:"; echo = 2;
408   ring  s = 0,(x,y,z),lp;
409   ideal i=x3-x2-x+1,xy-y;
410   ideal i1=x-1;
411   ideal i2=x-1;
412   ideal i3=y,x2-2x+1;
413   ideal i4=y,x-1;
414   ideal i5=y,x+1;
415   ideal i6=y,x+1;
416   list pr=i1,i2,i3,i4,i5,i6;
417   testPrimary(pr,i);
418   pr[5]=y+1,x+1;
419   testPrimary(pr,i);
420}
421
422////////////////////////////////////////////////////////////////////////////////
423proc printPrimary( list l, list #)
424USAGE:   printPrimary(l) l list;
425RETURN:  nothing
426NOTE:   
427EXAMPLE: example printPrimary; shows an example
428{
429   if(size(#)>0)
430   {
431      "                                            ";
432      "  The primary decomposition of the ideal    ";
433      #[1];
434      "                                            ";
435      "  is:                                       ";
436      "                                            ";
437   }
438   int k;
439   for (k=1;k<=size(l)/2;k=k+1)
440   {
441      "                                            ";
442      "primary ideal:                              ";
443      l[2*k-1];
444      "                                            ";
445      "associated prime ideal                      ";
446      l[2*k];
447      "                                            ";
448   }
449}
450example
451{ "EXAMPLE:"; echo = 2;
452}
453
454////////////////////////////////////////////////////////////////////////////////
455
456
457proc randomLast(int b)
458USAGE:   randomLast
459RETURN:  ideal = maxideal(1) but the last variable exchanged by
460         a sum of it with a linear random combination of the other
461         variables
462NOTE:   
463EXAMPLE: example randomLast; shows an example
464{
465
466  ideal i=maxideal(1);
467  int k=size(i);
468  i[k]=0;
469  i=randomid(i,size(i),b);
470  ideal ires=maxideal(1);
471  ires[k]=i[1]+var(k);
472  return(ires);
473}
474example
475{ "EXAMPLE:"; echo = 2;
476   ring  r = 0,(x,y,z),lp;
477   ideal i = randomLast(10);
478   i;
479}
480
481
482////////////////////////////////////////////////////////////////////////////////
483
484
485proc primaryTest (ideal i, poly p)
486USAGE:   primaryTest(i,p); i ideal p poly
487RETURN:  ideal = radical of i, if i is primary in general position,
488         zerodimensional and the radical of i intersected with K[z]
489         is (p), z the smallest variable with respect to the lexico-
490         graphical ordering, and 0 else
491NOTE:    It is necessary that i is a standardbasis with respect to
492         the lexicographical ordering and the first element in i is
493         a power of p.
494EXAMPLE: example primaryTest; shows an example
495{
496   int m=1;
497   int n=nvars(basering);
498   int e;
499   poly t;
500   ideal h;
501
502   //the first generator of the prim ideal for the result
503   ideal prm=p;
504   attrib(prm,"isSB",1);
505
506   while (n>1)
507   {
508      n=n-1;
509      m=m+1;
510
511      //search for i[m] which has a power of var(n) as leading term
512      if (n==1)
513      {
514         m=size(i);
515      }
516      else
517      {
518         while (lead(i[m])/var(n-1)==0)
519        {
520            m=m+1;
521         }
522         m=m-1;
523      }
524      //check whether i[m] =(c*var(n)+h)^e modulo prm for some
525      //h in K[var(n+1),...,var(nvars(basering))], c in K
526      //if not (0) is returned, else var(n)+h is added to prm
527     
528      e=deg(lead(i[m]));
529      t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
530      i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
531
532      if (reduce(i[m]-t^e,prm) !=0)
533      {
534        return(ideal(0));
535      }
536      h=interred(t);
537      t=h[1];
538
539      prm = prm,t;
540      attrib(prm,"isSB",1);
541   }
542   return(prm);
543}
544example
545{ "EXAMPLE:"; echo=2;
546   ring  r = (0,a,b),(x,y,z),lp;
547   poly  p = z^2+1;
548   ideal i = p^2,(a*y-z^3)^3,(b*x-yz+z4)^4;
549   ideal pr= primaryTest(i,p);
550   pr;
551   i = p^2,(y-z3)^3,(x-yz+z4)^4+1;
552   pr= primaryTest(i,p);
553   pr;
554   ring s=(0,e),(d,c,b,a,y,x,g,f),lp;
555   ideal i=f,g,x4,y,a,b3,c,d;
556   poly p=f;
557   ideal pr= primaryTest(i,p);
558   pr;
559}
560
561
562///////////////////////////////////////////////////////////////////////////////
563proc splitPrimary(list l,ideal ser,int @wr,list sact)
564{
565   int i,j,k,s,r,w;
566   list keepresult,act,keepprime;
567   poly @f;
568   int sl=size(l);
569
570   for(i=1;i<=sl/2;i++)
571   {
572      if(sact[2][i]>1)
573      {
574         keepprime[i]=l[2*i-1]+ideal(sact[1][i]);
575      }
576      else
577      {
578         keepprime[i]=l[2*i-1];
579      }
580   }
581   i=0;
582   while(i<size(l)/2)
583   {
584      i++;
585      if((size(ser)>0)&&(size(reduce(ser,l[2*i-1]))==0))
586      {
587         l[2*i-1]=ideal(1);
588         l[2*i]=ideal(1);
589         continue;
590      }
591
592
593      if(size(l[2*i])==0)
594      {
595         if(homog(l[2*i-1])==1)
596         {
597            l[2*i]=maxideal(1);
598            continue;
599         }
600         j=0;
601         if(i<=sl/2)
602         {
603            j=1;
604         }
605         while(j<size(l[2*i-1]))
606         {
607            j++;
608            act=factor(l[2*i-1][j]);
609            r=size(act[1]);
610            attrib(l[2*i-1],"isSB",1);
611            if((r==1)&&(vdim(l[2*i-1])==deg(l[2*i-1][j])))
612            {
613              l[2*i]=std(l[2*i-1],act[1][1]);
614              break;
615            }
616            if((r==1)&&(act[2][1]>1))
617            {
618               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
619               if(homog(keepprime[i])==1)
620               {
621                  l[2*i]=maxideal(1);
622                  break;
623               }
624            }
625            if(gcdTest(act[1])==1)
626            {
627               for(k=2;k<=r;k++)
628               {
629                  keepprime[size(l)/2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
630               }
631               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
632               for(k=1;k<=r;k++)
633               {
634                  if(@wr==0)
635                  {
636                     keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]);
637                  }
638                  else
639                  {
640                     keepresult[k]=std(l[2*i-1],act[1][k]);
641                  }
642               }
643               l[2*i-1]=keepresult[1];
644               if(vdim(keepresult[1])==deg(act[1][1]))
645               {
646                  l[2*i]=keepresult[1];
647               }
648               if((homog(keepresult[1])==1)||(homog(keepprime[i])==1))
649               {
650                  l[2*i]=maxideal(1);
651               }
652               s=size(l)-2;
653               for(k=2;k<=r;k++)
654               {
655                  l[s+2*k-1]=keepresult[k];
656                  keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k]));
657                  if(vdim(keepresult[k])==deg(act[1][k]))
658                  {
659                     l[s+2*k]=keepresult[k];
660                  }
661                  else
662                  {
663                     l[s+2*k]=ideal(0);
664                  }
665                  if((homog(keepresult[k])==1)||(homog(keepprime[s/2+k])==1))
666                  {
667                    l[s+2*k]=maxideal(1);
668                  }
669               }
670               i--;
671               break;             
672            }
673            if(r>=2)
674            {
675               s=size(l);
676               @f=act[1][1];
677               act=sat1(l[2*i-1],act[1][1]);
678               if(deg(act[1][1])>0)
679               {
680                  l[s+1]=std(l[2*i-1],act[2]);
681                  if(homog(l[s+1])==1)
682                  {
683                     l[s+2]=maxideal(1);
684                  }
685                  else
686                  {
687                     l[s+2]=ideal(0);
688                  }
689                  keepprime[s/2+1]=interred(keepprime[i]+ideal(@f));
690                  if(homog(keepprime[s/2+1])==1)
691                  {
692                     l[s+2]=maxideal(1);
693                  }
694                  keepprime[i]=act[1];               
695                  l[2*i-1]=act[1];
696                  attrib(l[2*i-1],"isSB",1);
697                  if(homog(l[2*i-1])==1)
698                  {
699                     l[2*i]=maxideal(1);
700                  }
701                 
702                  i--;
703                  break;
704               }
705            }
706         }
707      }
708   }
709   if(sl==size(l))
710   {
711      return(l);
712   }
713   for(i=1;i<=size(l)/2;i++)
714   {
715      if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
716      {
717         keepprime[i]=std(keepprime[i]); 
718         if(homog(keepprime[i])==1)
719         {     
720             l[2*i]=maxideal(1);
721         }
722         else
723         {
724            act=zero_decomp(keepprime[i],ideal(0),@wr,1);
725            if(size(act)==2)
726            {
727               l[2*i]=act[2];
728            }
729         }
730      }
731   }
732   return(l);
733}
734example
735{ "EXAMPLE:"; echo=2;
736   LIB "primdec.lib";
737   ring  r = 32003,(x,y,z),lp;
738   ideal i1=x*(x+1),yz,(z+1)*(z-1);
739   ideal i2=xy,yz,(x-2)*(x+3);
740   list l=i1,ideal(0),i2,ideal(0),i2,ideal(1);
741   list l1=splitPrimary(l,ideal(0),0);
742   l1;
743}
744
745////////////////////////////////////////////////////////////////////////////////
746
747proc zero_decomp (ideal j,ideal ser,int @wr,list #)
748USAGE:   zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
749         (@wr=0 for primary decomposition, @wr=1 for computaion of associated
750         primes)
751RETURN:  list = list of primary ideals and their radicals (at even positions
752         in the list) if the input is zero-dimensional and a standardbases
753         with respect to lex-ordering
754         If ser!=(0) and ser is contained in j or if j is not zero-dimen-
755         sional then ideal(1),ideal(1) is returned
756NOTE:    Algorithm of Gianni, Traeger, Zacharias
757EXAMPLE: example zero_decomp; shows an example
758{
759  def   @P = basering;
760  int nva = nvars(basering);
761  int @k,@s,@n,@k1;
762  list primary,lres,act,@lh,@wh;
763  map phi,psi;
764  ideal jmap,helpprim,@qh,@qht;
765  intvec @vh,@hilb;
766  string @ri;
767  poly @f;
768 
769  if (dim(j)>0)
770  {
771     primary[1]=ideal(1);
772     primary[2]=ideal(1);
773     return(primary);
774  }
775  j=interred(j);
776  attrib(j,"isSB",1);
777  if(vdim(j)==deg(j[1]))
778  {
779     if((size(ser)>0)&&(size(reduce(ser,j))==0))
780     {
781          primary[1]=ideal(1);
782          primary[2]=ideal(1);
783          return(primary);
784     }
785     act=factor(j[1]);
786     for(@k=1;@k<=size(act[1]);@k++)
787     {
788       @qh=j;
789       if(@wr==0)
790       {
791          @qh[1]=act[1][@k]^act[2][@k];
792       }
793       else
794       {
795          @qh[1]=act[1][@k];   
796       }
797       primary[2*@k-1]=interred(@qh);
798       @qh=j;
799       @qh[1]=act[1][@k];
800       primary[2*@k]=interred(@qh);
801     }
802        return(primary);
803  }
804
805  if(homog(j)==1)
806  {
807     primary[1]=j;
808     if((size(ser)>0)&&(size(reduce(ser,j))==0))
809     {
810          primary[1]=ideal(1);
811          primary[2]=ideal(1);
812          return(primary);
813     }
814     if(dim(j)==-1)
815     {
816        primary[1]=ideal(1);
817        primary[2]=ideal(1);
818     }
819     else
820     {
821        primary[2]=maxideal(1);
822     }
823     return(primary);
824  }
825 
826//the first element in the standardbase is factorized
827  if(deg(j[1])>0)
828  {
829    act=factor(j[1]);
830    testFactor(act,j[1]);
831  }
832  else
833  {
834     primary[1]=ideal(1);
835     primary[2]=ideal(1);
836     return(primary);
837  }
838
839//withe the factors new ideals (hopefully the primary decomposition)
840//are created
841
842  if(size(act[1])>1)
843  {
844     if(size(#)>1)
845     {
846        primary[1]=ideal(1);
847        primary[2]=ideal(1);
848        primary[3]=ideal(1);
849        primary[4]=ideal(1);
850        return(primary);
851     }
852     for(@k=1;@k<=size(act[1]);@k++)
853     {
854        if(@wr==0)
855        {
856           primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
857        }
858        else
859        {
860           primary[2*@k-1]=std(j,act[1][@k]);
861        }
862        if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
863        {
864           primary[2*@k]   = primary[2*@k-1];
865        }
866        else
867        {
868           primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
869        }
870     }
871  }
872  else
873  {
874     primary[1]=j;
875     if((size(#)>0)&&(act[2][1]>1))
876     {
877        act[2]=1;
878        primary[1]=std(primary[1],act[1][1]);
879     }
880     
881     if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
882     {
883        primary[2]=primary[1];
884     }
885     else
886     {
887        primary[2]=primaryTest(primary[1],act[1][1]);
888     }
889  }
890  if(size(#)==0)
891  {
892     primary=splitPrimary(primary,ser,@wr,act);
893  }
894
895//test whether all ideals in the decomposition are primary and
896//in general position
897//if not after a random coordinate transformation of the last
898//variable the corresponding ideal is decomposed again.
899
900  @k=0;
901  while(@k<(size(primary)/2))
902  {
903    @k++;
904    if (size(primary[2*@k])==0)
905    {
906       jmap=randomLast(100);
907       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
908       {
909          if(deg(lead(primary[2*@k-1][@n]))==1)
910          {
911             jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0);
912          }
913       }
914       phi=@P,jmap;
915       jmap[nva]=-(jmap[nva]-2*var(nva));
916       psi=@P,jmap;   
917       @qht=primary[2*@k-1];
918       @qh=phi(@qht);
919       if(npars(@P)>0)
920       {
921          @ri= "ring @Phelp ="
922                  +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;";
923       }
924       else
925       {
926          @ri= "ring @Phelp ="
927                  +string(char(@P))+",("+varstr(@P)+",@t),dp;";
928       }
929       execute(@ri);
930       ideal @qh=homog(imap(@P,@qht),@t);
931
932       ideal @qh1=std(@qh);
933       @hilb=hilb(@qh1,1);
934       @ri= "ring @Phelp1 ="
935                  +string(char(@P))+",("+varstr(@Phelp)+"),lp;";
936       execute(@ri);
937       ideal @qh=homog(imap(@P,@qh),@t);
938       kill @Phelp;
939       @qh=std(@qh,@hilb);
940       @qh=subst(@qh,@t,1);
941       setring @P;
942       @qh=imap(@Phelp1,@qh);
943       kill @Phelp1;
944       @qh=clearSB(@qh);
945       attrib(@qh,"isSB",1);     
946
947       @lh=zero_decomp (@qh,psi(ser),@wr);
948
949       kill lres;
950       list lres;
951       if(size(@lh)==2)
952       {
953          helpprim=@lh[2];
954          lres[1]=primary[2*@k-1];
955          lres[2]=psi(helpprim);
956          if(size(reduce(lres[2],lres[1]))==0)
957          {
958             primary[2*@k]=primary[2*@k-1];
959             continue;
960          }
961       }
962       else
963       {
964          act=factor(@qh[1]);
965          if(2*size(act[1])==size(@lh))
966          {
967             for(@n=1;@n<=size(act[1]);@n++)
968             {
969                @f=act[1][@n]^act[2][@n];
970                lres[2*@n-1]=interred(primary[2*@k-1]+psi(@f));
971                helpprim=@lh[2*@n];
972                lres[2*@n]=psi(helpprim);
973             }
974          }
975          else
976          {
977             lres=psi(@lh);
978          }
979       }
980       if(npars(@P)>0)
981       {
982          @ri= "ring @Phelp ="
983                  +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;";
984       }
985       else
986       {
987          @ri= "ring @Phelp ="
988                  +string(char(@P))+",("+varstr(@P)+",@t),dp;";
989       }
990       execute(@ri);
991       list @lvec;
992       list @lr=imap(@P,lres);
993       ideal @lr1;
994       
995       if(size(@lr)==2)
996       {
997          @lr[2]=homog(@lr[2],@t);
998          @lr1=std(@lr[2]);
999          @lvec[2]=hilb(@lr1,1);     
1000       }
1001       else
1002       {
1003          for(@n=1;@n<=size(@lr)/2;@n++)
1004          {
1005             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1006             {
1007                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1008                @lr1=std(@lr[2*@n-1]);
1009                @lvec[2*@n-1]=hilb(@lr1,1);
1010                @lvec[2*@n]=@lvec[2*@n-1];
1011             }
1012             else
1013             {
1014                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
1015                @lr1=std(@lr[2*@n-1]);
1016                @lvec[2*@n-1]=hilb(@lr1,1);
1017                @lr[2*@n]=homog(@lr[2*@n],@t);
1018                @lr1=std(@lr[2*@n]);
1019                @lvec[2*@n]=hilb(@lr1,1);
1020
1021             }
1022         }
1023       }
1024       @ri= "ring @Phelp1 ="
1025                  +string(char(@P))+",("+varstr(@Phelp)+"),lp;";
1026       execute(@ri);
1027       list @lr=imap(@Phelp,@lr);
1028
1029       kill @Phelp;
1030       if(size(@lr)==2)
1031       {
1032          @lr[2]=std(@lr[2],@lvec[2]);
1033          @lr[2]=subst(@lr[2],@t,1);
1034     
1035       }
1036       else
1037       {
1038          for(@n=1;@n<=size(@lr)/2;@n++)
1039          {
1040             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
1041             {
1042                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1043                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1044                @lr[2*@n]=@lr[2*@n-1];
1045                attrib(@lr[2*@n],"isSB",1);
1046             }
1047             else
1048             {
1049                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
1050                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
1051                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
1052                @lr[2*@n]=subst(@lr[2*@n],@t,1);
1053             }
1054          }
1055       }
1056       kill @lvec;
1057       setring @P;
1058       lres=imap(@Phelp1,@lr);
1059       kill @Phelp1;
1060       for(@n=1;@n<=size(lres);@n++)
1061       {
1062          lres[@n]=clearSB(lres[@n]);
1063          attrib(lres[@n],"isSB",1);
1064       }
1065       
1066       primary[2*@k-1]=lres[1];
1067       primary[2*@k]=lres[2];
1068       @s=size(primary)/2;
1069       for(@n=1;@n<=size(lres)/2-1;@n++)
1070       {
1071         primary[2*@s+2*@n-1]=lres[2*@n+1];
1072         primary[2*@s+2*@n]=lres[2*@n+2];
1073       }
1074       @k--;
1075     }
1076  }
1077  return(primary);
1078}
1079example
1080{ "EXAMPLE:"; echo = 2;
1081   ring  r = 0,(x,y,z),lp;
1082   poly  p = z2+1;
1083   poly  q = z4+2;
1084   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1085   i=std(i);
1086   list  pr= zero_decomp(i,ideal(0),0);
1087   pr;
1088}
1089
1090////////////////////////////////////////////////////////////////////////////////
1091
1092proc ggt (ideal i)
1093USAGE:   ggt(i); i list of polynomials
1094RETURN:  poly = ggt(i[1],...,i[size(i)])
1095NOTE:   
1096EXAMPLE: example ggt; shows an example
1097{
1098  int k;
1099  poly p=i[1];
1100  if(deg(p)==0)
1101  {
1102    return(1);
1103  }
1104
1105 
1106  for (k=2;k<=size(i);k++)
1107  {
1108     if(deg(i[k])==0)
1109     {
1110        return(1)
1111     }
1112     p=GCD(p,i[k]);
1113     if(deg(p)==0)
1114     {
1115        return(1);
1116     }
1117  }
1118  return(p);
1119}
1120example
1121{ "EXAMPLE:"; echo = 2;
1122   ring  r = 0,(x,y,z),lp;
1123   poly  p = (x+y)*(y+z);
1124   poly  q = (z4+2)*(y+z);
1125   ideal l=p,q;
1126   poly  pr= ggt(l);
1127   pr;
1128}
1129///////////////////////////////////////////////////////////////////////////////
1130proc gcdTest(ideal act)
1131{
1132  int i,j;
1133  if(size(act)<=1)
1134  {
1135     return(0);
1136  }
1137  for (i=1;i<=size(act)-1;i++)
1138  {
1139     for(j=i+1;j<=size(act);j++)
1140     {
1141        if(deg(std(ideal(act[i],act[j]))[1])>0)
1142        {
1143           return(0);
1144        }
1145     }
1146  }
1147  return(1);
1148}
1149
1150///////////////////////////////////////////////////////////////////////////////
1151proc coeffLcm(ideal h)
1152{
1153   string @pa=parstr(basering);
1154   if(size(@pa)==0)
1155   {
1156      return(lcm(h));
1157   }
1158   def bsr= basering;
1159   string @id=string(h);
1160   execute "ring @r=0,("+@pa+","+varstr(bsr)+"),dp;";
1161   execute "ideal @i="+@id+";";
1162   poly @p=lcm(@i);
1163   string @ps=string(@p);
1164   setring bsr;
1165   execute "poly @p="+@ps+";";
1166   return(@p);
1167}
1168example
1169{
1170   "EXAMPLE:"; echo = 2;
1171   ring  r =( 0,a,b),(x,y,z),lp;
1172   poly  p = (a+b)*(y-z);
1173   poly  q = (a+b)*(y+z);
1174   ideal l=p,q;
1175   poly  pr= coeffLcm(l);
1176   pr;
1177}
1178
1179///////////////////////////////////////////////////////////////////////////////
1180
1181proc lcm (ideal i)
1182USAGE:   lcm(i); i list of polynomials
1183RETURN:  poly = lcm(i[1],...,i[size(i)])
1184NOTE:   
1185EXAMPLE: example lcm; shows an example
1186{
1187  int k,j;
1188   poly p,q;
1189  i=simplify(i,10);
1190  for(j=1;j<=size(i);j++)
1191  {
1192    if(deg(i[j])>0)
1193    {
1194      p=i[j];
1195      break;
1196    }
1197  }
1198  if(deg(p)==-1)
1199  {
1200    return(1);
1201  }
1202  for (k=j+1;k<=size(i);k++)
1203  {
1204     if(deg(i[k])!=0)
1205     {
1206        q=GCD(p,i[k]);
1207        if(deg(q)==0)
1208        {
1209           p=p*i[k];
1210        }
1211        else
1212        {
1213           p=p/q;
1214           p=p*i[k];
1215        }
1216     }
1217   }
1218  return(p);
1219}
1220example
1221{ "EXAMPLE:"; echo = 2;
1222   ring  r = 0,(x,y,z),lp;
1223   poly  p = (x+y)*(y+z);
1224   poly  q = (z4+2)*(y+z);
1225   ideal l=p,q;
1226   poly  pr= lcm(l);
1227   pr;
1228   l=1,-1,p,1,-1,q,1;
1229   pr=lcm(l);
1230   pr;
1231}
1232
1233///////////////////////////////////////////////////////////////////////////////
1234proc clearSB (ideal i,list #)
1235USAGE:   clearSB(i); i ideal which is SB ordered by monomial ordering
1236RETURN:  ideal = minimal SB
1237NOTE:   
1238EXAMPLE: example clearSB; shows an example
1239{
1240  int k,j;
1241  poly m;
1242  int c=size(i);
1243 
1244  if(size(#)==0)
1245  {
1246    for(j=1;j<c;j++)
1247    {
1248      if(deg(i[j])==0)
1249      {
1250        i=ideal(1);
1251        return(i);
1252      }   
1253      if(deg(i[j])>0)
1254      {
1255        m=lead(i[j]);
1256        for(k=j+1;k<=c;k++)
1257        {
1258           if(size(lead(i[k])/m)>0)
1259           {
1260             i[k]=0;
1261           }
1262        }
1263      }
1264    }
1265  }
1266  else
1267  {
1268    j=0;
1269    while(j<c-1)
1270    {
1271      j++;
1272      if(deg(i[j])==0)
1273      {
1274        i=ideal(1);
1275        return(i);
1276      }   
1277      if(deg(i[j])>0)
1278      {
1279        m=lead(i[j]);
1280        for(k=j+1;k<=c;k++)
1281        {
1282           if(size(lead(i[k])/m)>0)
1283           {
1284             if((leadexp(m)!=leadexp(i[k]))||(#[j]<=#[k]))
1285             {
1286                i[k]=0;
1287             }
1288             else
1289             {
1290                i[j]=0;
1291                break; 
1292             }
1293           }
1294        }
1295      }
1296    }
1297  }
1298  return(simplify(i,2));
1299}
1300example
1301{ "EXAMPLE:"; echo = 2;
1302   LIB "primdec.lib";
1303   ring  r = (0,a,b),(x,y,z),dp;
1304   ideal i=ax2+y,a2x+y,bx;
1305   list l=1,2,1;
1306   ideal j=clearSB(i,l);
1307   j;
1308}
1309
1310///////////////////////////////////////////////////////////////////////////////
1311
1312proc independSet (ideal j)
1313USAGE:   independentSet(i); i ideal
1314RETURN:  list = new varstring with the independent set at the end,
1315                ordstring with the corresponding block ordering,
1316                the integer where the independent set starts in the varstring
1317NOTE:   
1318EXAMPLE: example independentSet; shows an example
1319{
1320   int n,k,di;
1321   list resu,hilf;
1322   string var1,var2;
1323   list v=system("indsetall",j,1);
1324   
1325   for(n=1;n<=size(v);n++)
1326   {
1327      di=0;
1328      var1="";
1329      var2="";
1330      for(k=1;k<=size(v[n]);k++)
1331      {
1332         if(v[n][k]!=0)
1333         {
1334            di++;
1335            var2=var2+"var("+string(k)+"),";
1336         }
1337         else
1338         {
1339            var1=var1+"var("+string(k)+"),";
1340         }
1341      }
1342      if(di>0)
1343      {
1344         var1=var1+var2;
1345         var1=var1[1..size(var1)-1];
1346         hilf[1]=var1;
1347         hilf[2]="lp";
1348         //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
1349         hilf[3]=di;
1350         resu[n]=hilf;
1351      }
1352      else
1353      {
1354         resu[n]=varstr(basering),ordstr(basering),0;
1355      }
1356   }
1357   return(resu);
1358}
1359example
1360{ "EXAMPLE:"; echo = 2;
1361   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1362   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1363   i=std(i);
1364   list  l=independSet(i);
1365   l;
1366   i=i,g;
1367   l=independSet(i);
1368   l;
1369
1370   ring s=0,(x,y,z),lp;
1371   ideal i=z,yx;
1372   list l=independSet(i);
1373   l;
1374
1375
1376}
1377///////////////////////////////////////////////////////////////////////////////
1378
1379proc maxIndependSet (ideal j)
1380USAGE:   maxIndependentSet(i); i ideal
1381RETURN:  list = new varstring with the maximal independent set at the end,
1382                ordstring with the corresponding block ordering,
1383                the integer where the independent set starts in the varstring
1384NOTE:   
1385EXAMPLE: example maxIndependentSet; shows an example
1386{
1387   int n,k,di;
1388   list resu,hilf;
1389   string var1,var2;
1390   list v=system("indsetall",j,0);
1391   
1392   for(n=1;n<=size(v);n++)
1393   {
1394      di=0;
1395      var1="";
1396      var2="";
1397      for(k=1;k<=size(v[n]);k++)
1398      {
1399         if(v[n][k]!=0)
1400         {
1401            di++;
1402            var2=var2+"var("+string(k)+"),";
1403         }
1404         else
1405         {
1406            var1=var1+"var("+string(k)+"),";
1407         }
1408      }
1409      if(di>0)
1410      {
1411         var1=var1+var2;
1412         var1=var1[1..size(var1)-1];
1413         hilf[1]=var1;
1414         hilf[2]="lp";
1415         hilf[3]=di;
1416         resu[n]=hilf;
1417      }
1418      else
1419      {
1420         resu[n]=varstr(basering),ordstr(basering),0;
1421      }
1422   }
1423   return(resu);
1424}
1425example
1426{ "EXAMPLE:"; echo = 2;
1427   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1428   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1429   i=std(i);
1430   list  l=maxIndependSet(i);
1431   l;
1432   i=i,g;
1433   l=maxIndependSet(i);
1434   l;
1435
1436   ring s=0,(x,y,z),lp;
1437   ideal i=z,yx;
1438   list l=maxIndependSet(i);
1439   l;
1440
1441
1442}
1443
1444///////////////////////////////////////////////////////////////////////////////
1445
1446proc prepareQuotientring (int nnp)
1447USAGE:   prepareQuotientring(nnp); nnp int
1448RETURN:  string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
1449NOTE:   
1450EXAMPLE: example independentSet; shows an example
1451{ 
1452  ideal @ih,@jh;
1453  int npar=npars(basering);
1454  int @n;
1455 
1456  string quotring= "ring quring = ("+charstr(basering);
1457  for(@n=nnp+1;@n<=nvars(basering);@n++)
1458  {
1459     quotring=quotring+",var("+string(@n)+")";
1460     @ih=@ih+var(@n);
1461  }
1462 
1463  quotring=quotring+"),(var(1)";
1464  @jh=@jh+var(1);
1465  for(@n=2;@n<=nnp;@n++)
1466  {
1467    quotring=quotring+",var("+string(@n)+")";
1468    @jh=@jh+var(@n);
1469  }
1470  quotring=quotring+"),lp;";
1471 
1472  return(quotring);
1473
1474}
1475example
1476{ "EXAMPLE:"; echo = 2;
1477   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
1478   def @Q=basering;
1479   list l= prepareQuotientring(3);
1480   l;
1481   execute l[1];
1482   execute l[2];
1483   basering;
1484   phi;
1485   setring @Q;
1486   
1487}
1488
1489///////////////////////////////////////////////////////////////////////////////
1490proc cleanPrimary(list l)
1491{
1492   int i,j;
1493   list lh;
1494   for(i=1;i<=size(l)/2;i++)
1495   {
1496      if(deg(l[2*i-1][1])>0)
1497      {
1498         j++;
1499         lh[j]=l[2*i-1];
1500         j++;
1501         lh[j]=l[2*i];
1502      }
1503   }
1504   return(lh);
1505}
1506///////////////////////////////////////////////////////////////////////////////
1507
1508proc minAssPrimes(ideal i, list #)
1509USAGE:   minAssPrimes(i); i ideal
1510         minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
1511RETURN:  list = the minimal associated prime ideals of i
1512NOTE:   
1513EXAMPLE: example minAssPrimes; shows an example
1514{
1515   def @P=basering;
1516   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
1517//   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;";
1518   ideal i=fetch(@P,i);
1519   if(size(#)==0)
1520   {
1521      int @wr;
1522      list tluser,@res;
1523      list primary=decomp(i,2);
1524
1525      @res[1]=primary;
1526
1527      tluser=union(@res);
1528      setring @P;
1529      list @res=imap(gnir,tluser);
1530      return(@res);
1531   }
1532   list @res,empty;
1533   option(redSB);
1534   list @pr=facstd(i);
1535   option(noredSB);
1536   int j,k,odim,ndim,count;
1537   attrib(@pr[1],"isSB",1);
1538   if(#[1]==77)
1539   {
1540     odim=dim(@pr[1]);
1541     count=1;
1542     intvec pos;
1543     pos[size(@pr)]=0;
1544     for(j=2;j<=size(@pr);j++)
1545     {
1546        attrib(@pr[j],"isSB",1);
1547        ndim=dim(@pr[j]);
1548        if(ndim>odim)
1549        {
1550           for(k=count;k<=j-1;k++)
1551           {
1552              pos[k]=1;
1553           }
1554           count=j;
1555           odim=ndim;
1556        }
1557        if(ndim<odim)     
1558        {
1559           pos[j]=1;
1560        }
1561     }
1562     for(j=1;j<=size(@pr);j++)
1563     {
1564        if(pos[j]!=1)
1565        {
1566            @res[j]=decomp(@pr[j],2);
1567        }
1568        else
1569        {
1570           @res[j]=empty;
1571        }
1572     }
1573   }
1574   else
1575   {
1576     for(j=1;j<=size(@pr);j++)
1577     {
1578       @res[j]=decomp(@pr[j],2);
1579     }     
1580   }
1581
1582   @res=union(@res);
1583   setring @P;
1584   list @res=imap(gnir,@res);
1585   return(@res);
1586}
1587example
1588{ "EXAMPLE:"; echo = 2;
1589   ring  r = 32003,(x,y,z),lp;
1590   poly  p = z2+1;
1591   poly  q = z4+2;
1592   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
1593   LIB "primaryDecomposition.lib";
1594   list pr= minAssPrimes(i);
1595   pr;
1596   pr= minAssPrimes(i,1);
1597}
1598
1599///////////////////////////////////////////////////////////////////////////////
1600
1601proc union(list l)
1602{
1603  int i,j,k;
1604  list @erg;
1605  i=0;
1606
1607  for(k=1;k<=size(l);k++)
1608  {
1609     for(j=1;j<=size(l[k])/2;j++)
1610     {
1611        if(deg(l[k][2*j][1])!=0)
1612        {
1613           i++;
1614           @erg[i]=l[k][2*j];
1615        }
1616     }
1617  }
1618
1619  list @wos;
1620  i=0;
1621  ideal i1,i2;
1622  while(i<size(@erg)-1)
1623  {
1624     i++;
1625     k=i+1;
1626     i1=lead(@erg[i]);
1627      attrib(i1,"isSB",1);
1628      attrib(@erg[i],"isSB",1);
1629
1630     while(k<=size(@erg))
1631     {
1632        if(deg(@erg[i][1])==0)
1633        {
1634           break;
1635        }
1636        i2=lead(@erg[k]);
1637        attrib(@erg[k],"isSB",1);
1638        attrib(i2,"isSB",1);
1639
1640        if(size(reduce(i1,i2,1))==0)
1641        {
1642           if(size(reduce(@erg[i],@erg[k]))==0)
1643           {
1644              @erg[k]=ideal(1);
1645              i2=ideal(1);
1646           }         
1647        }
1648        if(size(reduce(i2,i1,1))==0)
1649        {
1650           if(size(reduce(@erg[k],@erg[i]))==0)
1651           {
1652              break;
1653           }         
1654        }
1655        k++;
1656        if(k>size(@erg))
1657        {
1658           @wos[size(@wos)+1]=@erg[i];
1659        }
1660     }
1661  }
1662  if(deg(@erg[size(@erg)][1])!=0)
1663  {
1664     @wos[size(@wos)+1]=@erg[size(@erg)];
1665  }
1666  return(@wos);
1667}
1668///////////////////////////////////////////////////////////////////////////////
1669proc radical(ideal i)
1670{
1671   list pr=minAssPrimes(i,1);
1672   int j;
1673   ideal k=pr[1];
1674   for(j=2;j<=size(pr);j++)
1675   {
1676      k=intersect(k,pr[j]);
1677   }
1678   return(k);
1679}
1680///////////////////////////////////////////////////////////////////////////////
1681proc decomp (ideal i,list #)
1682USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
1683         decomp(i,1);        (for the minimal associated primes) )         
1684RETURN:  list = list of primary ideals and their associated primes
1685         (at even positions in the list)
1686         (resp. a list of the minimal associated primes)
1687NOTE:    Algorithm of Gianni, Traeger, Zacharias
1688EXAMPLE: example decomp; shows an example
1689{
1690  def  @P = basering;
1691  list primary,indep;
1692  intvec @vh,isat;
1693  int @wr,@k,@n,@m,@n1,@n2,@n3,homo;
1694  ideal peek=i;
1695  ideal ser,tras;
1696 
1697  int @aa=timer;
1698 
1699  homo=homog(i);
1700  if(size(#)>0)
1701  {
1702     if((#[1]==1)||(#[1]==2))
1703     {
1704        @wr=#[1];
1705        if(size(#)>1)
1706        {
1707          peek=#[2];
1708          ser=#[3];
1709        }
1710      }
1711      else
1712      {
1713         peek=#[1];
1714         ser=#[2];
1715      }
1716  }
1717
1718  if(homo==1)
1719  {
1720     tras=std(i);
1721     if(dim(tras)==0)
1722     {
1723        primary[1]=tras;
1724        primary[2]=maxideal(1);
1725        if(@wr>0)
1726        {
1727           list l;
1728           l[1]=maxideal(1);
1729           l[2]=maxideal(1);
1730           return(l);
1731        }
1732        return(primary);
1733     }
1734     intvec @hilb=hilb(tras,1);
1735  }
1736 
1737  //----------------------------------------------------------------
1738  //i is the zero-ideal
1739  //----------------------------------------------------------------
1740 
1741  if(size(i)==0)
1742  {
1743     primary=i,i;
1744     return(primary);
1745  }
1746 
1747  //----------------------------------------------------------------
1748  //pass to the lexicographical ordering and compute a standardbasis
1749  //----------------------------------------------------------------
1750
1751  execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
1752 
1753  option(redSB);
1754  ideal ser=fetch(@P,ser);
1755  ideal peek=std(fetch(@P,peek));
1756  homo=homog(peek);
1757   
1758  if(homo==1)
1759  {
1760     if(ordstr(@P)[1,2]!="lp")
1761     {
1762        ideal @j=std(fetch(@P,i),@hilb);
1763     }
1764     else
1765     {
1766        ideal @j=fetch(@P,tras);
1767        attrib(@j,"isSB",1);
1768     }
1769  }
1770  else
1771  {
1772     ideal @j=std(fetch(@P,i));
1773  }
1774
1775  //----------------------------------------------------------------
1776  //j is the ring
1777  //----------------------------------------------------------------
1778
1779  if (dim(@j)==-1)
1780  {
1781     setring @P;
1782     option(noredSB);
1783     return(ideal(0));
1784  }
1785 
1786  //----------------------------------------------------------------
1787  //  the case of one variable
1788  //----------------------------------------------------------------
1789
1790  if(nvars(basering)==1)
1791  {
1792     list fac=factor(@j[1]);
1793     list gprimary;
1794     for(@k=1;@k<=size(fac[1]);@k++)
1795     {
1796        if(@wr==0)
1797        {
1798           gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
1799           gprimary[2*@k]=ideal(fac[1][@k]);
1800        }
1801        else
1802        {
1803           gprimary[2*@k-1]=ideal(fac[1][@k]);
1804           gprimary[2*@k]=ideal(fac[1][@k]);
1805        }
1806     }
1807     setring @P;
1808     option(noredSB);
1809     primary=fetch(gnir,gprimary);
1810
1811     return(primary);
1812  }
1813
1814 //------------------------------------------------------------------
1815 //the zero-dimensional case
1816 //------------------------------------------------------------------
1817
1818  if (dim(@j)==0)
1819  {
1820      list gprimary= zero_decomp(@j,ser,@wr);
1821
1822      setring @P;
1823      option(noredSB);
1824      primary=fetch(gnir,gprimary);
1825      if(size(ser)>0)
1826      {
1827         primary=cleanPrimary(primary);
1828      }
1829      return(primary);
1830   }
1831
1832
1833  //------------------------------------------------------------------
1834  //search for a maximal independent set indep,i.e.
1835  //look for subring such that the intersection with the ideal is zero
1836  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
1837  //indep[1] is the new varstring and indep[2] the string for the block-ordering
1838  //------------------------------------------------------------------
1839
1840  poly @gs,@gh,@p;
1841  string @va,quotring;
1842  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
1843  ideal @h;
1844  int jdim=dim(@j);
1845  list fett;
1846  int lauf,di; 
1847
1848  if(@wr!=1)
1849  {
1850     allindep=independSet(@j);
1851     for(@m=1;@m<=size(allindep);@m++)
1852     {
1853        if(allindep[@m][3]==jdim)
1854        {
1855           di++;
1856           indep[di]=allindep[@m];
1857        }
1858        else
1859        {
1860           lauf++;
1861           restindep[lauf]=allindep[@m];
1862        }
1863     }
1864   }
1865   else
1866   {
1867      indep=maxIndependSet(@j);
1868   }
1869 
1870  ideal jkeep=@j;
1871
1872  if(ordstr(@P)[1]=="w")
1873  {
1874     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");";
1875  }
1876  else
1877  {
1878     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),dp;";
1879  }
1880  ideal jwork=std(imap(gnir,@j));
1881  poly @p,@q;
1882  ideal @h,fac;
1883  di=dim(jwork);
1884  setring gnir;
1885  for(@m=1;@m<=size(indep);@m++)
1886  {
1887     isat=0;
1888     @n2=0;
1889     if((indep[@m][1]==varstr(basering))&&(@m==1))
1890     //this is the good case, nothing to do, just to have the same notations
1891     //change the ring
1892     {
1893        execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
1894                              +ordstr(basering)+");";
1895        ideal @j=fetch(gnir,@j);
1896        attrib(@j,"isSB",1);
1897     }
1898     else
1899     {
1900        @va=string(maxideal(1));
1901        execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
1902                              +indep[@m][2]+");";
1903        execute "map phi=gnir,"+@va+";";
1904        if(homo==1)
1905        {
1906           ideal @j=std(phi(@j),@hilb);
1907        }
1908        else
1909        {
1910          ideal @j=std(phi(@j));
1911        }
1912      }
1913     if((deg(@j[1])==0)||(dim(@j)<jdim))
1914     {
1915       setring gnir;
1916       break;
1917     }
1918     for (lauf=1;lauf<=size(@j);lauf++)
1919     {
1920         fett[lauf]=size(@j[lauf]);
1921     }
1922     //------------------------------------------------------------------------------------
1923     //we have now the following situation:
1924     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
1925     //to this quotientring, j is their still a standardbasis, the
1926     //leading coefficients of the polynomials  there (polynomials in
1927     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
1928     //we need their ggt, gh, because of the following:
1929     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
1930     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
1931     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
1932 
1933     //------------------------------------------------------------------------------------
1934
1935     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
1936     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
1937     //-------------------------------------------------------------------------------------
1938
1939     quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
1940
1941     //---------------------------------------------------------------------
1942     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
1943     //---------------------------------------------------------------------
1944
1945     execute quotring;
1946
1947    // @j considered in the quotientring
1948     ideal @j=imap(gnir1,@j);
1949     ideal ser=imap(gnir,ser);
1950
1951     kill gnir1;
1952 
1953     //j is a standardbasis in the quotientring but usually not minimal
1954     //here it becomes minimal
1955
1956     @j=clearSB(@j,fett);
1957     attrib(@j,"isSB",1);
1958
1959     //we need later ggt(h[1],...)=gh for saturation
1960     ideal @h;
1961     if(deg(@j[1])>0)
1962     {
1963        for(@n=1;@n<=size(@j);@n++)
1964        {
1965           @h[@n]=leadcoef(@j[@n]);
1966        }
1967        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
1968   
1969        list uprimary= zero_decomp(@j,ser,@wr);
1970
1971     }
1972     else
1973     {
1974       list uprimary;
1975       uprimary[1]=ideal(1);
1976       uprimary[2]=ideal(1);
1977     }
1978
1979     //we need the intersection of the ideals in the list quprimary with the
1980     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
1981     //but fi polynomials, then the intersection of q with the polynomialring
1982     //is the saturation of the ideal generated by f1,...,fr with respect to
1983     //h which is the lcm of the leading coefficients of the fi considered in the
1984     //quotientring: this is coded in saturn
1985
1986     list saturn;
1987     ideal hpl;
1988 
1989     for(@n=1;@n<=size(uprimary);@n++)
1990     {
1991        hpl=0;
1992        for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
1993        {
1994           hpl=hpl,leadcoef(uprimary[@n][@n1]);
1995        }
1996        saturn[@n]=hpl;
1997     }
1998     //--------------------------------------------------------------------
1999     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2000     //back to the polynomialring
2001     //---------------------------------------------------------------------
2002     setring gnir;
2003       
2004     collectprimary=imap(quring,uprimary);
2005     lsau=imap(quring,saturn);
2006     @h=imap(quring,@h); 
2007
2008     kill quring;
2009
2010
2011     @n2=size(quprimary);
2012     @n3=@n2;
2013   
2014     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2015     {
2016        if(deg(collectprimary[2*@n1][1])>0)
2017        {
2018           @n2++;
2019           quprimary[@n2]=collectprimary[2*@n1-1];
2020           lnew[@n2]=lsau[2*@n1-1];
2021           @n2++;
2022           lnew[@n2]=lsau[2*@n1];
2023           quprimary[@n2]=collectprimary[2*@n1];
2024        }
2025     } 
2026
2027     //here the intersection with the polynomialring
2028     //mentioned above is really computed
2029
2030    for(@n=@n3/2+1;@n<=@n2/2;@n++)
2031     {
2032        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
2033        {
2034           quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2035           quprimary[2*@n]=quprimary[2*@n-1];
2036        }
2037        else
2038        {
2039           if(@wr==0)
2040           {
2041              quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2042           }
2043           quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
2044        }
2045     }
2046     if(size(@h)>0)
2047     {
2048           //---------------------------------------------------------------
2049           //we change to @Phelp to have the ordering dp for saturation
2050           //---------------------------------------------------------------
2051        setring @Phelp;
2052        @h=imap(gnir,@h);
2053        if(@wr!=1)
2054//        if(@wr==0)
2055        {
2056           @q=minSat(jwork,@h)[2];
2057        }
2058        else
2059        {
2060            fac=ideal(0);
2061            for(lauf=1;lauf<=ncols(@h);lauf++)
2062           {
2063              if(deg(@h[lauf])>0)
2064              {
2065                 fac=fac+factorize(@h[lauf],1);
2066              }
2067           }
2068           fac=simplify(fac,4);
2069           @q=1;
2070           for(lauf=1;lauf<=size(fac);lauf++)
2071           {
2072             @q=@q*fac[lauf];
2073           }
2074        }
2075        jwork=std(jwork,@q);
2076        if(dim(jwork)<di)
2077        {
2078           setring gnir;
2079           @j=imap(@Phelp,jwork);
2080           break;
2081        }
2082        if(homo==1)
2083        {
2084              @hilb=hilb(jwork,1);
2085        }
2086       
2087        setring gnir;
2088        @j=imap(@Phelp,jwork);
2089     }   
2090  }
2091  if((size(quprimary)==0)&&(@wr>0))
2092  {
2093     @j=ideal(1);
2094     quprimary[1]=ideal(1);
2095     quprimary[2]=ideal(1);
2096  }
2097  //---------------------------------------------------------------
2098  //notice that j=sat(j,gh) intersected with (j,gh^n)
2099  //we finished with sat(j,gh) and have to start with (j,gh^n)
2100  //---------------------------------------------------------------
2101  if((deg(@j[1])!=0)&&(@wr!=1))
2102  {
2103     int uq=size(quprimary);
2104     if(uq>0)
2105     {
2106        if(@wr==0)
2107        {
2108           ideal htest=quprimary[1];
2109
2110           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
2111           {
2112              htest=intersect(htest,quprimary[2*@n1-1]);
2113           }
2114        }
2115        else
2116        {
2117           ideal htest=quprimary[2];
2118
2119           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
2120           {
2121              htest=intersect(htest,quprimary[2*@n1]);
2122           }
2123        }
2124        if(size(ser)>0)
2125        {
2126           htest=intersect(htest,ser);
2127        }
2128        ser=std(htest);
2129     }
2130     //we are not ready yet
2131     if (specialIdealsEqual(ser,peek)!=1)
2132     {
2133        for(@m=1;@m<=size(restindep);@m++)
2134        {
2135           isat=0;
2136           @n2=0;
2137           if(restindep[@m][1]==varstr(basering))
2138           //this is the good case, nothing to do, just to have the same notations
2139           //change the ring
2140           {
2141              execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
2142                              +ordstr(basering)+");";
2143              ideal @j=fetch(gnir,jkeep);
2144              attrib(@j,"isSB",1);
2145           }
2146           else
2147           {
2148              @va=string(maxideal(1));
2149              execute "ring gnir1 = ("+charstr(basering)+"),("+restindep[@m][1]+"),("
2150                              +restindep[@m][2]+");";
2151              execute "map phi=gnir,"+@va+";";
2152              if(homo==1)
2153              {
2154                 ideal @j=std(phi(jkeep),@hilb);
2155              }
2156              else
2157              {
2158                ideal @j=std(phi(jkeep));
2159              }
2160           }
2161           
2162           for (lauf=1;lauf<=size(@j);lauf++)
2163           {
2164              fett[lauf]=size(@j[lauf]);
2165           }
2166           //------------------------------------------------------------------------------------
2167           //we have now the following situation:
2168           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
2169           //to this quotientring, j is their still a standardbasis, the
2170           //leading coefficients of the polynomials  there (polynomials in
2171           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
2172           //we need their ggt, gh, because of the following:
2173           //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
2174           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
2175           //on the other hand j=(j,gh^n) intersected with (j:gh^n)
2176 
2177           //------------------------------------------------------------------------------------
2178
2179           //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
2180           //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
2181           //-------------------------------------------------------------------------------------
2182
2183           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
2184
2185           //---------------------------------------------------------------------
2186           //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2187           //---------------------------------------------------------------------
2188           
2189           execute quotring;
2190
2191           // @j considered in the quotientring
2192           ideal @j=imap(gnir1,@j);
2193           ideal ser=imap(gnir,ser);
2194
2195           kill gnir1;
2196 
2197           //j is a standardbasis in the quotientring but usually not minimal
2198           //here it becomes minimal
2199           @j=clearSB(@j,fett);
2200           attrib(@j,"isSB",1);
2201
2202           //we need later ggt(h[1],...)=gh for saturation
2203           ideal @h;
2204
2205           for(@n=1;@n<=size(@j);@n++)
2206           {
2207              @h[@n]=leadcoef(@j[@n]);
2208           }
2209           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
2210   
2211            list uprimary= zero_decomp(@j,ser,@wr);
2212       
2213           //we need the intersection of the ideals in the list quprimary with the
2214           //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
2215           //but fi polynomials, then the intersection of q with the polynomialring
2216           //is the saturation of the ideal generated by f1,...,fr with respect to
2217           //h which is the lcm of the leading coefficients of the fi considered in the
2218           //quotientring: this is coded in saturn
2219
2220           list saturn;
2221           ideal hpl;
2222 
2223           for(@n=1;@n<=size(uprimary);@n++)
2224           {
2225              hpl=0;
2226              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
2227              {
2228                 hpl=hpl,leadcoef(uprimary[@n][@n1]);
2229              }
2230               saturn[@n]=hpl;
2231           }
2232           //--------------------------------------------------------------------
2233           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
2234           //back to the polynomialring
2235           //---------------------------------------------------------------------
2236           setring gnir;
2237       
2238           collectprimary=imap(quring,uprimary);
2239           lsau=imap(quring,saturn);
2240           @h=imap(quring,@h); 
2241
2242           kill quring;
2243
2244
2245           @n2=size(quprimary);
2246           @n3=@n2;
2247   
2248           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
2249           {
2250              if(deg(collectprimary[2*@n1][1])>0)
2251              {
2252                 @n2++;
2253                 quprimary[@n2]=collectprimary[2*@n1-1];
2254                 lnew[@n2]=lsau[2*@n1-1];
2255                 @n2++;
2256                 lnew[@n2]=lsau[2*@n1];
2257                 quprimary[@n2]=collectprimary[2*@n1];
2258              }
2259           } 
2260
2261           //here the intersection with the polynomialring
2262           //mentioned above is really computed
2263
2264           for(@n=@n3/2+1;@n<=@n2/2;@n++)
2265           {
2266              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
2267              {
2268                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2269                 quprimary[2*@n]=quprimary[2*@n-1];
2270              }
2271              else
2272              {
2273                 if(@wr==0)
2274                 {
2275                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
2276                 }
2277                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
2278             }
2279             if(@wr==0)
2280             {
2281                ser=std(intersect(ser,quprimary[2*@n-1]));
2282             }
2283             else
2284             {
2285                ser=std(intersect(ser,quprimary[2*@n]));
2286             }
2287           }
2288        }
2289        //we are not ready yet
2290        if (specialIdealsEqual(ser,peek)!=1)
2291        {
2292           if(@wr>0)
2293           {
2294              htprimary=decomp(@j,@wr,peek,ser);
2295           }
2296           else
2297           {
2298              htprimary=decomp(@j,peek,ser);
2299           }   
2300           // here we collect now both results primary(sat(j,gh))
2301           // and primary(j,gh^n)
2302   
2303           @n=size(quprimary);
2304           for (@k=1;@k<=size(htprimary);@k++)
2305           {
2306              quprimary[@n+@k]=htprimary[@k];
2307           }
2308        }
2309     }
2310   }
2311  //------------------------------------------------------------
2312  //back to the ring we started with
2313  //the final result: primary
2314  //------------------------------------------------------------
2315  setring @P;
2316  primary=imap(gnir,quprimary);
2317
2318  option(noredSB);
2319  return(primary);
2320}
2321
2322
2323example
2324{ "EXAMPLE:"; echo = 2;
2325   ring  r = 32003,(x,y,z),lp;
2326   poly  p = z2+1;
2327   poly  q = z4+2;
2328   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
2329   LIB "primdec.lib";
2330   list pr= decomp(i);
2331   pr;
2332   testPrimary( pr, i);
2333}
Note: See TracBrowser for help on using the repository browser.