singular-data 1:4.1.0-p3+ds-2build1 / usr / share / singular / LIB / modwalk.lib

///////////////////////////////////////////////////////////////////////////////
version="version modwalk.lib 4.0.0.0 Jun_2013 "; // $Id: fb19d7ae95f2dfd3c0d4ce218805b7f15c730b9b $
category = "Commutative Algebra";
info="
LIBRARY:  modwalk.lib      Groebner basis convertion

AUTHORS:  S. Oberfranz    oberfran@mathematik.uni-kl.de

OVERVIEW:

  A library for converting Groebner bases of an ideal in the polynomial
  ring over the rational numbers using modular methods. The procedures are
  inspired by the following paper:
  Elizabeth A. Arnold: Modular algorithms for computing Groebner bases.
  Journal of Symbolic Computation 35, 403-419 (2003).

PROCEDURES:

modWalk(I,#);                   standard basis conversion of I by Groebner Walk using modular methods
modrWalk(I,radius,#);           standard basis conversion of I by Random Walk using modular methods
modfWalk(I,#);                  standard basis conversion of I by Fractal Walk using modular methods
modfrWalk(I,radius,#);          standard basis conversion of I by Random Fractal Walk using modular methods

KEYWORDS: walk, groebner;Groebnerwalk
SEE ALSO: grwalk_lib;swalk_lib;rwalk_lib
";

LIB "rwalk.lib";
LIB "grwalk.lib";
LIB "modular.lib";

proc modWalk(ideal I, list #)
"USAGE:   modWalk(I, [, v, w]); I ideal, v intvec or string, w intvec
          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
          respect to dp or Dp, respectively.
          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
          the output will be a standard basis with respect to lp.
          If no optional argument is given, I is assumed to be a standard basis with respect to dp
          and a standard basis with respect to lp will be computed.
RETURN:   a standard basis of I
NOTE:     The procedure computes a standard basis of I (over the rational
          numbers) by using modular methods.
SEE ALSO: modular
EXAMPLE:  example modWalk; shows an example"
{
    /* save options */
    intvec opt = option(get);
    option(redSB);

    /* call modular() */
    if (size(#) > 0) {
        I = modular("gwalk", list(I,#), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }
    else {
        I = modular("gwalk", list(I), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }

    /* return the result */
    attrib(I, "isSB", 1);
    option(set, opt);
    return(I);
}
example
{
    "EXAMPLE:";
    echo = 2;
    ring R1 = 0, (x,y,z,t), dp;
    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
    I = std(I);
    ring R2 = 0, (x,y,z,t), lp;
    ideal I = fetch(R1, I);
    ideal J = modWalk(I);
    J;
    ring S1 = 0, (a,b,c,d), Dp;
    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
    I = std(I);
    ring S2 = 0, (c,d,b,a), lp;
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a lp-Groebner basis.
    ideal J = modWalk(I,"Dp");
    J;
    intvec w = 3,2,1,2;
    ring S3 = 0, (c,d,b,a), (a(w),lp);
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a (a(w),lp)-Groebner basis.
    ideal J = modWalk(I,"Dp",w);
    J;
}

proc modrWalk(ideal I, int radius, list #)
"USAGE:   modrWalk(I, radius[, v, w]);
          I ideal, radius int, pertdeg int, v intvec or string, w intvec
          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
          respect to dp or Dp, respectively.
          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
          the output will be a standard basis with respect to lp.
          If no optional argument is given, I is assumed to be a standard basis with respect to dp
          and a standard basis with respect to lp will be computed.
RETURN:   a standard basis of I
NOTE:     The procedure computes a standard basis of I (over the rational
          numbers) by using modular methods.
SEE ALSO: modular
EXAMPLE:  example modrWalk; shows an example"
{
    /* save options */
    intvec opt = option(get);
    option(redSB);

    /* call modular() */
    if (size(#) > 0) {
        I = modular("rwalk", list(I,radius,1,#), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }
    else {
        I = modular("rwalk", list(I,radius,1), primeTest_std, deleteUnluckyPrimes_std,
            pTest_std,finalTest_std);
    }

    /* return the result */
    attrib(I, "isSB", 1);
    option(set, opt);
    return(I);
}
example
{
    "EXAMPLE:";
    echo = 2;
    ring R1 = 0, (x,y,z,t), dp;
    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
    I = std(I);
    ring R2 = 0, (x,y,z,t), lp;
    ideal I = fetch(R1, I);
    int radius = 2;
    ideal J = modrWalk(I,radius);
    J;
    ring S1 = 0, (a,b,c,d), Dp;
    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
    I = std(I);
    ring S2 = 0, (c,d,b,a), lp;
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a lp-Groebner basis.
    ideal J = modrWalk(I,radius,"Dp");
    J;
    intvec w = 3,2,1,2;
    ring S3 = 0, (c,d,b,a), (a(w),lp);
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a (a(w),lp)-Groebner basis.
    ideal J = modrWalk(I,radius,"Dp",w);
    J;
}

proc modfWalk(ideal I, list #)
"USAGE:   modfWalk(I, [, v, w]); I ideal, v intvec or string, w intvec
          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
          respect to dp or Dp, respectively.
          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
          the output will be a standard basis with respect to lp.
          If no optional argument is given, I is assumed to be a standard basis with respect to dp
          and a standard basis with respect to lp will be computed.
RETURN:   a standard basis of I
NOTE:     The procedure computes a standard basis of I (over the rational
          numbers) by using modular methods.
SEE ALSO: modular
EXAMPLE:  example modfWalk; shows an example"
{
    /* save options */
    intvec opt = option(get);
    option(redSB);

    /* call modular() */
    if (size(#) > 0) {
        I = modular("fwalk", list(I,#), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }
    else {
        I = modular("fwalk", list(I), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }

    /* return the result */
    attrib(I, "isSB", 1);
    option(set, opt);
    return(I);
}
example
{
    "EXAMPLE:";
    echo = 2;
    ring R1 = 0, (x,y,z,t), dp;
    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
    I = std(I);
    ring R2 = 0, (x,y,z,t), lp;
    ideal I = fetch(R1, I);
    ideal J = modfWalk(I);
    J;
    ring S1 = 0, (a,b,c,d), Dp;
    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
    I = std(I);
    ring S2 = 0, (c,d,b,a), lp;
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a lp-Groebner basis.
    ideal J = modfWalk(I,"Dp");
    J;
    intvec w = 3,2,1,2;
    ring S3 = 0, (c,d,b,a), (a(w),lp);
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a (a(w),lp)-Groebner basis.
    ideal J = modfWalk(I,"Dp",w);
    J;
}

proc modfrWalk(ideal I, int radius, list #)
"USAGE:   modfrWalk(I, radius [, v, w]); I ideal, radius int, v intvec or string, w intvec
          If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp).
          If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with
          respect to dp or Dp, respectively.
          If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise,
          the output will be a standard basis with respect to lp.
          If no optional argument is given, I is assumed to be a standard basis with respect to dp
          and a standard basis with respect to lp will be computed.
RETURN:   a standard basis of I
NOTE:     The procedure computes a standard basis of I (over the rational
          numbers) by using modular methods.
SEE ALSO: modular
EXAMPLE:  example modfrWalk; shows an example"
{
    /* save options */
    intvec opt = option(get);
    option(redSB);

    /* call modular() */
    if (size(#) > 0) {
        I = modular("frandwalk", list(I,radius,#), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }
    else {
        I = modular("frandwalk", list(I,radius), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }

    /* return the result */
    attrib(I, "isSB", 1);
    option(set, opt);
    return(I);
}
example
{
    "EXAMPLE:";
    echo = 2;
    ring R1 = 0, (x,y,z,t), dp;
    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
    I = std(I);
    ring R2 = 0, (x,y,z,t), lp;
    ideal I = fetch(R1, I);
    int radius = 2;
    ideal J = modfrWalk(I,radius);
    J;
    ring S1 = 0, (a,b,c,d), Dp;
    ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c;
    I = std(I);
    ring S2 = 0, (c,d,b,a), lp;
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a lp-Groebner basis.
    ideal J = modfrWalk(I,radius,"Dp");
    J;
    intvec w = 3,2,1,2;
    ring S3 = 0, (c,d,b,a), (a(w),lp);
    ideal I = fetch(S1,I);
    // I is assumed to be a Dp-Groebner basis.
    // We compute a (a(w),lp)-Groebner basis.
    ideal J = modfrWalk(I,radius,"Dp",w);
    J;
}

/* test if the prime p is suitable for the input, i.e. it does not divide
 * the numerator or denominator of any of the coefficients */
static proc primeTest_std(int p, alias list args)
{
    /* erase zero generators */
    ideal I = simplify(args[1], 2);

    /* clear denominators and count the terms */
    ideal J;
    ideal K;
    int n = ncols(I);
    intvec sizes;
    number cnt;
    int i;
    for(i = n; i > 0; i--) {
        J[i] = cleardenom(I[i]);
        cnt = leadcoef(J[i])/leadcoef(I[i]);
        K[i] = numerator(cnt)*var(1)+denominator(cnt);
    }
    sizes = size(J[1..n]);

    /* change to characteristic p */
    def br = basering;
    list lbr = ringlist(br);
    if (typeof(lbr[1]) == "int") {
        lbr[1] = p;
    }
    else {
        lbr[1][1] = p;
    }
    def rp = ring(lbr);
    setring(rp);
    ideal Jp = fetch(br, J);
    ideal Kp = fetch(br, K);

    /* test if any coefficient is missing */
    if (intvec(size(Kp[1..n])) != 2:n) {
        setring(br);
        return(0);
    }
    if (intvec(size(Jp[1..n])) != sizes) {
        setring(br);
        return(0);
    }
    setring(br);
    return(1);
}

/* find entries in modresults which come from unlucky primes.
 * For this, sort the entries into categories depending on their leading
 * ideal and return the indices in all but the biggest category. */
static proc deleteUnluckyPrimes_std(alias list modresults)
{
    int size_modresults = size(modresults);

    /* sort results into categories.
     * each category is represented by three entries:
     * - the corresponding leading ideal
     * - the number of elements
     * - the indices of the elements
     */
    list cat;
    int size_cat;
    ideal L;
    int i;
    int j;
    for (i = 1; i <= size_modresults; i++) {
        L = lead(modresults[i]);
        attrib(L, "isSB", 1);
        for (j = 1; j <= size_cat; j++) {
            if (size(L) == size(cat[j][1])
                && size(reduce(L, cat[j][1])) == 0
                && size(reduce(cat[j][1], L)) == 0) {
                cat[j][2] = cat[j][2]+1;
                cat[j][3][cat[j][2]] = i;
                break;
            }
        }
        if (j > size_cat) {
            size_cat++;
            cat[size_cat] = list();
            cat[size_cat][1] = L;
            cat[size_cat][2] = 1;
            cat[size_cat][3] = list(i);
        }
    }

    /* find the biggest categories */
    int cat_max = 1;
    int max = cat[1][2];
    for (i = 2; i <= size_cat; i++) {
        if (cat[i][2] > max) {
            cat_max = i;
            max = cat[i][2];
        }
    }

    /* return all other indices */
    list unluckyIndices;
    for (i = 1; i <= size_cat; i++) {
        if (i != cat_max) {
            unluckyIndices = unluckyIndices + cat[i][3];
        }
    }
    return(unluckyIndices);
}

/* test if 'command' applied to 'args' in characteristic p is the same as
   'result' mapped to characteristic p */
static proc pTest_std(string command, alias list args, alias ideal result,
    int p)
{
    /* change to characteristic p */
    def br = basering;
    list lbr = ringlist(br);
    if (typeof(lbr[1]) == "int") {
        lbr[1] = p;
    }
    else {
        lbr[1][1] = p;
    }
    def rp = ring(lbr);
    setring(rp);
    ideal Ip = fetch(br, args)[1];
    list Arg = fetch(br, args);
    string exstr;
    ideal Gp = fetch(br, result);
    attrib(Gp, "isSB", 1);

    /* test if Ip is in Gp */
    int i;
    for (i = ncols(Ip); i > 0; i--) {
        if (reduce(Ip[i], Gp, 1) != 0) {
            setring(br);
            return(0);
        }
    }

    /* compute command(args) */
    exstr = "Ip = "+command+" (Ip";

    for(i=2; i<=size(Arg); i++) {
      exstr = exstr+",Arg["+string(eval(i))+"]";
      }
    exstr = exstr+");";

    execute(exstr);

    /* test if Gp is in Ip */
    for (i = ncols(Gp); i > 0; i--) {
        if (reduce(Gp[i], Ip, 1) != 0) {
            setring(br);
            return(0);
        }
    }

    setring(br);
    return(1);
}

/* test if 'result' is a GB of the input ideal */
static proc finalTest_std(string command, alias list args, ideal result)
{
    /* test if args[1] is in result */
    attrib(result, "isSB", 1);
    int i;
    for (i = ncols(args[1]); i > 0; i--) {
        if (reduce(args[1][i], result, 1) != 0) {
            return(0);
        }
    }

   /* test if result is in args[1].                      */
   /* args[1] is given by a Groebner basis. Thus we may  */
   /* reduce the result with respect to args[1].         */
   int n=nvars(basering);
   string ord_str = "dp";

   for(i=2; i<=size(args); i++)
   {
     if(typeof(args[i]) == "list") {
       if(typeof(args[i][1]) == "intvec") {
         ord_str = "(a("+string(args[i][1])+"),lp("+string(n) + "),C)";
         break;
       }
       if(typeof(args[i][1]) == "string") {
         if(args[i][1] == "Dp") {
           ord_str = "Dp";
         }
         break;
       }
     }
   }
   ideal xI = args[1];
   ring xR = basering;
   execute("ring yR = ("+charstr(xR)+"),("+varstr(xR)+"),"+ord_str+";");
   ideal yI = fetch(xR,xI);
   ideal yresult = fetch(xR,result);
   attrib(yI, "isSB", 1);
   for(i=size(yresult); i>0; i--)
   {
     if(reduce(yresult[i],yI) != 0)
     {
       return(0);
     }
   }
   setring xR;
   kill yR;

   /* test if result is a Groebner basis */
    link l1="ssi:fork";
    open(l1);
    link l2="ssi:fork";
    open(l2);
    list l=list(l1,l2);
    write(l1,quote(TestSBred(result)));
    write(l2,quote(TestSBstd(result)));
    i=waitfirst(l);
    if(i==1) {
      i=read(l1);
      }
    else {
      i=read(l2);
      }
    close(l1);
    close(l2);
    return(i);
}

/* return 1, if I_reduce is _not_ in G_reduce,
 *        0, otherwise
 * (same as size(reduce(I_reduce, G_reduce))).
 * Uses parallelization. */
static proc reduce_parallel(def I_reduce, def G_reduce)
{
    exportto(Modwalk, I_reduce);
    exportto(Modwalk, G_reduce);
    int size_I = ncols(I_reduce);
    int chunks = Modular::par_range(size_I);
    intvec range;
    int i;
    for (i = chunks; i > 0; i--) {
        range = Modular::par_range(size_I, i);
        task t(i) = "Modwalk::reduce_task", list(range);
    }
    startTasks(t(1..chunks));
    waitAllTasks(t(1..chunks));
    int result = 0;
    for (i = chunks; i > 0; i--) {
        if (getResult(t(i))) {
            result = 1;
            break;
        }
    }
    kill I_reduce;
    kill G_reduce;
    return(result);
}

/* compute a chunk of reductions for reduce_parallel */
static proc reduce_task(intvec range)
{
    int result = 0;
    int i;
    for (i = range[1]; i <= range[2]; i++) {
        if (reduce(I_reduce[i], G_reduce, 1) != 0) {
            result = 1;
            break;
        }
    }
    return(result);
}

/* test if result is a GB with std*/
static proc TestSBstd(ideal result)
{
  ideal G = std(result);
  if(reduce_parallel(G,result)) {
    return(0);
    }
  return(1);
}

/* test if result is a GB by reducing s-polynomials*/
static proc TestSBred(ideal result)
{
  int i,j;
  for(i=1; i<=size(result); i++) {
    for(j=i; j<=size(result); j++) {
      if(reduce(sPolynomial(result[i],result[j]),result)!=0) {
        return(0);
        }
      }
    }
  return(1);
}

/* compute s-polynomial of f and g */
static proc sPolynomial(poly f,poly g)
{
  int i;
  poly lcmp = 1;

  intvec lexpf = leadexp(f);
  intvec lexpg = leadexp(g);

  for(i=1; i<=nvars(basering); i++) {
    if(lexpf[i]>=lexpg[i]) {
      lcmp=lcmp*var(i)**lexpf[i];
      }
    else {
      lcmp=lcmp*var(i)**lexpg[i];
      }
   }

  poly fmult=lcmp/leadmonom(f);
  poly gmult=lcmp/leadmonom(g);
  poly result=leadcoef(g)*fmult*f-leadcoef(f)*gmult*g;

  return(result);
}