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

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