source: git/Singular/LIB/modstd.lib @ 5abb0e7

jengelh-datetimespielwiese
Last change on this file since 5abb0e7 was 5abb0e7, checked in by Andreas Steenpass <steenpass@…>, 11 years ago
fix: erase zero generators in modstd::primeTest()
  • Property mode set to 100644
File size: 38.4 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   I = simplify(I, 2);   // erase zero generators
389
390   int i,j;
391   poly f;
392   number cnt;
393   for(i = 1; i <= size(I); i++)
394   {
395      f = cleardenom(I[i]);
396      if(f == 0) { return(0); }
397      cnt = leadcoef(I[i])/leadcoef(f);
398      if((numerator(cnt) mod p) == 0) { return(0); }
399      if((denominator(cnt) mod p) == 0) { return(0); }
400      for(j = size(f); j > 0; j--)
401      {
402         if((leadcoef(f[j]) mod p) == 0) { return(0); }
403      }
404   }
405   return(1);
406}
407
408////////////////////////////////////////////////////////////////////////////////
409
410proc primeList(ideal I, int n, list #)
411"USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
412         n integer
413RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
414         < L[size(L)] union with L) such that none of these primes divides any
415         coefficient occuring in I
416NOTE:    The number of cores to use can be defined by ncores, default is 1.
417EXAMPLE: example primList; shows an example
418"
419{
420   intvec L;
421   int i,p;
422   int ncores = 1;
423   
424//-----------------  Initialize optional parameter ncores  ---------------------
425   if(size(#) > 0)
426   {
427      if(size(#) == 1)
428      {
429         if(typeof(#[1]) == "int")
430         {
431            ncores = #[1];
432            # = list();
433         }
434      }
435      else
436      {
437         ncores = #[2];
438      }
439   }
440   
441   if(size(#) == 0)
442   {
443      p = 2147483647;
444      while(!primeTest(I,p))
445      {
446         p = prime(p-1);
447         if(p == 2) { ERROR("no more primes"); }
448      }
449      L[1] = p;
450   }
451   else
452   {
453      L = #[1];
454      p = prime(L[size(L)]-1);
455      while(!primeTest(I,p))
456      {
457         p = prime(p-1);
458         if(p == 2) { ERROR("no more primes"); }
459      }
460      L[size(L)+1] = p;
461   }
462   if(p == 2) { ERROR("no more primes"); }
463   if(ncores == 1)
464   {
465      for(i = 2; i <= n; i++)
466      {
467         p = prime(p-1);
468         while(!primeTest(I,p))
469         {
470            p = prime(p-1);
471            if(p == 2) { ERROR("no more primes"); }
472         }
473         L[size(L)+1] = p;
474      }
475   }
476   else
477   {
478      int neededSize = size(L)+n-1;;
479      list parallelResults;
480      list arguments;
481      int neededPrimes = neededSize-size(L);
482      while(neededPrimes > 0)
483      {
484         arguments = list();
485         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
486            *ncores; i > 0; i--)
487         {
488            p = prime(p-1);
489            if(p == 2) { ERROR("no more primes"); }
490            arguments[i] = list("I", p);
491         }
492         parallelResults = parallelWaitAll("primeTest", arguments,
493            list(list(list(ncores))));
494         for(i = size(arguments); i > 0; i--)
495         {
496            if(parallelResults[i])
497            {
498               L[size(L)+1] = arguments[i][2];
499            }
500         }
501         neededPrimes = neededSize-size(L);
502      }
503      if(size(L) > neededSize)
504      {
505         L = L[1..neededSize];
506      }
507   }
508   return(L);
509}
510example
511{ "EXAMPLE:"; echo = 2;
512   ring r = 0,(x,y,z),dp;
513   ideal I = 2147483647x+y, z-181;
514   intvec L = primeList(I,10);
515   size(L);
516   L[1];
517   L[size(L)];
518   L = primeList(I,5,L);
519   size(L);
520   L[size(L)];
521}
522
523////////////////////////////////////////////////////////////////////////////////
524
525static proc liftstd1(ideal I)
526{
527   def R = basering;
528   list rl = ringlist(R);
529   list ordl = rl[3];
530
531   int i;
532   for(i = 1; i <= size(ordl); i++)
533   {
534      if((ordl[i][1] == "C") || (ordl[i][1] == "c"))
535      {
536         ordl = delete(ordl, i);
537         break;
538      }
539   }
540
541   ordl = insert(ordl, list("c", 0));
542   rl[3] = ordl;
543   def newR = ring(rl);
544   setring newR;
545   ideal I = imap(R,I);
546
547   intvec opt = option(get);
548   option(none);
549   option(prompt);
550
551   module M;
552   for(i = 1; i <= size(I); i++)
553   {
554      M = M + module(I[i]*gen(1) + gen(i+1));
555      M = M + module(gen(i+1));
556   }
557
558   module sM = std(M);
559
560   ideal sI;
561   if(attrib(R,"global"))
562   {
563      for(i = size(I)+1; i <= size(sM); i++)
564      {
565         sI[size(sI)+1] = sM[i][1];
566      }
567      matrix T = submat(sM,2..nrows(sM),size(I)+1..ncols(sM));
568   }
569   else
570   {
571      //"==========================================================";
572      //"WARNING: Algorithm is not applicable if ordering is mixed.";
573      //"==========================================================";
574      for(i = 1; i <= size(sM)-size(I); i++)
575      {
576         sI[size(sI)+1] = sM[i][1];
577      }
578      matrix T = submat(sM,2..nrows(sM),1..ncols(sM)-size(I));
579   }
580
581   setring R;
582   option(set, opt);
583   return(imap(newR,sI),imap(newR,T));
584}
585example
586{ "EXAMPLE:"; echo = 2;
587   ring R = 0,(x,y,z),dp;
588   poly f = x3+y7+z2+xyz;
589   ideal i = jacob(f);
590   matrix T;
591   ideal sm = liftstd(i,T);
592   sm;
593   print(T);
594   matrix(sm) - matrix(i)*T;
595
596
597   ring S = 32003, x(1..5), lp;
598   ideal I = cyclic(5);
599   ideal sI;
600   matrix T;
601   sI,T = liftstd1(I);
602   matrix(sI) - matrix(I)*T;
603}
604
605////////////////////////////////////////////////////////////////////////////////
606
607proc modpStd(ideal I, int p, int variant, list #)
608"USAGE:  modpStd(I,p,variant,#); I ideal, p integer, variant integer
609ASSUME:  If size(#) > 0, then #[1] is an intvec describing the Hilbert series.
610RETURN:  ideal - a standard basis of I mod p, integer - p
611NOTE:    The procedure computes a standard basis of the ideal I modulo p and
612         fetches the result to the basering. If size(#) > 0 the Hilbert driven
613         standard basis computation std(.,#[1]) is used instead of groebner.
614         The standard basis computation modulo p does also vary depending on the
615         integer variant, namely
616@*       - variant = 1: std(.,#[1]) resp. groebner,
617@*       - variant = 2: groebner,
618@*       - variant = 3: homog. - std(.,#[1]) resp. groebner - dehomog.,
619@*       - variant = 4: fglm.
620EXAMPLE: example modpStd; shows an example
621"
622{
623   def R0 = basering;
624   list rl = ringlist(R0);
625   rl[1] = p;
626   def @r = ring(rl);
627   setring @r;
628   ideal i = fetch(R0,I);
629
630   option(redSB);
631
632   if(variant == 1)
633   {
634      if(size(#) > 0)
635      {
636         i = std(i, #[1]);
637      }
638      else
639      {
640         i = groebner(i);
641      }
642   }
643
644   if(variant == 2)
645   {
646      i = groebner(i);
647   }
648
649   if(variant == 3)
650   {
651      list rl = ringlist(@r);
652      int nvar@r = nvars(@r);
653
654      int k;
655      intvec w;
656      for(k = 1; k <= nvar@r; k++)
657      {
658         w[k] = deg(var(k));
659      }
660      w[nvar@r + 1] = 1;
661
662      rl[2][nvar@r + 1] = "homvar";
663      rl[3][2][2] = w;
664
665      def HomR = ring(rl);
666      setring HomR;
667      ideal i = imap(@r, i);
668      i = homog(i, homvar);
669
670      if(size(#) > 0)
671      {
672         if(w == 1)
673         {
674            i = std(i, #[1]);
675         }
676         else
677         {
678            i = std(i, #[1], w);
679         }
680      }
681      else
682      {
683         i = groebner(i);
684      }
685
686      i = subst(i, homvar, 1);
687      i = simplify(i, 34);
688
689      setring @r;
690      i = imap(HomR, i);
691      i = interred(i);
692      kill HomR;
693   }
694   
695   if(variant == 4)
696   {
697      def R1 = changeord("dp");
698      setring R1;
699      ideal i = fetch(@r,i);
700      i = std(i);
701      setring @r;
702      i = fglm(R1,i);
703   }
704
705   setring R0;
706   return(list(fetch(@r,i),p));
707}
708example
709{ "EXAMPLE:"; echo = 2;
710   ring r = 0, x(1..4), dp;
711   ideal I = cyclic(4);
712   int p = 181;
713   list P = modpStd(I,p,5);
714   P;
715
716   int q = 32003;
717   list Q = modpStd(I,q,2);
718   Q;
719}
720
721////////////////////////////// main procedures /////////////////////////////////
722
723proc modStd(ideal I, list #)
724"USAGE:  modStd(I); I ideal
725ASSUME:  If size(#) > 0, then # contains either 1, 2 or 4 integers such that
726@*       - #[1] is the number of available processors for the computation,
727@*       - #[2] is an optional parameter for the exactness of the computation,
728                if #[2] = 1, the procedure computes a standard basis for sure,
729@*       - #[3] is the number of primes until the first lifting,
730@*       - #[4] is the constant number of primes between two liftings until
731           the computation stops.
732RETURN:  a standard basis of I if no warning appears;
733NOTE:    The procedure computes a standard basis of I (over the rational
734         numbers) by using  modular methods.
735         By default the procedure computes a standard basis of I for sure, but
736         if the optional parameter #[2] = 0, it computes a standard basis of I
737         with high probability.
738         The procedure distinguishes between different variants for the standard
739         basis computation in positive characteristic depending on the ordering
740         of the basering, the parameter #[2] and if the ideal I is homogeneous.
741@*       - variant = 1, if I is homogeneous,
742@*       - variant = 2, if I is not homogeneous, 1-block-ordering,
743@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
744                        > 1 block),
745@*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
746EXAMPLE: example modStd; shows an example
747"
748{
749   int TT = timer;
750   int RT = rtimer;
751
752   def R0 = basering;
753   list rl = ringlist(R0);
754   if((npars(R0) > 0) || (rl[1] > 0))
755   {
756      ERROR("Characteristic of basering should be zero, basering should
757             have no parameters.");
758   }
759
760   int index = 1;
761   int i,k,c;
762   int j = 1;
763   int pTest, sizeTest;
764   int en = 2134567879;
765   int an = 1000000000;
766   bigint N;
767
768//--------------------  Initialize optional parameters  ------------------------
769   if(size(#) > 0)
770   {
771      if(size(#) == 1)
772      {
773         int n1 = #[1];
774         int exactness = 1;
775         if(n1 >= 10)
776         {
777            int n2 = n1 + 1;
778            int n3 = n1;
779         }
780         else
781         {
782            int n2 = 10;
783            int n3 = 10;
784         }
785      }
786      if(size(#) == 2)
787      {
788         int n1 = #[1];
789         int exactness = #[2];
790         if(n1 >= 10)
791         {
792            int n2 = n1 + 1;
793            int n3 = n1;
794         }
795         else
796         {
797            int n2 = 10;
798            int n3 = 10;
799         }
800      }
801      if(size(#) == 4)
802      {
803         int n1 = #[1];
804         int exactness = #[2];
805         if(n1 >= #[3])
806         {
807            int n2 = n1 + 1;
808         }
809         else
810         {
811            int n2 = #[3];
812         }
813         if(n1 >= #[4])
814         {
815            int n3 = n1;
816         }
817         else
818         {
819            int n3 = #[4];
820         }
821      }
822   }
823   else
824   {
825      int n1 = 1;
826      int exactness = 1;
827      int n2 = 10;
828      int n3 = 10;
829   }
830
831   if(printlevel >= 10)
832   {
833      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
834       +", exactness = "+string(exactness);
835   }
836
837//-------------------------  Save current options  -----------------------------
838   intvec opt = option(get);
839
840   option(redSB);
841
842//--------------------  Initialize the list of primes  -------------------------
843   int tt = timer;
844   int rt = rtimer;
845   intvec L = primeList(I,n2,n1);
846   if(printlevel >= 10)
847   {
848      "CPU-time for primeList: "+string(timer-tt)+" seconds.";
849      "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
850   }
851   L[5] = prime(random(an,en));
852
853//---------------------  Decide which variant to take  -------------------------
854   int variant;
855   int h = homog(I);
856
857   tt = timer;
858   rt = rtimer;
859
860   if(!mixedTest())
861   {
862      if(h)
863      {
864         variant = 1;
865         if(printlevel >= 10) { "variant = 1"; }
866         
867         rl[1] = L[5];
868         def @r = ring(rl);
869         setring @r;
870         def @s = changeord("dp");
871         setring @s;
872         ideal I = std(fetch(R0,I));
873         intvec hi = hilb(I,1);
874         setring R0;
875         kill @r,@s;
876      }
877      else
878      {
879         string ordstr_R0 = ordstr(R0);
880         int neg = 1 - attrib(R0,"global");
881
882         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
883         {
884            variant = 2;
885            if(printlevel >= 10) { "variant = 2"; }
886         }
887         else
888         {
889            string order;
890            if(system("nblocks") <= 2)
891            {
892               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
893                                       + find(ordstr_R0, "rp") <= 0)
894               {
895                  order = "simple";
896               }
897            }
898
899            if((order == "simple") || (size(rl) > 4))
900            {
901               variant = 2;
902               if(printlevel >= 10) { "variant = 2"; }
903            }
904            else
905            {
906               rl[1] = L[5];
907               def @r = ring(rl);
908               setring @r;
909               
910               def @s = changeord("dp");
911               setring @s;
912               ideal I = std(fetch(R0,I));
913               if(dim(I) == 0)
914               {
915                  variant = 4;
916                  if(printlevel >= 10) { "variant = 4"; }
917               }
918               else
919               {
920                  variant = 3;
921                  if(printlevel >= 10) { "variant = 3"; }
922                 
923                  int nvar@r = nvars(@r);
924                  intvec w;
925                  for(i = 1; i <= nvar@r; i++)
926                  {
927                     w[i] = deg(var(i));
928                  }
929                  w[nvar@r + 1] = 1;
930
931                  list hiRi = hilbRing(fetch(R0,I),w);
932                  intvec W = hiRi[2];
933                  @s = hiRi[1];
934                  setring @s;
935
936                  Id(1) = std(Id(1));
937                  intvec hi = hilb(Id(1), 1, W);
938               }
939               
940               setring R0;
941               kill @r,@s;
942            }
943         }
944      }
945   }
946   else
947   {
948      if(exactness == 1) { return(groebner(I)); }
949      if(h)
950      {
951         variant = 1;
952         if(printlevel >= 10) { "variant = 1"; }
953         rl[1] = L[5];
954         def @r = ring(rl);
955         setring @r;
956         def @s = changeord("dp");
957         setring @s;
958         ideal I = std(fetch(R0,I));
959         intvec hi = hilb(I,1);
960         setring R0;
961         kill @r,@s;
962      }
963      else
964      {
965         string ordstr_R0 = ordstr(R0);
966         int neg = 1 - attrib(R0,"global");
967
968         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
969         {
970            variant = 2;
971            if(printlevel >= 10) { "variant = 2"; }
972         }
973         else
974         {
975            string order;
976            if(system("nblocks") <= 2)
977            {
978               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
979                                       + find(ordstr_R0, "rp") <= 0)
980               {
981                  order = "simple";
982               }
983            }
984
985            if((order == "simple") || (size(rl) > 4))
986            {
987               variant = 2;
988               if(printlevel >= 10) { "variant = 2"; }
989            }
990            else
991            {
992               variant = 3;
993               if(printlevel >= 10) { "variant = 3"; }
994
995               rl[1] = L[5];
996               def @r = ring(rl);
997               setring @r;
998               int nvar@r = nvars(@r);
999               intvec w;
1000               for(i = 1; i <= nvar@r; i++)
1001               {
1002                  w[i] = deg(var(i));
1003               }
1004               w[nvar@r + 1] = 1;
1005
1006               list hiRi = hilbRing(fetch(R0,I),w);
1007               intvec W = hiRi[2];
1008               def @s = hiRi[1];
1009               setring @s;
1010
1011               Id(1) = std(Id(1));
1012               intvec hi = hilb(Id(1), 1, W);
1013
1014               setring R0;
1015               kill @r,@s;
1016            }
1017         }
1018      }
1019   }
1020
1021   list P,T1,T2,T3,LL;
1022
1023   ideal J,K,H;
1024
1025//-----  If there is more than one processor available, we parallelize the  ----
1026//-----  main standard basis computations in positive characteristic        ----
1027
1028   if(n1 > 1)
1029   {
1030      ideal I_for_fork = I;
1031      export(I_for_fork);           // I available for each link
1032
1033//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
1034//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
1035
1036      for(i = 1; i <= n1; i++)
1037      {
1038         //link l(i) = "MPtcp:fork";
1039         link l(i) = "ssi:fork";
1040         open(l(i));
1041         if((variant == 1) || (variant == 3))
1042         {
1043            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1044                                                  eval(variant), eval(hi))));
1045         }
1046         if((variant == 2) || (variant == 4))
1047         {
1048            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1049                                                  eval(variant))));
1050         }
1051      }
1052
1053      int t = timer;
1054      if((variant == 1) || (variant == 3))
1055      {
1056         P = modpStd(I_for_fork, L[1], variant, hi);
1057      }
1058      if((variant == 2) || (variant == 4))
1059      {
1060         P = modpStd(I_for_fork, L[1], variant);
1061      }
1062      t = timer - t;
1063      if(t > 60) { t = 60; }
1064      int i_sleep = system("sh", "sleep "+string(t));
1065      T1[1] = P[1];
1066      T2[1] = bigint(P[2]);
1067      index++;
1068
1069      j = j + n1 + 1;
1070   }
1071
1072//--------------  Main standard basis computations in positive  ----------------
1073//----------------------  characteristic start here  ---------------------------
1074
1075   list arguments_farey, results_farey;
1076
1077   while(1)
1078   {
1079      tt = timer; rt = rtimer;
1080     
1081      if(printlevel >= 10) { "size(L) = "+string(size(L)); }
1082         
1083      if(n1 > 1)
1084      {
1085         while(j <= size(L) + 1)
1086         {
1087            for(i = 1; i <= n1; i++)
1088            {
1089               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
1090               if(status(l(i), "read", "ready"))
1091               {
1092                  //--- read the result from l(i) ---
1093                  P = read(l(i));
1094                  T1[index] = P[1];
1095                  T2[index] = bigint(P[2]);
1096                  index++;
1097
1098                  if(j <= size(L))
1099                  {
1100                     if((variant == 1) || (variant == 3))
1101                     {
1102                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
1103                                                  eval(variant), eval(hi))));
1104                        j++;
1105                     }
1106                     if((variant == 2) || (variant == 4))
1107                     {
1108                        write(l(i), quote(modpStd(I_for_fork,
1109                                                  eval(L[j]), eval(variant))));
1110                        j++;
1111                     }
1112                  }
1113                  else
1114                  {
1115                     k++;
1116                     close(l(i));
1117                  }
1118               }
1119            }
1120            //--- k describes the number of closed links ---
1121            if(k == n1)
1122            {
1123               j++;
1124            }
1125            i_sleep = system("sh", "sleep "+string(t));
1126         }
1127      }
1128      else
1129      {
1130         while(j <= size(L))
1131         {
1132            if((variant == 1) || (variant == 3))
1133            {
1134               P = modpStd(I, L[j], variant, hi);
1135            }
1136            if((variant == 2) || (variant == 4))
1137            {
1138               P = modpStd(I, L[j], variant);
1139            }
1140
1141            T1[index] = P[1];
1142            T2[index] = bigint(P[2]);
1143            index++;
1144            j++;
1145         }
1146      }
1147
1148      if(printlevel >= 10)
1149      {
1150         "CPU-time for computing list is "+string(timer - tt)+" seconds.";
1151         "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
1152      }
1153
1154//------------------------  Delete unlucky primes  -----------------------------
1155//-------------  unlucky if and only if the leading ideal is wrong  ------------
1156
1157      LL = deleteUnluckyPrimes(T1,T2,h);
1158      T1 = LL[1];
1159      T2 = LL[2];
1160
1161//-------------------  Now all leading ideals are the same  --------------------
1162//-------------------  Lift results to basering via farey  ---------------------
1163
1164      tt = timer; rt = rtimer;
1165      N = T2[1];
1166      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
1167      H = chinrem(T1,T2);
1168      if(n1 == 1)
1169      {
1170         J = farey(H,N);
1171      }
1172      else
1173      {
1174         for(i = ncols(H); i > 0; i--)
1175         {
1176            arguments_farey[i] = list(ideal(H[i]), N);
1177         }
1178         results_farey = parallelWaitAll("farey", arguments_farey,
1179                                         list(list(list(n1))));
1180         for(i = ncols(H); i > 0; i--)
1181         {
1182            J[i] = results_farey[i][1];
1183         }
1184      }
1185      if(printlevel >= 10)
1186      {
1187         "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
1188         "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
1189      }
1190
1191//----------------  Test if we already have a standard basis of I --------------
1192
1193      tt = timer; rt = rtimer;
1194      if((variant == 1) || (variant == 3))
1195      {
1196         pTest = pTestSB(I,J,L,variant,hi);
1197      }
1198      if((variant == 2) || (variant == 4))
1199      {
1200         pTest = pTestSB(I,J,L,variant);
1201      }
1202
1203      if(printlevel >= 10)
1204      {
1205         "CPU-time for pTest is "+string(timer - tt)+" seconds.";
1206         "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
1207      }
1208
1209      if(pTest)
1210      {
1211         if(printlevel >= 10)
1212         {
1213            "CPU-time for computation without final tests is "
1214            +string(timer - TT)+" seconds.";
1215            "Real-time for computation without final tests is "
1216            +string(rtimer - RT)+" seconds.";
1217         }
1218
1219         attrib(J,"isSB",1);
1220         
1221         if(exactness == 0)
1222         {
1223            option(set, opt);
1224            if(n1 > 1) { kill I_for_fork; }
1225            return(J);
1226         }
1227         
1228         if(exactness == 1)
1229         {
1230            tt = timer; rt = rtimer;
1231            sizeTest = 1 - isIncluded(I,J,n1);
1232
1233            if(printlevel >= 10)
1234            {
1235               "CPU-time for checking if I subset <G> is "
1236               +string(timer - tt)+" seconds.";
1237               "Real-time for checking if I subset <G> is "
1238               +string(rtimer - rt)+" seconds.";
1239            }
1240
1241            if(sizeTest == 0)
1242            {
1243               tt = timer; rt = rtimer;
1244               K = std(J);
1245
1246               if(printlevel >= 10)
1247               {
1248                  "CPU-time for last std-computation is "
1249                  +string(timer - tt)+" seconds.";
1250                  "Real-time for last std-computation is "
1251                  +string(rtimer - rt)+" seconds.";
1252               }
1253
1254               if(size(reduce(K,J)) == 0)
1255               {
1256                  option(set, opt);
1257                  if(n1 > 1) { kill I_for_fork; }
1258                  return(J);
1259               }
1260            }
1261         }
1262      }
1263
1264//--------------  We do not already have a standard basis of I  ----------------
1265//-----------  Therefore do the main computation for more primes  --------------
1266
1267      T1 = H;
1268      T2 = N;
1269      index = 2;
1270
1271      j = size(L) + 1;
1272      tt = timer; rt = rtimer;
1273      L = primeList(I,n3,L,n1);
1274      if(printlevel >= 10)
1275      {
1276         "CPU-time for primeList: "+string(timer-tt)+" seconds.";
1277         "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
1278      }
1279
1280      if(n1 > 1)
1281      {
1282         for(i = 1; i <= n1; i++)
1283         {
1284            open(l(i));
1285            if((variant == 1) || (variant == 3))
1286            {
1287               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1288                                                     eval(variant), eval(hi))));
1289            }
1290            if((variant == 2) || (variant == 4))
1291            {
1292               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1293                                                     eval(variant))));
1294            }
1295         }
1296         j = j + n1;
1297         k = 0;
1298      }
1299   }
1300}
1301example
1302{ "EXAMPLE:"; echo = 2;
1303   ring R1 = 0,(x,y,z,t),dp;
1304   ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
1305   ideal J = modStd(I);
1306   J;
1307   I = homog(I,t);
1308   J = modStd(I);
1309   J;
1310
1311   ring R2 = 0,(x,y,z),ds;
1312   ideal I = jacob(x5+y6+z7+xyz);
1313   ideal J1 = modStd(I,1,0);
1314   J1;
1315
1316   ring R3 = 0,x(1..4),lp;
1317   ideal I = cyclic(4);
1318   ideal J1 = modStd(I,1);
1319   ideal J2 = modStd(I,1,0);
1320   size(reduce(J1,J2));
1321   size(reduce(J2,J1));
1322}
1323
1324////////////////////////////////////////////////////////////////////////////////
1325
1326proc modS(ideal I, list L, list #)
1327"USAGE:  modS(I,L); I ideal, L intvec of primes
1328         if size(#)>0 std is used instead of groebner
1329RETURN:  an ideal which is with high probability a standard basis
1330NOTE:    This procedure is designed for fast experiments.
1331         It is not tested whether the result is a standard basis.
1332         It is not tested whether the result generates I.
1333EXAMPLE: example modS; shows an example
1334"
1335{
1336   int j;
1337   bigint N = 1;
1338   def R0 = basering;
1339   ideal J;
1340   list T;
1341   list rl = ringlist(R0);
1342   if((npars(R0)>0) || (rl[1]>0))
1343   {
1344      ERROR("Characteristic of basering should be zero.");
1345   }
1346   for(j = 1; j <= size(L); j++)
1347   {
1348      N = N*L[j];
1349      rl[1] = L[j];
1350      def @r = ring(rl);
1351      setring @r;
1352      ideal I = fetch(R0,I);
1353      if(size(#) > 0)
1354      {
1355         I = std(I);
1356      }
1357      else
1358      {
1359         I = groebner(I);
1360      }
1361      setring R0;
1362      T[j] = fetch(@r,I);
1363      kill @r;
1364   }
1365   L = deleteUnluckyPrimes(T,L,homog(I));
1366   // unlucky if and only if the leading ideal is wrong
1367   J = farey(chinrem(L[1],L[2]),N);
1368   attrib(J,"isSB",1);
1369   return(J);
1370}
1371example
1372{ "EXAMPLE:"; echo = 2;
1373   list L = 3,5,11,13,181,32003;
1374   ring r = 0,(x,y,z,t),dp;
1375   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1376   I = homog(I,t);
1377   ideal J = modS(I,L);
1378   J;
1379}
1380
1381////////////////////////////////////////////////////////////////////////////////
1382
1383proc modHenselStd(ideal I, list #)
1384"USAGE:  modHenselStd(I);
1385RETURN:  a standard basis of I;
1386NOTE:    The procedure computes a standard basis of I (over the rational
1387         numbers) by using  modular computations and Hensellifting.
1388         For further experiments see procedure modS.
1389EXAMPLE: example modHenselStd; shows an example
1390"
1391{
1392   int i,j;
1393
1394   bigint p = 2134567879;
1395   if(size(#)!=0) { p=#[1]; }
1396   while(!primeTest(I,p))
1397   {
1398      p = prime(random(2000000000,2134567879));
1399   }
1400
1401   def R = basering;
1402   module F,PrevG,PrevZ,Z2;
1403   ideal testG,testG1,G1,G2,G3,Gp;
1404   list L = p;
1405   list rl = ringlist(R);
1406   rl[1] = int(p);
1407
1408   def S = ring(rl);
1409   setring S;
1410   intvec opt = option(get);
1411   option(redSB);
1412   module Z,M,Z2;
1413   ideal I = imap(R,I);
1414   ideal Gp,G1,G2,G3;
1415   Gp,Z = liftstd1(I);
1416   attrib(Gp,"isSB",1);
1417   module ZZ = syz(I);
1418   attrib(ZZ,"isSB",1);
1419   Z = reduce(Z,ZZ);
1420
1421   setring R;
1422   Gp = imap(S,Gp);
1423   PrevZ = imap(S,Z);
1424   PrevG = module(Gp);
1425   F = module(I);
1426   testG = farey(Gp,p);
1427   attrib(testG,"isSB",1);
1428   while(1)
1429   {
1430      i++;
1431      G1 = ideal(1/(p^i) * sum(F*PrevZ,(-1)*PrevG));
1432      setring S;
1433      G1 = imap(R,G1);
1434      G2 = reduce(G1,Gp);
1435      G3 = sum(G1,(-1)*G2);
1436      M = lift(Gp,G3);
1437      Z2 = (-1)*Z*M;
1438
1439      setring R;
1440      G2 = imap(S,G2);
1441      Z2 = imap(S,Z2);
1442      PrevG = sum(PrevG, module(p^i*G2));
1443      PrevZ = sum(PrevZ, multiply(poly(p^i),Z2));
1444      testG1 = farey(ideal(PrevG),p^(i+1));
1445      attrib(testG1,"isSB",1);
1446      if(size(reduce(testG1,testG)) == 0)
1447      {
1448         if(size(reduce(I,testG1)) == 0)       // I is in testG1
1449         {
1450            if(pTestSB(I,testG1,L,2))
1451            {
1452               G3 = std(testG1);               // testG1 is SB
1453               if(size(reduce(G3,testG1)) == 0)
1454               {
1455                  option(set, opt);
1456                  return(G3);
1457               }
1458            }
1459         }
1460      }
1461      testG = testG1;
1462      attrib(testG,"isSB",1);
1463   }
1464}
1465example
1466{ "EXAMPLE:"; echo = 2;
1467   ring r = 0,(x,y,z),dp;
1468   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1469   ideal J = modHenselStd(I);
1470   J;
1471}
1472
1473////////////////////////////////////////////////////////////////////////////////
1474
1475static proc sum(list #)
1476{
1477   if(typeof(#[1])=="ideal")
1478   {
1479      ideal M;
1480   }
1481   else
1482   {
1483      module M;
1484   }
1485
1486   int i;
1487   for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; }
1488   return(M);
1489}
1490
1491////////////////////////////////////////////////////////////////////////////////
1492
1493static proc multiply(poly p, list #)
1494{
1495   if(typeof(#[1])=="ideal")
1496   {
1497      ideal M;
1498   }
1499   else
1500   {
1501      module M;
1502   }
1503
1504   int i;
1505   for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; }
1506   return(M);
1507}
1508
1509
1510////////////////////////////// further examples ////////////////////////////////
1511
1512/*
1513ring r = 0, (x,y,z), lp;
1514poly s1 = 5x3y2z+3y3x2z+7xy2z2;
1515poly s2 = 3xy2z2+x5+11y2z2;
1516poly s3 = 4xyz+7x3+12y3+1;
1517poly s4 = 3x3-4y3+yz2;
1518ideal i =  s1, s2, s3, s4;
1519
1520ring r = 0, (x,y,z), lp;
1521poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
1522poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
1523poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
1524poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
1525ideal i =  s1, s2, s3, s4;
1526
1527ring r = 0, (x,y,z), lp;
1528poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
1529poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
1530poly s3 = 8x3 + 12y3 + xz2 + 3;
1531poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
1532ideal i = s1, s2, s3, s4;
1533
1534int n = 6;
1535ring r = 0,(x(1..n)),lp;
1536ideal i = cyclic(n);
1537ring s = 0, (x(1..n),t), lp;
1538ideal i = imap(r,i);
1539i = homog(i,t);
1540
1541ring r = 0, (x(1..4),s), (dp(4),dp);
1542poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
1543poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
1544poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
1545poly 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);
1546poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
1547ideal i =  s1, s2, s3, s4, s5;
1548
1549ring r = 0, (x,y,z), ds;
1550int a = 16;
1551int b = 15;
1552int c = 4;
1553int t = 1;
1554poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1555         +x^(c-2)*y^c*(y2+t*x)^2;
1556ideal i = jacob(f);
1557
1558ring r = 0, (x,y,z), ds;
1559int a = 25;
1560int b = 25;
1561int c = 5;
1562int t = 1;
1563poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1564         +x^(c-2)*y^c*(y2+t*x)^2;
1565ideal i = jacob(f),f;
1566
1567ring r = 0, (x,y,z), ds;
1568int a = 10;
1569poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
1570ideal i = jacob(f);
1571
1572ring r = 0, (x,y,z), ds;
1573int a = 6;
1574int b = 8;
1575int c = 10;
1576int alpha = 5;
1577int beta = 5;
1578int t = 1;
1579poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
1580         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
1581ideal i = jacob(f);
1582*/
1583
Note: See TracBrowser for help on using the repository browser.