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
Line 
1////////////////////////////////////////////////////////////////////////////////
2version="$Id: modstd.lib 14375 2011-08-23 09:29:47Z steidel $";
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
33proc mixedTest()
34"USAGE:  mixedTest();
35RETURN:  1 if ordering of basering is mixed, 0 else
36EXAMPLE: example mixedTest(); shows an example
37"
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}
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}
63
64////////////////////////////////////////////////////////////////////////////////
65
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,
77@*       0 if there is an element f in I which does not reduce to 0 w.r.t. J.
78EXAMPLE: example isIncluded; shows an example
79"
80{
81   def R = basering;
82   setring R;
83
84   attrib(J,"isSB",1);
85   int i,j,k;
86
87   if(size(#) > 0)
88   {
89      int n = #[1];
90      if(n >= ncols(I)) { n = ncols(I); }
91      if(n > 1)
92      {
93         for(i = 1; i <= n - 1; i++)
94         {
95            //link l(i) = "MPtcp:fork";
96            link l(i) = "ssi:fork";
97            open(l(i));
98
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"))
122               {
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
225   if(printlevel >= 11)
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
245      if(printlevel >= 11)
246      {
247         "isIncluded(J,K) takes "+string(timer - t)+" seconds";
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
268RETURN:  lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #),
269         lT ideal
270NOTE:    - if ho = 1, the polynomials in T are homogeneous, else ho = 0,
271@*       - lT is prevalent, i.e. the most appearing leading ideal in T
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      }
293      if(c > size(T) div 2){ lT = cK; }
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         {
313            if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; }
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
377proc primeTest(def II, bigint p)
378{
379   if(typeof(II) == "string")
380   {
381      execute("ideal I = "+II+";");
382   }
383   else
384   {
385      ideal I = II;
386   }
387
388   int i,j;
389   poly f;
390   number cnt;
391   for(i = 1; i <= size(I); i++)
392   {
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--)
399      {
400         if((leadcoef(f[j]) mod p) == 0) { return(0); }
401      }
402   }
403   return(1);
404}
405
406////////////////////////////////////////////////////////////////////////////////
407
408proc primeList(ideal I, int n, list #)
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
412         < L[size(L)] union with L) such that none of these primes divides any
413         coefficient occuring in I
414NOTE:    The number of cores to use can be defined by ncores, default is 1.
415EXAMPLE: example primList; shows an example
416"
417{
418   intvec L;
419   int i,p;
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   
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"); }
461   if(ncores == 1)
462   {
463      for(i = 2; i <= n; i++)
464      {
465         p = prime(p-1);
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];
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
523static proc liftstd1(ideal I)
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
545   intvec opt = option(get);
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   {
569      //"==========================================================";
570      //"WARNING: Algorithm is not applicable if ordering is mixed.";
571      //"==========================================================";
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;
580   option(set, opt);
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
595   ring S = 32003, x(1..5), lp;
596   ideal I = cyclic(5);
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.
608RETURN:  ideal - a standard basis of I mod p, integer - p
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
614@*       - variant = 1: std(.,#[1]) resp. groebner,
615@*       - variant = 2: groebner,
616@*       - variant = 3: homog. - std(.,#[1]) resp. groebner - dehomog.,
617@*       - variant = 4: fglm.
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
630   if(variant == 1)
631   {
632      if(size(#) > 0)
633      {
634         i = std(i, #[1]);
635      }
636      else
637      {
638         i = groebner(i);
639      }
640   }
641
642   if(variant == 2)
643   {
644      i = groebner(i);
645   }
646
647   if(variant == 3)
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   }
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   }
702
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
724@*       - #[1] is the number of available processors for the computation,
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
732         numbers) by using  modular methods.
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
735         with high probability.
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.
739@*       - variant = 1, if I is homogeneous,
740@*       - variant = 2, if I is not homogeneous, 1-block-ordering,
741@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
742                        > 1 block),
743@*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
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   {
754      ERROR("Characteristic of basering should be zero, basering should
755             have no parameters.");
756   }
757
758   int index = 1;
759   int i,k,c;
760   int j = 1;
761   int pTest, sizeTest;
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];
772         int exactness = 1;
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         }
783      }
784      if(size(#) == 2)
785      {
786         int n1 = #[1];
787         int exactness = #[2];
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         }
798      }
799      if(size(#) == 4)
800      {
801         int n1 = #[1];
802         int exactness = #[2];
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         }
819      }
820   }
821   else
822   {
823      int n1 = 1;
824      int exactness = 1;
825      int n2 = 10;
826      int n3 = 10;
827   }
828
829   if(printlevel >= 10)
830   {
831      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
832       +", exactness = "+string(exactness);
833   }
834
835//-------------------------  Save current options  -----------------------------
836   intvec opt = option(get);
837
838   option(redSB);
839
840//--------------------  Initialize the list of primes  -------------------------
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   }
849   L[5] = prime(random(an,en));
850
851//---------------------  Decide which variant to take  -------------------------
852   int variant;
853   int h = homog(I);
854
855   tt = timer;
856   rt = rtimer;
857
858   if(!mixedTest())
859   {
860      if(h)
861      {
862         variant = 1;
863         if(printlevel >= 10) { "variant = 1"; }
864         
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         {
882            variant = 2;
883            if(printlevel >= 10) { "variant = 2"; }
884         }
885         else
886         {
887            string order;
888            if(system("nblocks") <= 2)
889            {
890               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
891                                       + find(ordstr_R0, "rp") <= 0)
892               {
893                  order = "simple";
894               }
895            }
896
897            if((order == "simple") || (size(rl) > 4))
898            {
899               variant = 2;
900               if(printlevel >= 10) { "variant = 2"; }
901            }
902            else
903            {
904               rl[1] = L[5];
905               def @r = ring(rl);
906               setring @r;
907               
908               def @s = changeord("dp");
909               setring @s;
910               ideal I = std(fetch(R0,I));
911               if(dim(I) == 0)
912               {
913                  variant = 4;
914                  if(printlevel >= 10) { "variant = 4"; }
915               }
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;
928
929                  list hiRi = hilbRing(fetch(R0,I),w);
930                  intvec W = hiRi[2];
931                  @s = hiRi[1];
932                  setring @s;
933
934                  Id(1) = std(Id(1));
935                  intvec hi = hilb(Id(1), 1, W);
936               }
937               
938               setring R0;
939               kill @r,@s;
940            }
941         }
942      }
943   }
944   else
945   {
946      if(exactness == 1) { return(groebner(I)); }
947      if(h)
948      {
949         variant = 1;
950         if(printlevel >= 10) { "variant = 1"; }
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
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            {
976               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
977                                       + find(ordstr_R0, "rp") <= 0)
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"; }
992
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   }
1018
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      {
1036         //link l(i) = "MPtcp:fork";
1037         link l(i) = "ssi:fork";
1038         open(l(i));
1039         if((variant == 1) || (variant == 3))
1040         {
1041            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1042                                                  eval(variant), eval(hi))));
1043         }
1044         if((variant == 2) || (variant == 4))
1045         {
1046            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1047                                                  eval(variant))));
1048         }
1049      }
1050
1051      int t = timer;
1052      if((variant == 1) || (variant == 3))
1053      {
1054         P = modpStd(I_for_fork, L[1], variant, hi);
1055      }
1056      if((variant == 2) || (variant == 4))
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
1073   list arguments_farey, results_farey;
1074
1075   while(1)
1076   {
1077      tt = timer; rt = rtimer;
1078     
1079      if(printlevel >= 10) { "size(L) = "+string(size(L)); }
1080         
1081      if(n1 > 1)
1082      {
1083         while(j <= size(L) + 1)
1084         {
1085            for(i = 1; i <= n1; i++)
1086            {
1087               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
1088               if(status(l(i), "read", "ready"))
1089               {
1090                  //--- read the result from l(i) ---
1091                  P = read(l(i));
1092                  T1[index] = P[1];
1093                  T2[index] = bigint(P[2]);
1094                  index++;
1095
1096                  if(j <= size(L))
1097                  {
1098                     if((variant == 1) || (variant == 3))
1099                     {
1100                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
1101                                                  eval(variant), eval(hi))));
1102                        j++;
1103                     }
1104                     if((variant == 2) || (variant == 4))
1105                     {
1106                        write(l(i), quote(modpStd(I_for_fork,
1107                                                  eval(L[j]), eval(variant))));
1108                        j++;
1109                     }
1110                  }
1111                  else
1112                  {
1113                     k++;
1114                     close(l(i));
1115                  }
1116               }
1117            }
1118            //--- k describes the number of closed links ---
1119            if(k == n1)
1120            {
1121               j++;
1122            }
1123            i_sleep = system("sh", "sleep "+string(t));
1124         }
1125      }
1126      else
1127      {
1128         while(j <= size(L))
1129         {
1130            if((variant == 1) || (variant == 3))
1131            {
1132               P = modpStd(I, L[j], variant, hi);
1133            }
1134            if((variant == 2) || (variant == 4))
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
1155      LL = deleteUnluckyPrimes(T1,T2,h);
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
1162      tt = timer; rt = rtimer;
1163      N = T2[1];
1164      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
1165      H = chinrem(T1,T2);
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      }
1188
1189//----------------  Test if we already have a standard basis of I --------------
1190
1191      tt = timer; rt = rtimer;
1192      if((variant == 1) || (variant == 3))
1193      {
1194         pTest = pTestSB(I,J,L,variant,hi);
1195      }
1196      if((variant == 2) || (variant == 4))
1197      {
1198         pTest = pTestSB(I,J,L,variant);
1199      }
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      }
1206
1207      if(pTest)
1208      {
1209         if(printlevel >= 10)
1210         {
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.";
1215         }
1216
1217         attrib(J,"isSB",1);
1218         
1219         if(exactness == 0)
1220         {
1221            option(set, opt);
1222            if(n1 > 1) { kill I_for_fork; }
1223            return(J);
1224         }
1225         
1226         if(exactness == 1)
1227         {
1228            tt = timer; rt = rtimer;
1229            sizeTest = 1 - isIncluded(I,J,n1);
1230
1231            if(printlevel >= 10)
1232            {
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.";
1237            }
1238
1239            if(sizeTest == 0)
1240            {
1241               tt = timer; rt = rtimer;
1242               K = std(J);
1243
1244               if(printlevel >= 10)
1245               {
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.";
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;
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      }
1277
1278      if(n1 > 1)
1279      {
1280         for(i = 1; i <= n1; i++)
1281         {
1282            open(l(i));
1283            if((variant == 1) || (variant == 3))
1284            {
1285               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1286                                                     eval(variant), eval(hi))));
1287            }
1288            if((variant == 2) || (variant == 4))
1289            {
1290               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1291                                                     eval(variant))));
1292            }
1293         }
1294         j = j + n1;
1295         k = 0;
1296      }
1297   }
1298}
1299example
1300{ "EXAMPLE:"; echo = 2;
1301   ring R1 = 0,(x,y,z,t),dp;
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
1309   ring R2 = 0,(x,y,z),ds;
1310   ideal I = jacob(x5+y6+z7+xyz);
1311   ideal J1 = modStd(I,1,0);
1312   J1;
1313
1314   ring R3 = 0,x(1..4),lp;
1315   ideal I = cyclic(4);
1316   ideal J1 = modStd(I,1);
1317   ideal J2 = modStd(I,1,0);
1318   size(reduce(J1,J2));
1319   size(reduce(J2,J1));
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);
1340   if((npars(R0)>0) || (rl[1]>0))
1341   {
1342      ERROR("Characteristic of basering should be zero.");
1343   }
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   }
1363   L = deleteUnluckyPrimes(T,L,homog(I));
1364   // unlucky if and only if the leading ideal is wrong
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;
1408   intvec opt = option(get);
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               {
1453                  option(set, opt);
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;
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;
1554ideal i = jacob(f);
1555
1556ring r = 0, (x,y,z), ds;
1557int a = 25;
1558int b = 25;
1559int c = 5;
1560int t = 1;
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;
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;
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;
1579ideal i = jacob(f);
1580*/
1581
Note: See TracBrowser for help on using the repository browser.