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

///////////////////////////////////////////////////////////////////////////////
version="version modstd.lib 4.0.0.0 May_2014 "; // $Id$
category="Commutative Algebra";
info="
LIBRARY:  modstd.lib      Groebner bases of ideals using modular methods

AUTHORS:  A. Hashemi      Amir.Hashemi@lip6.fr
          G. Pfister      pfister@mathematik.uni-kl.de
          H. Schoenemann  hannes@mathematik.uni-kl.de
          A. Steenpass    steenpass@mathematik.uni-kl.de
          S. Steidel      steidel@mathematik.uni-kl.de

OVERVIEW:
  A library for computing Groebner bases of ideals in the polynomial ring over
  the rational numbers using modular methods.

REFERENCES:
  E. A. Arnold: Modular algorithms for computing Groebner bases.
  J. Symb. Comp. 35, 403-419 (2003).

  N. Idrees, G. Pfister, S. Steidel: Parallelization of Modular Algorithms.
  J. Symb. Comp. 46, 672-684 (2011).

PROCEDURES:
  modStd(I);    standard basis of I using modular methods
";

LIB "poly.lib";
LIB "modular.lib";

proc modStd(ideal I, list #)
"USAGE:   modStd(I[, exactness]); I ideal, exactness int
RETURN:   a standard basis of I
NOTE:     The procedure computes a standard basis of I (over the rational
          numbers) by using modular methods.
       @* An optional parameter 'exactness' can be provided.
          If exactness = 1, the procedure computes a standard basis of I for
          sure; if exactness = 0, it computes a standard basis of I
          with high probability.
SEE ALSO: modular
EXAMPLE:  example modStd; shows an example"
{
    /* read optional parameter */
    int exactness = 1;
    if (size(#) > 0) {
        /* For compatibility, we only test size(#) > 4. This can be changed to
         * size(#) > 1 in the future. */
        if (size(#) > 4 || typeof(#[1]) != "int") {
            ERROR("wrong optional parameter");
        }
        exactness = #[1];
    }

    /* save options */
    intvec opt = option(get);
    option(redSB);

    /* choose the right command */
    string command = "groebner";
    if (npars(basering) > 0) {
        command = "Modstd::groebner_norm";
    }

    /* call modular() */
    if (exactness) {
        I = modular(command, list(I), primeTest_std,
            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
    }
    else {
        I = modular(command, list(I), primeTest_std,
            deleteUnluckyPrimes_std, pTest_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;
    ideal J = modStd(I);
    J;
    I = homog(I, t);
    J = modStd(I);
    J;

    ring R2 = 0, (x,y,z), ds;
    ideal I = jacob(x5+y6+z7+xyz);
    ideal J = modStd(I, 0);
    J;

    ring R3 = 0, x(1..4), lp;
    ideal I = cyclic(4);
    ideal J1 = modStd(I, 1);   // default
    ideal J2 = modStd(I, 0);
    size(reduce(J1, J2));
    size(reduce(J2, J1));
}

/* compute a normalized GB via groebner() */
static proc groebner_norm(ideal I)
{
    I = simplify(groebner(I), 1);
    attrib(I, "isSB", 1);
    return(I);
}

/* 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, list args, 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];
    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) */
    execute("Ip = "+command+"(Ip);");

    /* 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 a GB */
    ideal G = std(result);
    if (reduce_parallel(G, result)) {
        return(0);
    }
    return(1);
}

/* 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(ideal I_reduce, ideal G_reduce)
{
    exportto(Modstd, I_reduce);
    exportto(Modstd, 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) = "Modstd::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);
}


////////////////////////////// further examples ////////////////////////////////

/*
ring r = 0, (x,y,z), lp;
poly s1 = 5x3y2z+3y3x2z+7xy2z2;
poly s2 = 3xy2z2+x5+11y2z2;
poly s3 = 4xyz+7x3+12y3+1;
poly s4 = 3x3-4y3+yz2;
ideal i =  s1, s2, s3, s4;

ring r = 0, (x,y,z), lp;
poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
ideal i =  s1, s2, s3, s4;

ring r = 0, (x,y,z), lp;
poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
poly s3 = 8x3 + 12y3 + xz2 + 3;
poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
ideal i = s1, s2, s3, s4;

int n = 6;
ring r = 0,(x(1..n)),lp;
ideal i = cyclic(n);
ring s = 0, (x(1..n),t), lp;
ideal i = imap(r,i);
i = homog(i,t);

ring r = 0, (x(1..4),s), (dp(4),dp);
poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
poly 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);
poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
ideal i =  s1, s2, s3, s4, s5;

ring r = 0, (x,y,z), ds;
int a = 16;
int b = 15;
int c = 4;
int t = 1;
poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
         +x^(c-2)*y^c*(y2+t*x)^2;
ideal i = jacob(f);

ring r = 0, (x,y,z), ds;
int a = 25;
int b = 25;
int c = 5;
int t = 1;
poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
         +x^(c-2)*y^c*(y2+t*x)^2;
ideal i = jacob(f),f;

ring r = 0, (x,y,z), ds;
int a = 10;
poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
ideal i = jacob(f);

ring r = 0, (x,y,z), ds;
int a = 6;
int b = 8;
int c = 10;
int alpha = 5;
int beta = 5;
int t = 1;
poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
ideal i = jacob(f);
*/