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

spielwiese
Last change on this file since fd0df1 was fd0df1, checked in by Stefan Steidel <steidel@…>, 12 years ago
Some optimization regarding the number of processes used im modStd, some minor documentation improvements. git-svn-id: file:///usr/local/Singular/svn/trunk@14273 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 36.6 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((numerator(leadcoef(I[i][j])) mod p) == 0) { return(0); }
371         if((denominator(leadcoef(I[i][j])) mod p) == 0) { return(0); }
372      }
373   }
374   return(1);
375}
376
377////////////////////////////////////////////////////////////////////////////////
378
379proc primeList(ideal I, int n, list #)
380"USAGE:  primeList(I,n); ( resp. primeList(I,n,L); ) I ideal, n integer
381RETURN:  the intvec of n greatest primes  <= 2147483647 (resp. n greatest primes
382         < L[size(L)] union with L) such that none of these primes divides any
383         coefficient occuring in I
384EXAMPLE: example primList; shows an example
385"
386{
387   intvec L;
388   int i,p;
389   if(size(#) == 0)
390   {
391      p = 2147483647;
392      while(!primeTest(I,p))
393      {
394         p = prime(p-1);
395         if(p == 2) { ERROR("no more primes"); }
396      }
397      L[1] = p;
398   }
399   else
400   {
401      L = #[1];
402      p = prime(L[size(L)]-1);
403      while(!primeTest(I,p))
404      {
405         p = prime(p-1);
406         if(p == 2) { ERROR("no more primes"); }
407      }
408      L[size(L)+1] = p;
409   }
410   if(p == 2) { ERROR("no more primes"); }
411   for(i = 2; i <= n; i++)
412   {
413      p = prime(p-1);
414      while(!primeTest(I,p))
415      {
416         p = prime(p-1);
417         if(p == 2) { ERROR("no more primes"); }
418      }
419      L[size(L)+1] = p;
420   }
421   return(L);
422}
423example
424{ "EXAMPLE:"; echo = 2;
425   ring r = 0,(x,y,z),dp;
426   ideal I = 2147483647x+y, z-181;
427   intvec L = primeList(I,10);
428   size(L);
429   L[1];
430   L[size(L)];
431   L = primeList(I,5,L);
432   size(L);
433   L[size(L)];
434}
435
436////////////////////////////////////////////////////////////////////////////////
437
438static proc liftstd1(ideal I)
439{
440   def R = basering;
441   list rl = ringlist(R);
442   list ordl = rl[3];
443
444   int i;
445   for(i = 1; i <= size(ordl); i++)
446   {
447      if((ordl[i][1] == "C") || (ordl[i][1] == "c"))
448      {
449         ordl = delete(ordl, i);
450         break;
451      }
452   }
453
454   ordl = insert(ordl, list("c", 0));
455   rl[3] = ordl;
456   def newR = ring(rl);
457   setring newR;
458   ideal I = imap(R,I);
459
460   option(none);
461   option(prompt);
462
463   module M;
464   for(i = 1; i <= size(I); i++)
465   {
466      M = M + module(I[i]*gen(1) + gen(i+1));
467      M = M + module(gen(i+1));
468   }
469
470   module sM = std(M);
471
472   ideal sI;
473   if(attrib(R,"global"))
474   {
475      for(i = size(I)+1; i <= size(sM); i++)
476      {
477         sI[size(sI)+1] = sM[i][1];
478      }
479      matrix T = submat(sM,2..nrows(sM),size(I)+1..ncols(sM));
480   }
481   else
482   {
483      //"==========================================================";
484      //"WARNING: Algorithm is not applicable if ordering is mixed.";
485      //"==========================================================";
486      for(i = 1; i <= size(sM)-size(I); i++)
487      {
488         sI[size(sI)+1] = sM[i][1];
489      }
490      matrix T = submat(sM,2..nrows(sM),1..ncols(sM)-size(I));
491   }
492
493   setring R;
494   return(imap(newR,sI),imap(newR,T));
495}
496example
497{ "EXAMPLE:"; echo = 2;
498   ring R = 0,(x,y,z),dp;
499   poly f = x3+y7+z2+xyz;
500   ideal i = jacob(f);
501   matrix T;
502   ideal sm = liftstd(i,T);
503   sm;
504   print(T);
505   matrix(sm) - matrix(i)*T;
506
507
508   ring S = 32003, x(1..5), lp;
509   ideal I = cyclic(5);
510   ideal sI;
511   matrix T;
512   sI,T = liftstd1(I);
513   matrix(sI) - matrix(I)*T;
514}
515
516////////////////////////////////////////////////////////////////////////////////
517
518proc modpStd(ideal I, int p, int variant, list #)
519"USAGE:  modpStd(I,p,variant,#); I ideal, p integer, variant integer
520ASSUME:  If size(#) > 0, then #[1] is an intvec describing the Hilbert series.
521RETURN:  ideal - a standard basis of I mod p, integer - p
522NOTE:    The procedure computes a standard basis of the ideal I modulo p and
523         fetches the result to the basering. If size(#) > 0 the Hilbert driven
524         standard basis computation std(.,#[1]) is used instead of groebner.
525         The standard basis computation modulo p does also vary depending on the
526         integer variant, namely
527@*       - variant = 1/4: std(.,#[1]) resp. groebner,
528@*       - variant = 2/5: groebner,
529@*       - variant = 3/6: homog. - std(.,#[1]) resp. groebner - dehomog..
530EXAMPLE: example modpStd; shows an example
531"
532{
533   def R0 = basering;
534   list rl = ringlist(R0);
535   rl[1] = p;
536   def @r = ring(rl);
537   setring @r;
538   ideal i = fetch(R0,I);
539
540   option(redSB);
541
542   int t = timer;
543   if((variant == 1) || (variant == 4))
544   {
545      if(size(#) > 0)
546      {
547         i = std(i, #[1]);
548      }
549      else
550      {
551         i = groebner(i);
552      }
553   }
554
555   if((variant == 2) || (variant == 5))
556   {
557      i = groebner(i);
558   }
559
560   if((variant == 3) || (variant == 6))
561   {
562      list rl = ringlist(@r);
563      int nvar@r = nvars(@r);
564
565      int k;
566      intvec w;
567      for(k = 1; k <= nvar@r; k++)
568      {
569         w[k] = deg(var(k));
570      }
571      w[nvar@r + 1] = 1;
572
573      rl[2][nvar@r + 1] = "homvar";
574      rl[3][2][2] = w;
575
576      def HomR = ring(rl);
577      setring HomR;
578      ideal i = imap(@r, i);
579      i = homog(i, homvar);
580
581      if(size(#) > 0)
582      {
583         if(w == 1)
584         {
585            i = std(i, #[1]);
586         }
587         else
588         {
589            i = std(i, #[1], w);
590         }
591      }
592      else
593      {
594         i = groebner(i);
595      }
596
597      t = timer;
598      i = subst(i, homvar, 1);
599      i = simplify(i, 34);
600
601      setring @r;
602      i = imap(HomR, i);
603      i = interred(i);
604      kill HomR;
605   }
606
607   setring R0;
608   return(list(fetch(@r,i),p));
609}
610example
611{ "EXAMPLE:"; echo = 2;
612   ring r = 0, x(1..4), dp;
613   ideal I = cyclic(4);
614   int p = 181;
615   list P = modpStd(I,p,5);
616   P;
617
618   int q = 32003;
619   list Q = modpStd(I,q,2);
620   Q;
621}
622
623////////////////////////////// main procedures /////////////////////////////////
624
625proc modStd(ideal I, list #)
626"USAGE:  modStd(I); I ideal
627ASSUME:  If size(#) > 0, then # contains either 1, 2 or 4 integers such that
628@*       - #[1] is the number of available processors for the computation,
629@*       - #[2] is an optional parameter for the exactness of the computation,
630                if #[2] = 1, the procedure computes a standard basis for sure,
631@*       - #[3] is the number of primes until the first lifting,
632@*       - #[4] is the constant number of primes between two liftings until
633           the computation stops.
634RETURN:  a standard basis of I if no warning appears;
635NOTE:    The procedure computes a standard basis of I (over the rational
636         numbers) by using  modular methods. If a warning appears then the
637         result is a standard basis containing I and with high probability
638         a standard basis of I.
639         By default the procedure computes a standard basis of I for sure, but
640         if the optional parameter #[2] = 0, it computes a standard basis of I
641         with high probability and a warning appears at the end.
642         The procedure distinguishes between different variants for the standard
643         basis computation in positive characteristic depending on the ordering
644         of the basering, the parameter #[2] and if the ideal I is homogeneous.
645@*       - variant = 1, if I is homogeneous and exactness = 0,
646@*       - variant = 2, if I is not homogeneous, 1-block-ordering and
647                        exactness = 0,
648@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
649                        > 1 block) and exactness = 0,
650@*       - variant = 4, if I is homogeneous and exactness = 1,
651@*       - variant = 5, if I is not homogeneous, 1-block-ordering and
652                        exactness = 1,
653@*       - variant = 6, if I is not homogeneous, complicated ordering (lp or
654                        > 1 block) and exactness = 1.
655EXAMPLE: example modStd; shows an example
656"
657{
658   int TT = timer;
659   int RT = rtimer;
660
661   def R0 = basering;
662   list rl = ringlist(R0);
663   if((npars(R0) > 0) || (rl[1] > 0))
664   {
665      ERROR("Characteristic of basering should be zero, basering should
666             have no parameters.");
667   }
668
669   int index = 1;
670   int i,k,c;
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         if(n1 >= 10)
685         {
686            int n2 = n1 + 1;
687            int n3 = n1;
688         }
689         else
690         {
691            int n2 = 10;
692            int n3 = 10;
693         }
694      }
695      if(size(#) == 2)
696      {
697         int n1 = #[1];
698         int exactness = #[2];
699         if(n1 >= 10)
700         {
701            int n2 = n1 + 1;
702            int n3 = n1;
703         }
704         else
705         {
706            int n2 = 10;
707            int n3 = 10;
708         }
709      }
710      if(size(#) == 4)
711      {
712         int n1 = #[1];
713         int exactness = #[2];
714         if(n1 >= #[3])
715         {
716            int n2 = n1 + 1;
717         }
718         else
719         {
720            int n2 = #[3];
721         }
722         if(n1 >= #[4])
723         {
724            int n3 = n1;
725         }
726         else
727         {
728            int n3 = #[4];
729         }
730      }
731   }
732   else
733   {
734      int n1 = 1;
735      int exactness = 1;
736      int n2 = 10;
737      int n3 = 10;
738   }
739
740   if(printlevel >= 10)
741   {
742      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3);
743   }
744
745//-------------------------  Save current options  -----------------------------
746   intvec opt = option(get);
747
748   option(redSB);
749
750//--------------------  Initialize the list of primes  -------------------------
751   intvec L = primeList(I,n2);
752   L[5] = prime(random(an,en));
753
754//---------------------  Decide which variant to take  -------------------------
755   int variant;
756   int h = homog(I);
757
758   int tt = timer;
759   int rt = rtimer;
760
761   if(!mixedTest())
762   {
763      if(h)
764      {
765         if(exactness == 0)
766         {
767            variant = 1;
768            if(printlevel >= 10) { "variant = 1"; }
769         }
770         if(exactness == 1)
771         {
772            variant = 4;
773            if(printlevel >= 10) { "variant = 4"; }
774         }
775         rl[1] = L[5];
776         def @r = ring(rl);
777         setring @r;
778         def @s = changeord("dp");
779         setring @s;
780         ideal I = std(fetch(R0,I));
781         intvec hi = hilb(I,1);
782         setring R0;
783         kill @r,@s;
784      }
785      else
786      {
787         string ordstr_R0 = ordstr(R0);
788         int neg = 1 - attrib(R0,"global");
789
790         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
791         {
792            if(exactness == 0)
793            {
794               variant = 2;
795               if(printlevel >= 10) { "variant = 2"; }
796            }
797            if(exactness == 1)
798            {
799               variant = 5;
800               if(printlevel >= 10) { "variant = 5"; }
801            }
802         }
803         else
804         {
805            string order;
806            if(system("nblocks") <= 2)
807            {
808               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
809                                       + find(ordstr_R0, "rp") <= 0)
810               {
811                  order = "simple";
812               }
813            }
814
815            if((order == "simple") || (size(rl) > 4))
816            {
817               if(exactness == 0)
818               {
819                  variant = 2;
820                  if(printlevel >= 10) { "variant = 2"; }
821               }
822               if(exactness == 1)
823               {
824                  variant = 5;
825                  if(printlevel >= 10) { "variant = 5"; }
826               }
827            }
828            else
829            {
830               if(exactness == 0)
831               {
832                  variant = 3;
833                  if(printlevel >= 10) { "variant = 3"; }
834               }
835               if(exactness == 1)
836               {
837                  variant = 6;
838                  if(printlevel >= 10) { "variant = 6"; }
839               }
840
841               rl[1] = L[5];
842               def @r = ring(rl);
843               setring @r;
844               int nvar@r = nvars(@r);
845               intvec w;
846               for(i = 1; i <= nvar@r; i++)
847               {
848                  w[i] = deg(var(i));
849               }
850               w[nvar@r + 1] = 1;
851
852               list hiRi = hilbRing(fetch(R0,I),w);
853               intvec W = hiRi[2];
854               def @s = hiRi[1];
855               setring @s;
856
857               Id(1) = std(Id(1));
858               intvec hi = hilb(Id(1), 1, W);
859
860               setring R0;
861               kill @r,@s;
862            }
863         }
864      }
865   }
866   else
867   {
868      if(exactness == 1) { return(groebner(I)); }
869      if(h)
870      {
871         variant = 1;
872         if(printlevel >= 10) { "variant = 1"; }
873         rl[1] = L[5];
874         def @r = ring(rl);
875         setring @r;
876         def @s = changeord("dp");
877         setring @s;
878         ideal I = std(fetch(R0,I));
879         intvec hi = hilb(I,1);
880         setring R0;
881         kill @r,@s;
882      }
883      else
884      {
885         string ordstr_R0 = ordstr(R0);
886         int neg = 1 - attrib(R0,"global");
887
888         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
889         {
890            variant = 2;
891            if(printlevel >= 10) { "variant = 2"; }
892         }
893         else
894         {
895            string order;
896            if(system("nblocks") <= 2)
897            {
898               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
899                                       + find(ordstr_R0, "rp") <= 0)
900               {
901                  order = "simple";
902               }
903            }
904
905            if((order == "simple") || (size(rl) > 4))
906            {
907               variant = 2;
908               if(printlevel >= 10) { "variant = 2"; }
909            }
910            else
911            {
912               variant = 3;
913               if(printlevel >= 10) { "variant = 3"; }
914
915               rl[1] = L[5];
916               def @r = ring(rl);
917               setring @r;
918               int nvar@r = nvars(@r);
919               intvec w;
920               for(i = 1; i <= nvar@r; i++)
921               {
922                  w[i] = deg(var(i));
923               }
924               w[nvar@r + 1] = 1;
925
926               list hiRi = hilbRing(fetch(R0,I),w);
927               intvec W = hiRi[2];
928               def @s = hiRi[1];
929               setring @s;
930
931               Id(1) = std(Id(1));
932               intvec hi = hilb(Id(1), 1, W);
933
934               setring R0;
935               kill @r,@s;
936            }
937         }
938      }
939   }
940
941   list P,T1,T2,T3,LL;
942
943   ideal J,K,H;
944
945//-----  If there is more than one processor available, we parallelize the  ----
946//-----  main standard basis computations in positive characteristic        ----
947
948   if(n1 > 1)
949   {
950      ideal I_for_fork = I;
951      export(I_for_fork);           // I available for each link
952
953//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
954//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
955
956      for(i = 1; i <= n1; i++)
957      {
958         //link l(i) = "MPtcp:fork";
959         link l(i) = "ssi:fork";
960         open(l(i));
961         if((variant == 1) || (variant == 3) ||
962            (variant == 4) || (variant == 6))
963         {
964            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
965                                                  eval(variant), eval(hi))));
966         }
967         if((variant == 2) || (variant == 5))
968         {
969            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
970                                                  eval(variant))));
971         }
972      }
973
974      int t = timer;
975      if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
976      {
977         P = modpStd(I_for_fork, L[1], variant, hi);
978      }
979      if((variant == 2) || (variant == 5))
980      {
981         P = modpStd(I_for_fork, L[1], variant);
982      }
983      t = timer - t;
984      if(t > 60) { t = 60; }
985      int i_sleep = system("sh", "sleep "+string(t));
986      T1[1] = P[1];
987      T2[1] = bigint(P[2]);
988      index++;
989
990      j = j + n1 + 1;
991   }
992
993//--------------  Main standard basis computations in positive  ----------------
994//----------------------  characteristic start here  ---------------------------
995
996   while(1)
997   {
998      tt = timer; rt = rtimer;
999     
1000      if(printlevel >= 10) { "size(L) = "+string(size(L)); }
1001         
1002      if(n1 > 1)
1003      {
1004         while(j <= size(L) + 1)
1005         {
1006            for(i = 1; i <= n1; i++)
1007            {
1008               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
1009               if(status(l(i), "read", "ready"))
1010               {
1011                  //--- read the result from l(i) ---
1012                  P = read(l(i));
1013                  T1[index] = P[1];
1014                  T2[index] = bigint(P[2]);
1015                  index++;
1016
1017                  if(j <= size(L))
1018                  {
1019                     if((variant == 1) || (variant == 3) ||
1020                        (variant == 4) || (variant == 6))
1021                     {
1022                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
1023                                                  eval(variant), eval(hi))));
1024                        j++;
1025                     }
1026                     if((variant == 2) || (variant == 5))
1027                     {
1028                        write(l(i), quote(modpStd(I_for_fork,
1029                                                  eval(L[j]), eval(variant))));
1030                        j++;
1031                     }
1032                  }
1033                  else
1034                  {
1035                     k++;
1036                     close(l(i));
1037                  }
1038               }
1039            }
1040            //--- k describes the number of closed links ---
1041            if(k == n1)
1042            {
1043               j++;
1044            }
1045            i_sleep = system("sh", "sleep "+string(t));
1046         }
1047      }
1048      else
1049      {
1050         while(j <= size(L))
1051         {
1052            if((variant == 1) || (variant == 3) ||
1053               (variant == 4) || (variant == 6))
1054            {
1055               P = modpStd(I, L[j], variant, hi);
1056            }
1057            if((variant == 2) || (variant == 5))
1058            {
1059               P = modpStd(I, L[j], variant);
1060            }
1061
1062            T1[index] = P[1];
1063            T2[index] = bigint(P[2]);
1064            index++;
1065            j++;
1066         }
1067      }
1068
1069      if(printlevel >= 10)
1070      {
1071         "CPU-time for computing list is "+string(timer - tt)+" seconds.";
1072         "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
1073      }
1074
1075//------------------------  Delete unlucky primes  -----------------------------
1076//-------------  unlucky if and only if the leading ideal is wrong  ------------
1077
1078      LL = deleteUnluckyPrimes(T1,T2,h);
1079      T1 = LL[1];
1080      T2 = LL[2];
1081
1082//-------------------  Now all leading ideals are the same  --------------------
1083//-------------------  Lift results to basering via farey  ---------------------
1084
1085      tt = timer;
1086      N = T2[1];
1087      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
1088      H = chinrem(T1,T2);
1089      J = farey(H,N);
1090      if(printlevel >= 10) { "Lifting-process takes "+string(timer - tt)
1091                             +" seconds"; }
1092
1093//----------------  Test if we already have a standard basis of I --------------
1094
1095      tt = timer; rt = rtimer;
1096      if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
1097      {
1098         pTest = pTestSB(I,J,L,variant,hi);
1099      }
1100      if((variant == 2) || (variant == 5))
1101      {
1102         pTest = pTestSB(I,J,L,variant);
1103      }
1104
1105      if(printlevel >= 10)
1106      {
1107         "CPU-time for pTest is "+string(timer - tt)+" seconds.";
1108         "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
1109      }
1110
1111      if(pTest)
1112      {
1113         if(printlevel >= 10)
1114         {
1115            "CPU-time for computation without final tests is "
1116            +string(timer - TT)+" seconds.";
1117            "Real-time for computation without final tests is "
1118            +string(rtimer - RT)+" seconds.";
1119         }
1120
1121         attrib(J,"isSB",1);
1122         tt = timer; rt = rtimer;
1123         sizeTest = 1 - isIncluded(I,J,n1);
1124
1125         if(printlevel >= 10)
1126         {
1127            "CPU-time for checking if I subset <G> is "
1128            +string(timer - tt)+" seconds.";
1129            "Real-time for checking if I subset <G> is "
1130            +string(rtimer - rt)+" seconds.";
1131         }
1132
1133         if(sizeTest == 0)
1134         {
1135            if((variant == 1) || (variant == 2) || (variant == 3))
1136            {
1137               "===========================================================";
1138               "WARNING: Output might not be a standard basis of the input.";
1139               "===========================================================";
1140               option(set, opt);
1141               if(n1 > 1) { kill I_for_fork; }
1142               return(J);
1143            }
1144            if((variant == 4) || (variant == 5) || (variant == 6))
1145            {
1146               tt = timer; rt = rtimer;
1147               K = std(J);
1148
1149               if(printlevel >= 10)
1150               {
1151                  "CPU-time for last std-computation is "
1152                  +string(timer - tt)+" seconds.";
1153                  "Real-time for last std-computation is "
1154                  +string(rtimer - rt)+" seconds.";
1155               }
1156
1157               if(size(reduce(K,J)) == 0)
1158               {
1159                  option(set, opt);
1160                  if(n1 > 1) { kill I_for_fork; }
1161                  return(J);
1162               }
1163            }
1164         }
1165      }
1166
1167//--------------  We do not already have a standard basis of I  ----------------
1168//-----------  Therefore do the main computation for more primes  --------------
1169
1170      T1 = H;
1171      T2 = N;
1172      index = 2;
1173
1174      j = size(L) + 1;
1175      L = primeList(I,n3,L);
1176
1177      if(n1 > 1)
1178      {
1179         for(i = 1; i <= n1; i++)
1180         {
1181            open(l(i));
1182            if((variant == 1) || (variant == 3) ||
1183               (variant == 4) || (variant == 6))
1184            {
1185               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1186                                                     eval(variant), eval(hi))));
1187            }
1188            if((variant == 2) || (variant == 5))
1189            {
1190               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
1191                                                     eval(variant))));
1192            }
1193         }
1194         j = j + n1;
1195         k = 0;
1196      }
1197   }
1198}
1199example
1200{ "EXAMPLE:"; echo = 2;
1201   ring R1 = 0,(x,y,z,t),dp;
1202   ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
1203   ideal J = modStd(I);
1204   J;
1205   I = homog(I,t);
1206   J = modStd(I);
1207   J;
1208
1209   ring R2 = 0,(x,y,z),ds;
1210   ideal I = jacob(x5+y6+z7+xyz);
1211   ideal J1 = modStd(I,1,0);
1212   J1;
1213
1214   ring R3 = 0,x(1..4),lp;
1215   ideal I = cyclic(4);
1216   ideal J1 = modStd(I,1);
1217   ideal J2 = modStd(I,1,0);
1218   size(reduce(J1,J2));
1219   size(reduce(J2,J1));
1220}
1221
1222////////////////////////////////////////////////////////////////////////////////
1223
1224proc modS(ideal I, list L, list #)
1225"USAGE:  modS(I,L); I ideal, L intvec of primes
1226         if size(#)>0 std is used instead of groebner
1227RETURN:  an ideal which is with high probability a standard basis
1228NOTE:    This procedure is designed for fast experiments.
1229         It is not tested whether the result is a standard basis.
1230         It is not tested whether the result generates I.
1231EXAMPLE: example modS; shows an example
1232"
1233{
1234   int j;
1235   bigint N = 1;
1236   def R0 = basering;
1237   ideal J;
1238   list T;
1239   list rl = ringlist(R0);
1240   if((npars(R0)>0) || (rl[1]>0))
1241   {
1242      ERROR("Characteristic of basering should be zero.");
1243   }
1244   for(j = 1; j <= size(L); j++)
1245   {
1246      N = N*L[j];
1247      rl[1] = L[j];
1248      def @r = ring(rl);
1249      setring @r;
1250      ideal I = fetch(R0,I);
1251      if(size(#) > 0)
1252      {
1253         I = std(I);
1254      }
1255      else
1256      {
1257         I = groebner(I);
1258      }
1259      setring R0;
1260      T[j] = fetch(@r,I);
1261      kill @r;
1262   }
1263   L = deleteUnluckyPrimes(T,L,homog(I));
1264   // unlucky if and only if the leading ideal is wrong
1265   J = farey(chinrem(L[1],L[2]),N);
1266   attrib(J,"isSB",1);
1267   return(J);
1268}
1269example
1270{ "EXAMPLE:"; echo = 2;
1271   list L = 3,5,11,13,181,32003;
1272   ring r = 0,(x,y,z,t),dp;
1273   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1274   I = homog(I,t);
1275   ideal J = modS(I,L);
1276   J;
1277}
1278
1279////////////////////////////////////////////////////////////////////////////////
1280
1281proc modHenselStd(ideal I, list #)
1282"USAGE:  modHenselStd(I);
1283RETURN:  a standard basis of I;
1284NOTE:    The procedure computes a standard basis of I (over the rational
1285         numbers) by using  modular computations and Hensellifting.
1286         For further experiments see procedure modS.
1287EXAMPLE: example modHenselStd; shows an example
1288"
1289{
1290   int i,j;
1291
1292   bigint p = 2134567879;
1293   if(size(#)!=0) { p=#[1]; }
1294   while(!primeTest(I,p))
1295   {
1296      p = prime(random(2000000000,2134567879));
1297   }
1298
1299   def R = basering;
1300   module F,PrevG,PrevZ,Z2;
1301   ideal testG,testG1,G1,G2,G3,Gp;
1302   list L = p;
1303   list rl = ringlist(R);
1304   rl[1] = int(p);
1305
1306   def S = ring(rl);
1307   setring S;
1308   option(redSB);
1309   module Z,M,Z2;
1310   ideal I = imap(R,I);
1311   ideal Gp,G1,G2,G3;
1312   Gp,Z = liftstd1(I);
1313   attrib(Gp,"isSB",1);
1314   module ZZ = syz(I);
1315   attrib(ZZ,"isSB",1);
1316   Z = reduce(Z,ZZ);
1317
1318   setring R;
1319   Gp = imap(S,Gp);
1320   PrevZ = imap(S,Z);
1321   PrevG = module(Gp);
1322   F = module(I);
1323   testG = farey(Gp,p);
1324   attrib(testG,"isSB",1);
1325   while(1)
1326   {
1327      i++;
1328      G1 = ideal(1/(p^i) * sum(F*PrevZ,(-1)*PrevG));
1329      setring S;
1330      G1 = imap(R,G1);
1331      G2 = reduce(G1,Gp);
1332      G3 = sum(G1,(-1)*G2);
1333      M = lift(Gp,G3);
1334      Z2 = (-1)*Z*M;
1335
1336      setring R;
1337      G2 = imap(S,G2);
1338      Z2 = imap(S,Z2);
1339      PrevG = sum(PrevG, module(p^i*G2));
1340      PrevZ = sum(PrevZ, multiply(poly(p^i),Z2));
1341      testG1 = farey(ideal(PrevG),p^(i+1));
1342      attrib(testG1,"isSB",1);
1343      if(size(reduce(testG1,testG)) == 0)
1344      {
1345         if(size(reduce(I,testG1)) == 0)       // I is in testG1
1346         {
1347            if(pTestSB(I,testG1,L,2))
1348            {
1349               G3 = std(testG1);               // testG1 is SB
1350               if(size(reduce(G3,testG1)) == 0)
1351               {
1352                  return(G3);
1353               }
1354            }
1355         }
1356      }
1357      testG = testG1;
1358      attrib(testG,"isSB",1);
1359   }
1360}
1361example
1362{ "EXAMPLE:"; echo = 2;
1363   ring r = 0,(x,y,z),dp;
1364   ideal I = 3x3+x2+1,11y5+y3+2,5z4+z2+4;
1365   ideal J = modHenselStd(I);
1366   J;
1367}
1368
1369////////////////////////////////////////////////////////////////////////////////
1370
1371static proc sum(list #)
1372{
1373   if(typeof(#[1])=="ideal")
1374   {
1375      ideal M;
1376   }
1377   else
1378   {
1379      module M;
1380   }
1381
1382   int i;
1383   for(i = 1; i <= ncols(#[1]); i++) { M[i] = #[1][i] + #[2][i]; }
1384   return(M);
1385}
1386
1387////////////////////////////////////////////////////////////////////////////////
1388
1389static proc multiply(poly p, list #)
1390{
1391   if(typeof(#[1])=="ideal")
1392   {
1393      ideal M;
1394   }
1395   else
1396   {
1397      module M;
1398   }
1399
1400   int i;
1401   for(i = 1; i <= ncols(#[1]); i++) { M[i] = p * #[1][i]; }
1402   return(M);
1403}
1404
1405
1406////////////////////////////// further examples ////////////////////////////////
1407
1408/*
1409ring r = 0, (x,y,z), lp;
1410poly s1 = 5x3y2z+3y3x2z+7xy2z2;
1411poly s2 = 3xy2z2+x5+11y2z2;
1412poly s3 = 4xyz+7x3+12y3+1;
1413poly s4 = 3x3-4y3+yz2;
1414ideal i =  s1, s2, s3, s4;
1415
1416ring r = 0, (x,y,z), lp;
1417poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
1418poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
1419poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
1420poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
1421ideal i =  s1, s2, s3, s4;
1422
1423ring r = 0, (x,y,z), lp;
1424poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
1425poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
1426poly s3 = 8x3 + 12y3 + xz2 + 3;
1427poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
1428ideal i = s1, s2, s3, s4;
1429
1430int n = 6;
1431ring r = 0,(x(1..n)),lp;
1432ideal i = cyclic(n);
1433ring s = 0, (x(1..n),t), lp;
1434ideal i = imap(r,i);
1435i = homog(i,t);
1436
1437ring r = 0, (x(1..4),s), (dp(4),dp);
1438poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
1439poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
1440poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
1441poly 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);
1442poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
1443ideal i =  s1, s2, s3, s4, s5;
1444
1445ring r = 0, (x,y,z), ds;
1446int a = 16;
1447int b = 15;
1448int c = 4;
1449int t = 1;
1450poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1451         +x^(c-2)*y^c*(y2+t*x)^2;
1452ideal i = jacob(f);
1453
1454ring r = 0, (x,y,z), ds;
1455int a = 25;
1456int b = 25;
1457int c = 5;
1458int t = 1;
1459poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
1460         +x^(c-2)*y^c*(y2+t*x)^2;
1461ideal i = jacob(f),f;
1462
1463ring r = 0, (x,y,z), ds;
1464int a = 10;
1465poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
1466ideal i = jacob(f);
1467
1468ring r = 0, (x,y,z), ds;
1469int a = 6;
1470int b = 8;
1471int c = 10;
1472int alpha = 5;
1473int beta = 5;
1474int t = 1;
1475poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
1476         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
1477ideal i = jacob(f);
1478*/
1479
Note: See TracBrowser for help on using the repository browser.