source: git/Singular/LIB/modstd.lib @ 1387b04

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