source: git/Singular/LIB/modwalk.lib @ 291c20

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