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

spielwiese
Last change on this file since db347cb was db347cb, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix from master: exceptional cases (algebra.lib), rename sum (modstd.lib)
  • Property mode set to 100644
File size: 38.4 KB
Line 
1////////////////////////////////////////////////////////////////////////////////
2version="$Id: modstd.lib 15103 2012-07-11 10:00:13Z motsak $";
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,2);
714   P;
715
716   ring r2 = 0, x(1..5), lp;
717   ideal I = cyclic(5);
718   int q = 32003;
719   list Q = modpStd(I,q,4);
720   Q;
721}
722
723////////////////////////////// main procedures /////////////////////////////////
724
725proc modStd(ideal I, list #)
726"USAGE:  modStd(I); I ideal
727ASSUME:  If size(#) > 0, then # contains either 1, 2 or 4 integers such that
728@*       - #[1] is the number of available processors for the computation,
729@*       - #[2] is an optional parameter for the exactness of the computation,
730                if #[2] = 1, the procedure computes a standard basis for sure,
731@*       - #[3] is the number of primes until the first lifting,
732@*       - #[4] is the constant number of primes between two liftings until
733           the computation stops.
734RETURN:  a standard basis of I if no warning appears;
735NOTE:    The procedure computes a standard basis of I (over the rational
736         numbers) by using modular methods.
737         By default the procedure computes a standard basis of I for sure, but
738         if the optional parameter #[2] = 0, it computes a standard basis of I
739         with high probability.
740         The procedure distinguishes between different variants for the standard
741         basis computation in positive characteristic depending on the ordering
742         of the basering, the parameter #[2] and if the ideal I is homogeneous.
743@*       - variant = 1, if I is homogeneous,
744@*       - variant = 2, if I is not homogeneous, 1-block-ordering,
745@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
746                        > 1 block),
747@*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
748EXAMPLE: example modStd; shows an example
749"
750{
751   int TT = timer;
752   int RT = rtimer;
753
754   def R0 = basering;
755   list rl = ringlist(R0);
756   if((npars(R0) > 0) || (rl[1] > 0))
757   {
758      ERROR("Characteristic of basering should be zero, basering should
759             have no parameters.");
760   }
761
762   int index = 1;
763   int i,k,c;
764   int j = 1;
765   int pTest, sizeTest;
766   int en = 2134567879;
767   int an = 1000000000;
768   bigint N;
769
770//--------------------  Initialize optional parameters  ------------------------
771   if(size(#) > 0)
772   {
773      if(size(#) == 1)
774      {
775         int n1 = #[1];
776         int exactness = 1;
777         if(n1 >= 10)
778         {
779            int n2 = n1 + 1;
780            int n3 = n1;
781         }
782         else
783         {
784            int n2 = 10;
785            int n3 = 10;
786         }
787      }
788      if(size(#) == 2)
789      {
790         int n1 = #[1];
791         int exactness = #[2];
792         if(n1 >= 10)
793         {
794            int n2 = n1 + 1;
795            int n3 = n1;
796         }
797         else
798         {
799            int n2 = 10;
800            int n3 = 10;
801         }
802      }
803      if(size(#) == 4)
804      {
805         int n1 = #[1];
806         int exactness = #[2];
807         if(n1 >= #[3])
808         {
809            int n2 = n1 + 1;
810         }
811         else
812         {
813            int n2 = #[3];
814         }
815         if(n1 >= #[4])
816         {
817            int n3 = n1;
818         }
819         else
820         {
821            int n3 = #[4];
822         }
823      }
824   }
825   else
826   {
827      int n1 = 1;
828      int exactness = 1;
829      int n2 = 10;
830      int n3 = 10;
831   }
832
833   if(printlevel >= 10)
834   {
835      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
836       +", exactness = "+string(exactness);
837   }
838
839//-------------------------  Save current options  -----------------------------
840   intvec opt = option(get);
841
842   option(redSB);
843
844//--------------------  Initialize the list of primes  -------------------------
845   int tt = timer;
846   int rt = rtimer;
847   intvec L = primeList(I,n2,n1);
848   if(printlevel >= 10)
849   {
850      "CPU-time for primeList: "+string(timer-tt)+" seconds.";
851      "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
852   }
853   L[5] = prime(random(an,en));
854
855//---------------------  Decide which variant to take  -------------------------
856   int variant;
857   int h = homog(I);
858
859   tt = timer;
860   rt = rtimer;
861
862   if(!mixedTest())
863   {
864      if(h)
865      {
866         variant = 1;
867         if(printlevel >= 10) { "variant = 1"; }
868
869         rl[1] = L[5];
870         def @r = ring(rl);
871         setring @r;
872         def @s = changeord("dp");
873         setring @s;
874         ideal I = std(fetch(R0,I));
875         intvec hi = hilb(I,1);
876         setring R0;
877         kill @r,@s;
878      }
879      else
880      {
881         string ordstr_R0 = ordstr(R0);
882         int neg = 1 - attrib(R0,"global");
883
884         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
885         {
886            variant = 2;
887            if(printlevel >= 10) { "variant = 2"; }
888         }
889         else
890         {
891            string order;
892            if(system("nblocks") <= 2)
893            {
894               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
895                                       + find(ordstr_R0, "rp") <= 0)
896               {
897                  order = "simple";
898               }
899            }
900
901            if((order == "simple") || (size(rl) > 4))
902            {
903               variant = 2;
904               if(printlevel >= 10) { "variant = 2"; }
905            }
906            else
907            {
908               rl[1] = L[5];
909               def @r = ring(rl);
910               setring @r;
911
912               def @s = changeord("dp");
913               setring @s;
914               ideal I = std(fetch(R0,I));
915               if(dim(I) == 0)
916               {
917                  variant = 4;
918                  if(printlevel >= 10) { "variant = 4"; }
919               }
920               else
921               {
922                  variant = 3;
923                  if(printlevel >= 10) { "variant = 3"; }
924
925                  int nvar@r = nvars(@r);
926                  intvec w;
927                  for(i = 1; i <= nvar@r; i++)
928                  {
929                     w[i] = deg(var(i));
930                  }
931                  w[nvar@r + 1] = 1;
932
933                  list hiRi = hilbRing(fetch(R0,I),w);
934                  intvec W = hiRi[2];
935                  @s = hiRi[1];
936                  setring @s;
937
938                  Id(1) = std(Id(1));
939                  intvec hi = hilb(Id(1), 1, W);
940               }
941
942               setring R0;
943               kill @r,@s;
944            }
945         }
946      }
947   }
948   else
949   {
950      if(exactness == 1) { return(groebner(I)); }
951      if(h)
952      {
953         variant = 1;
954         if(printlevel >= 10) { "variant = 1"; }
955         rl[1] = L[5];
956         def @r = ring(rl);
957         setring @r;
958         def @s = changeord("dp");
959         setring @s;
960         ideal I = std(fetch(R0,I));
961         intvec hi = hilb(I,1);
962         setring R0;
963         kill @r,@s;
964      }
965      else
966      {
967         string ordstr_R0 = ordstr(R0);
968         int neg = 1 - attrib(R0,"global");
969
970         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
971         {
972            variant = 2;
973            if(printlevel >= 10) { "variant = 2"; }
974         }
975         else
976         {
977            string order;
978            if(system("nblocks") <= 2)
979            {
980               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
981                                       + find(ordstr_R0, "rp") <= 0)
982               {
983                  order = "simple";
984               }
985            }
986
987            if((order == "simple") || (size(rl) > 4))
988            {
989               variant = 2;
990               if(printlevel >= 10) { "variant = 2"; }
991            }
992            else
993            {
994               variant = 3;
995               if(printlevel >= 10) { "variant = 3"; }
996
997               rl[1] = L[5];
998               def @r = ring(rl);
999               setring @r;
1000               int nvar@r = nvars(@r);
1001               intvec w;
1002               for(i = 1; i <= nvar@r; i++)
1003               {
1004                  w[i] = deg(var(i));
1005               }
1006               w[nvar@r + 1] = 1;
1007
1008               list hiRi = hilbRing(fetch(R0,I),w);
1009               intvec W = hiRi[2];
1010               def @s = hiRi[1];
1011               setring @s;
1012
1013               Id(1) = std(Id(1));
1014               intvec hi = hilb(Id(1), 1, W);
1015
1016               setring R0;
1017               kill @r,@s;
1018            }
1019         }
1020      }
1021   }
1022
1023   list P,T1,T2,T3,LL;
1024
1025   ideal J,K,H;
1026
1027//-----  If there is more than one processor available, we parallelize the  ----
1028//-----  main standard basis computations in positive characteristic        ----
1029
1030   if(n1 > 1)
1031   {
1032      ideal I_for_fork = I;
1033      export(I_for_fork);           // I available for each link
1034
1035//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
1036//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
1037
1038      for(i = 1; i <= n1; i++)
1039      {
1040         //link l(i) = "MPtcp:fork";
1041         link l(i) = "ssi:fork";
1042         open(l(i));
1043         if((variant == 1) || (variant == 3))
1044         {
1045            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1046                                                  eval(variant), eval(hi))));
1047         }
1048         if((variant == 2) || (variant == 4))
1049         {
1050            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
1051                                                  eval(variant))));
1052         }
1053      }
1054
1055      int t = timer;
1056      if((variant == 1) || (variant == 3))
1057      {
1058         P = modpStd(I_for_fork, L[1], variant, hi);
1059      }
1060      if((variant == 2) || (variant == 4))
1061      {
1062         P = modpStd(I_for_fork, L[1], variant);
1063      }
1064      t = timer - t;
1065      if(t > 60) { t = 60; }
1066      int i_sleep = system("sh", "sleep "+string(t));
1067      T1[1] = P[1];
1068      T2[1] = bigint(P[2]);
1069      index++;
1070
1071      j = j + n1 + 1;
1072   }
1073
1074//--------------  Main standard basis computations in positive  ----------------
1075//----------------------  characteristic start here  ---------------------------
1076
1077   list arguments_farey, results_farey;
1078
1079   while(1)
1080   {
1081      tt = timer; rt = rtimer;
1082
1083      if(printlevel >= 10) { "size(L) = "+string(size(L)); }
1084
1085      if(n1 > 1)
1086      {
1087         while(j <= size(L) + 1)
1088         {
1089            for(i = 1; i <= n1; i++)
1090            {
1091               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
1092               if(status(l(i), "read", "ready"))
1093               {
1094                  //--- read the result from l(i) ---
1095                  P = read(l(i));
1096                  T1[index] = P[1];
1097                  T2[index] = bigint(P[2]);
1098                  index++;
1099
1100                  if(j <= size(L))
1101                  {
1102                     if((variant == 1) || (variant == 3))
1103                     {
1104                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
1105                                                  eval(variant), eval(hi))));
1106                        j++;
1107                     }
1108                     if((variant == 2) || (variant == 4))
1109                     {
1110                        write(l(i), quote(modpStd(I_for_fork,
1111                                                  eval(L[j]), eval(variant))));
1112                        j++;
1113                     }
1114                  }
1115                  else
1116                  {
1117                     k++;
1118                     close(l(i));
1119                  }
1120               }
1121            }
1122            //--- k describes the number of closed links ---
1123            if(k == n1)
1124            {
1125               j++;
1126            }
1127            i_sleep = system("sh", "sleep "+string(t));
1128         }
1129      }
1130      else
1131      {
1132         while(j <= size(L))
1133         {
1134            if((variant == 1) || (variant == 3))
1135            {
1136               P = modpStd(I, L[j], variant, hi);
1137            }
1138            if((variant == 2) || (variant == 4))
1139            {
1140               P = modpStd(I, L[j], variant);
1141            }
1142
1143            T1[index] = P[1];
1144            T2[index] = bigint(P[2]);
1145            index++;
1146            j++;
1147         }
1148      }
1149
1150      if(printlevel >= 10)
1151      {
1152         "CPU-time for computing list is "+string(timer - tt)+" seconds.";
1153         "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
1154      }
1155
1156//------------------------  Delete unlucky primes  -----------------------------
1157//-------------  unlucky if and only if the leading ideal is wrong  ------------
1158
1159      LL = deleteUnluckyPrimes(T1,T2,h);
1160      T1 = LL[1];
1161      T2 = LL[2];
1162
1163//-------------------  Now all leading ideals are the same  --------------------
1164//-------------------  Lift results to basering via farey  ---------------------
1165
1166      tt = timer; rt = rtimer;
1167      N = T2[1];
1168      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
1169      H = chinrem(T1,T2);
1170      if(n1 == 1)
1171      {
1172         J = farey(H,N);
1173      }
1174      else
1175      {
1176         for(i = ncols(H); i > 0; i--)
1177         {
1178            arguments_farey[i] = list(ideal(H[i]), N);
1179         }
1180         results_farey = parallelWaitAll("farey", arguments_farey,
1181                                         list(list(list(n1))));
1182         for(i = ncols(H); i > 0; i--)
1183         {
1184            J[i] = results_farey[i][1];
1185         }
1186      }
1187      if(printlevel >= 10)
1188      {
1189         "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
1190         "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
1191      }
1192
1193//----------------  Test if we already have a standard basis of I --------------
1194
1195      tt = timer; rt = rtimer;
1196      if((variant == 1) || (variant == 3))
1197      {
1198         pTest = pTestSB(I,J,L,variant,hi);
1199      }
1200      if((variant == 2) || (variant == 4))
1201      {
1202         pTest = pTestSB(I,J,L,variant);
1203      }
1204
1205      if(printlevel >= 10)
1206      {
1207         "CPU-time for pTest is "+string(timer - tt)+" seconds.";
1208         "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
1209      }
1210
1211      if(pTest)
1212      {
1213         if(printlevel >= 10)
1214         {
1215            "CPU-time for computation without final tests is "
1216            +string(timer - TT)+" seconds.";
1217            "Real-time for computation without final tests is "
1218            +string(rtimer - RT)+" seconds.";
1219         }
1220
1221         attrib(J,"isSB",1);
1222
1223         if(exactness == 0)
1224         {
1225            option(set, opt);
1226            if(n1 > 1) { kill I_for_fork; }
1227            return(J);
1228         }
1229
1230         if(exactness == 1)
1231         {
1232            tt = timer; rt = rtimer;
1233            sizeTest = 1 - isIncluded(I,J,n1);
1234
1235            if(printlevel >= 10)
1236            {
1237               "CPU-time for checking if I subset <G> is "
1238               +string(timer - tt)+" seconds.";
1239               "Real-time for checking if I subset <G> is "
1240               +string(rtimer - rt)+" seconds.";
1241            }
1242
1243            if(sizeTest == 0)
1244            {
1245               tt = timer; rt = rtimer;
1246               K = std(J);
1247
1248               if(printlevel >= 10)
1249               {
1250                  "CPU-time for last std-computation is "
1251                  +string(timer - tt)+" seconds.";
1252                  "Real-time for last std-computation is "
1253                  +string(rtimer - rt)+" seconds.";
1254               }
1255
1256               if(size(reduce(K,J)) == 0)
1257               {
1258                  option(set, opt);
1259                  if(n1 > 1) { kill I_for_fork; }
1260                  return(J);
1261               }
1262            }
1263         }
1264      }
1265
1266//--------------  We do not already have a standard basis of I  ----------------
1267//-----------  Therefore do the main computation for more primes  --------------
1268
1269      T1 = H;
1270      T2 = N;
1271      index = 2;
1272
1273      j = size(L) + 1;
1274      tt = timer; rt = rtimer;
1275      L = primeList(I,n3,L,n1);
1276      if(printlevel >= 10)
1277      {
1278         "CPU-time for primeList: "+string(timer-tt)+" seconds.";
1279         "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
1280      }
1281
1282      if(n1 > 1)
1283      {
1284         for(i = 1; i <= n1; i++)
1285         {
1286            open(l(i));
1287            if((variant == 1) || (variant == 3))
1288            {
1289               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1290                                                     eval(variant), eval(hi))));
1291            }
1292            if((variant == 2) || (variant == 4))
1293            {
1294               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1295                                                     eval(variant))));
1296            }
1297         }
1298         j = j + n1;
1299         k = 0;
1300      }
1301   }
1302}
1303example
1304{ "EXAMPLE:"; echo = 2;
1305   ring R1 = 0,(x,y,z,t),dp;
1306   ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
1307   ideal J = modStd(I);
1308   J;
1309   I = homog(I,t);
1310   J = modStd(I);
1311   J;
1312
1313   ring R2 = 0,(x,y,z),ds;
1314   ideal I = jacob(x5+y6+z7+xyz);
1315   ideal J1 = modStd(I,1,0);
1316   J1;
1317
1318   ring R3 = 0,x(1..4),lp;
1319   ideal I = cyclic(4);
1320   ideal J1 = modStd(I,1);
1321   ideal J2 = modStd(I,1,0);
1322   size(reduce(J1,J2));
1323   size(reduce(J2,J1));
1324}
1325
1326////////////////////////////////////////////////////////////////////////////////
1327
1328proc modS(ideal I, list L, list #)
1329"USAGE:  modS(I,L); I ideal, L intvec of primes
1330         if size(#)>0 std is used instead of groebner
1331RETURN:  an ideal which is with high probability a standard basis
1332NOTE:    This procedure is designed for fast experiments.
1333         It is not tested whether the result is a standard basis.
1334         It is not tested whether the result generates I.
1335EXAMPLE: example modS; shows an example
1336"
1337{
1338   int j;
1339   bigint N = 1;
1340   def R0 = basering;
1341   ideal J;
1342   list T;
1343   list rl = ringlist(R0);
1344   if((npars(R0)>0) || (rl[1]>0))
1345   {
1346      ERROR("Characteristic of basering should be zero.");
1347   }
1348   for(j = 1; j <= size(L); j++)
1349   {
1350      N = N*L[j];
1351      rl[1] = L[j];
1352      def @r = ring(rl);
1353      setring @r;
1354      ideal I = fetch(R0,I);
1355      if(size(#) > 0)
1356      {
1357         I = std(I);
1358      }
1359      else
1360      {
1361         I = groebner(I);
1362      }
1363      setring R0;
1364      T[j] = fetch(@r,I);
1365      kill @r;
1366   }
1367   L = deleteUnluckyPrimes(T,L,homog(I));
1368   // unlucky if and only if the leading ideal is wrong
1369   J = farey(chinrem(L[1],L[2]),N);
1370   attrib(J,"isSB",1);
1371   return(J);
1372}
1373example
1374{ "EXAMPLE:"; echo = 2;
1375   list L = 3,5,11,13,181,32003;
1376   ring r = 0,(x,y,z,t),dp;
1377   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1378   I = homog(I,t);
1379   ideal J = modS(I,L);
1380   J;
1381}
1382
1383////////////////////////////////////////////////////////////////////////////////
1384
1385proc modHenselStd(ideal I, list #)
1386"USAGE:  modHenselStd(I);
1387RETURN:  a standard basis of I;
1388NOTE:    The procedure computes a standard basis of I (over the rational
1389         numbers) by using  modular computations and Hensellifting.
1390         For further experiments see procedure modS.
1391EXAMPLE: example modHenselStd; shows an example
1392"
1393{
1394   int i,j;
1395
1396   bigint p = 2134567879;
1397   if(size(#)!=0) { p=#[1]; }
1398   while(!primeTest(I,p))
1399   {
1400      p = prime(random(2000000000,2134567879));
1401   }
1402
1403   def R = basering;
1404   module F,PrevG,PrevZ,Z2;
1405   ideal testG,testG1,G1,G2,G3,Gp;
1406   list L = p;
1407   list rl = ringlist(R);
1408   rl[1] = int(p);
1409
1410   def S = ring(rl);
1411   setring S;
1412   intvec opt = option(get);
1413   option(redSB);
1414   module Z,M,Z2;
1415   ideal I = imap(R,I);
1416   ideal Gp,G1,G2,G3;
1417   Gp,Z = liftstd1(I);
1418   attrib(Gp,"isSB",1);
1419   module ZZ = syz(I);
1420   attrib(ZZ,"isSB",1);
1421   Z = reduce(Z,ZZ);
1422
1423   setring R;
1424   Gp = imap(S,Gp);
1425   PrevZ = imap(S,Z);
1426   PrevG = module(Gp);
1427   F = module(I);
1428   testG = farey(Gp,p);
1429   attrib(testG,"isSB",1);
1430   while(1)
1431   {
1432      i++;
1433      G1 = ideal(1/(p^i) * sum_id(F*PrevZ,(-1)*PrevG));
1434      setring S;
1435      G1 = imap(R,G1);
1436      G2 = reduce(G1,Gp);
1437      G3 = sum_id(G1,(-1)*G2);
1438      M = lift(Gp,G3);
1439      Z2 = (-1)*Z*M;
1440
1441      setring R;
1442      G2 = imap(S,G2);
1443      Z2 = imap(S,Z2);
1444      PrevG = sum_id(PrevG, module(p^i*G2));
1445      PrevZ = sum_id(PrevZ, multiply(poly(p^i),Z2));
1446      testG1 = farey(ideal(PrevG),p^(i+1));
1447      attrib(testG1,"isSB",1);
1448      if(size(reduce(testG1,testG)) == 0)
1449      {
1450         if(size(reduce(I,testG1)) == 0)       // I is in testG1
1451         {
1452            if(pTestSB(I,testG1,L,2))
1453            {
1454               G3 = std(testG1);               // testG1 is SB
1455               if(size(reduce(G3,testG1)) == 0)
1456               {
1457                  option(set, opt);
1458                  return(G3);
1459               }
1460            }
1461         }
1462      }
1463      testG = testG1;
1464      attrib(testG,"isSB",1);
1465   }
1466}
1467example
1468{ "EXAMPLE:"; echo = 2;
1469   ring r = 0,(x,y,z),dp;
1470   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1471   ideal J = modHenselStd(I);
1472   J;
1473}
1474
1475////////////////////////////////////////////////////////////////////////////////
1476
1477static proc sum_id(list #)
1478{
1479   if(typeof(#[1])=="ideal")
1480   {
1481      ideal M;
1482   }
1483   else
1484   {
1485      module M;
1486   }
1487
1488   int i;
1489   for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; }
1490   return(M);
1491}
1492
1493////////////////////////////////////////////////////////////////////////////////
1494
1495static proc multiply(poly p, list #)
1496{
1497   if(typeof(#[1])=="ideal")
1498   {
1499      ideal M;
1500   }
1501   else
1502   {
1503      module M;
1504   }
1505
1506   int i;
1507   for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; }
1508   return(M);
1509}
1510
1511
1512////////////////////////////// further examples ////////////////////////////////
1513
1514/*
1515ring r = 0, (x,y,z), lp;
1516poly s1 = 5x3y2z+3y3x2z+7xy2z2;
1517poly s2 = 3xy2z2+x5+11y2z2;
1518poly s3 = 4xyz+7x3+12y3+1;
1519poly s4 = 3x3-4y3+yz2;
1520ideal i =  s1, s2, s3, s4;
1521
1522ring r = 0, (x,y,z), lp;
1523poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
1524poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
1525poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
1526poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
1527ideal i =  s1, s2, s3, s4;
1528
1529ring r = 0, (x,y,z), lp;
1530poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
1531poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
1532poly s3 = 8x3 + 12y3 + xz2 + 3;
1533poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
1534ideal i = s1, s2, s3, s4;
1535
1536int n = 6;
1537ring r = 0,(x(1..n)),lp;
1538ideal i = cyclic(n);
1539ring s = 0, (x(1..n),t), lp;
1540ideal i = imap(r,i);
1541i = homog(i,t);
1542
1543ring r = 0, (x(1..4),s), (dp(4),dp);
1544poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
1545poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
1546poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
1547poly 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);
1548poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
1549ideal i =  s1, s2, s3, s4, s5;
1550
1551ring r = 0, (x,y,z), ds;
1552int a = 16;
1553int b = 15;
1554int c = 4;
1555int t = 1;
1556poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1557         +x^(c-2)*y^c*(y2+t*x)^2;
1558ideal i = jacob(f);
1559
1560ring r = 0, (x,y,z), ds;
1561int a = 25;
1562int b = 25;
1563int c = 5;
1564int t = 1;
1565poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1566         +x^(c-2)*y^c*(y2+t*x)^2;
1567ideal i = jacob(f),f;
1568
1569ring r = 0, (x,y,z), ds;
1570int a = 10;
1571poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
1572ideal i = jacob(f);
1573
1574ring r = 0, (x,y,z), ds;
1575int a = 6;
1576int b = 8;
1577int c = 10;
1578int alpha = 5;
1579int beta = 5;
1580int t = 1;
1581poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
1582         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
1583ideal i = jacob(f);
1584*/
1585
Note: See TracBrowser for help on using the repository browser.