source: git/Singular/LIB/modstd.lib @ a7050d

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