source: git/Singular/LIB/symodstd.lib @ b0732eb

spielwiese
Last change on this file since b0732eb was b0732eb, checked in by Hans Schoenemann <hannes@…>, 12 years ago
add: new libs from master
  • Property mode set to 100644
File size: 41.9 KB
Line 
1////////////////////////////////////////////////////////////////////////////////
2category = "Commutative Algebra";
3info="
4LIBRARY:  symodstd.lib    Procedures for computing Groebner basis of ideals
5                          being invariant under certain variable permutations.
6
7AUTHOR:   Stefan Steidel, steidel@mathematik.uni-kl.de
8
9OVERVIEW:
10
11   A library for computing the Groebner basis of an ideal in the polynomial
12   ring over the rational numbers, that is invariant under certain permutations
13   of the variables, using the symmetry and modular methods.
14   More precisely let I = <f1,...,fr> be an ideal in Q[x(1),...,x(n)] and
15   sigma a permutation of order k in Sym(n) such that sigma(I) = I.
16   We assume that sigma({f1,...,fr}) = {f1,...,fr}. This can always be obtained
17   by adding sigma(fi) to {f1,...,fr}.
18   To compute a standard basis of I we apply a modification of the modular
19   version of the standard basis algorithm (improving the calculations in
20   positive characteristic). Therefore we only allow primes p such that p-1 is
21   divisible by k. This guarantees the existance of a k-th primitive root of
22   unity in Z/pZ.
23
24PROCEDURES:
25 genSymId(I,sigma);      compute ideal J such that sigma(J) = J and J includes I
26 isSymmetric(I,sigma);   check if I is invariant under permutation sigma
27 primRoot(p,k);          int describing a k-th primitive root of unity in Z/pZ
28 eigenvalues(I,sigma);   list of eigenvalues of generators of I under sigma
29 symmStd(I,sigma);       standard basis of I using invariance of I under sigma
30 syModStd(I,sigma);      SB of I using modular methods and sigma(I) = I
31";
32
33LIB "brnoeth.lib";
34LIB "modstd.lib";
35LIB "parallel.lib";
36
37////////////////////////////////////////////////////////////////////////////////
38
39static proc mixedTest()
40"USAGE:  mixedTest();
41RETURN:  1 if ordering of basering is mixed, 0 else
42EXAMPLE: example mixedTest(); shows an example
43"
44{
45   int i,p,m;
46   for(i = 1; i <= nvars(basering); i++)
47   {
48      if(var(i) > 1)
49      {
50         p++;
51      }
52      else
53      {
54         m++;
55      }
56   }
57   if((p > 0) && (m > 0)) { return(1); }
58   return(0);
59}
60example
61{ "EXAMPLE:"; echo = 2;
62   ring R1 = 0,(x,y,z),dp;
63   mixedTest();
64   ring R2 = 31,(x(1..4),y(1..3)),(ds(4),lp(3));
65   mixedTest();
66   ring R3 = 181,x(1..9),(dp(5),lp(4));
67   mixedTest();
68}
69
70////////////////////////////////////////////////////////////////////////////////
71
72proc genSymId(ideal I, intvec sigma)
73"USAGE:  genSymId(I,sigma); I ideal, sigma intvec
74ASSUME:  size(sigma) = nvars(basering =: n
75RETURN:  ideal J such that sigma(J) = J and J includes I
76NOTE:    sigma is a permutation of the variables of the basering, i.e.
77@*       sigma: var(i) ----> var(sigma[i]), 1 <= i <= n.
78EXAMPLE: example genSymId; shows an example
79"
80{
81   if(nvars(basering) != size(sigma))
82   {
83      ERROR("The input is no permutation of the ring-variables!!");
84   }
85
86   int i;
87
88   ideal perm;
89   for(i = 1; i <= size(sigma); i++)
90   {
91      perm[size(perm)+1] = var(sigma[i]);
92   }
93
94   map f = basering, perm;
95   ideal J = I;
96   ideal helpJ = I;
97   for(i = 1; i <= order(sigma) - 1; i++)
98   {
99      helpJ = f(helpJ);
100      J = J, helpJ;
101   }
102
103   return(simplify(simplify(J,4),2));
104}
105example
106{ "EXAMPLE:"; echo = 2;
107   ring R = 0,(u,v,w,x,y),dp;
108   intvec pi = 2,3,4,5,1;
109   ideal I = u2v + x3y - w2;
110   genSymId(I,pi);
111}
112
113////////////////////////////////////////////////////////////////////////////////
114
115proc isSymmetric(ideal I, intvec sigma)
116"USAGE:  isSymmetric(I,sigma); I ideal, sigma intvec
117ASSUME:  size(sigma) = nvars(basering) =: n
118RETURN:  1, if the set of generators of I is invariant under sigma;
119@*       0, if the set of generators of I is not invariant under sigma
120NOTE:    sigma is a permutation of the variables of the basering, i.e.
121@*       sigma: var(i) ----> var(sigma[i]), 1 <= i <= n.
122EXAMPLE: example isSymmetric; shows an example
123"
124{
125   if(nvars(basering) != size(sigma))
126   {
127      ERROR("The input is no permutation of the ring-variables!!");
128   }
129
130   int i, j;
131
132   list L;
133   for(i = 1; i <= size(I); i++)
134   {
135      L[size(L)+1] = I[i];
136   }
137
138   ideal perm;
139   for(i = 1; i <= size(sigma); i++)
140   {
141      perm[size(perm)+1] = var(sigma[i]);
142   }
143
144   map f = basering, perm;
145   ideal J = f(I);
146
147   poly g;
148   while(size(L) > 0)
149   {
150      j = size(L);
151      g = L[1];
152
153      for(i = 1; i <= size(J); i++)
154      {
155         if(g - J[i] == 0)
156         {
157            L = delete(L, 1);
158            break;
159         }
160      }
161
162      if(j == size(L))
163      {
164         return(0);
165      }
166   }
167
168   return(1);
169}
170example
171{ "EXAMPLE:"; echo = 2;
172   ring R = 0,x(1..5),dp;
173   ideal I = cyclic(5);
174   intvec pi = 2,3,4,5,1;
175   isSymmetric(I,pi);
176   intvec tau = 2,5,1,4,3;
177   isSymmetric(I,tau);
178}
179
180////////////////////////////////////////////////////////////////////////////////
181
182static proc permute(intvec v, intvec P)
183{
184// permute the intvec v according to the permutation given by P
185
186   int s = size(v);
187   int n = size(P);
188   int i;
189   if(s < n)
190   {
191      for(i = s+1; i <= n; i = i+1)
192      {
193         v[i] = 0;
194      }
195      s = size(v);
196   }
197   intvec auxv = v;
198   for(i = 1; i <= n; i = i+1)
199   {
200      auxv[i] = v[P[i]];
201   }
202
203   return(auxv);
204}
205
206////////////////////////////////////////////////////////////////////////////////
207
208static proc order(intvec sigma)
209{
210// compute the order of sigma in Sym({1,...,n}) with n := size(sigma)
211
212   int ORD = 1;
213   intvec id = 1..size(sigma);
214   intvec tau = sigma;
215
216   while(tau != id)
217   {
218      tau = permute(tau, sigma);
219      ORD = ORD + 1;
220   }
221
222   return(ORD);
223}
224
225////////////////////////////////////////////////////////////////////////////////
226
227static proc modExpo(bigint x, bigint a, bigint n)
228{
229// compute x^a mod n
230
231   bigint z = 1;
232
233   while(a != 0)
234   {
235      while((a mod 2) == 0)
236      {
237         a = a div 2;
238         x = x^2 mod n;
239      }
240
241      a = a - 1;
242      z = (z*x) mod n;
243   }
244
245   return(z);
246}
247
248////////////////////////////////////////////////////////////////////////////////
249
250proc primRoot(int p, int k)
251"USAGE:  primRoot(p,k); p,k integers
252ASSUME:  p is a prime and k divides p-1.
253RETURN:  int: a k-th primitive root of unity in Z/pZ
254EXAMPLE: example primRoot; shows an example
255"
256{
257   if(k == 2)
258   {
259      return(-1);
260   }
261
262   if(p == 0)
263   {
264      return(0);
265   }
266
267   int i, j;
268
269   if(((p-1) mod k) != 0)
270   {
271      ERROR("There is no "+string(k)+"-th primitive root of unity "
272             +"in Z/"+string(p)+"Z.");
273      return(0);
274   }
275
276   list PF = primefactors(p-1)[1];
277
278   bigint a;
279
280   for(i = 2; i <= p-1; i++)
281   {
282      a = i;
283
284      for(j = 1; j <= size(PF); j++)
285      {
286         if(modExpo(a, (p-1) div PF[j], p) == 1)
287         {
288            break;
289         }
290      }
291
292      if(j == size(PF)+1)
293      {
294         a = modExpo(a, (p-1) div k, p);
295
296         string str = "int xi = "+string(a);
297         execute(str);
298
299         return(xi);
300      }
301   }
302}
303example
304{ "EXAMPLE:";  echo = 2;
305   primRoot(181,10);
306   ring R = 2147482801, x, lp;
307   number a = primRoot(2147482801,5);
308   a;
309   a^2;
310   a^3;
311   a^4;
312   a^5;
313}
314
315////////////////////////////////////////////////////////////////////////////////
316
317static proc permMat(intvec sigma, list #)
318{
319// compute an intmat such that i-th row describes sigma^i
320
321   int i;
322   int n = size(sigma);
323
324   if(size(#) == 0)
325   {
326      int ORD = order(sigma);
327   }
328   else
329   {
330      int ORD = #[1];
331   }
332
333   intmat P[ORD][n];
334   intvec sigmai = sigma;
335   for(i = 1; i <= ORD; i++)
336   {
337      P[i,1..n] = sigmai;
338      sigmai = permute(sigmai, sigma);
339   }
340
341   return(P);
342}
343
344////////////////////////////////////////////////////////////////////////////////
345
346static proc genTransId(intvec sigma, list #)
347{
348// list L of two ideals such that
349// - L[1] describes the transformation and
350// - L[2] describes the retransformation (inverse mapping).
351// ORD is the order of sigma in Sym({1,...,n}) with n := size(sigma) and
352// sigma: {1,...,n} ---> {1,...,n}: sigma(j) = sigma[j]. Since sigma is a
353// permutation of variables it induces an automorphism phi of the
354// basering, more precisely a linear variable transformation which is
355// generated by this procedure. In terms it holds:
356//
357//             phi : basering  --------->  basering
358//                    var(i)     |---->     L[1][i]
359//                    L[2][i]    <----|     var(i)
360
361   int n = nvars(basering);
362
363   if(n != size(sigma))
364   {
365      ERROR("The input is no permutation of the ring-variables!!");
366   }
367
368   int i, j, k;
369
370   if(size(#) == 0)
371   {
372      int CHAR = char(basering);
373      int ORD = order(sigma);
374
375      if((((CHAR - 1) mod ORD) != 0) && (CHAR > 0))
376      {
377         ERROR("Basering of characteristic "+string(CHAR)+" has no "
378                +string(ORD)+"-th primitive root of unity!!");
379      }
380
381      if((CHAR == 0) && (ORD > 2))
382      {
383         "========================================
384          ========================================";
385         "If basering really has a "+string(ORD)+"-th "
386          +"primitive root of unity then insert it as input!!";
387         "========================================
388          ========================================";
389         return(list());
390      }
391      else
392      {
393         int xi = primRoot(CHAR, ORD);
394         number a = xi;
395      }
396   }
397   else
398   {
399      int ORD = #[1];
400      number a = #[2];
401   }
402
403   intmat PERM = permMat(sigma,ORD);
404
405   ideal TR, RETR;
406   poly s_trans;
407   matrix C;
408
409//-------------- retransformation ideal RETR is generated here -----------------
410   for(i = 1; i <= n; i++)
411   {
412      for(j = 0; j < ORD; j++)
413      {
414         for(k = 1; k <= ORD; k++)
415         {
416            // for each variable var(i):
417            // s_trans^(j) = sum_{k=1}^{ORD} a^(k*j)*sigma^k(var(i))
418            // for j = 0,...,ORD-1
419            s_trans = s_trans + a^(k*j)*var(PERM[k,i]);
420         }
421         RETR = RETR + simplify(s_trans, 1);
422         s_trans = 0;
423      }
424   }
425
426//---------------- transformation ideal TR is generated here -------------------
427   for(i = 1; i <= n; i++)
428   {
429      for(j = 1; j <= size(RETR); j++)
430      {
431         C = coeffs(RETR[j], var(i));
432         if(nrows(C) > 1)
433         {
434            // var(j) = RETR[j] = sum_{i in J} c_ij*var(i), J subset {1,...,n},
435            // and therefore var(i) = (sum_{j} s(j)/c_ij)/#(summands)
436            s_trans = s_trans + var(j)/(C[nrows(C),1]);
437         }
438      }
439
440      TR = TR + s_trans/number(size(s_trans));
441      s_trans = 0;
442   }
443
444   return(list(TR,RETR));
445}
446
447////////////////////////////////////////////////////////////////////////////////
448
449proc eigenvalues(ideal I, intvec sigma)
450"USAGE:  eigenvalues(I,sigma); I ideal, sigma intvec
451ASSUME:  size(sigma) = nvars(basering) =: n
452RETURN:  list of eigenvalues of generators of I under permutation sigma
453NOTE:    sigma is a permutation of the variables of the basering, i.e.
454         sigma: var(i) ----> var(sigma[i]), 1 <= i <= n.
455EXAMPLE: example eigenvalues; shows an example
456"
457{
458   int i, j;
459
460   def A = basering;
461   int n = nvars(A);
462
463   poly ev;
464   list EV;
465   poly s, help_var;
466   matrix C1, C2;
467
468   ideal perm;
469   for(i = 1; i <= n; i++)
470   {
471      perm[size(perm)+1] = var(sigma[i]);
472   }
473
474   map f = A, perm;
475
476   for(i = 1; i <= size(I); i++)
477   {
478//-------------- s is the image of I[i] under permutation sigma ----------------
479      s = I[i];
480      s = f(s);
481
482      for(j = 1; j <= n; j++)
483      {
484         C1 = coeffs(I[i], var(j));
485         C2 = coeffs(s, var(j));
486         if(nrows(C1) > 1)
487         {
488            ev = C2[nrows(C2),1]/C1[nrows(C1),1];
489
490//------ Furthermore check that I[i] is eigenvector of permutation sigma. ------
491            if(s == ev*I[i])
492            {
493               break;
494            }
495            else
496            {
497               ERROR("I["+string(i)+"] is no eigenvector "
498                      +"of permutation sigma!!");
499            }
500         }
501      }
502
503      EV[size(EV)+1] = ev;
504   }
505
506   return(EV);
507}
508example
509{ "EXAMPLE:";  echo = 2;
510   ring R = 11, x(1..5), dp;
511   poly p1 = x(1)+x(2)+x(3)+x(4)+x(5);
512   poly p2 = x(1)+4*x(2)+5*x(3)-2*x(4)+3*x(5);
513   poly p3 = x(1)+5*x(2)+3*x(3)+4*x(4)-2*x(5);
514   poly p4 = x(1)-2*x(2)+4*x(3)+3*x(4)+5*x(5);
515   poly p5 = x(1)+3*x(2)-2*x(3)+5*x(4)+4*x(5);
516   ideal I = p1,p2,p3,p4,p5;
517   intvec tau = 2,3,4,5,1;
518   eigenvalues(I,tau);
519}
520
521////////////////////////////////////////////////////////////////////////////////
522
523proc symmStd(ideal I, intvec sigma, list #)
524"USAGE:  symmStd(I,sigma,#); I ideal, sigma intvec
525ASSUME:  size(sigma) = nvars(basering) =: n, basering has an order(sigma)-th
526         primitive root of unity a (if char(basering) > 0) and sigma(I) = I
527RETURN:  ideal, a standard basis of I
528NOTE:    Assuming that the ideal I is invariant under the variable permutation
529         sigma and the basering has an order(sigma)-th primitive root of unity
530         the procedure uses linear transformation of variables in order to
531         improve standard basis computation.
532         If char(basering) = 0 all computations are done in the polynomial ring
533         over the smallest field extension that has an order(sigma)-th primitive
534         root of unity.
535EXAMPLE: example symmStd; shows an example
536"
537{
538   if((nvars(basering) != size(sigma)) || (!isSymmetric(I,sigma)))
539   {
540      ERROR("The input is no permutation of the ring-variables!!");
541   }
542
543   option(redSB);
544
545   def R = basering;
546   int CHAR = char(R);
547   int n = nvars(R);
548
549   int t;
550
551//-------- (1) Compute the order of variable permutation sigma. ----------------
552   int ORD = order(sigma);
553   if((((CHAR - 1) mod ORD) != 0) && (CHAR > 0))
554   {
555      ERROR("Basering of characteristic "+string(CHAR)+" has no "
556             +string(ORD)+"-th primitive root of unity!!");
557   }
558
559//-------- (2) Compute the order(sigma)-th primitive root of unity -------------
560//--------     in basering or move to ring extension.              -------------
561   if((CHAR == 0) && (ORD > 2))
562   {
563      def S = changechar("0,p");
564      setring S;
565      minpoly = rootofUnity(ORD);
566      ideal I = imap(R, I);
567      number a = p;
568   }
569   else
570   {
571      int xi = primRoot(CHAR, ORD);
572      number a = xi;
573   }
574
575//--------- (3) Define the linear transformation of variables with -------------
576//---------     respect to sigma.                                  -------------
577   list L = genTransId(sigma,ORD,a);
578   ideal TR = L[1];
579   ideal RETR = L[2];
580
581//--------- (4) Compute the eigenvalues of the "new" variables of --------------
582//---------     sigma after transformation.                       --------------
583   list EV = eigenvalues(RETR, sigma);
584
585//--------- (5) Transformation of the input-ideal is done here. ----------------
586   map f = basering, TR;
587   t = timer;
588   ideal I_trans = f(I);
589   if(printlevel >= 11) { "Transformation: "+string(timer - t)+" seconds"; }
590
591//--------- (6) Compute a standard basis of the transformed ideal. -------------
592   t = timer;
593   if(size(#) > 0) { ideal sI_trans = std(I_trans, #[1]); }
594   else { ideal sI_trans = std(I_trans); }
595   if(printlevel >= 11) { "1st Groebner basis: "+string(timer - t)+" seconds"; }
596
597//--------- (7) Retransformation is done here. ---------------------------------
598   map g = basering, RETR;
599   t = timer;
600   ideal I_retrans = g(sI_trans);
601   if(printlevel >= 11) { "Reverse Transformation: "+string(timer - t)
602                           +" seconds"; }
603
604//--------- (8) Compute a standard basis of the retransformaed ideal -----------
605//---------     which is then a standard basis of the input-ideal.   -----------
606   t = timer;
607   ideal sI_retrans = std(I_retrans);
608   if(printlevel >= 11) { "2nd Groebner basis: "+string(timer - t)+" seconds"; }
609
610   if((CHAR == 0) && (ORD > 2))
611   {
612      setring R;
613      ideal sI_retrans = fetch(S, sI_retrans);
614      return(sI_retrans);
615   }
616   else
617   {
618      return(sI_retrans);
619   }
620}
621example
622{ "EXAMPLE:";  echo = 2;
623   ring R = 0, x(1..4), dp;
624   ideal I = cyclic(4);
625   I;
626   intvec pi = 2,3,4,1;
627   ideal sI = symmStd(I,pi);
628   sI;
629
630   ring S = 31, (x,y,z), dp;
631   ideal J;
632   J[1] = xy-y2+xz;
633   J[2] = xy+yz-z2;
634   J[3] = -x2+xz+yz;
635   intvec tau = 3,1,2;
636   ideal sJ = symmStd(J,tau);
637   sJ;
638}
639
640////////////////////////////////////////////////////////////////////////////////
641
642proc divPrimeTest(def II, bigint p, int k)
643{
644   if((p - 1) mod k != 0) { return(0); }
645
646   if(typeof(II) == "string")
647   {
648      execute("ideal I = "+II+";");
649   }
650   else
651   {
652      ideal I = II;
653   }
654
655   int i,j;
656   poly f;
657   number cnt;
658   for(i = 1; i <= size(I); i++)
659   {
660      f = cleardenom(I[i]);
661      if(f == 0) { return(0); }
662      cnt = leadcoef(I[i])/leadcoef(f);
663      if((numerator(cnt) mod p) == 0) { return(0); }
664      if((denominator(cnt) mod p) == 0) { return(0); }
665      for(j = size(f); j > 0; j--)
666      {
667         if((leadcoef(f[j]) mod p) == 0) { return(0); }
668      }
669   }
670   return(1);
671}
672
673////////////////////////////////////////////////////////////////////////////////
674
675proc divPrimeList(int k, ideal I, int n, list #)
676{
677// the intvec of n greatest primes p <= 2147483647 (resp. n greatest primes
678// < L[size(L)] union with L) such that each (p-1) is divisible by k, and none
679// of these primes divides any coefficient occuring in I
680// --> similar to procedure primeList in modstd.lib
681
682   intvec L;
683   int i,p;
684   int ncores = 1;
685
686//-----------------  Initialize optional parameter ncores  ---------------------
687   if(size(#) > 0)
688   {
689      if(size(#) == 1)
690      {
691         if(typeof(#[1]) == "int")
692         {
693            ncores = #[1];
694            # = list();
695         }
696      }
697      else
698      {
699         ncores = #[2];
700      }
701   }
702
703   if(size(#) == 0)
704   {
705      p = 2147483647;
706      while(!divPrimeTest(I,p,k))
707      {
708         p = prime(p-1);
709         if(p == 2) { ERROR("No more primes."); }
710      }
711      L[1] = p;
712   }
713   else
714   {
715      L = #[1];
716      p = prime(L[size(L)]-1);
717      while(!divPrimeTest(I,p,k))
718      {
719         p = prime(p-1);
720         if(p == 2) { ERROR("No more primes."); }
721      }
722      L[size(L)+1] = p;
723   }
724   if(p == 2) { ERROR("No more primes."); }
725   if(ncores == 1)
726   {
727      for(i = 2; i <= n; i++)
728      {
729         p = prime(p-1);
730         while(!divPrimeTest(I,p,k))
731         {
732            p = prime(p-1);
733            if(p == 2) { ERROR("no more primes"); }
734         }
735         L[size(L)+1] = p;
736      }
737   }
738   else
739   {
740      int neededSize = size(L)+n-1;;
741      list parallelResults;
742      list arguments;
743      int neededPrimes = neededSize-size(L);
744      while(neededPrimes > 0)
745      {
746         arguments = list();
747         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
748            *ncores; i > 0; i--)
749         {
750            p = prime(p-1);
751            if(p == 2) { ERROR("no more primes"); }
752            arguments[i] = list("I", p, k);
753         }
754         parallelResults = parallelWaitAll("divPrimeTest", arguments,
755            list(list(list(ncores))));
756         for(i = size(arguments); i > 0; i--)
757         {
758            if(parallelResults[i])
759            {
760               L[size(L)+1] = arguments[i][2];
761            }
762         }
763         neededPrimes = neededSize-size(L);
764      }
765      if(size(L) > neededSize)
766      {
767         L = L[1..neededSize];
768      }
769   }
770   return(L);
771}
772example
773{ "EXAMPLE:"; echo = 2;
774   ring r = 0,(x,y,z),dp;
775   ideal I = 2147483647x+y, z-181;
776   intvec L = divPrimeList(4,I,10,10);
777   size(L);
778   L[1];
779   L[size(L)];
780   L = divPrimeList(4,I,5,L,5);
781   size(L);
782   L[size(L)];
783}
784
785////////////////////////////////////////////////////////////////////////////////
786
787proc spTestSB(ideal I, ideal J, list L, intvec sigma, int variant, list #)
788"USAGE:  spTestSB(I,J,L,sigma,variant,#); I,J ideals, L intvec of primes,
789                                          sigma intvec, variant integer
790RETURN:  1 (resp. 0) if for a randomly chosen prime p, that is not in L and
791         divisible by the order of sigma, J mod p is (resp. is not) a standard
792         basis of I mod p
793EXAMPLE: example spTestSB; shows an example
794"
795{
796   int i,j,k,p;
797   int ORD = order(sigma);
798   def R = basering;
799   list r = ringlist(R);
800
801   while(!j)
802   {
803      j = 1;
804      while(((p - 1) mod ORD) != 0)
805      {
806         p = prime(random(1000000000,2134567879));
807         if(p == 2){ ERROR("no more primes"); }
808      }
809      for(i = 1; i <= size(L); i++)
810      {
811         if(p == L[i]){ j = 0; break }
812      }
813      if(j)
814      {
815         for(i = 1; i <= ncols(I); i++)
816         {
817            for(k = 2; k <= size(I[i]); k++)
818            {
819               if((denominator(leadcoef(I[i][k])) mod p) == 0){ j = 0; break; }
820            }
821            if(!j){ break; }
822         }
823      }
824      if(j)
825      {
826         if(!primeTest(I,p)) { j = 0; }
827      }
828   }
829   r[1] = p;
830   def @R = ring(r);
831   setring @R;
832   ideal I = imap(R,I);
833   ideal J = imap(R,J);
834   attrib(J,"isSB",1);
835
836   int t = timer;
837   j = 1;
838   if(isIncluded(I,J) == 0){ j = 0; }
839
840   if(printlevel >= 11)
841   {
842      "isIncluded(I,J) takes "+string(timer - t)+" seconds";
843      "j = "+string(j);
844   }
845
846   t = timer;
847   if(j)
848   {
849      if(size(#) > 0)
850      {
851         ideal K = smpStd(I,sigma,p,variant,#[1])[1];
852      }
853      else
854      {
855         ideal K = smpStd(I,sigma,p,variant)[1];
856      }
857      t = timer;
858      if(isIncluded(J,K) == 0){ j = 0; }
859
860      if(printlevel >= 11)
861      {
862         "isIncluded(J,K) takes "+string(timer - t)+" seconds";
863         "j = "+string(j);
864      }
865   }
866   setring R;
867   return(j);
868}
869example
870{ "EXAMPLE:"; echo = 2;
871   intvec L = 2,3,5;
872   ring r = 0,(x,y,z),dp;
873   ideal I = x+1,y+1;
874   intvec sigma = 2,1,3;
875   ideal J = x+1,y;
876   spTestSB(I,J,L,sigma,2);
877   spTestSB(J,I,L,sigma,2);
878}
879
880////////////////////////////////////////////////////////////////////////////////
881
882proc smpStd(ideal I, intvec sigma, int p, int variant, list #)
883"USAGE:  smpStd(I,sigma,p,#); I ideal, sigma intvec, p integer, variant integer
884ASSUME:  If size(#) > 0, then #[1] is an intvec describing the Hilbert series.
885RETURN:  ideal - a standard basis of I mod p, integer - p
886NOTE:    The procedure computes a standard basis of the ideal I modulo p and
887         fetches the result to the basering. If size(#) > 0 the Hilbert driven
888         standard basis computation symmStd(.,.,#[1]) is used in symmStd.
889         The standard basis computation modulo p does also vary depending on the
890         integer variant, namely
891@*       - variant = 1: symmStd(.,.,#[1]) resp. symmStd,
892@*       - variant = 2: symmStd,
893@*       - variant = 3: homog. - symmStd(.,.,#[1]) resp. symmStd - dehomog.,
894@*       - variant = 4: fglm.
895EXAMPLE: example smpStd; shows an example
896"
897{
898   def R0 = basering;
899   list rl = ringlist(R0);
900   rl[1] = p;
901   def @r = ring(rl);
902   setring @r;
903   ideal i = fetch(R0,I);
904
905   option(redSB);
906
907   if(variant == 1)
908   {
909      if(size(#) > 0)
910      {
911         i = symmStd(i, sigma, #[1]);
912      }
913      else
914      {
915         i = symmStd(i, sigma);
916      }
917   }
918
919   if(variant == 2)
920   {
921      i = symmStd(i, sigma);
922   }
923
924   if(variant == 3)
925   {
926      list rl = ringlist(@r);
927      int nvar@r = nvars(@r);
928
929      int k;
930      intvec w;
931      for(k = 1; k <= nvar@r; k++)
932      {
933         w[k] = deg(var(k));
934      }
935      w[nvar@r + 1] = 1;
936
937      rl[2][nvar@r + 1] = "homvar";
938      rl[3][2][2] = w;
939
940      def HomR = ring(rl);
941      setring HomR;
942      ideal i = imap(@r, i);
943      i = homog(i, homvar);
944      intvec tau = sigma, size(sigma)+1;
945
946      if(size(#) > 0)
947      {
948         if(w == 1)
949         {
950            i = symmStd(i, tau, #[1]);
951         }
952         else
953         {
954            i = symmStd(i, tau, #[1], w);
955         }
956      }
957      else
958      {
959         i = symmStd(i, tau);
960      }
961
962      i = subst(i, homvar, 1);
963      i = simplify(i, 34);
964
965      setring @r;
966      i = imap(HomR, i);
967      i = interred(i);
968      kill HomR;
969   }
970
971   if(variant == 4)
972   {
973      def R1 = changeord("dp");
974      setring R1;
975      ideal i = fetch(@r,i);
976      i = symmStd(i, sigma);
977      setring @r;
978      i = fglm(R1,i);
979   }
980
981   setring R0;
982   return(list(fetch(@r,i),p));
983}
984example
985{ "EXAMPLE:"; echo = 2;
986   ring r1 = 0, x(1..4), dp;
987   ideal I = cyclic(4);
988   intvec sigma = 2,3,4,1;
989   int p = 181;
990   list P = smpStd(I,sigma,p,2);
991   P;
992
993   ring r2 = 0, x(1..5), lp;
994   ideal I = cyclic(5);
995   intvec tau = 2,3,4,5,1;
996   int q = 31981;
997   list Q = smpStd(I,tau,q,4);
998   Q;
999}
1000
1001////////////////////////////////////////////////////////////////////////////////
1002
1003proc syModStd(ideal I, intvec sigma, list #)
1004"USAGE:  syModStd(I,sigma); I ideal, sigma intvec
1005ASSUME:  size(sigma) = nvars(basering) and sigma(I) = I. If size(#) > 0, then
1006         # contains either 1, 2 or 4 integers such that
1007@*       - #[1] is the number of available processors for the computation,
1008@*       - #[2] is an optional parameter for the exactness of the computation,
1009                if #[2] = 1, the procedure computes a standard basis for sure,
1010@*       - #[3] is the number of primes until the first lifting,
1011@*       - #[4] is the constant number of primes between two liftings until
1012           the computation stops.
1013RETURN:  ideal, a standard basis of I if no warning appears;
1014NOTE:    The procedure computes a standard basis of the ideal I (over the
1015         rational numbers) by using modular methods and the fact that I is
1016         invariant under the variable permutation sigma.
1017         By default the procedure computes a standard basis of I for sure, but
1018         if the optional parameter #[2] = 0, it computes a standard basis of I
1019         with high probability.
1020         The procedure distinguishes between different variants for the standard
1021         basis computation in positive characteristic depending on the ordering
1022         of the basering, the parameter #[2] and if the ideal I is homogeneous.
1023@*       - variant = 1, if I is homogeneous,
1024@*       - variant = 2, if I is not homogeneous, 1-block-ordering,
1025@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
1026                        > 1 block),
1027@*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
1028EXAMPLE: example syModStd; shows an example
1029"
1030{
1031   if((nvars(basering) != size(sigma)) || (!isSymmetric(I,sigma)))
1032   {
1033      ERROR("The input is no permutation of the ring-variables!!");
1034   }
1035
1036   int TT = timer;
1037   int RT = rtimer;
1038
1039   def R0 = basering;
1040   list rl = ringlist(R0);
1041   if((npars(R0) > 0) || (rl[1] > 0))
1042   {
1043      ERROR("Characteristic of basering should be zero, basering should
1044             have no parameters.");
1045   }
1046
1047   int index = 1;
1048   int i,k,c;
1049   int j = 1;
1050   int pTest, sizeTest;
1051   int en = 2134567879;
1052   int an = 1000000000;
1053   bigint N;
1054   int ORD = order(sigma);
1055
1056//--------------------  Initialize optional parameters  ------------------------
1057   if(size(#) > 0)
1058   {
1059      if(size(#) == 1)
1060      {
1061         int n1 = #[1];
1062         int exactness = 1;
1063         if(n1 >= 10)
1064         {
1065            int n2 = n1 + 1;
1066            int n3 = n1;
1067         }
1068         else
1069         {
1070            int n2 = 10;
1071            int n3 = 10;
1072         }
1073      }
1074      if(size(#) == 2)
1075      {
1076         int n1 = #[1];
1077         int exactness = #[2];
1078         if(n1 >= 10)
1079         {
1080            int n2 = n1 + 1;
1081            int n3 = n1;
1082         }
1083         else
1084         {
1085            int n2 = 10;
1086            int n3 = 10;
1087         }
1088      }
1089      if(size(#) == 4)
1090      {
1091         int n1 = #[1];
1092         int exactness = #[2];
1093         if(n1 >= #[3])
1094         {
1095            int n2 = n1 + 1;
1096         }
1097         else
1098         {
1099            int n2 = #[3];
1100         }
1101         if(n1 >= #[4])
1102         {
1103            int n3 = n1;
1104         }
1105         else
1106         {
1107            int n3 = #[4];
1108         }
1109      }
1110   }
1111   else
1112   {
1113      int n1 = 1;
1114      int exactness = 1;
1115      int n2 = 10;
1116      int n3 = 10;
1117   }
1118
1119   if(printlevel >= 10)
1120   {
1121      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
1122       +", exactness = "+string(exactness);
1123   }
1124
1125//-------------------------- Save current options ------------------------------
1126   intvec opt = option(get);
1127
1128   option(redSB);
1129
1130//--------------------  Initialize the list of primes  -------------------------
1131   int tt = timer;
1132   int rt = rtimer;
1133   intvec L = divPrimeList(ORD,I,n2,n1);
1134   if(printlevel >= 10)
1135   {
1136      "CPU-time for divPrimeList: "+string(timer-tt)+" seconds.";
1137      "Real-time for divPrimeList: "+string(rtimer-rt)+" seconds.";
1138   }
1139
1140//---------------------  Decide which variant to take  -------------------------
1141   int variant;
1142   int h = homog(I);
1143
1144   tt = timer;
1145   rt = rtimer;
1146
1147   if(!mixedTest())
1148   {
1149      if(h)
1150      {
1151         variant = 1;
1152         if(printlevel >= 10) { "variant = 1"; }
1153
1154         rl[1] = L[5];
1155         def @r = ring(rl);
1156         setring @r;
1157         def @s = changeord("dp");
1158         setring @s;
1159         ideal I = std(fetch(R0,I));
1160         intvec hi = hilb(I,1);
1161         setring R0;
1162         kill @r,@s;
1163      }
1164      else
1165      {
1166         string ordstr_R0 = ordstr(R0);
1167         int neg = 1 - attrib(R0,"global");
1168
1169         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
1170         {
1171            variant = 2;
1172            if(printlevel >= 10) { "variant = 2"; }
1173         }
1174         else
1175         {
1176            string order;
1177            if(system("nblocks") <= 2)
1178            {
1179               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
1180                                       + find(ordstr_R0, "rp") <= 0)
1181               {
1182                  order = "simple";
1183               }
1184            }
1185
1186            if((order == "simple") || (size(rl) > 4))
1187            {
1188               variant = 2;
1189               if(printlevel >= 10) { "variant = 2"; }
1190            }
1191            else
1192            {
1193               rl[1] = L[5];
1194               def @r = ring(rl);
1195               setring @r;
1196
1197               def @s = changeord("dp");
1198               setring @s;
1199               ideal I = std(fetch(R0,I));
1200               if(dim(I) == 0)
1201               {
1202                  variant = 4;
1203                  if(printlevel >= 10) { "variant = 4"; }
1204               }
1205               else
1206               {
1207                  variant = 3;
1208                  if(printlevel >= 10) { "variant = 3"; }
1209
1210                  int nvar@r = nvars(@r);
1211                  intvec w;
1212                  for(i = 1; i <= nvar@r; i++)
1213                  {
1214                     w[i] = deg(var(i));
1215                  }
1216                  w[nvar@r + 1] = 1;
1217
1218                  list hiRi = hilbRing(fetch(R0,I),w);
1219                  intvec W = hiRi[2];
1220                  @s = hiRi[1];
1221                  setring @s;
1222                  intvec tau = sigma, size(sigma)+1;
1223
1224                  Id(1) = symmStd(Id(1),tau);
1225                  intvec hi = hilb(Id(1), 1, W);
1226               }
1227
1228               setring R0;
1229               kill @r,@s;
1230            }
1231         }
1232      }
1233   }
1234   else
1235   {
1236      if(exactness == 1) { return(groebner(I)); }
1237      if(h)
1238      {
1239         variant = 1;
1240         if(printlevel >= 10) { "variant = 1"; }
1241         rl[1] = L[5];
1242         def @r = ring(rl);
1243         setring @r;
1244         def @s = changeord("dp");
1245         setring @s;
1246         ideal I = std(fetch(R0,I));
1247         intvec hi = hilb(I,1);
1248         setring R0;
1249         kill @r,@s;
1250      }
1251      else
1252      {
1253         string ordstr_R0 = ordstr(R0);
1254         int neg = 1 - attrib(R0,"global");
1255
1256         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
1257         {
1258            variant = 2;
1259            if(printlevel >= 10) { "variant = 2"; }
1260         }
1261         else
1262         {
1263            string order;
1264            if(system("nblocks") <= 2)
1265            {
1266               if(find(ordstr_R0, "M") + find(ordstr_R0, "lp")
1267                                       + find(ordstr_R0, "rp") <= 0)
1268               {
1269                  order = "simple";
1270               }
1271            }
1272
1273            if((order == "simple") || (size(rl) > 4))
1274            {
1275               variant = 2;
1276               if(printlevel >= 10) { "variant = 2"; }
1277            }
1278            else
1279            {
1280               variant = 3;
1281               if(printlevel >= 10) { "variant = 3"; }
1282
1283               rl[1] = L[5];
1284               def @r = ring(rl);
1285               setring @r;
1286               int nvar@r = nvars(@r);
1287               intvec w;
1288               for(i = 1; i <= nvar@r; i++)
1289               {
1290                  w[i] = deg(var(i));
1291               }
1292               w[nvar@r + 1] = 1;
1293
1294               list hiRi = hilbRing(fetch(R0,I),w);
1295               intvec W = hiRi[2];
1296               def @s = hiRi[1];
1297               setring @s;
1298               intvec tau = sigma, size(sigma)+1;
1299
1300               Id(1) = symmStd(Id(1),tau);
1301               intvec hi = hilb(Id(1), 1, W);
1302
1303               setring R0;
1304               kill @r,@s;
1305            }
1306         }
1307      }
1308   }
1309
1310   list P,T1,T2,T3,LL;
1311
1312   ideal J,K,H;
1313
1314//-----  If there is more than one processor available, we parallelize the  ----
1315//-----  main standard basis computations in positive characteristic        ----
1316
1317   if(n1 > 1)
1318   {
1319      ideal I_for_fork = I;
1320      export(I_for_fork);           // I available for each link
1321
1322//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
1323//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
1324
1325      for(i = 1; i <= n1; i++)
1326      {
1327         //link l(i) = "MPtcp:fork";
1328         link l(i) = "ssi:fork";
1329         open(l(i));
1330         if((variant == 1) || (variant == 3))
1331         {
1332            write(l(i), quote(smpStd(I_for_fork, eval(sigma), eval(L[i + 1]),
1333                                                 eval(variant), eval(hi))));
1334         }
1335         if((variant == 2) || (variant == 4))
1336         {
1337            write(l(i), quote(smpStd(I_for_fork, eval(sigma), eval(L[i + 1]),
1338                                                 eval(variant))));
1339         }
1340      }
1341
1342      int t = timer;
1343      if((variant == 1) || (variant == 3))
1344      {
1345         P = smpStd(I_for_fork, sigma, L[1], variant, hi);
1346      }
1347      if((variant == 2) || (variant == 4))
1348      {
1349         P = smpStd(I_for_fork, sigma, L[1], variant);
1350      }
1351      t = timer - t;
1352      if(t > 60) { t = 60; }
1353      int i_sleep = system("sh", "sleep "+string(t));
1354      T1[1] = P[1];
1355      T2[1] = bigint(P[2]);
1356      index++;
1357
1358      j = j + n1 + 1;
1359   }
1360
1361//--------------  Main standard basis computations in positive  ----------------
1362//----------------------  characteristic start here  ---------------------------
1363
1364   list arguments_farey, results_farey;
1365
1366   while(1)
1367   {
1368      tt = timer; rt = rtimer;
1369
1370      if(printlevel >= 10) { "size(L) = "+string(size(L)); }
1371
1372      if(n1 > 1)
1373      {
1374         while(j <= size(L) + 1)
1375         {
1376            for(i = 1; i <= n1; i++)
1377            {
1378               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
1379               if(status(l(i), "read", "ready"))
1380               {
1381                  //--- read the result from l(i) ---
1382                  P = read(l(i));
1383                  T1[index] = P[1];
1384                  T2[index] = bigint(P[2]);
1385                  index++;
1386
1387                  if(j <= size(L))
1388                  {
1389                     if((variant == 1) || (variant == 3))
1390                     {
1391                        write(l(i), quote(smpStd(I_for_fork, eval(sigma),
1392                                        eval(L[j]), eval(variant), eval(hi))));
1393                        j++;
1394                     }
1395                     if((variant == 2) || (variant == 4))
1396                     {
1397                        write(l(i), quote(smpStd(I_for_fork, eval(sigma),
1398                                        eval(L[j]), eval(variant))));
1399                        j++;
1400                     }
1401                  }
1402                  else
1403                  {
1404                     k++;
1405                     close(l(i));
1406                  }
1407               }
1408            }
1409            //--- k describes the number of closed links ---
1410            if(k == n1)
1411            {
1412               j++;
1413            }
1414            i_sleep = system("sh", "sleep "+string(t));
1415         }
1416      }
1417      else
1418      {
1419         while(j <= size(L))
1420         {
1421            if((variant == 1) || (variant == 3))
1422            {
1423               P = smpStd(I, sigma, L[j], variant, hi);
1424            }
1425            if((variant == 2) || (variant == 4))
1426            {
1427               P = smpStd(I, sigma, L[j], variant);
1428            }
1429
1430            T1[index] = P[1];
1431            T2[index] = bigint(P[2]);
1432            index++;
1433            j++;
1434         }
1435      }
1436
1437      if(printlevel >= 10)
1438      {
1439         "CPU-time for computing list is "+string(timer - tt)+" seconds.";
1440         "Real-time for computing list is "+string(rtimer - rt)+" seconds.";
1441      }
1442
1443//------------------------  Delete unlucky primes  -----------------------------
1444//-------------  unlucky if and only if the leading ideal is wrong  ------------
1445
1446      LL = deleteUnluckyPrimes(T1,T2,h);
1447      T1 = LL[1];
1448      T2 = LL[2];
1449
1450//-------------------  Now all leading ideals are the same  --------------------
1451//-------------------  Lift results to basering via farey  ---------------------
1452
1453      tt = timer; rt = rtimer;
1454      N = T2[1];
1455      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
1456      H = chinrem(T1,T2);
1457      if(n1 == 1)
1458      {
1459         J = farey(H,N);
1460      }
1461      else
1462      {
1463         for(i = ncols(H); i > 0; i--)
1464         {
1465            arguments_farey[i] = list(ideal(H[i]), N);
1466         }
1467         results_farey = parallelWaitAll("farey", arguments_farey,
1468                                         list(list(list(n1))));
1469         for(i = ncols(H); i > 0; i--)
1470         {
1471            J[i] = results_farey[i][1];
1472         }
1473      }
1474      if(printlevel >= 10)
1475      {
1476         "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
1477         "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
1478      }
1479
1480//----------------  Test if we already have a standard basis of I --------------
1481
1482      tt = timer; rt = rtimer;
1483      if((variant == 1) || (variant == 3))
1484      {
1485         pTest = spTestSB(I,J,L,sigma,variant,hi);
1486      }
1487      if((variant == 2) || (variant == 4))
1488      {
1489         pTest = spTestSB(I,J,L,sigma,variant);
1490      }
1491
1492      if(printlevel >= 10)
1493      {
1494         "CPU-time for pTest is "+string(timer - tt)+" seconds.";
1495         "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
1496      }
1497
1498      if(pTest)
1499      {
1500         if(printlevel >= 10)
1501         {
1502            "CPU-time for computation without final tests is "
1503            +string(timer - TT)+" seconds.";
1504            "Real-time for computation without final tests is "
1505            +string(rtimer - RT)+" seconds.";
1506         }
1507
1508         attrib(J,"isSB",1);
1509
1510         if(exactness == 0)
1511         {
1512            option(set, opt);
1513            if(n1 > 1) { kill I_for_fork; }
1514            return(J);
1515         }
1516
1517         if(exactness == 1)
1518         {
1519            tt = timer; rt = rtimer;
1520            sizeTest = 1 - isIncluded(I,J,n1);
1521
1522            if(printlevel >= 10)
1523            {
1524               "CPU-time for checking if I subset <G> is "
1525               +string(timer - tt)+" seconds.";
1526               "Real-time for checking if I subset <G> is "
1527               +string(rtimer - rt)+" seconds.";
1528            }
1529
1530            if(sizeTest == 0)
1531            {
1532               tt = timer; rt = rtimer;
1533               K = std(J);
1534
1535               if(printlevel >= 10)
1536               {
1537                  "CPU-time for last std-computation is "
1538                  +string(timer - tt)+" seconds.";
1539                  "Real-time for last std-computation is "
1540                  +string(rtimer - rt)+" seconds.";
1541               }
1542
1543               if(size(reduce(K,J)) == 0)
1544               {
1545                  option(set, opt);
1546                  if(n1 > 1) { kill I_for_fork; }
1547                  return(J);
1548               }
1549            }
1550         }
1551      }
1552
1553//--------------  We do not already have a standard basis of I  ----------------
1554//-----------  Therefore do the main computation for more primes  --------------
1555
1556      T1 = H;
1557      T2 = N;
1558      index = 2;
1559
1560      j = size(L) + 1;
1561      tt = timer; rt = rtimer;
1562      L = divPrimeList(ORD,I,n3,L,n1);
1563      if(printlevel >= 10)
1564      {
1565         "CPU-time for divPrimeList: "+string(timer-tt)+" seconds.";
1566         "Real-time for divPrimeList: "+string(rtimer-rt)+" seconds.";
1567      }
1568
1569      if(n1 > 1)
1570      {
1571         for(i = 1; i <= n1; i++)
1572         {
1573            open(l(i));
1574            if((variant == 1) || (variant == 3))
1575            {
1576               write(l(i), quote(smpStd(I_for_fork, eval(sigma), eval(L[j+i-1]),
1577                                                    eval(variant), eval(hi))));
1578            }
1579            if((variant == 2) || (variant == 4))
1580            {
1581               write(l(i), quote(smpStd(I_for_fork, eval(sigma), eval(L[j+i-1]),
1582                                                    eval(variant))));
1583            }
1584         }
1585         j = j + n1;
1586         k = 0;
1587      }
1588   }
1589}
1590example
1591{ "EXAMPLE:"; echo = 2;
1592   ring R1 = 0, (x,y,z), dp;
1593   ideal I;
1594   I[1] = -2xyz4+xz5+xz;
1595   I[2] = -2xyz4+yz5+yz;
1596   intvec sigma = 2,1,3;
1597   ideal sI = syModStd(I,sigma);
1598   sI;
1599
1600   ring R2 = 0, x(1..4), dp;
1601   ideal I = cyclic(4);
1602   I;
1603   intvec pi = 2,3,4,1;
1604   ideal sJ1 = syModStd(I,pi,1);
1605   ideal sJ2 = syModStd(I,pi,1,0);
1606   size(reduce(sJ1,sJ2));
1607   size(reduce(sJ2,sJ1));
1608}
Note: See TracBrowser for help on using the repository browser.