source: git/Singular/LIB/modstd.lib @ 91c7251

spielwiese
Last change on this file since 91c7251 was 230493, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: syntax in modstd.lib
  • Property mode set to 100644
File size: 39.6 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version modstd.lib 4.0.0.0 Dec_2013 "; // $Id$
3category = "Commutative Algebra";
4info="
5LIBRARY:  modstd.lib      Groebner basis of ideals
6
7AUTHORS:  A. Hashemi      Amir.Hashemi@lip6.fr
8@*        G. Pfister      pfister@mathematik.uni-kl.de
9@*        H. Schoenemann  hannes@mathematik.uni-kl.de
10@*        A. Steenpass    steenpass@mathematik.uni-kl.de
11@*        S. Steidel      steidel@mathematik.uni-kl.de
12
13OVERVIEW:
14
15  A library for computing the Groebner basis of an ideal in the polynomial
16  ring over the rational numbers using modular methods. The procedures are
17  inspired by the following paper:
18  Elizabeth A. Arnold: Modular algorithms for computing Groebner bases.
19  Journal of Symbolic Computation 35, 403-419 (2003).
20
21PROCEDURES:
22 modStd(I);        standard basis of I using modular methods (chinese remainder)
23 modS(I,L);        liftings to Q of standard bases of I mod p for p in L
24 modHenselStd(I);  standard basis of I using modular methods (hensel lifting)
25";
26
27LIB "poly.lib";
28LIB "ring.lib";
29LIB "parallel.lib";
30
31////////////////////////////////////////////////////////////////////////////////
32
33static proc mod_init()
34{
35   newstruct("idealPrimeTest", "ideal Ideal");
36}
37
38////////////////////////////////////////////////////////////////////////////////
39
40proc mixedTest()
41"USAGE:  mixedTest();
42RETURN:  1 if ordering of basering is mixed, 0 else
43EXAMPLE: example mixedTest(); shows an example
44"
45{
46   int i,p,m;
47   for(i = 1; i <= nvars(basering); i++)
48   {
49      if(var(i) > 1)
50      {
51         p++;
52      }
53      else
54      {
55         m++;
56      }
57   }
58   if((p > 0) && (m > 0)) { return(1); }
59   return(0);
60}
61example
62{ "EXAMPLE:"; echo = 2;
63   ring R1 = 0,(x,y,z),dp;
64   mixedTest();
65   ring R2 = 31,(x(1..4),y(1..3)),(ds(4),lp(3));
66   mixedTest();
67   ring R3 = 181,x(1..9),(dp(5),lp(4));
68   mixedTest();
69}
70
71////////////////////////////////////////////////////////////////////////////////
72
73static proc redFork(ideal I, ideal J, int n)
74{
75   attrib(J,"isSB",1);
76   return(reduce(I,J,1));
77}
78
79////////////////////////////////////////////////////////////////////////////////
80
81proc isIncluded(ideal I, ideal J, list #)
82"USAGE:  isIncluded(I,J); I,J ideals
83RETURN:  1 if J includes I,
84@*       0 if there is an element f in I which does not reduce to 0 w.r.t. J.
85EXAMPLE: example isIncluded; shows an example
86"
87{
88   def R = basering;
89   setring R;
90
91   attrib(J,"isSB",1);
92   int i,j,k;
93
94   if(size(#) > 0)
95   {
96      int n = #[1];
97      if(n >= ncols(I)) { n = ncols(I); }
98      if(n > 1)
99      {
100         for(i = 1; i <= n - 1; i++)
101         {
102            //link l(i) = "MPtcp:fork";
103            link l(i) = "ssi:fork";
104            open(l(i));
105
106            write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1)));
107         }
108
109         int t = timer;
110         if(reduce(I[ncols(I)], J, 1) != 0)
111         {
112            for(i = 1; i <= n - 1; i++)
113            {
114               close(l(i));
115            }
116            return(0);
117         }
118         t = timer - t;
119         if(t > 60) { t = 60; }
120         int i_sleep = system("sh", "sleep "+string(t));
121
122         j = ncols(I) - n;
123
124         while(j >= 0)
125         {
126            for(i = 1; i <= n - 1; i++)
127            {
128               if(status(l(i), "read", "ready"))
129               {
130                  if(read(l(i)) != 0)
131                  {
132                     for(i = 1; i <= n - 1; i++)
133                     {
134                        close(l(i));
135                     }
136                     return(0);
137                  }
138                  else
139                  {
140                     if(j >= 1)
141                     {
142                        write(l(i), quote(redFork(eval(I[j]), eval(J), 1)));
143                        j--;
144                     }
145                     else
146                     {
147                        k++;
148                        close(l(i));
149                     }
150                  }
151               }
152            }
153            if(k == n - 1)
154            {
155               j--;
156            }
157            i_sleep = system("sh", "sleep "+string(t));
158         }
159         return(1);
160      }
161   }
162
163   for(i = ncols(I); i >= 1; i--)
164   {
165      if(reduce(I[i],J,1) != 0){ return(0); }
166   }
167   return(1);
168}
169example
170{ "EXAMPLE:"; echo = 2;
171   ring r=0,(x,y,z),dp;
172   ideal I = x+1,x+y+1;
173   ideal J = x+1,y;
174   isIncluded(I,J);
175   isIncluded(J,I);
176   isIncluded(I,J,4);
177
178   ring R = 0, x(1..5), dp;
179   ideal I1 = cyclic(4);
180   ideal I2 = I1,x(5)^2;
181   isIncluded(I1,I2,4);
182}
183
184////////////////////////////////////////////////////////////////////////////////
185
186proc pTestSB(ideal I, ideal J, list L, int variant, list #)
187"USAGE:  pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
188RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
189         J mod p is (resp. is not) a standard basis of I mod p
190EXAMPLE: example pTestSB; shows an example
191"
192{
193   int i,j,k,p;
194   def R = basering;
195   list r = ringlist(R);
196
197   while(!j)
198   {
199      j = 1;
200      p = prime(random(1000000000,2134567879));
201      for(i = 1; i <= size(L); i++)
202      {
203         if(p == L[i]) { j = 0; break; }
204      }
205      if(j)
206      {
207         for(i = 1; i <= ncols(I); i++)
208         {
209            for(k = 2; k <= size(I[i]); k++)
210            {
211               if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; }
212            }
213            if(!j){ break; }
214         }
215      }
216      if(j)
217      {
218         if(!primeTest(I,p)) { j = 0; }
219      }
220   }
221   r[1] = p;
222   def @R = ring(r);
223   setring @R;
224   ideal I = imap(R,I);
225   ideal J = imap(R,J);
226   attrib(J,"isSB",1);
227
228   int t = timer;
229   j = 1;
230   if(isIncluded(I,J) == 0) { j = 0; }
231
232   if(printlevel >= 11)
233   {
234      "isIncluded(I,J) takes "+string(timer - t)+" seconds";
235      "j = "+string(j);
236   }
237
238   t = timer;
239   if(j)
240   {
241      if(size(#) > 0)
242      {
243         ideal K = modpStd(I,p,variant,#[1])[1];
244      }
245      else
246      {
247         ideal K = groebner(I);
248      }
249      t = timer;
250      if(isIncluded(J,K) == 0) { j = 0; }
251
252      if(printlevel >= 11)
253      {
254         "isIncluded(J,K) takes "+string(timer - t)+" seconds";
255         "j = "+string(j);
256      }
257   }
258   setring R;
259   return(j);
260}
261example
262{ "EXAMPLE:"; echo = 2;
263   intvec L = 2,3,5;
264   ring r = 0,(x,y,z),dp;
265   ideal I = x+1,x+y+1;
266   ideal J = x+1,y;
267   pTestSB(I,I,L,2);
268   pTestSB(I,J,L,2);
269}
270
271////////////////////////////////////////////////////////////////////////////////
272
273proc deleteUnluckyPrimes(list T, list L, int ho, list #)
274"USAGE:  deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer
275RETURN:  lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #),
276         lT ideal
277NOTE:    - if ho = 1, the polynomials in T are homogeneous, else ho = 0,
278@*       - lT is prevalent, i.e. the most appearing leading ideal in T
279EXAMPLE: example deleteUnluckyPrimes; shows an example
280"
281{
282   ho = ((ho)||(ord_test(basering) == -1));
283   int j,k,c;
284   intvec hl,hc;
285   ideal cT,lT,cK;
286   lT = lead(T[size(T)]);
287   attrib(lT,"isSB",1);
288   if(!ho)
289   {
290      for(j = 1; j < size(T); j++)
291      {
292         cT = lead(T[j]);
293         attrib(cT,"isSB",1);
294         if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
295         {
296            cK = cT;
297            c++;
298         }
299      }
300      if(c > size(T) div 2){ lT = cK; }
301   }
302   else
303   {
304      hl = hilb(lT,1);
305      for(j = 1; j < size(T); j++)
306      {
307         cT = lead(T[j]);
308         attrib(cT,"isSB",1);
309         hc = hilb(cT,1);
310         if(hl == hc)
311         {
312            for(k = 1; k <= size(lT); k++)
313            {
314               if(lT[k] < cT[k]) { lT = cT; c++; break; }
315               if(lT[k] > cT[k]) { c++; break; }
316            }
317         }
318         else
319         {
320            if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; }
321         }
322      }
323   }
324
325   int addList;
326   if(size(#) > 0) { list M = #; addList = 1; }
327   j = 1;
328   attrib(lT,"isSB",1);
329   while((j <= size(T))&&(c > 0))
330   {
331      cT = lead(T[j]);
332      attrib(cT,"isSB",1);
333      if((size(reduce(cT,lT)) != 0)||(size(reduce(lT,cT)) != 0))
334      {
335         T = delete(T,j);
336         if(j == 1)
337         {
338            L = L[2..size(L)];
339            if(addList == 1) { M = M[2..size(M)]; }
340         }
341         else
342         {
343            if(j == size(L))
344            {
345               L = L[1..size(L)-1];
346               if(addList == 1) { M = M[1..size(M)-1]; }
347            }
348            else
349            {
350               L = L[1..j-1],L[j+1..size(L)];
351               if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; }
352            }
353         }
354         j--;
355      }
356      j++;
357   }
358
359   for(j = 1; j <= size(L); j++)
360   {
361      L[j] = bigint(L[j]);
362   }
363
364   if(addList == 0) { return(list(T,L,lT)); }
365   if(addList == 1) { return(list(T,L,M,lT)); }
366}
367example
368{ "EXAMPLE:"; echo = 2;
369   list L = 2,3,5,7,11;
370   ring r = 0,(y,x),Dp;
371   ideal I1 = 2y2x,y6;
372   ideal I2 = yx2,y3x,x5,y6;
373   ideal I3 = y2x,x3y,x5,y6;
374   ideal I4 = y2x,11x3y,x5;
375   ideal I5 = y2x,yx3,x5,7y6;
376   list T = I1,I2,I3,I4,I5;
377   deleteUnluckyPrimes(T,L,1);
378   list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5);
379   deleteUnluckyPrimes(T,L,1,P);
380}
381
382////////////////////////////////////////////////////////////////////////////////
383
384proc primeTest(def II, bigint p)
385{
386   if(typeof(II) == "string")
387   {
388      ideal I = `II`.Ideal;
389   }
390   else
391   {
392      ideal I = II;
393   }
394
395   I = simplify(I, 2);   // erase zero generators
396
397   int i,j;
398   poly f;
399   number cnt;
400   for(i = 1; i <= size(I); i++)
401   {
402      f = cleardenom(I[i]);
403      if(f == 0) { return(0); }
404      cnt = leadcoef(I[i])/leadcoef(f);
405      if((numerator(cnt) mod p) == 0) { return(0); }
406      if((denominator(cnt) mod p) == 0) { return(0); }
407      for(j = size(f); j > 0; j--)
408      {
409         if((leadcoef(f[j]) mod p) == 0) { return(0); }
410      }
411   }
412   return(1);
413}
414
415////////////////////////////////////////////////////////////////////////////////
416
417proc primeList(ideal I, int n, list #)
418"USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
419         n integer
420RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
421         < L[size(L)] union with L) such that none of these primes divides any
422         coefficient occuring in I
423NOTE:    The number of cores to use can be defined by ncores, default is 1.
424EXAMPLE: example primeList; shows an example
425"
426{
427   intvec L;
428   int i,p;
429   int ncores = 1;
430
431//-----------------  Initialize optional parameter ncores  ---------------------
432   if(size(#) > 0)
433   {
434      if(size(#) == 1)
435      {
436         if(typeof(#[1]) == "int")
437         {
438            ncores = #[1];
439            # = list();
440         }
441      }
442      else
443      {
444         ncores = #[2];
445      }
446   }
447
448   if(size(#) == 0)
449   {
450      p = 2147483647;
451      while(!primeTest(I,p))
452      {
453         p = prime(p-1);
454         if(p == 2) { ERROR("no more primes"); }
455      }
456      L[1] = p;
457   }
458   else
459   {
460      L = #[1];
461      p = prime(L[size(L)]-1);
462      while(!primeTest(I,p))
463      {
464         p = prime(p-1);
465         if(p == 2) { ERROR("no more primes"); }
466      }
467      L[size(L)+1] = p;
468   }
469   if(p == 2) { ERROR("no more primes"); }
470   if(ncores == 1)
471   {
472      for(i = 2; i <= n; i++)
473      {
474         p = prime(p-1);
475         while(!primeTest(I,p))
476         {
477            p = prime(p-1);
478            if(p == 2) { ERROR("no more primes"); }
479         }
480         L[size(L)+1] = p;
481      }
482   }
483   else
484   {
485      int neededSize = size(L)+n-1;;
486      list parallelResults;
487      list arguments;
488      int neededPrimes = neededSize-size(L);
489      idealPrimeTest Id;
490      Id.Ideal = I;
491      export(Id);
492      while(neededPrimes > 0)
493      {
494         arguments = list();
495         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
496            *ncores; i > 0; i--)
497         {
498            p = prime(p-1);
499            if(p == 2) { ERROR("no more primes"); }
500            arguments[i] = list("Id", p);
501         }
502         parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores);
503         for(i = size(arguments); i > 0; i--)
504         {
505            if(parallelResults[i])
506            {
507               L[size(L)+1] = arguments[i][2];
508            }
509         }
510         neededPrimes = neededSize-size(L);
511      }
512      kill Id;
513      if(size(L) > neededSize)
514      {
515         L = L[1..neededSize];
516      }
517   }
518   return(L);
519}
520example
521{ "EXAMPLE:"; echo = 2;
522   ring r = 0,(x,y,z),dp;
523   ideal I = 2147483647x+y, z-181;
524   intvec L = primeList(I,10);
525   size(L);
526   L[1];
527   L[size(L)];
528   L = primeList(I,5,L);
529   size(L);
530   L[size(L)];
531}
532
533////////////////////////////////////////////////////////////////////////////////
534
535static proc liftstd1(ideal I)
536{
537   def R = basering;
538   list rl = ringlist(R);
539   list ordl = rl[3];
540
541   int i;
542   for(i = 1; i <= size(ordl); i++)
543   {
544      if((ordl[i][1] == "C") || (ordl[i][1] == "c"))
545      {
546         ordl = delete(ordl, i);
547         break;
548      }
549   }
550
551   ordl = insert(ordl, list("c", 0));
552   rl[3] = ordl;
553   def newR = ring(rl);
554   setring newR;
555   ideal I = imap(R,I);
556
557   intvec opt = option(get);
558   option(none);
559   option(prompt);
560
561   module M;
562   for(i = 1; i <= size(I); i++)
563   {
564      M = M + module(I[i]*gen(1) + gen(i+1));
565      M = M + module(gen(i+1));
566   }
567
568   module sM = std(M);
569
570   ideal sI;
571   if(attrib(R,"global"))
572   {
573      for(i = size(I)+1; i <= size(sM); i++)
574      {
575         sI[size(sI)+1] = sM[i][1];
576      }
577      matrix T = submat(sM,2..nrows(sM),size(I)+1..ncols(sM));
578   }
579   else
580   {
581      //"==========================================================";
582      //"WARNING: Algorithm is not applicable if ordering is mixed.";
583      //"==========================================================";
584      for(i = 1; i <= size(sM)-size(I); i++)
585      {
586         sI[size(sI)+1] = sM[i][1];
587      }
588      matrix T = submat(sM,2..nrows(sM),1..ncols(sM)-size(I));
589   }
590
591   setring R;
592   option(set, opt);
593   return(imap(newR,sI),imap(newR,T));
594}
595example
596{ "EXAMPLE:"; echo = 2;
597   ring R = 0,(x,y,z),dp;
598   poly f = x3+y7+z2+xyz;
599   ideal i = jacob(f);
600   matrix T;
601   ideal sm = liftstd(i,T);
602   sm;
603   print(T);
604   matrix(sm) - matrix(i)*T;
605
606
607   ring S = 32003, x(1..5), lp;
608   ideal I = cyclic(5);
609   ideal sI;
610   matrix T;
611   sI,T = liftstd1(I);
612   matrix(sI) - matrix(I)*T;
613}
614
615////////////////////////////////////////////////////////////////////////////////
616
617proc modpStd(ideal I, int p, int variant, list #)
618"USAGE:  modpStd(I,p,variant,#); I ideal, p integer, variant integer
619ASSUME:  If size(#) > 0, then #[1] is an intvec describing the Hilbert series.
620RETURN:  ideal - a standard basis of I mod p, integer - p
621NOTE:    The procedure computes a standard basis of the ideal I modulo p and
622         fetches the result to the basering. If size(#) > 0 the Hilbert driven
623         standard basis computation std(.,#[1]) is used instead of groebner.
624         The standard basis computation modulo p does also vary depending on the
625         integer variant, namely
626@*       - variant = 1: std(.,#[1]) resp. groebner,
627@*       - variant = 2: groebner,
628@*       - variant = 3: homog. - std(.,#[1]) resp. groebner - dehomog.,
629@*       - variant = 4: fglm.
630EXAMPLE: example modpStd; shows an example
631"
632{
633   def R0 = basering;
634   list rl = ringlist(R0);
635   rl[1] = p;
636   def @r = ring(rl);
637   setring @r;
638   ideal i = fetch(R0,I);
639
640   option(redSB);
641
642   if(variant == 1)
643   {
644      if(size(#) > 0)
645      {
646         i = std(i, #[1]);
647      }
648      else
649      {
650         i = groebner(i);
651      }
652   }
653
654   if(variant == 2)
655   {
656      i = groebner(i);
657   }
658
659   if(variant == 3)
660   {
661      list rl = ringlist(@r);
662      int nvar@r = nvars(@r);
663
664      int k;
665      intvec w;
666      for(k = 1; k <= nvar@r; k++)
667      {
668         w[k] = deg(var(k));
669      }
670      w[nvar@r + 1] = 1;
671
672      rl[2][nvar@r + 1] = "homvar";
673      rl[3][2][2] = w;
674
675      def HomR = ring(rl);
676      setring HomR;
677      ideal i = imap(@r, i);
678      i = homog(i, homvar);
679
680      if(size(#) > 0)
681      {
682         if(w == 1)
683         {
684            i = std(i, #[1]);
685         }
686         else
687         {
688            i = std(i, #[1], w);
689         }
690      }
691      else
692      {
693         i = groebner(i);
694      }
695
696      i = subst(i, homvar, 1);
697      i = simplify(i, 34);
698
699      setring @r;
700      i = imap(HomR, i);
701      i = interred(i);
702      kill HomR;
703   }
704
705   if(variant == 4)
706   {
707      def R1 = changeord(list(list("dp",1:nvars(basering))));
708      setring R1;
709      ideal i = fetch(@r,i);
710      i = std(i);
711      setring @r;
712      i = fglm(R1,i);
713   }
714
715   setring R0;
716   return(list(fetch(@r,i),p));
717}
718example
719{ "EXAMPLE:"; echo = 2;
720   ring r = 0, x(1..4), dp;
721   ideal I = cyclic(4);
722   int p = 181;
723   list P = modpStd(I,p,2);
724   P;
725
726   ring r2 = 0, x(1..5), lp;
727   ideal I = cyclic(5);
728   int q = 32003;
729   list Q = modpStd(I,q,4);
730   Q;
731}
732
733static proc algeModStd(ideal i,list #)
734{
735//reduces modStd over algebraic extensions to the one over a polynomial ring
736   list L=#;
737   def R=basering;
738   int n=nvars(R);
739   list rl=ringlist(R);
740   poly p=minpoly;
741   rl[2][n+1]=rl[1][2][1];
742   rl[1]=rl[1][1];
743   rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
744   rl[3][size(rl[3])-1]=list("dp",1);
745   def S=ring(rl);
746   setring S;
747   poly p=imap(R,p);
748   ideal i=imap(R,i);
749   i=i,p;
750   ideal j=modStd(i,L);
751   if(j[1]==p)
752   {
753      j[1]=0;
754   }
755   j=simplify(j,2);
756   setring R;
757   ideal j=imap(S,j);
758   return(j);
759}
760////////////////////////////// main procedures /////////////////////////////////
761
762proc modStd(ideal I, list #)
763"USAGE:  modStd(I); I ideal
764ASSUME:  If size(#) > 0, then # contains either 1, 2 or 4 integers such that
765@*       - #[1] is the number of available processors for the computation,
766@*       - #[2] is an optional parameter for the exactness of the computation,
767                if #[2] = 1, the procedure computes a standard basis for sure,
768@*       - #[3] is the number of primes until the first lifting,
769@*       - #[4] is the constant number of primes between two liftings until
770           the computation stops.
771RETURN:  a standard basis of I if no warning appears;
772NOTE:    The procedure computes a standard basis of I (over the rational
773         numbers) by using modular methods.
774         By default the procedure computes a standard basis of I for sure, but
775         if the optional parameter #[2] = 0, it computes a standard basis of I
776         with high probability.
777         The procedure distinguishes between different variants for the standard
778         basis computation in positive characteristic depending on the ordering
779         of the basering, the parameter #[2] and if the ideal I is homogeneous.
780@*       - variant = 1, if I is homogeneous,
781@*       - variant = 2, if I is not homogeneous, 1-block-ordering,
782@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
783                        > 1 block),
784@*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
785EXAMPLE: example modStd; shows an example
786"
787{
788   int TT = timer;
789   int RT = rtimer;
790
791   def R0 = basering;
792   list rl = ringlist(R0);
793
794   int algebraic;
795   if(size(#)>0)
796   {
797      if(#[1]<=0)
798      {
799         algebraic=1;
800         #[1]=-#[1];
801         if(#[1]==0){list LK;#=LK;}
802      }
803   }
804
805   if((npars(R0) > 0) || (rl[1][1] > 0))
806   {
807      if(npars(R0)==1)
808      {
809         if(minpoly!=0)
810         {
811            list LM=#;
812            if(size(LM)==0){LM[1]=0;}
813            LM[1]=-LM[1];
814            return(algeModStd(I,LM));
815         }
816      }
817
818      ERROR("Characteristic of basering should be zero, basering should
819             have no parameters.");
820   }
821
822   int index = 1;
823   int i,k,c;
824   int j = 1;
825   int pTest, sizeTest;
826   int en = 2134567879;
827   int an = 1000000000;
828   bigint N;
829
830//--------------------  Initialize optional parameters  ------------------------
831   if(size(#) > 0)
832   {
833      if(size(#) == 1)
834      {
835         int n1 = #[1];
836         int exactness = 1;
837         if(n1 >= 10)
838         {
839            int n2 = n1 + 1;
840            int n3 = n1;
841         }
842         else
843         {
844            int n2 = 10;
845            int n3 = 10;
846         }
847      }
848      if(size(#) == 2)
849      {
850         int n1 = #[1];
851         int exactness = #[2];
852         if(n1 >= 10)
853         {
854            int n2 = n1 + 1;
855            int n3 = n1;
856         }
857         else
858         {
859            int n2 = 10;
860            int n3 = 10;
861         }
862      }
863      if(size(#) == 4)
864      {
865         int n1 = #[1];
866         int exactness = #[2];
867         if(n1 >= #[3])
868         {
869            int n2 = n1 + 1;
870         }
871         else
872         {
873            int n2 = #[3];
874         }
875         if(n1 >= #[4])
876         {
877            int n3 = n1;
878         }
879         else
880         {
881            int n3 = #[4];
882         }
883      }
884   }
885   else
886   {
887      int n1 = 1;
888      int exactness = 1;
889      int n2 = 10;
890      int n3 = 10;
891   }
892
893   if(printlevel >= 10)
894   {
895      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
896       +", exactness = "+string(exactness);
897   }
898
899//-------------------------  Save current options  -----------------------------
900   intvec opt = option(get);
901
902   option(redSB);
903
904//--------------------  Initialize the list of primes  -------------------------
905   int tt = timer;
906   int rt = rtimer;
907   intvec L = primeList(I,n2,n1);
908   if(printlevel >= 10)
909   {
910      "CPU-time for primeList: "+string(timer-tt)+" seconds.";
911      "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
912   }
913   L[5] = prime(random(an,en));
914
915//---------------------  Decide which variant to take  -------------------------
916   int variant;
917   int h = homog(I);
918
919   tt = timer;
920   rt = rtimer;
921
922   if(!mixedTest())
923   {
924      if(h)
925      {
926         variant = 1;
927         if(printlevel >= 10) { "variant = 1"; }
928
929         rl[1] = L[5];
930         def @r = ring(rl);
931         setring @r;
932         def @s = changeord(list(list("dp",1:nvars(basering))));
933         setring @s;
934         ideal I = std(fetch(R0,I));
935         intvec hi = hilb(I,1);
936         setring R0;
937         kill @r,@s;
938      }
939      else
940      {
941         string ordstr_R0 = ordstr(R0);
942         int neg = 1 - attrib(R0,"global");
943
944         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
945         {
946            variant = 2;
947            if(printlevel >= 10) { "variant = 2"; }
948         }
949         else
950         {
951            string order;
952            if(system("nblocks") <= 2)
953            {
954               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
955                                       + find(ordstr_R0, "rp") <= 0)
956               {
957                  order = "simple";
958               }
959            }
960
961            if((order == "simple") || (size(rl) > 4))
962            {
963               variant = 2;
964               if(printlevel >= 10) { "variant = 2"; }
965            }
966            else
967            {
968               rl[1] = L[5];
969               def @r = ring(rl);
970               setring @r;
971
972               def @s = changeord(list(list("dp",1:nvars(basering))));
973               setring @s;
974               ideal I = std(fetch(R0,I));
975               if(dim(I) == 0)
976               {
977                  variant = 4;
978                  if(printlevel >= 10) { "variant = 4"; }
979               }
980               else
981               {
982                  variant = 3;
983                  if(printlevel >= 10) { "variant = 3"; }
984
985                  int nvar@r = nvars(@r);
986                  intvec w;
987                  for(i = 1; i <= nvar@r; i++)
988                  {
989                     w[i] = deg(var(i));
990                  }
991                  w[nvar@r + 1] = 1;
992
993                  list hiRi = hilbRing(fetch(R0,I),w);
994                  intvec W = hiRi[2];
995                  @s = hiRi[1];
996                  setring @s;
997
998                  Id(1) = std(Id(1));
999                  intvec hi = hilb(Id(1), 1, W);
1000               }
1001
1002               setring R0;
1003               kill @r,@s;
1004            }
1005         }
1006      }
1007   }
1008   else
1009   {
1010      if(exactness == 1) { return(groebner(I)); }
1011      if(h)
1012      {
1013         variant = 1;
1014         if(printlevel >= 10) { "variant = 1"; }
1015         rl[1] = L[5];
1016         def @r = ring(rl);
1017         setring @r;
1018         def @s = changeord(list(list("dp",1:nvars(basering))));
1019         setring @s;
1020         ideal I = std(fetch(R0,I));
1021         intvec hi = hilb(I,1);
1022         setring R0;
1023         kill @r,@s;
1024      }
1025      else
1026      {
1027         string ordstr_R0 = ordstr(R0);
1028         int neg = 1 - attrib(R0,"global");
1029
1030         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
1031         {
1032            variant = 2;
1033            if(printlevel >= 10) { "variant = 2"; }
1034         }
1035         else
1036         {
1037            string order;
1038            if(system("nblocks") <= 2)
1039            {
1040               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
1041                                       + find(ordstr_R0, "rp") <= 0)
1042               {
1043                  order = "simple";
1044               }
1045            }
1046
1047            if((order == "simple") || (size(rl) > 4))
1048            {
1049               variant = 2;
1050               if(printlevel >= 10) { "variant = 2"; }
1051            }
1052            else
1053            {
1054               variant = 3;
1055               if(printlevel >= 10) { "variant = 3"; }
1056
1057               rl[1] = L[5];
1058               def @r = ring(rl);
1059               setring @r;
1060               int nvar@r = nvars(@r);
1061               intvec w;
1062               for(i = 1; i <= nvar@r; i++)
1063               {
1064                  w[i] = deg(var(i));
1065               }
1066               w[nvar@r + 1] = 1;
1067
1068               list hiRi = hilbRing(fetch(R0,I),w);
1069               intvec W = hiRi[2];
1070               def @s = hiRi[1];
1071               setring @s;
1072
1073               Id(1) = std(Id(1));
1074               intvec hi = hilb(Id(1), 1, W);
1075
1076               setring R0;
1077               kill @r,@s;
1078            }
1079         }
1080      }
1081   }
1082   if(algebraic){variant=2;}
1083
1084   list P,T1,T2,T3,LL;
1085
1086   ideal J,K,H;
1087
1088//-----  If there is more than one processor available, we parallelize the  ----
1089//-----  main standard basis computations in positive characteristic        ----
1090
1091   if(n1 > 1)
1092   {
1093      ideal I_for_fork = I;
1094      export(I_for_fork);           // I available for each link
1095
1096//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
1097//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
1098
1099      for(i = 1; i <= n1; i++)
1100      {
1101         //link l(i) = "MPtcp:fork";
1102         link l(i) = "ssi:fork";
1103         open(l(i));
1104         if((variant == 1) || (variant == 3))
1105         {
1106            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1107                                                  eval(variant), eval(hi))));
1108         }
1109         if((variant == 2) || (variant == 4))
1110         {
1111            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1112                                                  eval(variant))));
1113         }
1114      }
1115
1116      int t = timer;
1117      if((variant == 1) || (variant == 3))
1118      {
1119         P = modpStd(I_for_fork, L[1], variant, hi);
1120      }
1121      if((variant == 2) || (variant == 4))
1122      {
1123         P = modpStd(I_for_fork, L[1], variant);
1124      }
1125      t = timer - t;
1126      if(t > 60) { t = 60; }
1127      int i_sleep = system("sh", "sleep "+string(t));
1128      T1[1] = P[1];
1129      T2[1] = bigint(P[2]);
1130      index++;
1131
1132      j = j + n1 + 1;
1133   }
1134
1135//--------------  Main standard basis computations in positive  ----------------
1136//----------------------  characteristic start here  ---------------------------
1137
1138   list arguments_farey, results_farey;
1139
1140   while(1)
1141   {
1142      tt = timer; rt = rtimer;
1143
1144      if(printlevel >= 10) { "size(L) = "+string(size(L)); }
1145
1146      if(n1 > 1)
1147      {
1148         while(j <= size(L) + 1)
1149         {
1150            for(i = 1; i <= n1; i++)
1151            {
1152               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
1153               if(status(l(i), "read", "ready"))
1154               {
1155                  //--- read the result from l(i) ---
1156                  P = read(l(i));
1157                  T1[index] = P[1];
1158                  T2[index] = bigint(P[2]);
1159                  index++;
1160
1161                  if(j <= size(L))
1162                  {
1163                     if((variant == 1) || (variant == 3))
1164                     {
1165                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
1166                                                  eval(variant), eval(hi))));
1167                        j++;
1168                     }
1169                     if((variant == 2) || (variant == 4))
1170                     {
1171                        write(l(i), quote(modpStd(I_for_fork,
1172                                                  eval(L[j]), eval(variant))));
1173                        j++;
1174                     }
1175                  }
1176                  else
1177                  {
1178                     k++;
1179                     close(l(i));
1180                  }
1181               }
1182            }
1183            //--- k describes the number of closed links ---
1184            if(k == n1)
1185            {
1186               j++;
1187            }
1188            i_sleep = system("sh", "sleep "+string(t));
1189         }
1190      }
1191      else
1192      {
1193         while(j <= size(L))
1194         {
1195            if((variant == 1) || (variant == 3))
1196            {
1197               P = modpStd(I, L[j], variant, hi);
1198            }
1199            if((variant == 2) || (variant == 4))
1200            {
1201               P = modpStd(I, L[j], variant);
1202            }
1203
1204            T1[index] = P[1];
1205            T2[index] = bigint(P[2]);
1206            index++;
1207            j++;
1208         }
1209      }
1210
1211      if(printlevel >= 10)
1212      {
1213         "CPU-time for computing list is "+string(timer - tt)+" seconds.";
1214         "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
1215      }
1216
1217//------------------------  Delete unlucky primes  -----------------------------
1218//-------------  unlucky if and only if the leading ideal is wrong  ------------
1219
1220      LL = deleteUnluckyPrimes(T1,T2,h);
1221      T1 = LL[1];
1222      T2 = LL[2];
1223
1224//-------------------  Now all leading ideals are the same  --------------------
1225//-------------------  Lift results to basering via farey  ---------------------
1226
1227      tt = timer; rt = rtimer;
1228      N = T2[1];
1229      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
1230      H = chinrem(T1,T2);
1231      if(n1 == 1)
1232      {
1233         J = farey(H,N);
1234      }
1235      else
1236      {
1237         for(i = ncols(H); i > 0; i--)
1238         {
1239            arguments_farey[i] = list(ideal(H[i]), N);
1240         }
1241         results_farey = parallelWaitAll("farey", arguments_farey, 0, n1);
1242         for(i = ncols(H); i > 0; i--)
1243         {
1244            J[i] = results_farey[i][1];
1245         }
1246      }
1247      if(printlevel >= 10)
1248      {
1249         "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
1250         "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
1251      }
1252
1253//----------------  Test if we already have a standard basis of I --------------
1254
1255      tt = timer; rt = rtimer;
1256      if((variant == 1) || (variant == 3))
1257      {
1258         pTest = pTestSB(I,J,L,variant,hi);
1259      }
1260      if((variant == 2) || (variant == 4))
1261      {
1262         pTest = pTestSB(I,J,L,variant);
1263      }
1264
1265      if(printlevel >= 10)
1266      {
1267         "CPU-time for pTest is "+string(timer - tt)+" seconds.";
1268         "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
1269      }
1270
1271      if(pTest)
1272      {
1273         if(printlevel >= 10)
1274         {
1275            "CPU-time for computation without final tests is "
1276            +string(timer - TT)+" seconds.";
1277            "Real-time for computation without final tests is "
1278            +string(rtimer - RT)+" seconds.";
1279         }
1280
1281         attrib(J,"isSB",1);
1282
1283         if(exactness == 0)
1284         {
1285            option(set, opt);
1286            if(n1 > 1) { kill I_for_fork; }
1287            return(J);
1288         }
1289
1290         if(exactness == 1)
1291         {
1292            tt = timer; rt = rtimer;
1293            sizeTest = 1 - isIncluded(I,J,n1);
1294
1295            if(printlevel >= 10)
1296            {
1297               "CPU-time for checking if I subset <G> is "
1298               +string(timer - tt)+" seconds.";
1299               "Real-time for checking if I subset <G> is "
1300               +string(rtimer - rt)+" seconds.";
1301            }
1302
1303            if(sizeTest == 0)
1304            {
1305               tt = timer; rt = rtimer;
1306               K = std(J);
1307
1308               if(printlevel >= 10)
1309               {
1310                  "CPU-time for last std-computation is "
1311                  +string(timer - tt)+" seconds.";
1312                  "Real-time for last std-computation is "
1313                  +string(rtimer - rt)+" seconds.";
1314               }
1315
1316               if(size(reduce(K,J)) == 0)
1317               {
1318                  option(set, opt);
1319                  if(n1 > 1) { kill I_for_fork; }
1320                  return(J);
1321               }
1322            }
1323         }
1324      }
1325
1326//--------------  We do not already have a standard basis of I  ----------------
1327//-----------  Therefore do the main computation for more primes  --------------
1328
1329      T1 = H;
1330      T2 = N;
1331      index = 2;
1332
1333      j = size(L) + 1;
1334      tt = timer; rt = rtimer;
1335      L = primeList(I,n3,L,n1);
1336      if(printlevel >= 10)
1337      {
1338         "CPU-time for primeList: "+string(timer-tt)+" seconds.";
1339         "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
1340      }
1341
1342      if(n1 > 1)
1343      {
1344         for(i = 1; i <= n1; i++)
1345         {
1346            open(l(i));
1347            if((variant == 1) || (variant == 3))
1348            {
1349               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1350                                                     eval(variant), eval(hi))));
1351            }
1352            if((variant == 2) || (variant == 4))
1353            {
1354               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1355                                                     eval(variant))));
1356            }
1357         }
1358         j = j + n1;
1359         k = 0;
1360      }
1361   }
1362}
1363example
1364{ "EXAMPLE:"; echo = 2;
1365   ring R1 = 0,(x,y,z,t),dp;
1366   ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
1367   ideal J = modStd(I);
1368   J;
1369   I = homog(I,t);
1370   J = modStd(I);
1371   J;
1372
1373   ring R2 = 0,(x,y,z),ds;
1374   ideal I = jacob(x5+y6+z7+xyz);
1375   ideal J1 = modStd(I,1,0);
1376   J1;
1377
1378   ring R3 = 0,x(1..4),lp;
1379   ideal I = cyclic(4);
1380   ideal J1 = modStd(I,1);
1381   ideal J2 = modStd(I,1,0);
1382   size(reduce(J1,J2));
1383   size(reduce(J2,J1));
1384}
1385
1386////////////////////////////////////////////////////////////////////////////////
1387
1388proc modS(ideal I, list L, list #)
1389"USAGE:  modS(I,L); I ideal, L intvec of primes
1390         if size(#)>0 std is used instead of groebner
1391RETURN:  an ideal which is with high probability a standard basis
1392NOTE:    This procedure is designed for fast experiments.
1393         It is not tested whether the result is a standard basis.
1394         It is not tested whether the result generates I.
1395EXAMPLE: example modS; shows an example
1396"
1397{
1398   int j;
1399   bigint N = 1;
1400   def R0 = basering;
1401   ideal J;
1402   list T;
1403   list rl = ringlist(R0);
1404   if((npars(R0)>0) || (rl[1]>0))
1405   {
1406      ERROR("Characteristic of basering should be zero.");
1407   }
1408   for(j = 1; j <= size(L); j++)
1409   {
1410      N = N*L[j];
1411      rl[1] = L[j];
1412      def @r = ring(rl);
1413      setring @r;
1414      ideal I = fetch(R0,I);
1415      if(size(#) > 0)
1416      {
1417         I = std(I);
1418      }
1419      else
1420      {
1421         I = groebner(I);
1422      }
1423      setring R0;
1424      T[j] = fetch(@r,I);
1425      kill @r;
1426   }
1427   L = deleteUnluckyPrimes(T,L,homog(I));
1428   // unlucky if and only if the leading ideal is wrong
1429   J = farey(chinrem(L[1],L[2]),N);
1430   attrib(J,"isSB",1);
1431   return(J);
1432}
1433example
1434{ "EXAMPLE:"; echo = 2;
1435   list L = 3,5,11,13,181,32003;
1436   ring r = 0,(x,y,z,t),dp;
1437   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1438   I = homog(I,t);
1439   ideal J = modS(I,L);
1440   J;
1441}
1442
1443////////////////////////////////////////////////////////////////////////////////
1444
1445proc modHenselStd(ideal I, list #)
1446"USAGE:  modHenselStd(I);
1447RETURN:  a standard basis of I;
1448NOTE:    The procedure computes a standard basis of I (over the rational
1449         numbers) by using  modular computations and Hensellifting.
1450         For further experiments see procedure modS.
1451EXAMPLE: example modHenselStd; shows an example
1452"
1453{
1454   int i,j;
1455
1456   bigint p = 2134567879;
1457   if(size(#)!=0) { p=#[1]; }
1458   while(!primeTest(I,p))
1459   {
1460      p = prime(random(2000000000,2134567879));
1461   }
1462
1463   def R = basering;
1464   module F,PrevG,PrevZ,Z2;
1465   ideal testG,testG1,G1,G2,G3,Gp;
1466   list L = p;
1467   list rl = ringlist(R);
1468   rl[1] = int(p);
1469
1470   def S = ring(rl);
1471   setring S;
1472   intvec opt = option(get);
1473   option(redSB);
1474   module Z,M,Z2;
1475   ideal I = imap(R,I);
1476   ideal Gp,G1,G2,G3;
1477   Gp,Z = liftstd1(I);
1478   attrib(Gp,"isSB",1);
1479   module ZZ = syz(I);
1480   attrib(ZZ,"isSB",1);
1481   Z = reduce(Z,ZZ);
1482
1483   setring R;
1484   Gp = imap(S,Gp);
1485   PrevZ = imap(S,Z);
1486   PrevG = module(Gp);
1487   F = module(I);
1488   testG = farey(Gp,p);
1489   attrib(testG,"isSB",1);
1490   while(1)
1491   {
1492      i++;
1493      G1 = ideal(1/(p^i) * sum_id(F*PrevZ,(-1)*PrevG));
1494      setring S;
1495      G1 = imap(R,G1);
1496      G2 = reduce(G1,Gp);
1497      G3 = sum_id(G1,(-1)*G2);
1498      M = lift(Gp,G3);
1499      Z2 = (-1)*Z*M;
1500
1501      setring R;
1502      G2 = imap(S,G2);
1503      Z2 = imap(S,Z2);
1504      PrevG = sum_id(PrevG, module(p^i*G2));
1505      PrevZ = sum_id(PrevZ, multiply(poly(p^i),Z2));
1506      testG1 = farey(ideal(PrevG),p^(i+1));
1507      attrib(testG1,"isSB",1);
1508      if(size(reduce(testG1,testG)) == 0)
1509      {
1510         if(size(reduce(I,testG1)) == 0)       // I is in testG1
1511         {
1512            if(pTestSB(I,testG1,L,2))
1513            {
1514               G3 = std(testG1);               // testG1 is SB
1515               if(size(reduce(G3,testG1)) == 0)
1516               {
1517                  option(set, opt);
1518                  return(G3);
1519               }
1520            }
1521         }
1522      }
1523      testG = testG1;
1524      attrib(testG,"isSB",1);
1525   }
1526}
1527example
1528{ "EXAMPLE:"; echo = 2;
1529   ring r = 0,(x,y,z),dp;
1530   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1531   ideal J = modHenselStd(I);
1532   J;
1533}
1534
1535////////////////////////////////////////////////////////////////////////////////
1536
1537static proc sum_id(list #)
1538{
1539   if(typeof(#[1])=="ideal")
1540   {
1541      ideal M;
1542   }
1543   else
1544   {
1545      module M;
1546   }
1547
1548   int i;
1549   for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; }
1550   return(M);
1551}
1552
1553////////////////////////////////////////////////////////////////////////////////
1554
1555static proc multiply(poly p, list #)
1556{
1557   if(typeof(#[1])=="ideal")
1558   {
1559      ideal M;
1560   }
1561   else
1562   {
1563      module M;
1564   }
1565
1566   int i;
1567   for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; }
1568   return(M);
1569}
1570
1571
1572////////////////////////////// further examples ////////////////////////////////
1573
1574/*
1575ring r = 0, (x,y,z), lp;
1576poly s1 = 5x3y2z+3y3x2z+7xy2z2;
1577poly s2 = 3xy2z2+x5+11y2z2;
1578poly s3 = 4xyz+7x3+12y3+1;
1579poly s4 = 3x3-4y3+yz2;
1580ideal i =  s1, s2, s3, s4;
1581
1582ring r = 0, (x,y,z), lp;
1583poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
1584poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
1585poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
1586poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
1587ideal i =  s1, s2, s3, s4;
1588
1589ring r = 0, (x,y,z), lp;
1590poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
1591poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
1592poly s3 = 8x3 + 12y3 + xz2 + 3;
1593poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
1594ideal i = s1, s2, s3, s4;
1595
1596int n = 6;
1597ring r = 0,(x(1..n)),lp;
1598ideal i = cyclic(n);
1599ring s = 0, (x(1..n),t), lp;
1600ideal i = imap(r,i);
1601i = homog(i,t);
1602
1603ring r = 0, (x(1..4),s), (dp(4),dp);
1604poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
1605poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
1606poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
1607poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4);
1608poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
1609ideal i =  s1, s2, s3, s4, s5;
1610
1611ring r = 0, (x,y,z), ds;
1612int a = 16;
1613int b = 15;
1614int c = 4;
1615int t = 1;
1616poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1617         +x^(c-2)*y^c*(y2+t*x)^2;
1618ideal i = jacob(f);
1619
1620ring r = 0, (x,y,z), ds;
1621int a = 25;
1622int b = 25;
1623int c = 5;
1624int t = 1;
1625poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1626         +x^(c-2)*y^c*(y2+t*x)^2;
1627ideal i = jacob(f),f;
1628
1629ring r = 0, (x,y,z), ds;
1630int a = 10;
1631poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
1632ideal i = jacob(f);
1633
1634ring r = 0, (x,y,z), ds;
1635int a = 6;
1636int b = 8;
1637int c = 10;
1638int alpha = 5;
1639int beta = 5;
1640int t = 1;
1641poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
1642         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
1643ideal i = jacob(f);
1644*/
1645
Note: See TracBrowser for help on using the repository browser.