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

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