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

spielwiese
Last change on this file since b4f7eb was b4f7eb, checked in by Stefan Steidel <steidel@…>, 12 years ago
Some printlevel improvements. git-svn-id: file:///usr/local/Singular/svn/trunk@14231 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 36.2 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)
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 >= 11)
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 >= 11)
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) div 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@*       - #[2] is an optional parameter for the exactness of the computation,
629                if #[2] = 1, the procedure computes a standard basis for sure,
630@*       - #[3] is the number of primes until the first lifting,
631@*       - #[4] is the constant number of primes between two liftings until
632           the computation stops.
633RETURN:  a standard basis of I if no warning appears;
634NOTE:    The procedure computes a standard basis of I (over the rational
635         numbers) by using  modular methods. If a warning appears then the
636         result is a standard basis containing I and with high probability
637         a standard basis of I.
638         By default the procedure computes a standard basis of I for sure, but
639         if the optional parameter #[2] = 0, it computes a standard basis of I
640         with high probability and a warning appears at the end.
641         The procedure distinguishes between different variants for the standard
642         basis computation in positive characteristic depending on the ordering
643         of the basering, the parameter #[2] and if the ideal I is homogeneous.
644@*       - variant = 1, if I is homogeneous and exactness = 0,
645@*       - variant = 2, if I is not homogeneous, 1-block-ordering and
646                        exactness = 0,
647@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
648                        > 1 block) and exactness = 0,
649@*       - variant = 4, if I is homogeneous and exactness = 1,
650@*       - variant = 5, if I is not homogeneous, 1-block-ordering and
651                        exactness = 1,
652@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
653                        > 1 block) and exactness = 1.
654EXAMPLE: example modStd; shows an example
655"
656{
657   int TT = timer;
658   int RT = rtimer;
659
660   def R0 = basering;
661   list rl = ringlist(R0);
662   if((npars(R0) > 0) || (rl[1] > 0))
663   {
664      ERROR("Characteristic of basering should be zero, basering should
665             have no parameters.");
666   }
667
668   int index = 1;
669   int i,k,c;
670   int pd = printlevel-voice+2;
671   int j = 1;
672   int pTest, sizeTest;
673   int en = 2134567879;
674   int an = 1000000000;
675   bigint N;
676
677//--------------------  Initialize optional parameters  ------------------------
678   if(size(#) > 0)
679   {
680      if(size(#) == 1)
681      {
682         int n1 = #[1];
683         int exactness = 1;
684         int n2 = 10;
685         int n3 = 10;
686      }
687      if(size(#) == 2)
688      {
689         int n1 = #[1];
690         int exactness = #[2];
691         int n2 = 10;
692         int n3 = 10;
693      }
694      if(size(#) == 4)
695      {
696         int n1 = #[1];
697         int exactness = #[2];
698         int n2 = #[3];
699         int n3 = #[4];
700      }
701   }
702   else
703   {
704      int n1 = 1;
705      int exactness = 1;
706      int n2 = 10;
707      int n3 = 10;
708   }
709
710//-------------------------  Save current options  -----------------------------
711   intvec opt = option(get);
712
713   option(redSB);
714
715//--------------------  Initialize the list of primes  -------------------------
716   intvec L = primeList(I,n2);
717   L[5] = prime(random(an,en));
718
719//---------------------  Decide which variant to take  -------------------------
720   int variant;
721   int h = homog(I);
722
723   int tt = timer;
724   int rt = rtimer;
725
726   if(!mixedTest())
727   {
728      if(h)
729      {
730         if(exactness == 0)
731         {
732            variant = 1;
733            if(printlevel >= 10) { "variant = 1"; }
734         }
735         if(exactness == 1)
736         {
737            variant = 4;
738            if(printlevel >= 10) { "variant = 4"; }
739         }
740         rl[1] = L[5];
741         def @r = ring(rl);
742         setring @r;
743         def @s = changeord("dp");
744         setring @s;
745         ideal I = std(fetch(R0,I));
746         intvec hi = hilb(I,1);
747         setring R0;
748         kill @r,@s;
749      }
750      else
751      {
752         string ordstr_R0 = ordstr(R0);
753         int neg = 1 - attrib(R0,"global");
754
755         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
756         {
757            if(exactness == 0)
758            {
759               variant = 2;
760               if(printlevel >= 10) { "variant = 2"; }
761            }
762            if(exactness == 1)
763            {
764               variant = 5;
765               if(printlevel >= 10) { "variant = 5"; }
766            }
767         }
768         else
769         {
770            string order;
771            if(system("nblocks") <= 2)
772            {
773               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
774                                       + find(ordstr_R0, "rp") <= 0)
775               {
776                  order = "simple";
777               }
778            }
779
780            if((order == "simple") || (size(rl) > 4))
781            {
782               if(exactness == 0)
783               {
784                  variant = 2;
785                  if(printlevel >= 10) { "variant = 2"; }
786               }
787               if(exactness == 1)
788               {
789                  variant = 5;
790                  if(printlevel >= 10) { "variant = 5"; }
791               }
792            }
793            else
794            {
795               if(exactness == 0)
796               {
797                  variant = 3;
798                  if(printlevel >= 10) { "variant = 3"; }
799               }
800               if(exactness == 1)
801               {
802                  variant = 6;
803                  if(printlevel >= 10) { "variant = 6"; }
804               }
805
806               rl[1] = L[5];
807               def @r = ring(rl);
808               setring @r;
809               int nvar@r = nvars(@r);
810               intvec w;
811               for(i = 1; i <= nvar@r; i++)
812               {
813                  w[i] = deg(var(i));
814               }
815               w[nvar@r + 1] = 1;
816
817               list hiRi = hilbRing(fetch(R0,I),w);
818               intvec W = hiRi[2];
819               def @s = hiRi[1];
820               setring @s;
821
822               Id(1) = std(Id(1));
823               intvec hi = hilb(Id(1), 1, W);
824
825               setring R0;
826               kill @r,@s;
827            }
828         }
829      }
830   }
831   else
832   {
833      if(exactness == 1) { return(groebner(I)); }
834      if(h)
835      {
836         variant = 1; if(printlevel >= 10) { "variant = 1"; }
837         rl[1] = L[5];
838         def @r = ring(rl);
839         setring @r;
840         def @s = changeord("dp");
841         setring @s;
842         ideal I = std(fetch(R0,I));
843         intvec hi = hilb(I,1);
844         setring R0;
845         kill @r,@s;
846      }
847      else
848      {
849         string ordstr_R0 = ordstr(R0);
850         int neg = 1 - attrib(R0,"global");
851
852         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
853         {
854            variant = 2;
855            if(printlevel >= 10) { "variant = 2"; }
856         }
857         else
858         {
859            string order;
860            if(system("nblocks") <= 2)
861            {
862               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
863                                       + find(ordstr_R0, "rp") <= 0)
864               {
865                  order = "simple";
866               }
867            }
868
869            if((order == "simple") || (size(rl) > 4))
870            {
871               variant = 2;
872               if(printlevel >= 10) { "variant = 2"; }
873            }
874            else
875            {
876               variant = 3;
877               if(printlevel >= 10) { "variant = 3"; }
878
879               rl[1] = L[5];
880               def @r = ring(rl);
881               setring @r;
882               int nvar@r = nvars(@r);
883               intvec w;
884               for(i = 1; i <= nvar@r; i++)
885               {
886                  w[i] = deg(var(i));
887               }
888               w[nvar@r + 1] = 1;
889
890               list hiRi = hilbRing(fetch(R0,I),w);
891               intvec W = hiRi[2];
892               def @s = hiRi[1];
893               setring @s;
894
895               Id(1) = std(Id(1));
896               intvec hi = hilb(Id(1), 1, W);
897
898               setring R0;
899               kill @r,@s;
900            }
901         }
902      }
903   }
904
905   list P,T1,T2,T3,LL;
906
907   ideal J,K,H;
908
909//-----  If there is more than one processor available, we parallelize the  ----
910//-----  main standard basis computations in positive characteristic        ----
911
912   if(n1 > 1)
913   {
914      ideal I_for_fork = I;
915      export(I_for_fork);           // I available for each link
916
917//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
918//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
919
920      for(i = 1; i <= n1; i++)
921      {
922         //link l(i) = "MPtcp:fork";
923         link l(i) = "ssi:fork";
924         open(l(i));
925         if((variant == 1) || (variant == 3) ||
926            (variant == 4) || (variant == 6))
927         {
928            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
929                                                  eval(variant), eval(hi))));
930         }
931         if((variant == 2) || (variant == 5))
932         {
933            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
934                                                  eval(variant))));
935         }
936      }
937
938      int t = timer;
939      if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
940      {
941         P = modpStd(I_for_fork, L[1], variant, hi);
942      }
943      if((variant == 2) || (variant == 5))
944      {
945         P = modpStd(I_for_fork, L[1], variant);
946      }
947      t = timer - t;
948      if(t > 60) { t = 60; }
949      int i_sleep = system("sh", "sleep "+string(t));
950      T1[1] = P[1];
951      T2[1] = bigint(P[2]);
952      index++;
953
954      j = j + n1 + 1;
955   }
956
957//--------------  Main standard basis computations in positive  ----------------
958//----------------------  characteristic start here  ---------------------------
959
960   while(1)
961   {
962      tt = timer; rt = rtimer;
963
964      if(n1 > 1)
965      {
966         if(printlevel >= 10) { "size(L) = "+string(size(L)); }
967         while(j <= size(L) + 1)
968         {
969            for(i = 1; i <= n1; i++)
970            {
971               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
972               if(status(l(i), "read", "ready"))
973               {
974                  //--- read the result from l(i) ---
975                  P = read(l(i));
976                  T1[index] = P[1];
977                  T2[index] = bigint(P[2]);
978                  index++;
979
980                  if(j <= size(L))
981                  {
982                     if((variant == 1) || (variant == 3) ||
983                        (variant == 4) || (variant == 6))
984                     {
985                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
986                                                  eval(variant), eval(hi))));
987                        j++;
988                     }
989                     if((variant == 2) || (variant == 5))
990                     {
991                        write(l(i), quote(modpStd(I_for_fork,
992                                                  eval(L[j]), eval(variant))));
993                        j++;
994                     }
995                  }
996                  else
997                  {
998                     k++;
999                     close(l(i));
1000                  }
1001               }
1002            }
1003            //--- k describes the number of closed links ---
1004            if(k == n1)
1005            {
1006               j++;
1007            }
1008            i_sleep = system("sh", "sleep "+string(t));
1009         }
1010      }
1011      else
1012      {
1013         while(j <= size(L))
1014         {
1015            if((variant == 1) || (variant == 3) ||
1016               (variant == 4) || (variant == 6))
1017            {
1018               P = modpStd(I, L[j], variant, hi);
1019            }
1020            if((variant == 2) || (variant == 5))
1021            {
1022               P = modpStd(I, L[j], variant);
1023            }
1024
1025            T1[index] = P[1];
1026            T2[index] = bigint(P[2]);
1027            index++;
1028            j++;
1029         }
1030      }
1031
1032      if(printlevel >= 10)
1033      {
1034         "CPU-time for computing list is "+string(timer - tt)+" seconds.";
1035         "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
1036      }
1037
1038//------------------------  Delete unlucky primes  -----------------------------
1039//-------------  unlucky if and only if the leading ideal is wrong  ------------
1040
1041      LL = deleteUnluckyPrimes(T1,T2,h);
1042      T1 = LL[1];
1043      T2 = LL[2];
1044
1045//-------------------  Now all leading ideals are the same  --------------------
1046//-------------------  Lift results to basering via farey  ---------------------
1047
1048      tt = timer;     
1049      N = T2[1];
1050      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
1051      H = chinrem(T1,T2);
1052      J = farey(H,N);
1053      if(printlevel >= 10) { "Lifting-process takes "+string(timer - tt)
1054                             +" seconds"; }
1055
1056//----------------  Test if we already have a standard basis of I --------------
1057
1058      tt = timer; rt = rtimer;
1059      if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
1060      {
1061         pTest = pTestSB(I,J,L,variant,hi);
1062      }
1063      if((variant == 2) || (variant == 5))
1064      {
1065         pTest = pTestSB(I,J,L,variant);
1066      }
1067
1068      if(printlevel >= 10)
1069      {
1070         "CPU-time for pTest is "+string(timer - tt)+" seconds.";
1071         "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
1072      }
1073
1074      if(pTest)
1075      {
1076         if(printlevel >= 10)
1077         {
1078            "CPU-time for computation without final tests is "
1079            +string(timer - TT)+" seconds.";
1080            "Real-time for computation without final tests is "
1081            +string(rtimer - RT)+" seconds.";
1082         }
1083
1084         attrib(J,"isSB",1);
1085         tt = timer; rt = rtimer;
1086         sizeTest = 1 - isIncluded(I,J,n1);
1087
1088         if(printlevel >= 10)
1089         {
1090            "CPU-time for checking if I subset <G> is "
1091            +string(timer - tt)+" seconds.";
1092            "Real-time for checking if I subset <G> is "
1093            +string(rtimer - rt)+" seconds.";
1094         }
1095
1096         if(sizeTest == 0)
1097         {
1098            if((variant == 1) || (variant == 2) || (variant == 3))
1099            {
1100               "===========================================================";
1101               "WARNING: Output might not be a standard basis of the input.";
1102               "===========================================================";
1103               option(set, opt);
1104               if(n1 > 1) { kill I_for_fork; }
1105               return(J);
1106            }
1107            if((variant == 4) || (variant == 5) || (variant == 6))
1108            {
1109               tt = timer; rt = rtimer;
1110               K = std(J);
1111
1112               if(printlevel >= 10)
1113               {
1114                  "CPU-time for last std-computation is "
1115                  +string(timer - tt)+" seconds.";
1116                  "Real-time for last std-computation is "
1117                  +string(rtimer - rt)+" seconds.";
1118               }
1119
1120               if(size(reduce(K,J)) == 0)
1121               {
1122                  option(set, opt);
1123                  if(n1 > 1) { kill I_for_fork; }
1124                  return(J);
1125               }
1126            }
1127            if(pd > 2) { "pTest o.k. but result wrong"; }
1128         }
1129         if(pd > 2) { "pTest o.k. but result wrong"; }
1130      }
1131
1132//--------------  We do not already have a standard basis of I  ----------------
1133//-----------  Therefore do the main computation for more primes  --------------
1134
1135      T1 = H;
1136      T2 = N;
1137      index = 2;
1138
1139      j = size(L) + 1;
1140      L = primeList(I,n3,L);
1141
1142      if(n1 > 1)
1143      {
1144         for(i = 1; i <= n1; i++)
1145         {
1146            open(l(i));
1147            if((variant == 1) || (variant == 3) ||
1148               (variant == 4) || (variant == 6))
1149            {
1150               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1151                                                     eval(variant), eval(hi))));
1152            }
1153            if((variant == 2) || (variant == 5))
1154            {
1155               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1156                                                     eval(variant))));
1157            }
1158         }
1159         j = j + n1;
1160         k = 0;
1161      }
1162   }
1163}
1164example
1165{ "EXAMPLE:"; echo = 2;
1166   ring R1 = 0,(x,y,z,t),dp;
1167   ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
1168   ideal J = modStd(I);
1169   J;
1170   I = homog(I,t);
1171   J = modStd(I);
1172   J;
1173
1174   ring R2 = 0,(x,y,z),ds;
1175   ideal I = jacob(x5+y6+z7+xyz);
1176   ideal J1 = modStd(I,1,0);
1177   J1;
1178
1179   ring R3 = 0,x(1..4),lp;
1180   ideal I = cyclic(4);
1181   ideal J1 = modStd(I,1);
1182   ideal J2 = modStd(I,1,0);
1183   size(reduce(J1,J2));
1184   size(reduce(J2,J1));
1185
1186   /*
1187   ring R4 = 0,x(1..4),wp(1,-1,-1,1);
1188   ideal I = cyclic(4);
1189   ideal J1 = modStd(I,1,0);
1190   */
1191}
1192
1193////////////////////////////////////////////////////////////////////////////////
1194
1195proc modS(ideal I, list L, list #)
1196"USAGE:  modS(I,L); I ideal, L intvec of primes
1197         if size(#)>0 std is used instead of groebner
1198RETURN:  an ideal which is with high probability a standard basis
1199NOTE:    This procedure is designed for fast experiments.
1200         It is not tested whether the result is a standard basis.
1201         It is not tested whether the result generates I.
1202EXAMPLE: example modS; shows an example
1203"
1204{
1205   int j;
1206   bigint N = 1;
1207   def R0 = basering;
1208   ideal J;
1209   list T;
1210   list rl = ringlist(R0);
1211   if((npars(R0)>0) || (rl[1]>0))
1212   {
1213      ERROR("Characteristic of basering should be zero.");
1214   }
1215   for(j = 1; j <= size(L); j++)
1216   {
1217      N = N*L[j];
1218      rl[1] = L[j];
1219      def @r = ring(rl);
1220      setring @r;
1221      ideal I = fetch(R0,I);
1222      if(size(#) > 0)
1223      {
1224         I = std(I);
1225      }
1226      else
1227      {
1228         I = groebner(I);
1229      }
1230      setring R0;
1231      T[j] = fetch(@r,I);
1232      kill @r;
1233   }
1234   L = deleteUnluckyPrimes(T,L,homog(I));
1235   // unlucky if and only if the leading ideal is wrong
1236   J = farey(chinrem(L[1],L[2]),N);
1237   attrib(J,"isSB",1);
1238   return(J);
1239}
1240example
1241{ "EXAMPLE:"; echo = 2;
1242   list L = 3,5,11,13,181,32003;
1243   ring r = 0,(x,y,z,t),dp;
1244   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1245   I = homog(I,t);
1246   ideal J = modS(I,L);
1247   J;
1248}
1249
1250////////////////////////////////////////////////////////////////////////////////
1251
1252proc modHenselStd(ideal I, list #)
1253"USAGE:  modHenselStd(I);
1254RETURN:  a standard basis of I;
1255NOTE:    The procedure computes a standard basis of I (over the rational
1256         numbers) by using  modular computations and Hensellifting.
1257         For further experiments see procedure modS.
1258EXAMPLE: example modHenselStd; shows an example
1259"
1260{
1261   int i,j;
1262
1263   bigint p = 2134567879;
1264   if(size(#)!=0) { p=#[1]; }
1265   while(!primeTest(I,p))
1266   {
1267      p = prime(random(2000000000,2134567879));
1268   }
1269
1270   def R = basering;
1271   module F,PrevG,PrevZ,Z2;
1272   ideal testG,testG1,G1,G2,G3,Gp;
1273   list L = p;
1274   list rl = ringlist(R);
1275   rl[1] = int(p);
1276
1277   def S = ring(rl);
1278   setring S;
1279   option(redSB);
1280   module Z,M,Z2;
1281   ideal I = imap(R,I);
1282   ideal Gp,G1,G2,G3;
1283   Gp,Z = liftstd1(I);
1284   attrib(Gp,"isSB",1);
1285   module ZZ = syz(I);
1286   attrib(ZZ,"isSB",1);
1287   Z = reduce(Z,ZZ);
1288
1289   setring R;
1290   Gp = imap(S,Gp);
1291   PrevZ = imap(S,Z);
1292   PrevG = module(Gp);
1293   F = module(I);
1294   testG = farey(Gp,p);
1295   attrib(testG,"isSB",1);
1296   while(1)
1297   {
1298      i++;
1299      G1 = ideal(1/(p^i) * sum(F*PrevZ,(-1)*PrevG));
1300      setring S;
1301      G1 = imap(R,G1);
1302      G2 = reduce(G1,Gp);
1303      G3 = sum(G1,(-1)*G2);
1304      M = lift(Gp,G3);
1305      Z2 = (-1)*Z*M;
1306
1307      setring R;
1308      G2 = imap(S,G2);
1309      Z2 = imap(S,Z2);
1310      PrevG = sum(PrevG, module(p^i*G2));
1311      PrevZ = sum(PrevZ, multiply(poly(p^i),Z2));
1312      testG1 = farey(ideal(PrevG),p^(i+1));
1313      attrib(testG1,"isSB",1);
1314      if(size(reduce(testG1,testG)) == 0)
1315      {
1316         if(size(reduce(I,testG1)) == 0)       // I is in testG1
1317         {
1318            if(pTestSB(I,testG1,L,2))
1319            {
1320               G3 = std(testG1);               // testG1 is SB
1321               if(size(reduce(G3,testG1)) == 0)
1322               {
1323                  return(G3);
1324               }
1325            }
1326         }
1327      }
1328      testG = testG1;
1329      attrib(testG,"isSB",1);
1330   }
1331}
1332example
1333{ "EXAMPLE:"; echo = 2;
1334   ring r = 0,(x,y,z),dp;
1335   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1336   ideal J = modHenselStd(I);
1337   J;
1338}
1339
1340////////////////////////////////////////////////////////////////////////////////
1341
1342static proc sum(list #)
1343{
1344   if(typeof(#[1])=="ideal")
1345   {
1346      ideal M;
1347   }
1348   else
1349   {
1350      module M;
1351   }
1352
1353   int i;
1354   for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; }
1355   return(M);
1356}
1357
1358////////////////////////////////////////////////////////////////////////////////
1359
1360static proc multiply(poly p, list #)
1361{
1362   if(typeof(#[1])=="ideal")
1363   {
1364      ideal M;
1365   }
1366   else
1367   {
1368      module M;
1369   }
1370
1371   int i;
1372   for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; }
1373   return(M);
1374}
1375
1376
1377////////////////////////////// further examples ////////////////////////////////
1378
1379/*
1380ring r = 0, (x,y,z), lp;
1381poly s1 = 5x3y2z+3y3x2z+7xy2z2;
1382poly s2 = 3xy2z2+x5+11y2z2;
1383poly s3 = 4xyz+7x3+12y3+1;
1384poly s4 = 3x3-4y3+yz2;
1385ideal i =  s1, s2, s3, s4;
1386
1387ring r = 0, (x,y,z), lp;
1388poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
1389poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
1390poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
1391poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
1392ideal i =  s1, s2, s3, s4;
1393
1394ring r = 0, (x,y,z), lp;
1395poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
1396poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
1397poly s3 = 8x3 + 12y3 + xz2 + 3;
1398poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
1399ideal i = s1, s2, s3, s4;
1400
1401int n = 6;
1402ring r = 0,(x(1..n)),lp;
1403ideal i = cyclic(n);
1404ring s = 0, (x(1..n),t), lp;
1405ideal i = imap(r,i);
1406i = homog(i,t);
1407
1408ring r = 0, (x(1..4),s), (dp(4),dp);
1409poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
1410poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
1411poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
1412poly 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);
1413poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
1414ideal i =  s1, s2, s3, s4, s5;
1415
1416ring r = 0, (x,y,z), ds;
1417int a = 16;
1418int b = 15;
1419int c = 4;
1420int t = 1;
1421poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1422         +x^(c-2)*y^c*(y2+t*x)^2;
1423ideal i = jacob(f);
1424
1425ring r = 0, (x,y,z), ds;
1426int a = 25;
1427int b = 25;
1428int c = 5;
1429int t = 1;
1430poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1431         +x^(c-2)*y^c*(y2+t*x)^2;
1432ideal i = jacob(f),f;
1433
1434ring r = 0, (x,y,z), ds;
1435int a = 10;
1436poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
1437ideal i = jacob(f);
1438
1439ring r = 0, (x,y,z), ds;
1440int a = 6;
1441int b = 8;
1442int c = 10;
1443int alpha = 5;
1444int beta = 5;
1445int t = 1;
1446poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
1447         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
1448ideal i = jacob(f);
1449*/
1450
Note: See TracBrowser for help on using the repository browser.