source: git/Singular/LIB/modstd.lib @ 78a42b1

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