source: git/Singular/LIB/symodstd.lib @ 380a17b

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