source: git/Singular/LIB/modstd.lib @ 6368d1

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