# source:git/Singular/LIB/primdecint.lib@6e8b02

spielwiese
Last change on this file since 6e8b02 was 6e8b02, checked in by Hans Schoenemann <hannes@…>, 9 years ago
fix: prepareQuotientring
• Property mode set to `100644`
File size: 43.8 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="version primdecint.lib 4.0.0.0 Jun_2013 "; // \$Id\$
3category = "Commutative Algebra";
4info="
5LIBRARY:  primdecint.lib   primary decomposition of an ideal in the polynomial
6                           ring over the integers
7
8AUTHORS:  G. Pfister       pfister@mathematik.uni-kl.de
10@*        S. Steidel       steidel@mathematik.uni-kl.de
11
12OVERVIEW:
13
14  A library for computing the primary decomposition of an ideal in the
15  polynomial ring over the integers, Z[x_1,...,x_n].
16  The first procedure 'primdecZ' can be used in parallel.
17
18  Reference: Pfister,   Sadiq, Steidel , \"An Algorithm for primary decomposition in polynomial rings over the integers\" , arXiv:1008.2074
19
20PROCEDURES:
21 primdecZ(I);       compute the primary decomposition of ideal I
22 primdecZM(I);      compute the primary decomposition of module I
23 minAssZ(I);        compute the minimal associated primes of I
25 heightZ(I);        compute the height of I
26 equidimZ(I);       compute the equidimensional part of I
27 intersectZ(I,J)    compute the intersection of I and J
28";
29
30LIB "primdec.lib";
31
32////////////////////////////////////////////////////////////////////////////////
33
34proc primdecZ(ideal I, list #)
35"USAGE:  primdecZ(I[, n]); I ideal, n integer (number of processors)
36NOTE:    If size(#) > 0, then #[1] is the number of available processors for
37         the computation.
38RETURN:  a list pr of primary ideals and their associated primes:
39@format
40   pr[i][1]   the i-th primary component,
41   pr[i][2]   the i-th prime component.
42@end format
43EXAMPLE: example primdecZ; shows an example
44"
45{
46   if(size(I)==0){return(list(list(ideal(0),ideal(0))));}
47
48//--------------------  Initialize optional parameters  ------------------------
49   if(size(#) > 0)
50   {
51      if(size(#) == 1)
52      {
53         int n = #[1];
54         ideal TES = 1;
55      }
56      if(size(#) == 2)
57      {
58         int n = #[1];
59         ideal TES = #[2];
60      }
61   }
62   else
63   {
64      int n = 1;
65      ideal TES = 1;
66   }
67
68
69   if(size(I)==1 && deg(I[1]) == 0)
70   {
71      ideal J = I;
72   }
73   else
74   {
75      ideal J = stdZ(I);
76   }
77
78   ideal K,N;
79   def R=basering;
80   number s;
81   list rl=ringlist(R);
83   list P,B,IS;
84   ideal Q,JJ;
85   ideal TQ=1;
86   if(deg(J[1])==0)
87   {
88      //=== I intersected with Z is not zero
89      list rp=rl;
90      rp[1]=0;
91      //=== q is generator of I intersect Z
93      def Rhelp=ring(rp);
94      setring Rhelp;
95      number q=imap(R,q);
96      //=== computes the primes occuring in a generator of I intersect Z
97      list L = primefactors(q);
98
99      list A;
100      A[1]= ideal(0);
101      ideal J = imap(R,J);
102
103      for(j=1;j<=size(L[2]);j++)
104      {
105         if(L[2][j] > 1){ ex = 1; break; }
106      }
107
108      if(printlevel >= 10)
109      {
110         "n = "+string(n);
111         "size(L[2]) = "+string(size(L[2]));
112      }
113
114      int RT = rtimer;
115      if((n > 1) && (n < size(L[2])))
116      {
117
118//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
119//-----  standard basis for the primes L[1][2],...,L[1][n + 1].        ---------
120
121         for(i = 1; i <= n; i++)
122         {
123            p=int(L[1][i + 1]);
124            nu=int(L[2][i + 1]);
127            open(l(i));
128            write(l(i), quote(modp(eval(J), eval(p), eval(nu))));
129         }
130
131         p = int(L[1][1]);
132         nu = int(L[2][1]);
133         int t = timer;
134         A[size(A)+1] = modp(J, p, nu);
135         t = timer - t;
136         if(t > 60) { t = 60; }
137         int i_sleep = system("sh", "sleep "+string(t));
138
139         j = n + 2;
140
141         while(j <= size(L[2]) + 1)
142         {
143            for(i = 1; i <= n; i++)
144            {
147               {
148                  //=== read the result from l(i)
150
151                  if(j <= size(L[2]))
152                  {
153                     p=int(L[1][j]);
154                     nu=int(L[2][j]);
155                     write(l(i), quote(modp(eval(J), eval(p), eval(nu))));
156                     j++;
157                  }
158                  else
159                  {
161                     close(l(i));
162                  }
163               }
164            }
167            {
168               j++;
169            }
170            i_sleep = system("sh", "sleep "+string(t));
171         }
172      }
173      else
174      {
175         for(j=1;j<=size(L[2]);j++)
176         {
177            p=int(L[1][j]);
178            nu=int(L[2][j]);
179            A[size(A)+1] = modp(J, p,nu );
180         }
181      }
182
183      setring R;
184      list A = imap(Rhelp,A);
185      if(printlevel >= 10)
186      {
187         "A is computed in "+string(rtimer - RT)+" seconds.";
188      }
189      for(i=2;i<=size(A);i++)
190      {
191      //=== computes for all p in L the minimal associated primes of
192      //=== IZ/p[variables]
193         p = int(A[i][2]);
194         if(printlevel >= 10)
195         {
196            "p = "+string(p);
197            RT = rtimer;
198         }
199         nu = int(A[i][3]);
200         //=== maximal power of p dividing q, generator of I intersect Z
201         s = p^nu;
202
203         rp[1] = p;
204         def S = ring(rp);
205         setring S;
206         ideal J = imap(R,J);
207         setring R;
208
209         if(nu>1)
210         {
211            //=== p is of multiplicity > 1 in q
212
213            B = A[i][1];
214            for(j=1;j<=size(B);j++)
215            {
216               //=== the minimal associated primes of I
217               K=B[j],p;
218               K=stdZ(K);
219               B[j]=K;
220            }
221            for(j=1;j<=size(B);j++)
222            {
223               K=B[j];
224               //=== compute maximal independent set for KZ/p[variables]
225
226               setring S;
227               J=imap(R,K);
228               J=simplify(J,2);
229               attrib(J,"isSB",1);
230               IS=Primdec::maxIndependSet(J);
231               setring R;
232               //=== computing the pseudo primary and extract it
233               N=J,s;
234               N=stdZ(N);
235               Q=extractZ(N,j,IS,B);
236               //=== test for useless primaries
237               if(size(reduce(TES,Q))>0)
238               {
239                  TQ=intersectZ(TQ,Q);
240                  //TQ=intersect(TQ,Q);
241                  P[size(P)+1]=list(Q,K);
242               }
243            }
244         }
245         else
246         {
247            //=== p is of multiplicity 1 in q we can compute the
248            //=== primary decomposition directly
249
250            B = A[i][1];
251            for(j=1;j<=size(B);j++)
252            {
253               K=B[j][2],p;
254               K=stdZ(K);
255               Q=B[j][1],p;
256               Q=stdZ(Q);
257               if(size(reduce(TES,Q))>0)
258               {
259                  //TQ=intersectZ(TQ,Q);
260                  P[size(P)+1]=list(Q,K);
261               }
262            }
263            if(ex)
264            {
265               JJ=imap(S,J);
266               JJ=JJ,p;
267               JJ=stdZ(JJ);
268               TQ=intersectZ(TQ,JJ);
269               //TQ=intersect(TQ,JJ);
270            }
271         }
272         kill S;
273         if(printlevel >= 10)
274         {
275            string(p)+" done in "+string(rtimer - RT)+" seconds.";
276         }
277      }
278
279      setring R;
280      if(!ex){return(P);}
281      J=stdZ(J);
282      TQ=intersectZ(TQ,TES);
283      //TQ=intersect(TQ,TES);
284      if(size(reduce(TQ,J))!=0)
285      {
286         //=== taking care about embedded components
287         K=stdZ(quotientZ(J,TQ));
288         ideal W=K;
289         m++;
290         while(size(reduce(intersectZ(W,TQ),J))!=0)
291         //while(size(reduce(intersect(W,TQ),J))!=0)
292         {
293            W=stdZ(I+specialPowerZ(K,m));
294            m++;
295         }
296         list E=primdecZ(W,n,TQ);
297         for(i=1;i<=size(E);i++)
298         {
299            P[size(P)+1]=E[i];
300         }
301      }
302      return(P);
303   }
304
305   //==== the ideal intersected with Z is zero
306   rl[1]=0;
307   def Rhelp=ring(rl);
308   setring Rhelp;
309   ideal J=imap(R,J);
310   J=std(J);
311   //=== the primary decomposition over Q which gives the primary
312   //=== decomposition of I:h for a suitable integer h
313   list pr=primdecGTZ(J);
314   for(i=1;i<=size(pr);i++)
315   {
316      pr[i]=list(std(pr[i][1]),std(pr[i][2]));
317   }
318   setring R;
319   list pr=imap(Rhelp,pr);
320   //=== intersection with Z[variables]
321   for(i=1;i<=size(pr);i++)
322   {
323      pr[i]=list(coefZ(pr[i][1])[1],coefZ(pr[i][2])[1]);
324   }
325   //=== find h in Z such that I is the intersection of I:h and <I,h>
326   //=== and I:h = IQ[variables] intersected with Z[varables]
327   list H =coefZ(J);
328   ideal Y=H[1];
329   int h=H[2];
330   J=J,h;
331   //=== call primary decomposition over Z for <I,h>
332   ideal testid = 1;
333   for(i = 1; i<=size(pr); i++)
334   {
335      testid = intersect(testid, pr[i][1]);
336   }
337   list M;
338   if(h!=1)
339   {
340       M=primdecZ(J,n,Y);
341       j=0;
342       //=== remove useless primary ideals
343       while(j<size(M))
344       {
345          j++;
346          M[j][1]=stdZ(M[j][1]);
347          testid = intersect(testid, M[j][1]);
348          if((reduce(testid, std(I))+0 == 0) && (reduce(I, std(testid))+0 == 0))
349          {
350            for(i=1;i<=j;i++)
351            {
352              pr[size(pr)+1]=M[i];
353            }
354            return(pr);
355          }
356          for(i=1;i<=size(pr);i++)
357          {
358             if(size(reduce(pr[i][1],M[j][1]))==0)
359             {
360                M=delete(M,j);
361                j--;
362                break;
363             }
364          }
365      }
366      for(i=1;i<=size(M);i++)
367      {
368         pr[size(pr)+1]=M[i];
369      }
370   }
371   return(pr);
372}
373example
374{ "EXAMPLE:";  echo = 2;
375   ring R=integer,(a,b,c,d),dp;
376   ideal I1=9,a,b;
377   ideal I2=3,c;
378   ideal I3=11,2a,7b;
379   ideal I4=13a2,17b4;
380   ideal I5=9c5,6d5;
381   ideal I6=17,a15,b15,c15,d15;
382   ideal I=intersectZ(I1,I2);
383   I=intersectZ(I,I3);
384   I=intersectZ(I,I4);
385   I=intersectZ(I,I5);
386   I=intersectZ(I,I6);
387   primdecZ(I);
388   ideal J=intersectZ(ideal(17,a),ideal(17,a2,b));
389   primdecZ(J);
390   ideal K=intersectZ(ideal(9,a+3),ideal(9,b+3));
391   primdecZ(K);
392}
393
394////////////////////////////////////////////////////////////////////////////////
395
396proc minAssZ(ideal I)
397"USAGE:  minAssZ(I); I ideal
398RETURN:  a list pr of associated primes:
399EXAMPLE: example minAssZ; shows an example
400"
401{
402   if(size(I)==0){return(list(ideal(0)));}
403   if(deg(I[1])==0)
404   {
405      ideal J=I;
406   }
407   else
408   {
409      ideal J=stdZ(I);
410   }
411   ideal K;
412   def R=basering;
413   list rl=ringlist(R);
414   int i,j,p,m;
415   list P,B;
416   if(deg(J[1])==0)
417   {
418      //=== I intersected with Z is not zero
419      list rp=rl;
420      rp[1]=0;
422      def Rhelp=ring(rp);
423      setring Rhelp;
424      number q=imap(R,q);
425      //=== computes the primes occuring in a generator of I intersect Z
426      //list L=PollardRho(q,5000,1);
427      list L=primefactors(q)[1];
428      for(i=1;i<=size(L);i++)
429      {
430      //=== computes for all p in L the minimal associated primes of
431      //=== IZ/p[variables]
432         p=int(L[i]);
433         setring R;
434         rp[1]=p;
435         def S=ring(rp);
436         setring S;
437         ideal J=imap(R,J);
438         list A=minAssGTZ(J);
439         setring R;
440         B=imap(S,A);
441         kill S;
442         for(j=1;j<=size(B);j++)
443         {
444            //=== the minimal associated primes of I
445            if(B[j][1]!=1)
446            {
447               K=B[j],p;
448               K=stdZ(K);
449               P[size(P)+1]=K;
450            }
451         }
452         setring Rhelp;
453      }
454      setring R;
455      return(P);
456   }
457   //==== the ideal intersected with Z is zero
458   rl[1]=0;
459   def Rhelp=ring(rl);
460   setring Rhelp;
461   ideal J=imap(R,J);
462   J=std(J);
463   //=== the primary decomposition over Q which gives the primary
464   //=== decomposition of I:h for a suitable integer h
465   list pr=minAssGTZ(J);
466   for(i=1;i<=size(pr);i++)
467   {
468      pr[i]=std(pr[i]);
469   }
470   setring R;
471   list pr=imap(Rhelp,pr);
472   //=== intersection with Z[variables]
473   for(i=1;i<=size(pr);i++)
474   {
475      pr[i]=coefZ(pr[i])[1];
476   }
477   //=== find h in Z such that I is the intersection of I:h and I,h
478   //=== and I:h =IQ[variables] intersected with Z[varables]
479   list H=coefZ(J);
480   int h=H[2];
481   J=J,h;
482   //=== call associated primes over Z for I,h
483   list M;
484   if(h!=1)
485   {
486       M=minAssZ(J);
487       //=== remove non-minimal primes
488       j=0;
489       while(j<size(M))
490       {
491          j++;
492          M[j]=stdZ(M[j]);
493          for(i=1;i<=size(pr);i++)
494          {
495             if(size(reduce(pr[i],M[j]))==0)
496             {
497                M=delete(M,j);
498                j--;
499                break;
500             }
501          }
502      }
503      for(i=1;i<=size(M);i++)
504      {
505         pr[size(pr)+1]=M[i];
506      }
507   }
508   return(pr);
509}
510example
511{ "EXAMPLE:";  echo = 2;
512   ring R=integer,(a,b,c,d),dp;
513   ideal I1=9,a,b;
514   ideal I2=3,c;
515   ideal I3=11,2a,7b;
516   ideal I4=13a2,17b4;
517   ideal I5=9c5,6d5;
518   ideal I6=17,a15,b15,c15,d15;
519   ideal I=intersectZ(I1,I2);
520   I=intersectZ(I,I3);
521   I=intersectZ(I,I4);
522   I=intersectZ(I,I5);
523   I=intersectZ(I,I6);
524   minAssZ(I);
525   ideal J=intersectZ(ideal(17,a),ideal(17,a2,b));
526   minAssZ(J);
527   ideal K=intersectZ(ideal(9,a+3),ideal(9,b+3));
528   minAssZ(K);
529}
530
531////////////////////////////////////////////////////////////////////////////////
532
533proc heightZ(ideal I)
534"USAGE:  heightZ(I); I ideal
535RETURN:  the height of the input ideal
536EXAMPLE: example heightZ; shows an example
537"
538{
539   if(size(I)==0){return(0);}
540   if(deg(I[1])==0)
541   {
542      ideal J=I;
543   }
544   else
545   {
546      ideal J=stdZ(I);
547   }
548   ideal K=1;
549   def R=basering;
550   list rl=ringlist(R);
551   int i,j,p,m;
552   list P;
553   ideal B;
554   if(deg(J[1])==0)
555   {
556      //=== I intersected with Z is not zero
557      m=nvars(R);
558      list rp=rl;
559      rp[1]=0;
561      def Rhelp=ring(rp);
562      setring Rhelp;
563      number q=imap(R,q);
564      //=== computes the primes occuring in a generator of I intersect Z
565      //list L=PollardRho(q,5000,1);
566      list L=primefactors(q)[1];
567      for(i=1;i<=size(L);i++)
568      {
569      //=== computes for all p in L the std of IZ/p[variables]
570         p=int(L[i]);
571         setring R;
572         rp[1]=p;
573         def S=ring(rp);
574         setring S;
575         ideal J=imap(R,J);
576         j=nvars(R)-dim(std(J));
577         if(j<m){m=j;}
578         setring Rhelp;
579         kill S;
580      }
581      setring R;
582      return(m+1);
583   }
584   //==== the ideal intersected with Z is zero
585   rl[1]=0;
586   def Rhelp=ring(rl);
587   setring Rhelp;
588   ideal J=imap(R,J);
589   J=std(J);
590   m=nvars(R)-dim(J);
591   //=== the height over Q
592   //=== of I:h for a suitable integer h
593   setring R;
594   //=== find h in Z such that I is the intersection of I:h and I,h
595   //=== and I:h =IQ[variables] intersected with Z[varables]
596   list H=coefZ(J);
597   int h=H[2];
598   J=J,h;
599   //=== call height over Z for I,h
600   if(h!=1)
601   {
602      j=heightZ(J);
603      if(j<m){m=j;}
604   }
605   return(m);
606}
607example
608{ "EXAMPLE:";  echo = 2;
609   ring R=integer,(a,b,c,d),dp;
610   ideal I1=9,a,b;
611   ideal I2=3,c;
612   ideal I3=11,2a,7b;
613   ideal I4=13a2,17b4;
614   ideal I5=9c5,6d5;
615   ideal I6=17,a15,b15,c15,d15;
616   ideal I=intersectZ(I1,I2);
617   I=intersectZ(I,I3);
618   I=intersectZ(I,I4);
619   I=intersectZ(I,I5);
620   I=intersectZ(I,I6);
621   heightZ(I);
622}
623
624////////////////////////////////////////////////////////////////////////////////
625
628RETURN:  the radcal of the input ideal
629EXAMPLE: example radicalZ; shows an example
630"
631{
632   if(size(I)==0){return(ideal(0));}
633   if(deg(I[1])==0)
634   {
635      ideal J=I;
636   }
637   else
638   {
639      ideal J=stdZ(I);
640   }
641   ideal K=1;
642   def R=basering;
643   list rl=ringlist(R);
644   int i,j,p,m;
645   list P;
646   ideal B;
647   if(deg(J[1])==0)
648   {
649      //=== I intersected with Z is not zero
650      list rp=rl;
651      rp[1]=0;
653      def Rhelp=ring(rp);
654      setring Rhelp;
655      number q=imap(R,q);
656      //=== computes the primes occuring in a generator of I intersect Z
657      //list L=PollardRho(q,5000,1);
658      list L=primefactors(q)[1];
659      for(i=1;i<=size(L);i++)
660      {
661      //=== computes for all p in L the radical of IZ/p[variables]
662         p=int(L[i]);
663         setring R;
664         rp[1]=p;
665         def S=ring(rp);
666         setring S;
667         ideal J=imap(R,J);
669         setring R;
670         B=imap(S,A);
671         kill S;
672         B=B,p;
673         B=stdZ(B);
674         K=stdZ(intersectZ(K,B));
675         //K=stdZ(intersect(K,B));
676         setring Rhelp;
677      }
678      setring R;
679      return(K);
680   }
681   //==== the ideal intersected with Z is zero
682   rl[1]=0;
683   def Rhelp=ring(rl);
684   setring Rhelp;
685   ideal J=imap(R,J);
686   J=std(J);
688   //=== of I:h for a suitable integer h
690   setring R;
691   K=imap(Rhelp,K);
692   //=== intersection with Z[variables]
693   K=coefZ(K)[1];
694   //=== find h in Z such that I is the intersection of I:h and I,h
695   //=== and I:h =IQ[variables] intersected with Z[varables]
696   list H=coefZ(J);
697   int h=H[2];
698   J=J,h;
699   //=== call radical over Z for I,h
700   if(h!=1)
701   {
703      K=intersectZ(K,M);
704      //K=intersect(K,M);
705   }
706   return(K);
707}
708example
709{ "EXAMPLE:";  echo = 2;
710   ring R=integer,(a,b,c,d),dp;
711   ideal I1=9,a,b;
712   ideal I2=3,c;
713   ideal I3=11,2a,7b;
714   ideal I4=13a2,17b4;
715   ideal I5=9c5,6d5;
716   ideal I6=17,a15,b15,c15,d15;
717   ideal I=intersectZ(I1,I2);
718   I=intersectZ(I,I3);
719   I=intersectZ(I,I4);
720   I=intersectZ(I,I5);
721   I=intersectZ(I,I6);
723   ideal J=intersectZ(ideal(17,a),ideal(17,a2,b));
725}
726
727////////////////////////////////////////////////////////////////////////////////
728
729proc equidimZ(ideal I)
730"USAGE:  equidimZ(I); I ideal
731RETURN:  the part of minimal height
732EXAMPLE: example equidimZ; shows an example
733"
734{
735   if(size(I)==0){return(ideal(0));}
736   if(deg(I[1])==0)
737   {
738      ideal J=I;
739   }
740   else
741   {
742      ideal J=stdZ(I);
743   }
744   int he=heightZ(J);
745   ideal K,N;
746   def R=basering;
747   bigint s;
748   list rl=ringlist(R);
749   int i,j,p,m,ex;
750   list P,IS,B;
751   ideal Q,JJ,E;
752   ideal TQ=1;
753   if(deg(J[1])==0)
754   {
755      //=== I intersected with Z is not zero
756      list rp=rl;
757      rp[1]=0;
758      //=== generator of I intersect Z
760      def Rhelp=ring(rp);
761      setring Rhelp;
762      number q=imap(R,q);
763      //=== computes the primes occuring in a generator of I intersect Z
764      //list L=PollardRho(q,5000,1);
765      list L=primefactors(q)[1];
766      list Le;
767      for(i=1;i<=size(L);i++)
768      {
769         L[i]=int(L[i]);
770         p=int(L[i]);
771         j=0;
772         s=bigint(q);
773         while((s mod p)==0)
774         {
775            j++;
776            s=s div p;
777         }
778         Le[i]=j;
779      }
780      for(i=1;i<=size(L);i++)
781      {
782      //=== computes for all p in L the minimal associated primes of
783      //=== IZ/p[variables]
784         p=int(L[i]);
785         j=Le[i];
786         setring R;
787         //=== maximal power of p dividing q, generator of I intersect Z
788         s=p^j;
789         rp[1]=p;
790         def S=ring(rp);
791         setring S;
792         ideal J=imap(R,J);
793         J=std(J);
794         if(nvars(R)-dim(J)+1==he)
795         {
796            if(j>1)
797            {
798               //=== p is of multiplicity >1 in q
799               list A=minAssGTZ(J);
800               j=0;
801               while(j<size(A))
802               {
803                  j++;
804                  if(dim(std(A[j]))!=nvars(R)-he+1)
805                  {
806                     A=delete(A,j);
807                     j--;
808                  }
809               }
810               setring R;
811               B=imap(S,A);
812               for(j=1;j<=size(B);j++)
813               {
814                  //=== the minimal associated primes of I
815                  K=B[j],p;
816                  K=stdZ(K);
817                  B[j]=K;
818               }
819               for(j=1;j<=size(B);j++)
820               {
821                  K=B[j];
822                  //=== compute maximal independent set for KZ/p[variables]
823                  setring S;
824                  J=imap(R,K);
825                  J=simplify(J,2);
826                  attrib(J,"isSB",1);
827                  IS=Primdec::maxIndependSet(J);
828                  setring R;
829                  //=== computing the pseudo primary and extract it
830                  N=J,s;
831                  N=stdZ(N);
832                  Q=extractZ(N,j,IS,B);
833                  TQ=intersectZ(TQ,Q);
834                  //TQ=intersect(TQ,Q);
835               }
836               setring Rhelp;
837            }
838            else
839            {
840               //=== p is of multiplicity 1 in q we can compute the
841               //=== equidimensional part directly
842               ideal E=equidimMax(J);
843               setring R;
844               E=imap(S,E);
845               E=E,p;
846               E=stdZ(E);
847               TQ=intersectZ(TQ,E);
848               //TQ=intersect(TQ,E);
849            }
850         }
851         kill S;
852         setring Rhelp;
853      }
854      setring R;
855      return(TQ);
856   }
857   //==== the ideal intersected with Z is zero
858   rl[1]=0;
859   def Rhelp=ring(rl);
860   setring Rhelp;
861   ideal J=imap(R,J);
862   J=std(J);
863   //=== the equidimensional part over Q which gives the equdimensional
864   //=== part of I:h for a suitable integer h
865   ideal E=1;
866   if(nvars(R)-he==dim(J))
867   {
868      E=std(equidimMax(J));
869   }
870   setring R;
871   E =imap(Rhelp,E);
872   //=== intersection with Z[variables]
873   E=coefZ(E)[1];
874   //=== find h in Z such that I is the intersection of I:h and I,h
875   //=== and I:h =IQ[variables] intersected with Z[varables]
876   int h =coefZ(J)[2];
877   J=J,h;
878   //=== call equidimensional part over Z for I,h
879   ideal M;
880   if(h!=1)
881   {
882       M=equidimZ(J);
883       if(he==heightZ(M))
884       {
885          E=intersectZ(M,E);
886          //E=intersect(M,E);
887       }
888   }
889   return(E);
890}
891example
892{ "EXAMPLE:";  echo = 2;
893   ring R=integer,(a,b,c,d),dp;
894   ideal I1=9,a,b;
895   ideal I2=3,c;
896   ideal I3=11,2a,7b;
897   ideal I4=13a2,17b4;
898   ideal I5=9c5,6d5;
899   ideal I6=17,a15,b15,c15,d15;
900   ideal I=intersectZ(I1,I2);
901   I=intersectZ(I,I3);
902   I=intersectZ(I,I4);
903   I=intersectZ(I,I5);
904   I=intersectZ(I,I6);
905   equidimZ(I);
906}
907
908////////////////////////////////////////////////////////////////////////////////
909
910proc intersectZ(ideal I, ideal J)
911"USAGE:  intersectZ(I,J); I,J ideals
912RETURN:  the intersection of the input ideals
913NOTE:    this is an alternative to intersect(I,J) over integers,
914         is faster for some examples and should be kept for debug purposes.
915EXAMPLE: example intersectZ; shows an example
916{
917   def R = basering;
918   execute("ring S=integer,( X(1..nvars(R)) ), ( dp(nvars(R)) );");
919   ideal I = fetch(R,I);
920   ideal J = fetch(R,J);
921   execute("ring St=integer, ( t, X(1..nvars(R))  ), ( dp(1), dp(nvars(R)) );");
922   ideal I = imap(S,I);
923   ideal J = imap(S,J);
924   ideal K =  var(1)*I+(1-var(1))*J;
925   K = stdZ(K);
926   int i;
927   ideal L;
928   for(i=1; i<=size(K); i++)
929   {
930      if ( lead(K[i])/var(1)==0 ) { L[size(L)+1] = K[i]; }
931   }
932   setring S;
933   ideal L = imap(St, L);
934   setring R;
935   ideal L = fetch(S, L);
936   return(L);
937}
938example
939{ "EXAMPLE:";  echo = 2;
940   ring R=integer,(a,b,c,d),dp;
941   ideal I1=9,a,b;
942   ideal I2=3,c;
943   ideal I3=11,2a,7b;
944   ideal I4=13a2,17b4;
945   ideal I5=9c5,6d5;
946   ideal I6=17,a15,b15,c15,d15;
947   ideal I=intersectZ(I1,I2); I;
948   I=intersectZ(I,I3); I;
949   I=intersectZ(I,I4); I;
950   I=intersectZ(I,I5); I;
951   I=intersectZ(I,I6); I;
952}
953
954////////////////////////////////////////////////////////////////////////////////
955
956static proc modp(ideal J, int p, int nu)
957{
958//=== computes the minimal associated primes (if nu > 1) resp. the primary
959//=== decomposition (else) of J in Z/p and maps the result back to the basering
960   def R = basering;
961   list rp = ringlist(R);
962   rp[1] = p;
963   def Rp = ring(rp);
964   setring Rp;
965   ideal J = imap(R,J);
966   int sizeA;
967   if(nu > 1)
968   {
969      //=== p is of multiplicity > 1 in q
970      list A = minAssGTZ(J);
971   }
972   else
973   {
974      list A = primdecGTZ(J);
975   }
976   sizeA = size(A);
977   setring R;
978   list A;
979   if (sizeA>0)   {  A = imap(Rp,A);    }
980   return(list(A,p,nu));
981}
982
983////////////////////////////////////////////////////////////////////////////////
984
985static proc coefPrimeZ(ideal I)
986{
987//=== computes the primes occuring in the product of the leading coefficients
988//=== of I
989   number h=1;
990   int i;
991   I = simplify(I,2); // del zero generators
992   for(i=1;i<=size(I);i++)
993   {
994      h=h*leadcoef(I[i]);  // besser machen (gleich zerlegen,
995                           // nicht ausmultiplizieren)
996   }
997   def R=basering;
998   ring Rhelp=0,x,dp;
999   number h=imap(R,h);
1000   //list L=PollardRho(h,5000,1);
1001   list L=primefactors(h)[1];
1002   for(i=1;i<=size(L);i++){L[i]=int(L[i]);}
1003   setring R;
1004   return(L);
1005}
1006
1007////////////////////////////////////////////////////////////////////////////////
1008
1009static proc coefZ(ideal I)
1010{
1011//=== assume IQ[variables]=<g_1,...,g_s>, Groebner basis, g_i in Z[variables]
1012//=== computes an integer h such that
1013//===   <g_1,...,g_s>Z[variables]:h^infinity = IQ[variables] intersected
1014//===                                          with Z[variables]
1015//=== returns a list with IQ[variables] intersected with Z[variables] and h
1016   int h=1;
1017   int i,e;
1018   ideal K=1;
1019   attrib(I,"isSB",1);
1020   list L=coefPrimeZ(I);
1021   if(size(L)==0){return(list(I,1));}
1022   int d=1;
1023   while(d!=0)
1024   {
1025      i++;
1026      K=quotientOneZ(I,L[i]);
1027      if(size(reduce(K,I))!=0)
1028      {
1029         h=h*L[i];
1030         I=stdZ(K);
1031         e=1;
1032      }
1033      if(i==size(L))
1034      {
1035         i=0;
1036         if(e)
1037         {
1038            e=0;
1039         }
1040         else
1041         {
1042           d=0;
1043         }
1044      }
1045   }
1046   if(h<0){h=-h;}
1047   return(list(K,h));
1048}
1049
1050////////////////////////////////////////////////////////////////////////////////
1051
1052static proc specialPowerZ(ideal I, int m)
1053{
1054//=== computes the ideal generated by the m-th power of the generators of I
1055   int i;
1056   for(i=1;i<=ncols(I);i++)  // use ncols to allow zero generators in I
1057   {
1058      I[i]=I[i]^m;
1059   }
1060   return(I);
1061}
1062
1063////////////////////////////////////////////////////////////////////////////////
1064
1065static proc separatorsZ(int j, list B)
1066{
1067//=== computes s such that s is not in B[j] but s is in B[i] for all i!=j
1068   int i,k;
1069   poly s=1;
1070   for(i=1;i<=size(B);i++)
1071   {
1072      if(i!=j)
1073      {
1074         for(k=1;k<=size(B[i]);k++)
1075         {
1076            if(reduce(B[i][k],B[j])!=0)
1077            {
1078               s=s*B[i][k];
1079               break;
1080            }
1081         }
1082      }
1083   }
1084   return(s);
1085}
1086
1087//////////////////////////////////////////////////////////////////////////////
1088static proc extractZ(ideal J, int j, list L, list B)
1089{
1090   //=== P is an associated prime of J, the corresponding primary ideal is
1091   //=== computed,
1092   //=== L is a list of maximal independent sets for P in Z/p[variables]
1093   def R=basering;
1094   ideal P=B[j];
1095
1096   //=== first compute a pseudo primary ideal I, radical of I is P
1097   //=== method of Eisenbud
1098   //ideal I=J+specialPowerZ(P,20);
1099
1100   //=== method of Shimoyama-Yokoyama
1101   poly s=separatorsZ(j,B);
1102   ideal I=satZ(J,s);
1103   //=== size(L)=0 means P is maximal ideal and I is primary
1104   if(size(L)>0)
1105   {
1106      if(L[1][3]!=0)
1107      {
1108         //=== if u in x is an independent set of L then we compute a Groebner
1109         //=== Basis in Z[u][x-u]
1110         execute("ring S=integer,("+L[1][1]+"),lp;");
1111         ideal I=imap(R,I);
1112         I=stdZ(I);
1113         list rl=ringlist(S);
1114         rl[1]=0;
1115         def Shelp =ring(rl);
1116         setring Shelp;
1117         ideal I=imap(S,I);
1118         I[1]=0;
1119         I=simplify(I,2);
1120         if(L[1][3]==nvars(basering))
1121         {
1122            list C;
1123            C[1]= ideal(0); //add dummy entry to tie C to current ring
1124            int i;
1125            ASSUME(1, size(I)==ncols(I) || size(I)==0);
1126            for(i=1;i<=size(I);i++)
1127            {
1128               C[i+1]=I[i];
1129            }
1130         }
1131         else
1132         {
1133            //=== this is our way to obtain the coefficients in Z[u] of the
1134            //=== leading terms of the Groebner basis above
1135            def quring=Primdec::prepareQuotientring(nvars(basering)-L[1][3],"lp");
1136            setring quring;
1137            ideal I=imap(Shelp,I);
1138            list C;
1139            C[1]= ideal(0); //add dummy entry to tie C to current ring
1140            int i;
1141            ASSUME(1, size(I)==ncols(I) || size(I)==0);
1142            for(i=1;i<=size(I);i++)
1143            {
1145            }
1146            setring Shelp;
1147            list C=imap(quring,C);
1148         }
1149         setring R;
1150         list C=imap(Shelp,C);
1151      }
1152      else
1153      {
1154         I=stdZ(I);
1155         list C;
1156         C[1]= ideal(0); //add dummy entry to tie C to current ring
1157         int i;
1158         ASSUME(1, size(I)==ncols(I) || size(I)==0);
1159         for(i=1;i<=size(I);i++)
1160         {
1161            C[i+1]=I[i];
1162         }
1163         list rl=ringlist(R);
1164         rl[1]=0;
1165         def Shelp =ring(rl);
1166      }
1167      poly h=1;
1168      for(i=2;i<=size(C);i++) // leave out first dummy entry in C.
1169      {
1170         if(deg(C[i])>0){h=h*C[i];}  // das muss noch besser gemacht werden,
1171                                     // nicht ausmultiplizieren!
1172      }
1173      setring Shelp;
1174      poly h=imap(R,h);
1175      ideal fac=factorize(h,1);
1176      setring R;
1177      ideal fac=imap(Shelp,fac);
1178      ASSUME(1, size(fac)==ncols(fac) || size(fac)==0);
1179      for(i=1;i<=size(fac);i++)
1180      {
1181         I=satZ(I,fac[i]);
1182      }
1183   }
1184   I=stdZ(I);
1185   return(I);
1186}
1187////////////////////////////////////////////////////////////////////////////////
1188
1189static proc normalizeZ(ideal I)
1190{
1191//=== if I[1]=q in Z, it replaces all other coeffs of polys in I by there value
1192//=== mod q, std should do this automatically and then this procedure should be
1193//=== removed
1194   if(deg(I[1])>0){return(I);}
1195   int i,j;
1196   number n;
1197   poly p;
1198   for(i=2;i<=ncols(I);i++)
1199   {
1200      j=1;
1201      while(j<=size(I[i]))
1202      {
1205         I[i]=I[i]-I[i][j]+p;
1206         if(p!=0){j++;}
1207      }
1208   }
1209   return(I);
1210}
1211
1212////////////////////////////////////////////////////////////////////////////////
1213
1214static proc satZ(ideal I,poly h)
1215{
1216//=== saturates I by h
1217   ideal J=quotientOneZ(I,h);
1218   while(size(reduce(J,stdZ(I)))!=0)
1219   {
1220      I=J;
1221      J=quotientOneZ(I,h);
1222      J=normalizeZ(J);
1223   }
1224   return(J);
1225}
1226
1227////////////////////////////////////////////////////////////////////////////////
1228
1229static proc quotientOneZ(ideal I, poly f)
1230{
1231//=== this is needed because quotient(I,f) does not work properly, should be
1232//=== replaced by quotient later
1233   if ( f==0 ) { return( ideal(1) ); }
1234   def R=basering;
1235   int i;
1236   ideal K=intersectZ(I,ideal(f));
1237   //ideal K=intersect(I,ideal(f));
1238   //=== K[i]/f; does not work in rings with integer! This should be replaced
1239   //=== later
1240   execute("ring Rhelp=0,("+varstr(R)+"),dp;");
1241   ideal K=imap(R,K);
1242   poly f=imap(R,f);
1243   ASSUME(1, ncols(K)==size(K) || size(K)==0 ); // postcondition intersectZ
1244   for(i=1;i<=ncols(K);i++)
1245   {
1246      K[i]=K[i]/f;
1247   }
1248   setring R;
1249   K=imap(Rhelp,K);
1250   return(K);
1251}
1252
1253////////////////////////////////////////////////////////////////////////////////
1254
1255static proc quotientZ(ideal I, ideal J)
1256{
1257//=== this is needed because quotient(I,J) does not work properly, should be
1258//=== replaced by quotient later
1259   if ( size(J)==0 ) { return( ideal(1) ); }
1260   J = simplify(J, 2); // del zero generators
1261   int i;
1262   ideal K=quotientOneZ(I,J[1]);
1263   // we had a bug here: either it is mandatory to kick out zero generators  of J(done) or to use ncols(J)
1264   for(i=2;i<=ncols(J);i++)
1265   {
1266      K=intersectZ(K,quotientOneZ(I,J[i]));
1267      //K=intersect(K,quotientOneZ(I,J[i]));
1268   }
1269   return(K);
1270}
1271
1272////////////////////////////////////////////////////////////////////////////////
1273
1274static proc reduceZ(poly f, ideal I)
1275{
1276//=== this is needed because reduce(f,I) does not work properly, should be
1277//=== replaced by reduce later
1278   if(f==0){return(f);}
1279   if (not attrib(I,"isSB") ) { print ("// ** I is no standard basis"); }
1280   def R=basering;
1281   execute("ring Rhelp=0,("+varstr(R)+"),dp;");
1282   ideal I=imap(R,I);
1283   poly f=imap(R,f);
1284   int i,j;
1285   poly m;
1286   number n;
1287   I = simplify(I, 2); // del zero gens (otherwise the following loop may fail with an error (div by 0))
1288   while(!i)
1289   {
1290      i=1;
1291      j=0;
1292      while(j<size(I))
1293      {
1294         j++;
1296         if(m!=0)
1297         {
1299            if(n==0)
1300            {
1302               if(f==0){setring R;return(poly(0));}
1303               i=0;
1304               break;
1305            }
1307            {
1309               i=0;
1310               break;
1311            }
1312         }
1313      }
1314   }
1315   setring R;
1316   f=imap(Rhelp,f);
1318}
1319
1320////////////////////////////////////////////////////////////////////////////////
1321
1322static proc stdZ(ideal I)
1323{
1324//=== this is needed because we want the leading coefficients to be positive
1325//=== otherwhise reduce gives wrong results! should be replaced later by std
1326   I=simplify(I,2);
1327   I=normalizeZ(I);  // why is this done before std() call?
1328                     // normalizeZ works only if I[1] == ( std(I) )[1]
1329   ideal J=std(I);
1330   int i;
1331   for(i=1;i<=size(J);i++)
1332   {
1334   }
1335   J=normalizeZ(J);
1336   attrib(J,"isSB",1);
1337   return(J);
1338}
1339
1340////////////////////////////////////////////////////////////////////////////////
1341
1342static proc testPrimaryZ(ideal I, list L)
1343{
1344//=== test whether I is the intersection of the primary ideals in L
1345   int i;
1346   ideal K=L[1][1];
1347   for(i=2;i<=size(L);i++)
1348   {
1349      K=intersectZ(K,L[i][1]);
1350      //K=intersect(K,L[i][1]);
1351   }
1352   i=size(reduce(K,stdZ(I)))+size(reduce(I,stdZ(K)));
1353   if(!i){return(1);}
1354   return(0);
1355}
1356
1357////////////////////////////////////////////////////////////////////////////////
1358
1359static proc pseudo_primdecZM(module N)
1360{
1361   ideal I=quotient(N,freemodule(nrows(N)));
1362   if(size(I)==0){return(list(list(N,I)));}
1363
1364   list B=minAssZ(I);
1365   list S,R,L;
1366   ideal K;
1367   if(size(B)==0){return(S);}
1368   for(int i=1;i<=size(B);i++)
1369   {
1370      S[i]=separatorsZ(i,B);
1371   }
1372   for(i=1;i<=size(B);i++)
1373   {
1374      L=sat(N,S[i]);
1375      K[i]=S[i]^L[2];
1376      R[i]=list(L[1],B[i]);
1377   }
1378   L=pseudo_primdecZM(N+K*freemodule(nrows(N)));
1379   for(i=1;i<=size(L);i++)
1380   {
1381      R[size(R)+1]=L[i];
1382   }
1383   return(R);
1384}
1385
1386
1387
1388static proc prepare_extractZM(list L)
1389{
1390   def R=basering;
1391   module N=L[1];
1392   ideal I=quotient(N,freemodule(nrows(N)));
1393   list B=primdecZ(I);
1394   list M;
1395   if(size(B)==1){return(M);}
1396   I=std(I);
1397   list rl=ringlist(R);
1398   if(deg(I[1])==0)
1399   {
1400      execute("int p="+string(I[1])+";");
1401      rl[1]=p;
1402   }
1403   else
1404   {
1405      rl[1]=0;
1406   }
1407   def Shelp =ring(rl);
1408   setring Shelp;
1409   ideal I=imap(R,I);
1410   I=std(I);
1411   M=Primdec::maxIndependSet(I);
1412   setring R;
1413   return(M);
1414}
1415
1416
1417static proc extractZM(list M, list L)
1418{
1419   //=== M is a list of a pseudo primary module and the corresponding prime
1420   //=== L is a list of maximal independent sets for P
1421   def R=basering;
1422   ideal P=M[2];
1423   module I=M[1];
1424   poly h=1;
1425
1426   //=== size(L)=0 means P is maximal ideal and I is primary
1427   if(size(L)>0)
1428   {
1429      if(L[1][3]!=0)
1430      {
1431         //=== if u in x is an independent set of L then we compute a Groebner
1432         //=== Basis in Z[u][x-u]
1433         execute("ring S=integer,("+L[1][1]+"),lp;");
1434         module I=imap(R,I);
1435         I=std(I);
1436         list rl=ringlist(S);
1437         rl[1]=0;
1438         def Shelp =ring(rl);
1439         setring Shelp;
1440         module I=imap(S,I);
1441         //=== this is our way to obtain the coefficients in Z[u] of the
1442         //=== leading terms of the Groebner basis above
1443         def quring=Primdec::prepareQuotientring(nvars(basering)-L[1][3],"lp");
1444         setring quring;
1445         module I=imap(Shelp,I);
1446         list C;
1447         int i;
1448         for(i=1;i<=size(I);i++)
1449         {
1451         }
1452         setring Shelp;
1453         list C=imap(quring,C);
1454         setring R;
1455         list C=imap(Shelp,C);
1456      }
1457      else
1458      {
1459         // this is the case that P=<p>, p prime
1460         I=std(I);
1462         list C;
1463         int i;
1464         for(i=1;i<=size(IC);i++)
1465         {
1466            C[i]=I[i];
1467         }
1468         list rl=ringlist(R);
1469         rl[1]=0;
1470         def Shelp =ring(rl);
1471      }
1472      for(i=1;i<=size(C);i++)
1473      {
1474         if(deg(C[i])>0){h=h*C[i];}  // das muss noch besser gemacht werden,
1475                                     // nicht ausmultiplizieren!
1476      }
1477      setring Shelp;
1478      poly h=imap(R,h);
1479      ideal fac=factorize(h,1);
1480      setring R;
1481      list II;
1482      h=1;
1483      ideal fac=imap(Shelp,fac);
1484      ASSUME(1, size(fac)==ncols(fac) || size(fac)==0);
1485      for(i=1;i<=size(fac);i++)
1486      {
1487         II=sat(I,fac[i]);
1488          I=II[1];
1489          h=h*fac[i]^II[2];
1490      }
1491   }
1492   I=std(I);
1493   return(list(I,h));
1494}
1495
1496
1497proc primdecZM(module N)
1498"USAGE:  primdecZM(N); N module
1499RETURN:  a list pr of primary modules and their associated primes:
1500@format
1501   pr[i][1]   the i-th primary component,
1502   pr[i][2]   the i-th prime component.
1503@end format
1504EXAMPLE: example primdecZM; shows an example
1505"
1506{
1507  list P,K,S;
1508  int i,j;
1509  list L=pseudo_primdecZM(N);
1510  list M,O;
1511  for(i=1;i<=size(L);i++)
1512  {
1513     if(size(L[i][2])!=0)
1514     {
1515        M=prepare_extractZM(L[i]);
1516        O=extractZM(L[i],M);
1517        P[size(P)+1]=list(O[1],L[i][2]);
1518        K[size(K)+1]=L[i][1]+O[2]*freemodule(nrows(L[i][1]));
1519     }
1520     else
1521     {
1522        P[size(P)+1]=L[i];
1523     }
1524  }
1525  for(j=1;j<=size(K);j++)
1526  {
1527     S=primdecZM(K[j]);
1528     for(i=1;i<=size(S);i++)
1529     {
1530        P[size(P)+1]=S[i];
1531     }
1532  }
1533  return(P);
1534}
1535example
1536{ "EXAMPLE:";  echo = 2;
1537   ring R=integer,(x,y),(c,lp);
1538   module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18x];
1539   primdecZM(N);
1540}
1541
1542
1543////////////////////////////////////////////////////////////////////////////////
1544
1545/*
1546Examples:
1547
1548//=== IQ[a,b,c,d,e,f,g] intersect Z[a,b,c,d,e,f,g] = I  (takes some time)
1549ring R1=integer,(a,b,c,d,e,f,g),dp;
1550ideal I=a2+2de+2cf+2bg+a,
1551        2ab+e2+2df+2cg+b,
1552        b2+2ac+2ef+2dg+c,
1554        c2+2bd+2ae+2fg+e,
1555        2cd+2be+2af+g2+f,
1556        d2+2ce+2bf+2ag+g;
1557
1558ring R2=integer,(a,b,c,d,e,f,g),dp;
1559ideal I=181*32003,
1560        a2+2de+2cf+2bg+a,
1561        2ab+e2+2df+2cg+b,
1562        b2+2ac+2ef+2dg+c,
1564        c2+2bd+2ae+2fg+e,
1565        2cd+2be+2af+g2+f,
1566        d2+2ce+2bf+2ag+g;
1567
1568ring R3=integer,(w,z,y,x),dp;
1569ideal I=xzw+(-y^2+y)*z^2,
1570        (-x^2+x)*w^2+yzw,
1571        ((y^4-2*y^3+y^2)*x-y^4+y^3)*z^3,
1572        y2z2w+(-y*4+2*y^3-y^2)*z3;
1573
1574ring R4=integer,(w,z,y,x),dp;
1575ideal I=-2*yxzw+(-yx-y^2+y)*z^2,
1576        xw^2-yz^2,
1577        (yx^2-(2*y^2+2*y)*x+y^3-2*y^2+y)*z^3,
1578        (-2*y^2+2*y)*z^2*w+(yx-3*y^2-y)*z^3;
1579
1580ring R5=integer,(x,y,z),dp;
1581ideal I=x2-y2-z2,
1582        xy-z2,
1583        y3+xz2-yz2+2z3+xy-z2,
1584        -y2z2+2z4+x2-y2+z2,
1585        y3z9+3y2z10+3yz11+z12-y2z2+2z4;
1586
1587ring R6=integer,(h, l, s, x, y, z),dp;  //takes some time
1588ideal I=hl-l2-4ls+hy,
1589        h2s-6ls3+h2z,
1590        xh2-l2s-h3;
1591
1592ring R7=integer,(x,y,z),dp;
1593ideal I=x2-y2-(z+2)^2,
1594        xy-(z+2)^2,
1595        y3+x*(z+2)^2-y*(z+2)^2+2*(z+2)^3+xy-(z+2)^2,
1596        -y^2*(z+2)^2+2*(z+2)^4+x2-y2+(z+2)^2,
1597        y3z9+3y2z10+3yz11+z12-y2z2+2z4;
1598
1599ring R8=integer,(x,y,z),dp;
1600ideal I=x2-y2-(z+2)^2,
1601        xy-(z+2)^2,
1602        y3+x*(z+2)^2-y*(z+2)^2+2*(z+2)^3+xy-(z+2)^2,
1603        -y^2*(z+2)^2+2*(z+2)^4+x2-y2+(z+2)^2,
1604        y3z9+3y2z10+3yz11+z12-y2z2+2z4;
1605
1606ring R9=integer,(w,z,y,x),dp;
1607ideal I=630,
1608        ((y^2-y)*x-y^3+y^2)*z^2,
1609        (x-y)*zw,
1610        (x-y^2)*zw+(-y^2+y)*z^2,
1611        (-x^2+x)*w^2+(-yx+y)*zw;
1612
1613ring R10=integer,(w,z,y,x),dp;
1614ideal I=1260,
1615        -yxzw+(-y^2+y)*z^2,
1616        (-x^2+x)*w^2-yxzw,
1617        ((-y^2+y)*x-y^3+2*y^2-y)*z^3,
1618        (y^2-y)*z^2*w+(-y^2+y)*z^2*w+(-y^2+y)*z^3;
1619
1620ring R11=integer,(w,z,y,x),dp;
1621ideal I=(4*y^2*x^2+(4*y^3+4*y^2-y)*x-y^2-y)*z^2,
1622        (x+y+1)*zw+(-4*y^2*x-4*y^3-4*y^2)*z^2,
1623        (-x-2*y^2 - 2*y - 1)*zw + (8*y^3*x + 8*y^4 + 8*y^3 + 2*y^2+y)*z^2,
1624        ((y^3 + y^2)*x - y^2 - y)*z^2,
1625        (y +1)*zw + (-y^3 -y^2)*z^2,
1626        (x + 1)*zw +(- y^2 -y)*z^2,
1627        (x^2 +x)*w^2 + (-yx - y)*zw;
1628
1629ring R12=integer,(w,z,y,x),dp;
1630ideal I=72,
1631        ((y^3 + y^2)*x - y^2 - y)*z^2,
1632        (y + 1)*zw + (-y^3 -y^2)*z^2,
1633        (x + 1)*zw + (-y^2 -y)*z^2, (x^2 + x)*w^2 + (-yx - y)*zw;
1634
1635ring R13=integer,(w,z,y,x),dp;
1636ideal I=(((12*y+8)*x^2 +(2*y+2)*x)*zw +((-15*y^2 -4*y)*x-4*y^2 -y)*z^2,
1637        -x*w^2 +((-12*y -8)*x+2*y)*zw +(15*y^2+4*y)*z^2,
1638        (81*y^4*x^2 +(-54*y^3 -12*y^2)*x-12*y^3 -3*y^2)*z^3,
1639        (-24*yx+6*y^2-6*y)*z^2*w + (-81*y^4*x + 81*y^3 + 24*y^2)*z^3,
1640        (48*x^2 + (-30*y + 12)*x - 6*y)*z^2*w + ((81*y^3 -54*y^2 -24*y)*x
1641        -21*y^2 -6*y)*z^3,
1642        (-96*yx-18*y^3 +18*y^2-24*y)*z^2*w +(243*y^5*x-243*y^4 +72*y^3
1643        +48*y^2)*z^3,
1644        6*y*z^2*w^2 +((576*y+384)*x^2 + (-81*y^3 -306*y^2 -168*y+96)*x+81*y^2
1645        -18*y)*z^3*w +((-720*y^2 - 192*y)*x + 450*y^3 - 60*y^2 - 48*y)*z^4);
1646
1647ring R14=integer,(x(1),x(2),x(3),x(4)),dp;
1648ideal I=181*49^2,
1649        x(4)^4,
1650        x(1)*x(4)^3,
1651        x(1)*x(2)*x(4)^2,
1652        x(2)^2*x(4)^2,
1653        x(2)^2*x(3)*x(4),
1654        x(1)*x(2)*x(3)*x(4),
1655        x(1)*x(3)^2*x(4),
1656        x(3)^3*x(4);
1657
1658
1659ring R15=integer,(x,y,z),dp;
1660ideal I=32003*181*64,
1661        ((z^2-z)*y^2 + (z^2 -z)*y)*x; (z*y^3 + z*y^2)*x,
1662        (y^4 - y^2)*x, (z^2 - z)*y*x^2, (y^3 - y^2)*x^2,
1663        (z^3 - z^2)*x^4 + (2*z^3 -2*z^2)*x^3 + (z^3 -z^2)*x^2,
1664        z*y^2*x^2, z*y*x^4 +z*y*x^3,
1665        2*y^2*x^4 +6*y^2*x^3 +6*y^2*x^2 + (y^3 +y^2)*x, z*x^5 + (z^2 +z)*x^4
1666        + (2*z^2 -z)*x^3 + (z^2 -z)*x^2,
1667        y*x^6 + 3*y*x^5 + 3*y*x^4 + y*x^3;
1668
1669
1670ring R16=integer,(x(1),x(2),x(3),x(4),x(5)),dp;
1671ideal I=x(5)^5,
1672        x(1)*x(5)^4,
1673        x(1)*x(2)*x(5)^3,
1674        x(2)^2*x(5)^3,
1675        x(2)^2*x(3)*x(5)^2,
1676        x(1)*x(2)*x(3)*x(5)^2,
1677        x(1)*x(3)^2*x(5)^2,
1678        x(3)^3*x(5)^2,
1679        x(3)^3*x(4)*x(5),
1680        x(1)*x(3)^2*x(4)*x(5),
1681        x(1)*x(2)*x(3)*x(4)*x(5),
1682        x(2)^2*x(3)*x(4)*x(5),
1683        x(2)^2*x(4)^2*x(5),
1684        x(1)*x(2)*x(4)^2*x(5),
1685        x(1)*x(4)^3*x(5),
1686        x(4)^4*x(5);
1687      I=intersectZ(I,ideal(64*181,x(1)^2));
1688
1689ring R17=integer,(x,y,z),dp;
1690ideal I=374,
1691        (z+2)^8-140z6+2622*(z+2)^4-1820*(z+2)^2+169,
1692        17y*(z+2)^4-374*y*(z+2)^2+221y+2z7-281z5+5240z3-3081z,
1693        204y2+136yz3-3128yz+z6-149z4+2739z2+117,
1694        17xz4-374xz2+221x+2z7-281z5+5240z3-3081z,
1695        136xy-136xz-136yz+2z6-281z4+5376z2-3081,
1696        204x2+136xz3-3128xz+z6-149z4+2739z2+117;
1697
1698ring R18=integer,(B,D,F,b,d,f),dp;
1699ideal I=6,
1700        (b-d)*(B-D)-2*F+2,
1701        (b-d)*(B+D-2*F)+2*(B-D),
1702        (b-d)^2-2*(b+d)+f+1,
1703        B^2*b^3-1,
1704        D^2*d^3-1,
1705        F^2*f^3-1;
1706
1707ring R19=integer,(a,b,c,d,e,f),dp;
1708ideal I=24,
1709        2*(f+2)*b+2ec+d2+a2+a,
1710        2*(f+2)*c+2ed+2ba+b,
1711        2*(f+2)*d+e2+2ca+c+b2,
1712        2*(f+2)*e+2da+d+2cb,
1713        (f+2)^2+2ea+e+2db+c2,
1714        2*(f+2)*a+f+2eb+2dc;
1715
1716ring R20=integer,(x,y,z,w,u),dp;
1717ideal I=24,
1718         2x2-2y2+2z2-2w2+2u2-1,
1719         2x3-2y3+2z3-2w3+2u3-1,
1720         2x4-2y4+2z4-2w4+2u4-1,
1721         2x5-2y5+2z5-2w5+2u5-1,
1722         2x6-2y6+2z6-2w6+2u6-1;
1723
1724ring R21=integer,(x,y,z,t,u,v,h),dp;
1725ideal I=66,
1726        2x2+2y2+2z2+2t2+2u2+v2-vh,
1727        xy+yz+2zt+2tu+2uv-uh,
1728        2xz+2yt+2zu+u2+2tv-th,
1729        2xt+2yu+2tu+2zv-zh,
1730        t2+2xv+2yv+2zv-yh,
1731        2x+2y+2z+2t+2u+v-h,
1732        x3+y3+z3+t3+u3+v3;
1733
1734ring R22=integer,(s,p,S,P,T,F,f),dp;
1735ideal I=35,
1736        2*T-S*s-2*F+2,
1737        8*F*p-4*p*S-2*F*s^2+S*s^2+4*T-2*S*s,
1738        -2*s-4*p+s^2+f+1,
1739        s*T^2-p*s*P-p*S*T-2,
1740        p^3*P^2-1,
1741        F^2*f^3-1;
1742
1743ring R=integer,(x,y),(c,lp);
1744module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18x];
1745
1746ring R=integer,(x,y),(c,lp);
1747module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,2xy-x],[x,0,-xy],[0,0,18];
1748
1749ring R=integer,(x,y),(c,lp);
1750module N=[-y,7,0],[2y3-y2],[3x,y2],[2y-y2,x],[4,5x3];
1751
1752ring r=integer,(x,y),(c,lp);
1753module N=[0,0,xy2-x2-xy],[0,y,x],[0,x,xy-x],[x,0,-xy],[5x,0,0];
1754
1755ring R2=integer,(a(1),a(2),a(3),b(1),b(2),b(3)),(c,lp);
1756module N=[a(1)*b(1),a(2)*b(1),a(3)*b(1)],[a(1)*b(2),a(2)*b(2),a(3)*b(2)],[a(1)*b(3),a(2)*b(3),a(3)*b(3)];
1757
1758ring R3=integer,(x,y,z),(c,lp);
1759module N=[y2+z2,xy,xz],[xy,x2+z2,yz],[xz,yz,x2+y2];
1760
1761ring R4=integer,(x,y,z,a,b,c),(c,lp);
1762module N=[x3y2z2c,x2y3z2c,x2y2z3c],[x3y2z2b,x2y3z2b,x2y2z3b],[x3y2z2a,x2y3z2a,x2y2z3a];
1763
1764*/
Note: See TracBrowser for help on using the repository browser.