source: git/Singular/LIB/modwalk.lib @ 8d8fef

spielwiese
Last change on this file since 8d8fef was 468f1c, checked in by Andreas Steenpass <steenpass@…>, 5 years ago
load ring.lib when using create_ring()
  • Property mode set to 100644
File size: 17.2 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version modwalk.lib 4.1.2.0 Feb_2019 "; // $Id$
3category = "Commutative Algebra";
4info="
5LIBRARY:  modwalk.lib      Groebner basis convertion
6
7AUTHORS:  S. Oberfranz    oberfran@mathematik.uni-kl.de
8
9OVERVIEW:
10
11  A library for converting Groebner bases of an ideal in the polynomial
12  ring over the rational numbers using modular methods. The procedures are
13  inspired by the following paper:
14  Elizabeth A. Arnold: Modular algorithms for computing Groebner bases.
15  Journal of Symbolic Computation 35, 403-419 (2003).
16
17PROCEDURES:
18
19modWalk(I,#);                   standard basis conversion of I by Groebner Walk using modular methods
20modrWalk(I,radius,#);           standard basis conversion of I by Random Walk using modular methods
21modfWalk(I,#);                  standard basis conversion of I by Fractal Walk using modular methods
22modfrWalk(I,radius,#);          standard basis conversion of I by Random Fractal Walk using modular methods
23
24KEYWORDS: walk, groebner;Groebnerwalk
25SEE ALSO: grwalk_lib;swalk_lib;rwalk_lib
26";
27
28LIB "rwalk.lib";
29LIB "grwalk.lib";
30LIB "modular.lib";
31LIB "ring.lib";
32
33proc modWalk(ideal I, list #)
34"USAGE:   modWalk(I, [, v, w]); I ideal, v intvec or string, w intvec
35          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
36          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
37          respect to dp or Dp, respectively.
38          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
39          the output will be a standard basis with respect to lp.
40          If no optional argument is given, I is assumed to be a standard basis with respect to dp
41          and a standard basis with respect to lp will be computed.
42RETURN:   a standard basis of I
43NOTE:     The procedure computes a standard basis of I (over the rational
44          numbers) by using modular methods.
45SEE ALSO: modular
46EXAMPLE:  example modWalk; shows an example"
47{
48    /* save options */
49    intvec opt = option(get);
50    option(redSB);
51
52    /* call modular() */
53    if (size(#) > 0) {
54        I = modular("gwalk", list(I,#), primeTest_std,
55            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
56    }
57    else {
58        I = modular("gwalk", list(I), primeTest_std,
59            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
60    }
61
62    /* return the result */
63    attrib(I, "isSB", 1);
64    option(set, opt);
65    return(I);
66}
67example
68{
69    "EXAMPLE:";
70    echo = 2;
71    ring R1 = 0, (x,y,z,t), dp;
72    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
73    I = std(I);
74    ring R2 = 0, (x,y,z,t), lp;
75    ideal I = fetch(R1, I);
76    ideal J = modWalk(I);
77    J;
78    ring S1 = 0, (a,b,c,d), Dp;
79    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
80    I = std(I);
81    ring S2 = 0, (c,d,b,a), lp;
82    ideal I = fetch(S1,I);
83    // I is assumed to be a Dp-Groebner basis.
84    // We compute a lp-Groebner basis.
85    ideal J = modWalk(I,"Dp");
86    J;
87    intvec w = 3,2,1,2;
88    ring S3 = 0, (c,d,b,a), (a(w),lp);
89    ideal I = fetch(S1,I);
90    // I is assumed to be a Dp-Groebner basis.
91    // We compute a (a(w),lp)-Groebner basis.
92    ideal J = modWalk(I,"Dp",w);
93    J;
94}
95
96proc modrWalk(ideal I, int radius, list #)
97"USAGE:   modrWalk(I, radius[, v, w]);
98          I ideal, radius int, pertdeg int, v intvec or string, w intvec
99          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
100          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
101          respect to dp or Dp, respectively.
102          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
103          the output will be a standard basis with respect to lp.
104          If no optional argument is given, I is assumed to be a standard basis with respect to dp
105          and a standard basis with respect to lp will be computed.
106RETURN:   a standard basis of I
107NOTE:     The procedure computes a standard basis of I (over the rational
108          numbers) by using modular methods.
109SEE ALSO: modular
110EXAMPLE:  example modrWalk; shows an example"
111{
112    /* save options */
113    intvec opt = option(get);
114    option(redSB);
115
116    /* call modular() */
117    if (size(#) > 0) {
118        I = modular("rwalk", list(I,radius,1,#), primeTest_std,
119            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
120    }
121    else {
122        I = modular("rwalk", list(I,radius,1), primeTest_std, deleteUnluckyPrimes_std,
123            pTest_std,finalTest_std);
124    }
125
126    /* return the result */
127    attrib(I, "isSB", 1);
128    option(set, opt);
129    return(I);
130}
131example
132{
133    "EXAMPLE:";
134    echo = 2;
135    ring R1 = 0, (x,y,z,t), dp;
136    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
137    I = std(I);
138    ring R2 = 0, (x,y,z,t), lp;
139    ideal I = fetch(R1, I);
140    int radius = 2;
141    ideal J = modrWalk(I,radius);
142    J;
143    ring S1 = 0, (a,b,c,d), Dp;
144    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
145    I = std(I);
146    ring S2 = 0, (c,d,b,a), lp;
147    ideal I = fetch(S1,I);
148    // I is assumed to be a Dp-Groebner basis.
149    // We compute a lp-Groebner basis.
150    ideal J = modrWalk(I,radius,"Dp");
151    J;
152    intvec w = 3,2,1,2;
153    ring S3 = 0, (c,d,b,a), (a(w),lp);
154    ideal I = fetch(S1,I);
155    // I is assumed to be a Dp-Groebner basis.
156    // We compute a (a(w),lp)-Groebner basis.
157    ideal J = modrWalk(I,radius,"Dp",w);
158    J;
159}
160
161proc modfWalk(ideal I, list #)
162"USAGE:   modfWalk(I, [, v, w]); I ideal, v intvec or string, w intvec
163          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
164          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
165          respect to dp or Dp, respectively.
166          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
167          the output will be a standard basis with respect to lp.
168          If no optional argument is given, I is assumed to be a standard basis with respect to dp
169          and a standard basis with respect to lp will be computed.
170RETURN:   a standard basis of I
171NOTE:     The procedure computes a standard basis of I (over the rational
172          numbers) by using modular methods.
173SEE ALSO: modular
174EXAMPLE:  example modfWalk; shows an example"
175{
176    /* save options */
177    intvec opt = option(get);
178    option(redSB);
179
180    /* call modular() */
181    if (size(#) > 0) {
182        I = modular("fwalk", list(I,#), primeTest_std,
183            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
184    }
185    else {
186        I = modular("fwalk", list(I), primeTest_std,
187            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
188    }
189
190    /* return the result */
191    attrib(I, "isSB", 1);
192    option(set, opt);
193    return(I);
194}
195example
196{
197    "EXAMPLE:";
198    echo = 2;
199    ring R1 = 0, (x,y,z,t), dp;
200    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
201    I = std(I);
202    ring R2 = 0, (x,y,z,t), lp;
203    ideal I = fetch(R1, I);
204    ideal J = modfWalk(I);
205    J;
206    ring S1 = 0, (a,b,c,d), Dp;
207    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
208    I = std(I);
209    ring S2 = 0, (c,d,b,a), lp;
210    ideal I = fetch(S1,I);
211    // I is assumed to be a Dp-Groebner basis.
212    // We compute a lp-Groebner basis.
213    ideal J = modfWalk(I,"Dp");
214    J;
215    intvec w = 3,2,1,2;
216    ring S3 = 0, (c,d,b,a), (a(w),lp);
217    ideal I = fetch(S1,I);
218    // I is assumed to be a Dp-Groebner basis.
219    // We compute a (a(w),lp)-Groebner basis.
220    ideal J = modfWalk(I,"Dp",w);
221    J;
222}
223
224proc modfrWalk(ideal I, int radius, list #)
225"USAGE:   modfrWalk(I, radius [, v, w]); I ideal, radius int, v intvec or string, w intvec
226          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
227          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
228          respect to dp or Dp, respectively.
229          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
230          the output will be a standard basis with respect to lp.
231          If no optional argument is given, I is assumed to be a standard basis with respect to dp
232          and a standard basis with respect to lp will be computed.
233RETURN:   a standard basis of I
234NOTE:     The procedure computes a standard basis of I (over the rational
235          numbers) by using modular methods.
236SEE ALSO: modular
237EXAMPLE:  example modfrWalk; shows an example"
238{
239    /* save options */
240    intvec opt = option(get);
241    option(redSB);
242
243    /* call modular() */
244    if (size(#) > 0) {
245        I = modular("frandwalk", list(I,radius,#), primeTest_std,
246            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
247    }
248    else {
249        I = modular("frandwalk", list(I,radius), primeTest_std,
250            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
251    }
252
253    /* return the result */
254    attrib(I, "isSB", 1);
255    option(set, opt);
256    return(I);
257}
258example
259{
260    "EXAMPLE:";
261    echo = 2;
262    ring R1 = 0, (x,y,z,t), dp;
263    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
264    I = std(I);
265    ring R2 = 0, (x,y,z,t), lp;
266    ideal I = fetch(R1, I);
267    int radius = 2;
268    ideal J = modfrWalk(I,radius);
269    J;
270    ring S1 = 0, (a,b,c,d), Dp;
271    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
272    I = std(I);
273    ring S2 = 0, (c,d,b,a), lp;
274    ideal I = fetch(S1,I);
275    // I is assumed to be a Dp-Groebner basis.
276    // We compute a lp-Groebner basis.
277    ideal J = modfrWalk(I,radius,"Dp");
278    J;
279    intvec w = 3,2,1,2;
280    ring S3 = 0, (c,d,b,a), (a(w),lp);
281    ideal I = fetch(S1,I);
282    // I is assumed to be a Dp-Groebner basis.
283    // We compute a (a(w),lp)-Groebner basis.
284    ideal J = modfrWalk(I,radius,"Dp",w);
285    J;
286}
287
288/* test if the prime p is suitable for the input, i.e. it does not divide
289 * the numerator or denominator of any of the coefficients */
290static proc primeTest_std(int p, alias list args)
291{
292    /* erase zero generators */
293    ideal I = simplify(args[1], 2);
294
295    /* clear denominators and count the terms */
296    ideal J;
297    ideal K;
298    int n = ncols(I);
299    intvec sizes;
300    number cnt;
301    int i;
302    for(i = n; i > 0; i--) {
303        J[i] = cleardenom(I[i]);
304        cnt = leadcoef(J[i])/leadcoef(I[i]);
305        K[i] = numerator(cnt)*var(1)+denominator(cnt);
306    }
307    sizes = size(J[1..n]);
308
309    /* change to characteristic p */
310    def br = basering;
311    list lbr = ringlist(br);
312    if (typeof(lbr[1]) == "int") {
313        lbr[1] = p;
314    }
315    else {
316        lbr[1][1] = p;
317    }
318    def rp = ring(lbr);
319    setring(rp);
320    ideal Jp = fetch(br, J);
321    ideal Kp = fetch(br, K);
322
323    /* test if any coefficient is missing */
324    if (intvec(size(Kp[1..n])) != 2:n) {
325        setring(br);
326        return(0);
327    }
328    if (intvec(size(Jp[1..n])) != sizes) {
329        setring(br);
330        return(0);
331    }
332    setring(br);
333    return(1);
334}
335
336/* find entries in modresults which come from unlucky primes.
337 * For this, sort the entries into categories depending on their leading
338 * ideal and return the indices in all but the biggest category. */
339static proc deleteUnluckyPrimes_std(alias list modresults)
340{
341    int size_modresults = size(modresults);
342
343    /* sort results into categories.
344     * each category is represented by three entries:
345     * - the corresponding leading ideal
346     * - the number of elements
347     * - the indices of the elements
348     */
349    list cat;
350    int size_cat;
351    ideal L;
352    int i;
353    int j;
354    for (i = 1; i <= size_modresults; i++) {
355        L = lead(modresults[i]);
356        attrib(L, "isSB", 1);
357        for (j = 1; j <= size_cat; j++) {
358            if (size(L) == size(cat[j][1])
359                && size(reduce(L, cat[j][1],5)) == 0
360                && size(reduce(cat[j][1], L,5)) == 0)
361            {
362                cat[j][2] = cat[j][2]+1;
363                cat[j][3][cat[j][2]] = i;
364                break;
365            }
366        }
367        if (j > size_cat) {
368            size_cat++;
369            cat[size_cat] = list();
370            cat[size_cat][1] = L;
371            cat[size_cat][2] = 1;
372            cat[size_cat][3] = list(i);
373        }
374    }
375
376    /* find the biggest categories */
377    int cat_max = 1;
378    int max = cat[1][2];
379    for (i = 2; i <= size_cat; i++) {
380        if (cat[i][2] > max) {
381            cat_max = i;
382            max = cat[i][2];
383        }
384    }
385
386    /* return all other indices */
387    list unluckyIndices;
388    for (i = 1; i <= size_cat; i++) {
389        if (i != cat_max) {
390            unluckyIndices = unluckyIndices + cat[i][3];
391        }
392    }
393    return(unluckyIndices);
394}
395
396/* test if 'command' applied to 'args' in characteristic p is the same as
397   'result' mapped to characteristic p */
398static proc pTest_std(string command, alias list args, alias ideal result,
399    int p)
400{
401    /* change to characteristic p */
402    def br = basering;
403    list lbr = ringlist(br);
404    if (typeof(lbr[1]) == "int") {
405        lbr[1] = p;
406    }
407    else {
408        lbr[1][1] = p;
409    }
410    def rp = ring(lbr);
411    setring(rp);
412    ideal Ip = fetch(br, args)[1];
413    list Arg = fetch(br, args);
414    string exstr;
415    ideal Gp = fetch(br, result);
416    attrib(Gp, "isSB", 1);
417
418    /* test if Ip is in Gp */
419    int i;
420    for (i = ncols(Ip); i > 0; i--) {
421        if (reduce(Ip[i], Gp, 1) != 0) {
422            setring(br);
423            return(0);
424        }
425    }
426
427    /* compute command(args) */
428    exstr = "Ip = "+command+" (Ip";
429
430    for(i=2; i<=size(Arg); i++) {
431      exstr = exstr+",Arg["+string(eval(i))+"]";
432      }
433    exstr = exstr+");";
434
435    execute(exstr);
436
437    /* test if Gp is in Ip */
438    for (i = ncols(Gp); i > 0; i--) {
439        if (reduce(Gp[i], Ip, 1) != 0) {
440            setring(br);
441            return(0);
442        }
443    }
444
445    setring(br);
446    return(1);
447}
448
449/* test if 'result' is a GB of the input ideal */
450static proc finalTest_std(string command, alias list args, ideal result)
451{
452    /* test if args[1] is in result */
453    attrib(result, "isSB", 1);
454    int i;
455    for (i = ncols(args[1]); i > 0; i--) {
456        if (reduce(args[1][i], result, 1) != 0) {
457            return(0);
458        }
459    }
460
461   /* test if result is in args[1].                      */
462   /* args[1] is given by a Groebner basis. Thus we may  */
463   /* reduce the result with respect to args[1].         */
464   int n=nvars(basering);
465   string ord_str = "dp";
466
467   for(i=2; i<=size(args); i++)
468   {
469     if(typeof(args[i]) == "list") {
470       if(typeof(args[i][1]) == "intvec") {
471         ord_str = "(a("+string(args[i][1])+"),lp("+string(n) + "),C)";
472         break;
473       }
474       if(typeof(args[i][1]) == "string") {
475         if(args[i][1] == "Dp") {
476           ord_str = "Dp";
477         }
478         break;
479       }
480     }
481   }
482   ideal xI = args[1];
483   ring xR = basering;
484   ring yR = create_ring(ringlist(xR)[1], "("+varstr(xR)+")", ord_str, "no_minpoly");
485   ideal yI = fetch(xR,xI);
486   ideal yresult = fetch(xR,result);
487   attrib(yI, "isSB", 1);
488   for(i=size(yresult); i>0; i--)
489   {
490     if(reduce(yresult[i],yI) != 0)
491     {
492       return(0);
493     }
494   }
495   setring xR;
496   kill yR;
497
498   /* test if result is a Groebner basis */
499    link l1="ssi:fork";
500    open(l1);
501    link l2="ssi:fork";
502    open(l2);
503    list l=list(l1,l2);
504    write(l1,quote(TestSBred(result)));
505    write(l2,quote(TestSBstd(result)));
506    i=waitfirst(l);
507    if(i==1) {
508      i=read(l1);
509      }
510    else {
511      i=read(l2);
512      }
513    close(l1);
514    close(l2);
515    return(i);
516}
517
518/* return 1, if I_reduce is _not_ in G_reduce,
519 *        0, otherwise
520 * (same as size(reduce(I_reduce, G_reduce))).
521 * Uses parallelization. */
522static proc reduce_parallel(def I_reduce, def G_reduce)
523{
524    exportto(Modwalk, I_reduce);
525    exportto(Modwalk, G_reduce);
526    int size_I = ncols(I_reduce);
527    int chunks = Modular::par_range(size_I);
528    intvec range;
529    int i;
530    for (i = chunks; i > 0; i--) {
531        range = Modular::par_range(size_I, i);
532        task t(i) = "Modwalk::reduce_task", list(range);
533    }
534    startTasks(t(1..chunks));
535    waitAllTasks(t(1..chunks));
536    int result = 0;
537    for (i = chunks; i > 0; i--) {
538        if (getResult(t(i))) {
539            result = 1;
540            break;
541        }
542    }
543    kill I_reduce;
544    kill G_reduce;
545    return(result);
546}
547
548/* compute a chunk of reductions for reduce_parallel */
549static proc reduce_task(intvec range)
550{
551    int result = 0;
552    int i;
553    for (i = range[1]; i <= range[2]; i++) {
554        if (reduce(I_reduce[i], G_reduce, 1) != 0) {
555            result = 1;
556            break;
557        }
558    }
559    return(result);
560}
561
562/* test if result is a GB with std*/
563static proc TestSBstd(ideal result)
564{
565  ideal G = std(result);
566  if(reduce_parallel(G,result)) {
567    return(0);
568    }
569  return(1);
570}
571
572/* test if result is a GB by reducing s-polynomials*/
573static proc TestSBred(ideal result)
574{
575  int i,j;
576  for(i=1; i<=size(result); i++) {
577    for(j=i; j<=size(result); j++) {
578      if(reduce(sPolynomial(result[i],result[j]),result)!=0) {
579        return(0);
580        }
581      }
582    }
583  return(1);
584}
585
586/* compute s-polynomial of f and g */
587static proc sPolynomial(poly f,poly g)
588{
589  int i;
590  poly lcmp = 1;
591
592  intvec lexpf = leadexp(f);
593  intvec lexpg = leadexp(g);
594
595  for(i=1; i<=nvars(basering); i++) {
596    if(lexpf[i]>=lexpg[i]) {
597      lcmp=lcmp*var(i)**lexpf[i];
598      }
599    else {
600      lcmp=lcmp*var(i)**lexpg[i];
601      }
602   }
603
604  poly fmult=lcmp/leadmonom(f);
605  poly gmult=lcmp/leadmonom(g);
606  poly result=leadcoef(g)*fmult*f-leadcoef(f)*gmult*g;
607
608  return(result);
609}
Note: See TracBrowser for help on using the repository browser.