modquotient.lib

////////////////////////////////////////////////////////////////////////////////
version=" ";
category="Commutative Algebra";
info=" ";

LIB "modular.lib";
LIB "modstd.lib";

proc modQuotient(ideal I, ideal J, list #)
{
    /* read optional parameter */
    int exactness = 1;
    if (size(#) > 0) {
        if (size(#) > 1 || typeof(#[1]) != "int") {
            ERROR("wrong optional parameter");
        }
        exactness = #[1];
    }

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

    /* compute standard basis of I */
    if (!attrib(I, "isSB")) {
        I = modStd(I);
    }

    /* set maximal prime number */
    int pmax = 2147483647;
    if (npars(basering) > 0) {
        pmax = 536870909;   // prime(2^29)
    }

    /* choose the right command */
    string command = "quotient";
    if (npars(basering) > 0) {
        command = "Modquotient::quotient_norm";
    }

    /* call modular() */
    ideal Q;
    if (exactness) {
        Q = modular(command, list(I, J), primeTest_quotient,
            Modstd::deleteUnluckyPrimes_std, pTest_quotient,
            finalTest_quotient, pmax);
    }
    else {
        Q = modular(command, list(I, J), primeTest_quotient,
            Modstd::deleteUnluckyPrimes_std, pTest_quotient, pmax);
    }

    /* return the result */
    attrib(Q, "isSB", 1);
    option(set, opt);
    return(Q);
}

static proc quotient_norm(ideal I, ideal J)
{
    ideal Q = simplify(quotient(I, J), 1);
    attrib(Q, "isSB", 1);
    return(Q);
}

static proc primeTest_quotient(int p, alias list args)
{
    if (!Modstd::primeTest_std(p, list(args[1]))) {
        return(0);
    }
    return(Modstd::primeTest_std(p, list(args[2])));
}

static proc pTest_quotient(string command, list args, ideal result, int 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 Jp = fetch(br, args)[2];
    ideal result = fetch(br, result);
    attrib(result, "isSB", 1);
    int i;
    execute("ideal Qp = "+command+"(Ip, Jp);");
    for (i = ncols(Qp); i > 0; i--) {
        if (reduce(Qp[i], result, 1) != 0) {
            setring(br);
            return(0);
        }
    }
    for (i = ncols(result); i > 0; i--) {
        if (reduce(result[i], Qp, 1) != 0) {
            setring(br);
            return(0);
        }
    }
    setring(br);
    return(1);
}

static proc finalTest_quotient(string command, list args, ideal result)
{
    ideal P = result*args[2];
    int i;
    for (i = ncols(P); i > 0; i--) {
        if (reduce(P[i], args[1], 1) != 0) {
            return(0);
        }
    }
    attrib(result, "isSB", 1);
    int res = quotient_test(args[1], args[2], result);
    return(res);
}

/* Tests if Q = (I:J). Assumes Q \subset (I:J) and that Q is a standard basis.
 */
static proc quotient_test(ideal I, ideal J, ideal Q)
{
    int i;
    int ncols_J = ncols(J);

    def br = basering;
    list rl = ringlist(br);
    for (i = size(rl[3]); i > 0; i--) {
        if (rl[3][i][1] == "c" || rl[3][i][1] == "C") {
            rl[3] = delete(rl[3], i);
        }
    }
    rl[3] = insert(rl[3], list("c", 0));
    def newR = ring(rl);
    setring(newR);
    option(noredThrough);   // this option is ring-dependent
    ideal I = fetch(br, I);
    ideal J = fetch(br, J);
    ideal Q = fetch(br, Q);
    attrib(Q, "isSB", 1);

    module M = module(transpose(J));
    M[1] = M[1] + gen(ncols_J+1);
    for (i = 1; i <= ncols_J; i++) {
        M = module(I*gen(i)) + M;
    }
    M = module(Q*gen(ncols_J+1)) + M;

    M = std(M);

    ideal QQ;
    for (i = ncols(M); i > 0; i--) {
        if (M[i][1..ncols_J] == 0) {
            QQ = M[i][ncols_J+1], QQ;
        }
    }
    QQ = simplify(QQ, 2);
    for (i = ncols(QQ); i > 0; i--) {
        if (reduce(QQ[i], Q, 1) != 0) {
            setring(br);
            return(0);
        }
    }
    setring(br);
    return(1);
}

proc modSat(def id, ideal j, int exactness)
"USAGE:  sat(id, j);  id ideal/module, j ideal
RETURN:  list of an ideal/module [1] and an integer [2]:
         [1] = saturation of id with respect to j (= union_(k=1...) of id:j^k)
         [2] = saturation exponent (= min( k | id:j^k = id:j^(k+1) ))
NOTE:    [1] is a standard basis in the basering
{
    int ii, kk;
    def i = id;
    id = modStd(id, exactness);
    while (ii <= size(i)) {
        i = modQuotient(id, j, exactness);
        for (ii = 1; ii <= size(i); ii++) {
            if (reduce(i[ii], id, 1) != 0) break;
        }
        id = modStd(i, exactness);
        kk++;
    }
    list L = id, kk-1;
    return(L);
}