source:git/Singular/LIB/symodstd.lib@4fb2ef5

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