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

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