source: git/Singular/LIB/modstd.lib @ 7f30e2

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