Parse negative numbers better

[libyakmo.git] / groebner.myr

use std

use t

use "traits"
use "Z"
use "Q"
use "polyring"

/*
   Algorithms for dealing with Gröbner bases. Everything is from Cox,
   Little, O'Shea "Ideals, Varieties, and Algorithms" (4th ed.), which
   is optimized for being readable, not for speed.
 */

pkg yakmo =
        const groebner_basis_of : (F : polynomial#[:] -> polynomial#[:])

        const positive_monomial_division : (a : monomial#, b : monomial# -> std.option(monomial))

        impl lcm_struct polynomial#
;;

/*
   As of this moment, I don't want to encourage using other monomial
   orderings. That's strictly a hack for lcm calculation via
   elimination.

   Also, this really wouldn't work with negative exponents.
 */
const lex_order = {a : monomial, b : monomial
        var ak = a.vars.len - 1
        var bk = b.vars.len - 1

        while ak >= 0 && bk >= 0
                match strnumcmp(a.vars[ak].variable, b.vars[bk].variable)
                | `std.Before: -> `std.Before
                | `std.After : -> `std.After
                | `std.Equal:
                        match t.cmp(a.vars[ak].power, b.vars[bk].power)
                        | `std.Before: -> `std.Before
                        | `std.After : -> `std.After
                        | `std.Equal:
                        ;;
                ;;

                ak--
                bk--
        ;;

        if ak >= 0
                -> `std.After
        ;;

        if bk >= 0
                -> `std.Before
        ;;

        -> `std.Equal
}

impl lcm_struct polynomial# =
        lcm = {f, g
                /*
                   Early case. If f or g are just constants, return the other, or maybe 1
                 */
                var f_trivial = f.terms.len == 1 && f.terms[0].vars.len == 0
                var g_trivial = g.terms.len == 1 && g.terms[0].vars.len == 0
                if f_trivial && g_trivial
                        -> rid()
                elif f_trivial
                        var ret = t.dup(g)
                        ensure_monic(ret)
                        -> ret
                elif g_trivial
                        var ret = t.dup(f)
                        ensure_monic(ret)
                        -> ret
                ;;

                /*
                   Slightly more involved case. Evaluate f and g with a
                   bunch of random constants, then see if they have an
                   lcm. If not, then just return f*g.
                 */
                if no_chance_of_lcm(f, g)
                        var ret = t.dup(f)
                        yakmo.rmul_ip(ret, g)
                        ensure_monic(ret)
                        -> ret
                ;;

                /*
                   Build t and (1-t) to multiply f and g. In particular,
                   we need t to be large enough that the standard
                   grevlex ordering we use has the same effect as lex
                   ordering: any monomial containing t is larger than
                   any monomial not containing t.
                 */
                var fake_var : byte[:] = make_var_larger_than(f, g)
                var vp_t : var_power = [
                        .variable = fake_var,
                        .power = rid(), // gadd(auto deg(f), auto deg(g)),
                ]
                var mon_t : monomial = [
                        .coeff = rid(),
                        .vars = std.sldup([ vp_t ][:])
                ]
                var mon_one : monomial = [
                        .coeff = rid(),
                        .vars = [][:]
                ]
                var ft = t.dup(f)
                mul_poly_mon_ip(ft, &mon_t)

                gneg_ip(mon_t.coeff)
                var poly_1mt : polynomial = [
                        .terms = [ mon_one, mon_t ][:]
                ]

                var g1mt = t.dup(g)
                rmul_ip(g1mt, &poly_1mt)
                auto mon_one
                auto mon_t // This kills fake_var as well

                /*
                   Force lex ordering. This will percolate all the way
                   through the Gröbner basis calculation, since we do it
                   all with dups.
                 */
                ft.monomial_cmp = `std.Some lex_order
                g1mt.monomial_cmp = `std.Some lex_order
                reducepoly(ft)
                reducepoly(g1mt)

                /* Now, compute a Gröbner basis for ⟨ ft, g(1-t) ⟩. */
                var groebner_basis = groebner_basis_of( [ ft, g1mt ][:] )

                /*
                   The promise is that, if we remove all the elements of
                   the groebner basis that contain fake_var, we'll be
                   left with at most one element, and it will be the
                   lcm.
                 */
                for var j = 0; j < groebner_basis.len; ++j
                        if contains_var(groebner_basis[j], fake_var)
                                __dispose__(groebner_basis[j])
                                std.sldel(&groebner_basis, j)
                                j--
                        ;;
                ;;
                if groebner_basis.len == 1
                        var ret = groebner_basis[0]
                        std.slfree(groebner_basis)
                        ret.monomial_cmp = `std.None
                        reducepoly(ret)
                        ensure_monic(ret)
                        -> ret
                elif groebner_basis.len == 0
                        std.slfree(groebner_basis)
                        var ret = rmul(f,g)
                        ensure_monic(ret)
                        -> ret
                ;;

                /* This should be impossible */
                var ret = groebner_basis[0]
                for var j = 1; j < groebner_basis.len; ++j
                        __dispose__(groebner_basis[j])
                ;;
                std.slfree(groebner_basis)
                ret.monomial_cmp = `std.None
                reducepoly(ret)
                ensure_monic(ret)
                -> ret
        }
;;

const no_chance_of_lcm = {f, g
        /* If f or g aren't entirely positive, the following fails */
        for var j = 0; j < f.terms.len; ++j
                if f.terms[j].coeff.p.sign < 0
                        -> false
                ;;
        ;;
        for var j = 0; j < g.terms.len; ++j
                if g.terms[j].coeff.p.sign < 0
                        -> false
                ;;
        ;;

        /*
           Now, assuming f and g are totally positive, suppose they factor into

               f = p1^e1 ... pn^en
               g = q1^d1 ... qm^dm

           All pi and qi are also positive, so evaluated on numbers
           (strictly greater than 1), they'll return results >= 2.
           Therefore, if f and g share a common factor pi = qj,
           evaluating f and g will result in numbers with a common
           factor >= 2.

           So, if we evaluate f and g on some random numbers and they
           don't share any common factors numerically, they can't share
           any abstractly.

           The number 5 was chosen after about three minutes of testing
           with one specific program. It should probably be calculated
           based on number of terms and variables.
         */
        for var j = 0; j < 5; ++j
                var rng : std.rng# = std.mksrng(j)
                var int_from_var : std.htab(byte[:], int)# = std.mkht()

                var fq : yakmo.Q# = eval_polynomial_by_seed(f, rng, int_from_var)
                var gq : yakmo.Q# = eval_polynomial_by_seed(g, rng, int_from_var)
                auto fq
                auto gq
                std.freerng(rng)

                var fzb : std.bigint# = std.bigdup(fq.p)
                std.bigmul(fzb, gq.q)
                var fz : yakmo.Z# = (fzb : yakmo.Z#)
                var gzb : std.bigint# = std.bigdup(gq.p)
                std.bigmul(gzb, fq.q)
                var gz : yakmo.Z# = (gzb : yakmo.Z#)
                auto fz
                auto gz

                /* cleanup */
                for (k, _) : std.byhtkeyvals(int_from_var)
                        std.slfree(k)
                ;;

                std.htfree(int_from_var)

                var fgz = yakmo.rmul(fz, gz)
                var fglcm = yakmo.lcm(fz, gz)
                auto fgz
                auto fglcm

                if std.eq(fgz, fglcm)
                        -> true
                ;;
        ;;

        -> false
}

/*
   TODO: all of this needs to be subsumed in some sort of "evaluable"
   trait. But traits are borked right now.
 */
const eval_polynomial_by_seed = { f : yakmo.polynomial#, rng : std.rng#, vars : std.htab(byte[:], int)#
        var ret : yakmo.Q# = yakmo.gid()

        for var j = 0; j < f.terms.len; ++j
                var temp : yakmo.Q# = eval_monomial_by_seed(f.terms[j], rng, vars)
                auto temp
                yakmo.gadd_ip(ret, temp)
        ;;

        -> ret
}

const eval_monomial_by_seed = { m : yakmo.monomial, rng : std.rng#, vars : std.htab(byte[:], int)#
        var ret : yakmo.Q# = t.dup(m.coeff)

        for var j = 0; j < m.vars.len; ++j
                var temp : yakmo.Q# = eval_var_power_by_seed(m.vars[j], rng, vars)
                auto temp
                yakmo.rmul_ip(ret, temp)
        ;;

        -> ret
}

const primes : int[:] = [
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
        61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
        131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193,
        197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269,
        271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
        353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431,
        433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499
][:]

const eval_var_power_by_seed = {vp : yakmo.var_power, rng : std.rng#, vars : std.htab(byte[:], int)#
        var var_name : byte[:] = vp.variable
        var var_value : int = 0
        var zret : std.bigint# = std.mkbigint(1)

        match std.htget(vars, var_name)
        | `std.None:
                var_value = primes[std.rngrand(rng, 0, primes.len)]
                std.htput(vars, std.sldup(var_name), var_value)
        | `std.Some z: var_value = z
        ;;

        var np : std.bigint# = std.bigdup((vp.power : std.bigint#))
        var double = std.mkbigint(var_value)
        while !std.bigiszero(np)
                if !std.bigiseven(np)
                        std.bigmul(zret, double)
                ;;

                std.bigmul(double, double)
                std.bigshri(np, 1)
        ;;

        std.bigfree(double)
        std.bigfree(np)
        var qret : yakmo.Q = [ .p = zret, .q = std.mkbigint(1) ]

        -> std.mk(qret)
}

const make_var_larger_than = {a, b
        var c : char = 'x'
        var n : int = 1
        for var j = 0; j < a.terms.len; ++j
                for var k = 0; k < a.terms[j].vars.len; ++k
                        var ct, st
                        (ct, st) = std.charstep(a.terms[j].vars[k].variable)
                        c = std.max(c, std.max(ct, ct + 1))
                        while c == ct
                                n++
                                (ct, st) = std.charstep(st)
                        ;;
                ;;
        ;;

        for var j = 0; j < b.terms.len; ++j
                for var k = 0; k < b.terms[j].vars.len; ++k
                        var ct, st
                        (ct, st) = std.charstep(b.terms[j].vars[k].variable)
                        c = std.max(c, std.max(ct, ct + 1))
                        while c == ct
                                n++
                                (ct, st) = std.charstep(st)
                        ;;
                ;;
        ;;

        var sb : std.strbuf# = std.mksb()
        while n > 0
                std.sbputc(sb, c)
                n--
        ;;

        -> std.sbfin(sb)
}

/*
   Return whether p contains a variable with given name, given the
   assumption that such a variable is the largest variable occuring in
   p, and that it only occurs with powers large enough to place it in
   the leading term.
 */
const contains_var = {p : polynomial#, name : byte[:]
        if p.terms.len == 0
                -> false
        ;;

        var mon : yakmo.monomial# = &p.terms[p.terms.len - 1]
        if mon.vars.len == 0
                -> false
        ;;

        -> std.eq(mon.vars[mon.vars.len - 1].variable, name)
}

/*
   Buchberger's Algorithm, specifically as described in [CLO] Chapter 2,
   Section 9, Theorem 9, with the incorporated reduction definition.
 */
const groebner_basis_of = {F : polynomial#[:]
        var G : polynomial#[:] = t.dup(F)

        /*
           Initial reduction step. Kill every monomial of any polynomial
           that appears as the leading term of any other polynomial.
           Repeat until done.
         */
        reduce_ideal(&G, false)
        remove_trivial(&G)

        /* Now run Buchberger's */
        var test_these : (std.size, std.size)[:] = [][:]
        for var j = 0; j < G.len; ++j
                for var k = 0; k < j; ++k
                        std.slpush(&test_these, (j, k))
                ;;
        ;;

        while test_these.len != 0
                var j, k
                (j, k) = std.slpop(&test_these)
                if !common_factors_in_LT(G[j], G[k])
                        continue
                ;;

                if proposition_8_criterion(G, test_these, j, k)
                        continue
                ;;

                var r = S_bar(G, G[j], G[k])
                if eq_gid(r)
                        __dispose__(r)
                else
                        std.slpush(&G, r)
                        reduce_ideal(&G, true)
                        for var i = 0; i < G.len - 1; ++i
                                std.slpush(&test_these, (i, G.len - 1))
                        ;;
                ;;
        ;;
        std.slfree(test_these)

        /*
           We now have a Gröbner basis. Make it slightly prettier by
           making all entries monic.
         */
        for var j = 0; j < G.len; ++j
                ensure_monic(G[j])
        ;;
        remove_trivial(&G)
        std.sort(G, t.cmp)

        /* Now, remove redundant entries */
        /*
           TODO: I'm pretty sure this is handled by the reduction step
           during the algorithm, but I'm too tired to check.
         */
        minimize_groebner_basis(&G)

        -> G
}

const ensure_monic = {p : polynomial#
        if p.terms.len == 0
                /* How did this get in here? */
                -> void
        ;;

        if eqQ(p.terms[p.terms.len - 1].coeff, 1, 1)
                -> void
        ;;

        match finv(p.terms[p.terms.len - 1].coeff)
        | `std.None:
        | `std.Some ci:
                for var k = 0; k < p.terms.len; ++k
                        rmul_ip(p.terms[k].coeff, ci)
                ;;
                __dispose__(ci)
        ;;
}

const reduce_ideal = {ps : polynomial#[:]#, only_check_end : bool
        /*
           Start at the end, because in the process of Buchberger's,
           that's where the new, interesting polynomials get slotted.
         */
        var test_these = [][:]
        if only_check_end
                std.slpush(&test_these, ps#.len - 1)
        else
                for var k = 0; k < ps#.len - 1; ++k
                        std.slpush(&test_these, k)
                ;;
        ;;
        while test_these.len != 0
                var j = std.slpop(&test_these)

                var p = ps#[j]

                if p.terms.len == 0
                        continue
                ;;

                /* Apply everything to p */
                for var k = 0; k < ps#.len; ++k
                        if k == j
                                continue
                        ;;

                        var q = ps#[k]

                        if q.terms.len == 0
                                continue
                        ;;
                        var ltq = &q.terms[q.terms.len - 1]

                        /*
                           We want to loop down through p, from the top.
                           If we ever change p, then we might add/remove
                           terms and change the indexing. But 1) we don't
                           want to check the top terms over and over,
                           and 2) as soon as we skip through a term, we
                           know it will never be affected by any future
                           action.
                         */
                        var good_terms = 1
                        for var l = p.terms.len - 1; l >= 0; --l
                                match positive_monomial_division(&p.terms[l], ltq)
                                | `std.None: good_terms++
                                | `std.Some r:
                                        auto r
                                        var summand = t.dup(q)
                                        mul_poly_mon_ip(summand, &r)
                                        gneg_ip(summand)
                                        auto summand
                                        gadd_ip(p, summand)
                                        l = p.terms.len - 1 // good_terms
                                        k = -1
                                ;;
                        ;;
                ;;

                if p.terms.len == 0
                        continue
                ;;

                /* Now apply p to everything else */
                var ltp : monomial# = &p.terms[p.terms.len - 1]
                for var k = 0; k < ps#.len; ++k
                        if k == j
                                continue
                        ;;

                        var q = ps#[k]

                        if q.terms.len == 0
                                continue
                        ;;

                        var good_terms = 1
                        for var l = q.terms.len - 1; l >= 0; --l
                                match positive_monomial_division(&q.terms[l], ltp)
                                | `std.None:
                                | `std.Some r:
                                        auto r
                                        var summand = t.dup(p)
                                        mul_poly_mon_ip(summand, &r)
                                        gneg_ip(summand)
                                        auto summand
                                        gadd_ip(q, summand)
                                        l = q.terms.len - 1 // good_terms
                                        ensure_contains(&test_these, k)
                                        k = -1
                                ;;
                        ;;
                ;;
        ;;
        std.slfree(test_these)
}

const ensure_contains = {s, e
        if e < 0
                /*
                   This can only happen if we're in the second loop for
                   a particular q, and in that case, the index for q has
                   already been added to the list.
                 */
                -> void
        ;;
        for var k = 0; k < s#.len; ++k
                if s#[k] == e
                        -> void
                ;;
        ;;
        std.slpush(s, e)
}

/* Weed out zeroes */
const remove_trivial = {ps : polynomial#[:]#
        for var j = 0; j < ps#.len; ++j
                if ps#[j].terms.len == 0
                        std.sldel(ps, j)
                        j--
                ;;
        ;;
}


/* True if LCM(LT(f), LT(g)) != LT(f) LT(g) */
const common_factors_in_LT = {f : polynomial#, g : polynomial#
        if f.terms.len == 0 || g.terms.len == 0
                -> false
        ;;
        var ltf = &f.terms[f.terms.len - 1]
        var ltg = &f.terms[f.terms.len - 1]

        for var j = 0; j < ltf.vars.len; ++j
                for var k = 0; k < ltg.vars.len; ++k
                        match strnumcmp(ltg.vars[k].variable, ltf.vars[j].variable)
                        | `std.Before:
                        | `std.Equal: -> true
                        | `std.After: break
                        ;;
                ;;
        ;;

        -> false
}

/*
   True if there is some l such that

    - l \notin { j, k }

    - both (j, l) and (k, l) are not in B (up to ordering)

    - LT(f_l) divides lcm(LT(f_j), LT(f_k))
 */
const proposition_8_criterion = {G : polynomial#[:], B : (std.size, std.size)[:], j : std.size, k : std.size
        var fj = G[j]
        var fk = G[k]
        if fj.terms.len == 0 || fk.terms.len == 0
                -> true
        ;;

        var ltj = &fj.terms[fj.terms.len - 1]
        var ltk = &fk.terms[fk.terms.len - 1]

        for var l = 0; l < G.len; ++l
                if l == j || l == k
                        continue
                ;;

                var already_in_b = false
                var jl1 = std.min(j, l)
                var jl2 = std.max(j, l)
                var kl1 = std.min(k, l)
                var kl2 = std.max(k, l)
                for (s, t) : B
                        if (s == jl1 && t == jl2) || (s == kl1 && t == kl2)
                                already_in_b = true
                                break
                        ;;
                ;;

                if already_in_b
                        continue
                ;;

                var fl = G[l]
                if fl.terms.len == 0
                        continue
                ;;

                var ltl = &fl.terms[fl.terms.len - 1]
                var this_l_works = true

                for vpl : ltl.vars
                        var this_var_in_lcm = false
                        for vpj : ltj.vars
                                match strnumcmp(vpj.variable, vpl.variable)
                                | `std.Before:
                                | `std.Equal:
                                        match t.cmp(vpj.power, vpl.power)
                                        | `std.Before:
                                        | _:
                                                this_var_in_lcm = true
                                                break
                                        ;;
                                | `std.After: break
                                ;;
                        ;;

                        if this_var_in_lcm
                                continue
                        ;;

                        for vpk : ltk.vars
                                match strnumcmp(vpk.variable, vpl.variable)
                                | `std.Before:
                                | `std.Equal:
                                        match t.cmp(vpk.power, vpl.power)
                                        | `std.Before:
                                        | _:
                                                this_var_in_lcm = true
                                                break
                                        ;;
                                | `std.After: break
                                ;;
                        ;;

                        if !this_var_in_lcm
                                this_l_works = false
                                break
                        ;;
                ;;

                if this_l_works
                        -> true
                ;;
        ;;

        -> false
}

/*
   Compute \overbar{S(f, g)}^{F}; the remainder of S(f, g) when reducing
   by everything in F.

   TODO: maybe break this out into a general polynomial div_maybe or div
   or reduce or something?
 */
const S_bar = {F, f, g
        var Sfg : polynomial# = S_polynomial(f, g)

        var again = true
        while again
                again = false
                for var j = 0; j < F.len; ++j
                        if F[j].terms.len == 0
                                continue
                        ;;

                        if Sfg.terms.len == 0
                                -> Sfg
                        ;;

                        match positive_monomial_division(&Sfg.terms[Sfg.terms.len - 1], &F[j].terms[F[j].terms.len - 1])
                        | `std.None:
                        | `std.Some q:
                                auto q
                                var summand = t.dup(F[j])
                                mul_poly_mon_ip(summand, &q)
                                gneg_ip(summand)
                                auto summand
                                gadd_ip(Sfg, summand)
                                again = true
                        ;;
                ;;
        ;;

        -> Sfg
}

/*
   Let X be the lcm of the leading terms of f and g (ignoring coefficients). Then

     S(f, g) := X/LT(f) * f - X/LT(g) * g
 */
const S_polynomial = {f : polynomial#, g : polynomial#
        var ltf : monomial# = &f.terms[f.terms.len - 1]
        var ltg : monomial# = &g.terms[g.terms.len - 1]
        var f_factor : monomial = [ .coeff = rid(), .vars = [][:] ]
        var g_factor : monomial = [ .coeff = rid(), .vars = [][:] ]

        /*
           First, construct X/LT(f) and X/LT(g). ltf.vars and ltg.vars
           are ordered, so we may walk through LT(f) and LT(g). A power
           of a variable which we can prove occurs only in LT(f) should
           be placed in X/LT(g), and vice-versa.
         */
        var fj = 0
        var gj = 0
        while fj < ltf.vars.len && gj < ltg.vars.len
                var fx = &ltf.vars[fj]
                var gx = &ltg.vars[gj]
                match strnumcmp(fx.variable, gx.variable)
                | `std.Before:
                        std.slpush(&g_factor.vars, t.dup(fx#))
                        fj++
                | `std.After:
                        std.slpush(&f_factor.vars, t.dup(gx#))
                        gj++
                | `std.Equal:
                        /* x^a and x^b. Add x^|a-b| to one side or the other */
                        var ngp : Z# = gneg(gx.power)
                        var p = gadd(fx.power, ngp)
                        __dispose__(ngp)
                        match cmp_zero(p)
                        | `std.Equal: __dispose__(p)
                        | `std.Before:
                                gneg_ip(p)
                                std.slpush(&f_factor.vars, [
                                        .variable = std.sldup(fx.variable),
                                        .power = p
                                ])
                        | `std.After:
                                std.slpush(&g_factor.vars, [
                                        .variable = std.sldup(gx.variable),
                                        .power = p
                                ])
                        ;;
                        fj++
                        gj++
                ;;
        ;;

        /* Clean up any extra, trailing terms */
        while fj < ltf.vars.len
                std.slpush(&g_factor.vars, t.dup(ltf.vars[fj]))
                fj++
        ;;

        while gj < ltg.vars.len
                std.slpush(&f_factor.vars, t.dup(ltg.vars[gj]))
                gj++
        ;;

        /* Now construct (X / LT(f)) * f - (X / LT(g)) * g */
        var f_term = t.dup(f)
        var g_term = auto t.dup(g)
        mul_poly_mon_ip(f_term, &f_factor)
        __dispose__(f_factor)
        mul_poly_mon_ip(g_term, &g_factor)
        __dispose__(g_factor)
        gneg_ip(g_term)
        gadd_ip(f_term, g_term)

        -> f_term
}

/*
   Assuming a and b have only non-negative exponents, try and build a/b,
   similarly with non-negative exponents.

   This is not a division_struct monomial because monomials in general
   are allowed to have negative exponents.
 */
const positive_monomial_division = { a : monomial#, b : monomial#
        var quotient : monomial = [ .coeff = t.dup(a.coeff), .vars = [][:] ]

        match finv(b.coeff)
        | `std.Some bci:
                rmul_ip(quotient.coeff, bci)
                __dispose__(bci)
        | `std.None:
                __dispose__(quotient)
                -> `std.None
        ;;

        var ak = 0
        var bk = 0
        while ak < a.vars.len && bk < b.vars.len
                var ax = &a.vars[ak]
                var bx = &b.vars[bk]
                match strnumcmp(ax.variable, bx.variable)
                | `std.Before:
                        std.slpush(&quotient.vars, t.dup(ax#))
                        ak++
                | `std.After:
                        __dispose__(quotient)
                        -> `std.None
                | `std.Equal:
                        match t.cmp(ax.power, bx.power)
                        | `std.Before:
                                __dispose__(quotient)
                                -> `std.None
                        | _:
                        ;;

                        var bp : Z# = t.dup(bx.power)
                        gneg_ip(bp)
                        auto bp
                        std.slpush(&quotient.vars, [
                                .variable = std.sldup(ax.variable),
                                .power = gadd(ax.power, bp)
                        ])
                        ak++
                        bk++
                ;;
        ;;

        if bk < b.vars.len
                __dispose__(quotient)
                -> `std.None
        ;;

        while ak < a.vars.len
                std.slpush(&quotient.vars, t.dup(a.vars[ak]))
                ak++
        ;;

        reducemonomial(&quotient)
        -> `std.Some quotient
}

const minimize_groebner_basis = { G : polynomial#[:]#
        /*
           Following [CLO] Chapter 2, Section 7, Lemma 3, try to remove
           any p in G such that LT(p) is generated by the leading terms
           of G - { p }.
         */
        std.sort(G#, t.cmp)
        for var j = G#.len - 1; j >= 0; --j
                var lt = t.dup(G#[j].terms[G#[j].terms.len - 1])
                for var k = 0; k < G#.len; ++k
                        if k == j
                                continue
                        ;;

                        match positive_monomial_division(&lt, &G#[k].terms[G#[k].terms.len - 1])
                        | `std.None:
                        | `std.Some q:
                                __dispose__(q)
                                __dispose__(G#[j])
                                std.sldel(G, j)
                                break
                        ;;
                ;;
        ;;
}