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

spielwiese
Last change on this file since a3c2eb6 was a3c2eb6, checked in by Hans Schoenemann <hannes@…>, 3 years ago
fix: typo
  • Property mode set to 100644
File size: 23.3 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version modstd.lib 4.2.0.0 Feb_2020 "; // $Id$
3category="Commutative Algebra";
4info="
5LIBRARY:  modstd.lib      Groebner bases of ideals/modules using modular methods
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  A library for computing Groebner bases of ideals/modules in the polynomial ring
15  over the rational numbers using modular methods.
16
17REFERENCES:
18  E. A. Arnold: Modular algorithms for computing Groebner bases.
19  J. Symb. Comp. 35, 403-419 (2003).
20
21  N. Idrees, G. Pfister, S. Steidel: Parallelization of Modular Algorithms.
22  J. Symb. Comp. 46, 672-684 (2011).
23
24PROCEDURES:
25  modStd(I);    standard basis of I using modular methods
26  modSyz(I);    syzygy module of I  using modular methods
27  modIntersect(I,J); intersection of I and J using modular methods
28";
29
30LIB "polylib.lib";
31LIB "modular.lib";
32
33proc modStd(def I, list #)
34"USAGE:   modStd(I[, exactness]); I ideal/module, exactness int
35RETURN:   a standard basis of I
36NOTE:     The procedure computes a standard basis of I (over the rational
37          numbers) by using modular methods.
38       @* An optional parameter 'exactness' can be provided.
39          If exactness = 1(default), the procedure computes a standard basis
40          of I for sure; if exactness = 0, it computes a standard basis of I
41          with high probability.
42SEE ALSO: modular
43EXAMPLE:  example modStd; shows an example"
44{
45    /* read optional parameter */
46    int exactness = 1;
47    if (size(#) > 0)
48    {
49        /* For compatibility, we only test size(#) > 4. This can be changed to
50         * size(#) > 1 in the future. */
51        if (size(#) > 4 || typeof(#[1]) != "int")
52        {
53            ERROR("wrong optional parameter");
54        }
55        exactness = #[1];
56    }
57
58    /* save options */
59    intvec opt = option(get);
60    option(redSB);
61
62    /* choose the right command */
63    string command = "groebner";
64    if (npars(basering) > 0) { command = "Modstd::groebner_norm"; }
65
66    /* call modular() */
67    if (exactness)
68    {
69      if(hasCommutativeVars(basering))
70      {
71        I = modular(command, list(I), primeTest_std,
72            deleteUnluckyPrimes_std, pTest_std, finalTest_std_comm);
73      }
74      else
75      {
76        I = modular(command, list(I), primeTest_std,
77            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
78      }
79    }
80    else
81    {
82        I = modular(command, list(I), primeTest_std,
83            deleteUnluckyPrimes_std, pTest_std);
84    }
85
86    /* return the result */
87    attrib(I, "isSB", 1);
88    option(set, opt);
89    return(I);
90}
91example
92{
93    "EXAMPLE:";
94    echo = 2;
95    ring R1 = 0, (x,y,z,t), dp;
96    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
97    ideal J = modStd(I);
98    J;
99    I = homog(I, t);
100    J = modStd(I);
101    J;
102
103    ring R2 = 0, (x,y,z), ds;
104    ideal I = jacob(x5+y6+z7+xyz);
105    ideal J = modStd(I, 0);
106    J;
107
108    ring R3 = 0, x(1..4), lp;
109    ideal I = cyclic(4);
110    ideal J1 = modStd(I, 1);   // default
111    ideal J2 = modStd(I, 0);
112    size(reduce(J1, J2));
113    size(reduce(J2, J1));
114}
115
116proc modSyz(def I)
117"USAGE:   modSyz(I); I ideal/module
118RETURN:   a generating set of syzygies of I
119NOTE:     The procedure computes a the syzygy module of I (over the rational
120          numbers) by using modular methods with high probability.
121          The property of being a syzygy is tested.
122SEE ALSO: modular
123EXAMPLE:  example modSyz; shows an example"
124{
125    /* save options */
126    intvec opt = option(get);
127    option(redSB);
128
129    /* choose the right command */
130    string command = "syz";
131
132    /* call modular() */
133    module M = modular(command, list(I), primeTest_std,
134         deleteUnluckyPrimes_std, pTest_syz);
135
136    /* return the result */
137    option(set, opt);
138    return(M);
139}
140example
141{
142    "EXAMPLE:"; echo = 2;
143    ring R1 = 0, (x,y,z,t), dp;
144    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
145    modSyz(I);
146    simplify(syz(I),1);
147}
148
149proc modIntersect(def I, def J)
150"USAGE:   modIntersect(I,J); I,J ideal/module
151RETURN:   a generating set of the intersection of I and J
152NOTE:     The procedure computes a the intersection of I and J
153          (over the rational numbers) by using modular methods
154          with high probability.
155          No additional tests are performed.
156SEE ALSO: modular
157EXAMPLE:  example modIntersect; shows an example"
158{
159    /* save options */
160    intvec opt = option(get);
161    option(redSB);
162
163    /* choose the right command */
164    string command = "intersect";
165
166    /* call modular() */
167    def M = modular(command, list(I,J), primeTest_std,
168         deleteUnluckyPrimes_std);
169
170    /* return the result */
171    option(set, opt);
172    return(M);
173}
174example
175{
176    "EXAMPLE:"; echo = 2;
177    ring R1 = 0, (x,y,z,t), dp;
178    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
179    ideal J = maxideal(2);
180    modIntersect(I,J);
181    simplify(intersect(I,J),1);
182}
183
184/* compute a normalized GB via groebner() */
185static proc groebner_norm(def I)
186{
187    I = simplify(groebner(I), 1);
188    attrib(I, "isSB", 1);
189    return(I);
190}
191
192/* test if the prime p is suitable for the input, i.e. it does not divide
193 * the numerator or denominator of any of the coefficients */
194static proc primeTest_std(int p, alias list args)
195{
196    /* erase zero generators */
197    def I = simplify(args[1], 2);
198
199    /* clear denominators and count the terms */
200    def J=I; // dummy assign, to get the type of I
201    ideal K;
202    int n = ncols(I);
203    intvec sizes;
204    number cnt;
205    int i;
206    for(i = n; i > 0; i--)
207    {
208        J[i] = cleardenom(I[i]);
209        cnt = leadcoef(J[i])/leadcoef(I[i]);
210        K[i] = numerator(cnt)*var(1)+denominator(cnt);
211    }
212    sizes = size(J[1..n]);
213
214    /* change to characteristic p */
215    def br = basering;
216    list lbr = ringlist(br);
217    if (typeof(lbr[1]) == "int") { lbr[1] = p; }
218    else { lbr[1][1] = p; }
219    def rp = ring(lbr);
220    setring(rp);
221    def Jp = fetch(br, J);
222    ideal Kp = fetch(br, K);
223
224    /* test if any coefficient is missing */
225    if (intvec(size(Kp[1..n])) != 2:n) { setring(br); return(0); }
226    if (intvec(size(Jp[1..n])) != sizes) { setring(br); return(0); }
227    setring(br);
228    return(1);
229}
230
231/* find entries in modresults which come from unlucky primes.
232 * For this, sort the entries into categories depending on their leading
233 * ideal and return the indices in all but the biggest category. */
234static proc deleteUnluckyPrimes_std(alias list modresults)
235{
236    int size_modresults = size(modresults);
237
238    /* sort results into categories.
239     * each category is represented by three entries:
240     * - the corresponding leading ideal
241     * - the number of elements
242     * - the indices of the elements
243     */
244    list cat;
245    int size_cat;
246    def L=modresults[1]; // dummy assign to get the type of L
247    int i;
248    int j;
249    for (i = 1; i <= size_modresults; i++)
250    {
251        L = lead(modresults[i]);
252        attrib(L, "isSB", 1);
253        for (j = 1; j <= size_cat; j++)
254        {
255            if (size(L) == size(cat[j][1])
256            && size(reduce(L, cat[j][1], 5)) == 0
257            && size(reduce(cat[j][1], L, 5)) == 0)
258            {
259                cat[j][2] = cat[j][2]+1;
260                cat[j][3][cat[j][2]] = i;
261                break;
262            }
263        }
264        if (j > size_cat)
265        {
266            size_cat++;
267            cat[size_cat] = list();
268            cat[size_cat][1] = L;
269            cat[size_cat][2] = 1;
270            cat[size_cat][3] = list(i);
271        }
272    }
273
274    /* find the biggest categories */
275    int cat_max = 1;
276    int max = cat[1][2];
277    for (i = 2; i <= size_cat; i++)
278    {
279        if (cat[i][2] > max)
280        {
281            cat_max = i;
282            max = cat[i][2];
283        }
284    }
285
286    /* return all other indices */
287    list unluckyIndices;
288    for (i = 1; i <= size_cat; i++)
289    {
290        if (i != cat_max) { unluckyIndices = unluckyIndices + cat[i][3]; }
291    }
292    return(unluckyIndices);
293}
294////////////////////////////////////////////////////////////////////////////////
295
296static proc cleardenomModule(def I)
297{
298  int t=ncols(I);
299  if(size(I)==0)
300  {
301      return(I);
302  }
303  else
304  {
305    for(int i=1;i<=t;i++)
306    {
307      I[i]=cleardenom(I[i]);
308    }
309  }
310  return(I);
311}
312
313static proc pTest_syz(string command, alias list args, alias def result, int p)
314{
315  module result_without_denom=cleardenomModule(result);
316  return(size(module(matrix(args[1])*matrix(result_without_denom)))==0);
317}
318
319/* test if 'command' applied to 'args' in characteristic p is the same as
320   'result' mapped to characteristic p */
321static proc pTest_std(string command, alias list args, alias def result,
322    int p)
323{
324    /* change to characteristic p */
325    def br = basering;
326    list lbr = ringlist(br);
327    if (typeof(lbr[1]) == "int") { lbr[1] = p; }
328    else { lbr[1][1] = p; }
329    def rp = ring(lbr);
330    setring(rp);
331    def Ip = fetch(br, args)[1];
332    def Gp = fetch(br, result);
333    attrib(Gp, "isSB", 1);
334
335    /* test if Ip is in Gp */
336    int i;
337    for (i = ncols(Ip); i > 0; i--)
338    {
339        if (reduce(Ip[i], Gp, 1) != 0)
340        {
341            setring(br);
342            return(0);
343        }
344    }
345
346    /* compute command(args) */
347    execute("Ip = "+command+"(Ip);");
348
349    /* test if Gp is in Ip */
350    for (i = ncols(Gp); i > 0; i--)
351    {
352        if (reduce(Gp[i], Ip, 1) != 0) { setring(br); return(0); }
353    }
354    setring(br);
355    return(1);
356}
357
358/* test if 'result' is a GB of the input ideal, commutative ring */
359static proc finalTest_std_comm(string command, alias list args, def result)
360{
361    /* test if args[1] is in result */
362    attrib(result, "isSB", 1);
363    int i;
364    for (i = ncols(args[1]); i > 0; i--)
365    {
366        if (reduce(args[1][i], result, 1) != 0) { return(0); }
367    }
368
369    /* test if result is a GB */
370    //def G = std(result);
371    //if (reduce_parallel(G, result)) { return(0); }
372    //return(1);
373    return(system("verifyGB",result));
374}
375
376/* test if 'result' is a GB of the input ideal, generic */
377static proc finalTest_std(string command, alias list args, def result)
378{
379    /* test if args[1] is in result */
380    attrib(result, "isSB", 1);
381    int i;
382    for (i = ncols(args[1]); i > 0; i--)
383    {
384        if (reduce(args[1][i], result, 1) != 0) { return(0); }
385    }
386
387    /* test if result is a GB */
388    def G = std(result);
389    if (reduce_parallel(G, result)) { return(0); }
390    return(1);
391}
392
393/* return 1, if I_reduce is _not_ in G_reduce,
394 *        0, otherwise
395 * (same as size(reduce(I_reduce, G_reduce))).
396 * Uses parallelization. */
397static proc reduce_parallel(def I_reduce, def G_reduce)
398{
399    exportto(Modstd, I_reduce);
400    exportto(Modstd, G_reduce);
401    int size_I = ncols(I_reduce);
402    int chunks = Modular::par_range(size_I);
403    intvec range;
404    int i;
405    for (i = chunks; i > 0; i--)
406    {
407        range = Modular::par_range(size_I, i);
408        task t(i) = "Modstd::reduce_task", list(range);
409    }
410    startTasks(t(1..chunks));
411    waitAllTasks(t(1..chunks));
412    int result = 0;
413    for (i = chunks; i > 0; i--)
414    {
415        if (getResult(t(i))) { result = 1; break; }
416    }
417    kill I_reduce;
418    kill G_reduce;
419    return(result);
420}
421
422/* compute a chunk of reductions for reduce_parallel */
423static proc reduce_task(intvec range)
424{
425    int result = 0;
426    int i;
427    for (i = range[1]; i <= range[2]; i++)
428    {
429        if (reduce(I_reduce[i], G_reduce, 1) != 0) { result = 1; break; }
430    }
431    return(result);
432}
433
434////////////////////////////////////////////////////////////////////////////////
435/*
436 * The following procedures are kept for backward compatibility with the old
437 * version of modstd.lib. As of now (May 2014), they are still needed in
438 * modnormal.lib, modwalk.lib, and symodstd.lib. They can be removed here as
439 * soon as they are not longer needed in these libraries.
440 */
441
442LIB "parallel.lib";
443
444static proc mod_init()
445{
446   newstruct("idealPrimeTest", "ideal Ideal");
447}
448
449static proc redFork(ideal I, ideal J, int n)
450{
451   attrib(J,"isSB",1);
452   return(reduce(I,J,1));
453}
454
455proc isIncluded(ideal I, ideal J, list #)
456"USAGE:  isIncluded(I,J); I,J ideals
457RETURN:  1 if J includes I,
458@*       0 if there is an element f in I which does not reduce to 0 w.r.t. J.
459EXAMPLE: example isIncluded; shows an example
460"
461{
462   def R = basering;
463   setring R;
464
465   attrib(J,"isSB",1);
466   int i,j,k;
467
468   if(size(#) > 0)
469   {
470      int n = #[1];
471      if(n >= ncols(I)) { n = ncols(I); }
472      if(n > 1)
473      {
474         for(i = 1; i <= n - 1; i++)
475         {
476            //link l(i) = "MPtcp:fork";
477            link l(i) = "ssi:fork";
478            open(l(i));
479
480            write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1)));
481         }
482
483         int t = timer;
484         if(reduce(I[ncols(I)], J, 1) != 0)
485         {
486            for(i = 1; i <= n - 1; i++)
487            {
488               close(l(i));
489            }
490            return(0);
491         }
492         t = timer - t;
493         if(t > 60) { t = 60; }
494         int i_sleep = system("sh", "sleep "+string(t));
495
496         j = ncols(I) - n;
497
498         while(j >= 0)
499         {
500            for(i = 1; i <= n - 1; i++)
501            {
502               if(status(l(i), "read", "ready"))
503               {
504                  if(read(l(i)) != 0)
505                  {
506                     for(i = 1; i <= n - 1; i++)
507                     {
508                        close(l(i));
509                     }
510                     return(0);
511                  }
512                  else
513                  {
514                     if(j >= 1)
515                     {
516                        write(l(i), quote(redFork(eval(I[j]), eval(J), 1)));
517                        j--;
518                     }
519                     else
520                     {
521                        k++;
522                        close(l(i));
523                     }
524                  }
525               }
526            }
527            if(k == n - 1)
528            {
529               j--;
530            }
531            i_sleep = system("sh", "sleep "+string(t));
532         }
533         return(1);
534      }
535   }
536
537   for(i = ncols(I); i >= 1; i--)
538   {
539      if(reduce(I[i],J,1) != 0){ return(0); }
540   }
541   return(1);
542}
543example
544{ "EXAMPLE:"; echo = 2;
545   ring r=0,(x,y,z),dp;
546   ideal I = x+1,x+y+1;
547   ideal J = x+1,y;
548   isIncluded(I,J);
549   isIncluded(J,I);
550   isIncluded(I,J,4);
551
552   ring R = 0, x(1..5), dp;
553   ideal I1 = cyclic(4);
554   ideal I2 = I1,x(5)^2;
555   isIncluded(I1,I2,4);
556}
557
558proc deleteUnluckyPrimes(list T, list L, int ho, list #)
559"USAGE:  deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer
560RETURN:  lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #),
561         lT ideal
562NOTE:    - if ho = 1, the polynomials in T are homogeneous, else ho = 0,
563@*       - lT is prevalent, i.e. the most appearing leading ideal in T
564EXAMPLE: example deleteUnluckyPrimes; shows an example
565"
566{
567   ho = ((ho)||(ord_test(basering) == -1));
568   int j,k,c;
569   intvec hl,hc;
570   ideal cT,lT,cK;
571   lT = lead(T[size(T)]);
572   attrib(lT,"isSB",1);
573   if(!ho)
574   {
575      for(j = 1; j < size(T); j++)
576      {
577         cT = lead(T[j]);
578         attrib(cT,"isSB",1);
579         if((size(reduce(cT,lT,5))!=0)||(size(reduce(lT,cT,5))!=0))
580         {
581            cK = cT;
582            c++;
583         }
584      }
585      if(c > size(T) div 2){ lT = cK; }
586   }
587   else
588   {
589      hl = hilb(lT,1);
590      for(j = 1; j < size(T); j++)
591      {
592         cT = lead(T[j]);
593         attrib(cT,"isSB",1);
594         hc = hilb(cT,1);
595         if(hl == hc)
596         {
597            for(k = 1; k <= size(lT); k++)
598            {
599               if(lT[k] < cT[k]) { lT = cT; c++; break; }
600               if(lT[k] > cT[k]) { c++; break; }
601            }
602         }
603         else
604         {
605            if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; }
606         }
607      }
608   }
609
610   int addList;
611   if(size(#) > 0) { list M = #; addList = 1; }
612   j = 1;
613   attrib(lT,"isSB",1);
614   while((j <= size(T))&&(c > 0))
615   {
616      cT = lead(T[j]);
617      attrib(cT,"isSB",1);
618      if((size(reduce(cT,lT,5)) != 0)||(size(reduce(lT,cT,5)) != 0))
619      {
620         T = delete(T,j);
621         if(j == 1)
622         {
623            L = L[2..size(L)];
624            if(addList == 1) { M = M[2..size(M)]; }
625         }
626         else
627         {
628            if(j == size(L))
629            {
630               L = L[1..size(L)-1];
631               if(addList == 1) { M = M[1..size(M)-1]; }
632            }
633            else
634            {
635               L = L[1..j-1],L[j+1..size(L)];
636               if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; }
637            }
638         }
639         j--;
640      }
641      j++;
642   }
643
644   for(j = 1; j <= size(L); j++)
645   {
646      L[j] = bigint(L[j]);
647   }
648
649   if(addList == 0) { return(list(T,L,lT)); }
650   if(addList == 1) { return(list(T,L,M,lT)); }
651}
652example
653{ "EXAMPLE:"; echo = 2;
654   list L = 2,3,5,7,11;
655   ring r = 0,(y,x),Dp;
656   ideal I1 = 2y2x,y6;
657   ideal I2 = yx2,y3x,x5,y6;
658   ideal I3 = y2x,x3y,x5,y6;
659   ideal I4 = y2x,11x3y,x5;
660   ideal I5 = y2x,yx3,x5,7y6;
661   list T = I1,I2,I3,I4,I5;
662   deleteUnluckyPrimes(T,L,1);
663   list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5);
664   deleteUnluckyPrimes(T,L,1,P);
665}
666
667proc primeTest(def II, bigint p)
668{
669   if(typeof(II) == "string")
670   {
671      ideal I = `II`.Ideal;
672   }
673   else
674   {
675      ideal I = II;
676   }
677
678   I = simplify(I, 2);   // erase zero generators
679
680   int i,j;
681   poly f;
682   number cnt;
683   for(i = 1; i <= size(I); i++)
684   {
685      f = cleardenom(I[i]);
686      if(f == 0) { return(0); }
687      cnt = leadcoef(I[i])/leadcoef(f);
688      if((bigint(numerator(cnt)) mod p) == 0) { return(0); }
689      if((bigint(denominator(cnt)) mod p) == 0) { return(0); }
690      for(j = size(f); j > 0; j--)
691      {
692         if((bigint(leadcoef(f[j])) mod p) == 0) { return(0); }
693      }
694   }
695   return(1);
696}
697
698proc primeList(ideal I, int n, list #)
699"USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
700         n integer
701RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
702         < L[size(L)] union with L) such that none of these primes divides any
703         coefficient occurring in I
704NOTE:    The number of cores to use can be defined by ncores, default is 1.
705EXAMPLE: example primeList; shows an example
706"
707{
708   intvec L;
709   int i,p;
710   int ncores = 1;
711
712//-----------------  Initialize optional parameter ncores  ---------------------
713   if(size(#) > 0)
714   {
715      if(size(#) == 1)
716      {
717         if(typeof(#[1]) == "int")
718         {
719            ncores = #[1];
720            # = list();
721         }
722      }
723      else
724      {
725         ncores = #[2];
726      }
727   }
728
729   if(size(#) == 0)
730   {
731      p = 2147483647;
732      while(!primeTest(I,p))
733      {
734         p = prime(p-1);
735         if(p == 2) { ERROR("no more primes"); }
736      }
737      L[1] = p;
738   }
739   else
740   {
741      L = #[1];
742      p = prime(L[size(L)]-1);
743      while(!primeTest(I,p))
744      {
745         p = prime(p-1);
746         if(p == 2) { ERROR("no more primes"); }
747      }
748      L[size(L)+1] = p;
749   }
750   if(p == 2) { ERROR("no more primes"); }
751   if(ncores == 1)
752   {
753      for(i = 2; i <= n; i++)
754      {
755         p = prime(p-1);
756         while(!primeTest(I,p))
757         {
758            p = prime(p-1);
759            if(p == 2) { ERROR("no more primes"); }
760         }
761         L[size(L)+1] = p;
762      }
763   }
764   else
765   {
766      int neededSize = size(L)+n-1;;
767      list parallelResults;
768      list arguments;
769      int neededPrimes = neededSize-size(L);
770      idealPrimeTest Id;
771      Id.Ideal = I;
772      export(Id);
773      while(neededPrimes > 0)
774      {
775         arguments = list();
776         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
777            *ncores; i > 0; i--)
778         {
779            p = prime(p-1);
780            if(p == 2) { ERROR("no more primes"); }
781            arguments[i] = list("Id", p);
782         }
783         parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores);
784         for(i = size(arguments); i > 0; i--)
785         {
786            if(parallelResults[i])
787            {
788               L[size(L)+1] = arguments[i][2];
789            }
790         }
791         neededPrimes = neededSize-size(L);
792      }
793      kill Id;
794      if(size(L) > neededSize)
795      {
796         L = L[1..neededSize];
797      }
798   }
799   return(L);
800}
801example
802{ "EXAMPLE:"; echo = 2;
803   ring r = 0,(x,y,z),dp;
804   ideal I = 2147483647x+y, z-181;
805   intvec L = primeList(I,10);
806   size(L);
807   L[1];
808   L[size(L)];
809   L = primeList(I,5,L);
810   size(L);
811   L[size(L)];
812}
813
814proc modStdL(def I, list #)
815"USAGE:   modStdL(I[, exactness]); I ideal/module, exactness int
816RETURN:   a standard basis of I
817NOTE:     The procedure computes a standard basis of I (over the rational
818          numbers) by using modular methods via an external Singular.
819       @* An optional parameter 'exactness' can be provided.
820          If exactness = 1(default), the procedure computes a standard basis
821          of I for sure; if exactness = 0, it computes a standard basis of I
822          with high probability.
823SEE ALSO: modular, modStd
824EXAMPLE:  example modStdL; shows an example"
825{
826  link l="ssi:tcp localhost:"+system("Singular");
827  write(l,quote(option(noloadLib))); // suppress "loaded..."
828  read(l); //dummy: return value of option
829  write(l,quote(load("modstd.lib","with"))); // load library
830  read(l); //dummy: return value of load
831  if (size(#)==0)
832  {
833    write(l,quote(modStd(eval(I))));
834  }
835  else
836  {
837    write(l,quote(modStd(eval(I),eval(#[1]))));
838  }
839  return(read(l));
840}
841example
842{
843    "EXAMPLE:";
844    echo = 2;
845    ring R1 = 0, (x,y,z,t), dp;
846    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
847    ideal J = modStdL(I);
848    J;
849    I = homog(I, t);
850    J = modStdL(I);
851    J;
852
853    ring R2 = 0, (x,y,z), ds;
854    ideal I = jacob(x5+y6+z7+xyz);
855    ideal J = modStdL(I, 0);
856    J;
857
858    ring R3 = 0, x(1..4), lp;
859    ideal I = cyclic(4);
860    ideal J1 = modStdL(I, 1);   // default
861    ideal J2 = modStdL(I, 0);
862    size(reduce(J1, J2));
863    size(reduce(J2, J1));
864}
865
866////////////////////////////// further examples ////////////////////////////////
867
868/*
869ring r = 0, (x,y,z), lp;
870poly s1 = 5x3y2z+3y3x2z+7xy2z2;
871poly s2 = 3xy2z2+x5+11y2z2;
872poly s3 = 4xyz+7x3+12y3+1;
873poly s4 = 3x3-4y3+yz2;
874ideal i =  s1, s2, s3, s4;
875
876ring r = 0, (x,y,z), lp;
877poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
878poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
879poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
880poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
881ideal i =  s1, s2, s3, s4;
882
883ring r = 0, (x,y,z), lp;
884poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
885poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
886poly s3 = 8x3 + 12y3 + xz2 + 3;
887poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
888ideal i = s1, s2, s3, s4;
889
890int n = 6;
891ring r = 0,(x(1..n)),lp;
892ideal i = cyclic(n);
893ring s = 0, (x(1..n),t), lp;
894ideal i = imap(r,i);
895i = homog(i,t);
896
897ring r = 0, (x(1..4),s), (dp(4),dp);
898poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
899poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
900poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
901poly 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);
902poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
903ideal i =  s1, s2, s3, s4, s5;
904
905ring r = 0, (x,y,z), ds;
906int a = 16;
907int b = 15;
908int c = 4;
909int t = 1;
910poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
911         +x^(c-2)*y^c*(y2+t*x)^2;
912ideal i = jacob(f);
913
914ring r = 0, (x,y,z), ds;
915int a = 25;
916int b = 25;
917int c = 5;
918int t = 1;
919poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
920         +x^(c-2)*y^c*(y2+t*x)^2;
921ideal i = jacob(f),f;
922
923ring r = 0, (x,y,z), ds;
924int a = 10;
925poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
926ideal i = jacob(f);
927
928ring r = 0, (x,y,z), ds;
929int a = 6;
930int b = 8;
931int c = 10;
932int alpha = 5;
933int beta = 5;
934int t = 1;
935poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
936         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
937ideal i = jacob(f);
938*/
939
Note: See TracBrowser for help on using the repository browser.