def critical_pair(f, g, ring):
    """
    Compute the critical pair corresponding to two labeled polynomials.

    A critical pair is a tuple (um, f, vm, g), where um and vm are
    terms such that um * f - vm * g is the S-polynomial of f and g (so,
    wlog assume um * f > vm * g).
    For performance sake, a critical pair is represented as a tuple
    (Sign(um * f), um, f, Sign(vm * g), vm, g), since um * f creates
    a new, relatively expensive object in memory, whereas Sign(um *
    f) and um are lightweight and f (in the tuple) is a reference to
    an already existing object in memory.
    """
    domain = ring.domain

    ltf = Polyn(f).LT
    ltg = Polyn(g).LT
    lt = (monomial_lcm(ltf[0], ltg[0]), domain.one)

    um = term_div(lt, ltf, domain)
    vm = term_div(lt, ltg, domain)

    # The full information is not needed (now), so only the product
    # with the leading term is considered:
    fr = lbp_mul_term(lbp(Sign(f), Polyn(f).leading_term, Num(f)), um)
    gr = lbp_mul_term(lbp(Sign(g), Polyn(g).leading_term, Num(g)), vm)

    # return in proper order, such that the S-polynomial is just
    # u_first * f_first - u_second * f_second:
    if lbp_cmp(fr, gr) == -1:
        return (Sign(gr), vm, g, Sign(fr), um, f)
    else:
        return (Sign(fr), um, f, Sign(gr), vm, g)

    def generate(R, P, G, B):
        while R:
            h = normal(F[R.pop()], G | P)

            if h is not None:
                k, LM = h

                G0 = set(g for g in G if monomial_div(sdp_LM(F[g], u), LM))
                P0 = set(p for p in P if monomial_div(sdp_LM(F[p], u), LM))

                G, P, R = G - G0, P - P0 | set([k]), R | G0 | P0

                for i, j in set(B):
                    if i in G0 or j in G0:
                        del B[(i, j)]

        G |= P

        for i in G:
            for j in P:
                if i == j:
                    continue

                if i < j:
                   k = (i, j)
                else:
                   k = (j, i)

                if k not in B:
                    B[k] = monomial_lcm(sdp_LM(F[i], u), sdp_LM(F[j], u))

        G = set([ normal(F[g], G - set([g]))[0] for g in G ])

        return R, P, G, B

def sdp_lcm(f, g, u, O, K):
    """
    Computes LCM of two polynomials in `K[X]`.

    The LCM is computed as the unique generater of the intersection
    of the two ideals generated by `f` and `g`. The approach is to
    compute a Groebner basis with respect to lexicographic ordering
    of `t*f` and `(1 - t)*g`, where `t` is an unrelated variable and
    then filtering out the solution that doesn't contain `t`.

    References
    ==========

    1. [Cox97]_

    """
    from sympy.polys.groebnertools import sdp_groebner

    if not f or not g:
        return []

    if sdp_term_p(f) and sdp_term_p(g):
        monom = monomial_lcm(sdp_LM(f, u), sdp_LM(g, u))

        fc, gc = sdp_LC(f, K), sdp_LC(g, K)

        if K.has_Field:
            coeff = K.one
        else:
            coeff = K.lcm(fc, gc)

        return [(monom, coeff)]

    if not K.has_Field:
        lcm = K.one
    else:
        fc, f = sdp_primitive(f, K)
        gc, g = sdp_primitive(g, K)

        lcm = K.lcm(fc, gc)

    f_terms = tuple( ((1,) + m,  c) for m, c in f )
    g_terms = tuple( ((0,) + m,  c) for m, c in g ) \
            + tuple( ((1,) + m, -c) for m, c in g )

    F = sdp_sort(f_terms, lex)
    G = sdp_sort(g_terms, lex)

    basis = sdp_groebner([F, G], u, lex, K)

    H = [ h for h in basis if sdp_indep_p(h, 0, u) ]

    if K.is_one(lcm):
        h = [ (m[1:], c)     for m, c in H[0] ]
    else:
        h = [ (m[1:], c * lcm) for m, c in H[0] ]

    return sdp_sort(h, O)

def groebner_lcm(f, g):
    """
    Computes LCM of two polynomials using Groebner bases.

    The LCM is computed as the unique generater of the intersection
    of the two ideals generated by `f` and `g`. The approach is to
    compute a Groebner basis with respect to lexicographic ordering
    of `t*f` and `(1 - t)*g`, where `t` is an unrelated variable and
    then filtering out the solution that doesn't contain `t`.

    References
    ==========

    1. [Cox97]_

    """
    assert f.ring == g.ring

    ring = f.ring
    domain = ring.domain

    if not f or not g:
        return ring.zero

    if len(f) <= 1 and len(g) <= 1:
        monom = monomial_lcm(f.LM, g.LM)
        coeff = domain.lcm(f.LC, g.LC)
        return ring.term_new(monom, coeff)

    fc, f = f.primitive()
    gc, g = g.primitive()

    lcm = domain.lcm(fc, gc)

    f_terms = [ ((1,) + monom, coeff) for monom, coeff in f.terms() ]
    g_terms = [ ((0,) + monom, coeff) for monom, coeff in g.terms() ] \
            + [ ((1,) + monom,-coeff) for monom, coeff in g.terms() ]

    t = Dummy("t")
    t_ring = ring.clone(symbols=(t,) + ring.symbols, order=lex)

    F = t_ring.from_terms(f_terms)
    G = t_ring.from_terms(g_terms)

    basis = groebner([F, G], t_ring)

    def is_independent(h, j):
        return all(not monom[j] for monom in h.monoms())

    H = [ h for h in basis if is_independent(h, 0) ]

    h_terms = [ (monom[1:], coeff*lcm) for monom, coeff in H[0].terms() ]
    h = ring.from_terms(h_terms)

    return h

def sdm_spoly(f, g, O, K, phantom=None):
    """
    Compute the generalized s-polynomial of ``f`` and ``g``.

    The ground field is assumed to be ``K``, and monomials ordered according to
    ``O``.

    This is invalid if either of ``f`` or ``g`` is zero.

    If the leading terms of `f` and `g` involve different basis elements of
    `F`, their s-poly is defined to be zero. Otherwise it is a certain linear
    combination of `f` and `g` in which the leading terms cancel.
    See [SCA, defn 2.3.6] for details.

    If ``phantom`` is not ``None``, it should be a pair of module elements on
    which to perform the same operation(s) as on ``f`` and ``g``. The in this
    case both results are returned.

    Examples
    ========

    >>> from sympy.polys.distributedmodules import sdm_spoly
    >>> from sympy.polys import QQ, lex
    >>> f = [((2, 1, 1), QQ(1)), ((1, 0, 1), QQ(1))]
    >>> g = [((2, 3, 0), QQ(1))]
    >>> h = [((1, 2, 3), QQ(1))]
    >>> sdm_spoly(f, h, lex, QQ)
    []
    >>> sdm_spoly(f, g, lex, QQ)
    [((1, 2, 1), 1/1)]
    """
    if not f or not g:
        return sdm_zero()
    LM1 = sdm_LM(f)
    LM2 = sdm_LM(g)
    if LM1[0] != LM2[0]:
        return sdm_zero()
    LM1 = LM1[1:]
    LM2 = LM2[1:]
    lcm = monomial_lcm(LM1, LM2)
    m1 = monomial_div(lcm, LM1)
    m2 = monomial_div(lcm, LM2)
    c = K.quo(-sdm_LC(f, K), sdm_LC(g, K))
    r1 = sdm_add(sdm_mul_term(f, (m1, K.one), O, K),
                 sdm_mul_term(g, (m2, c), O, K), O, K)
    if phantom is None:
        return r1
    r2 = sdm_add(sdm_mul_term(phantom[0], (m1, K.one), O, K),
                 sdm_mul_term(phantom[1], (m2, c), O, K), O, K)
    return r1, r2

def spoly(p1, p2):
    """
    Compute LCM(LM(p1), LM(p2))/LM(p1)*p1 - LCM(LM(p1), LM(p2))/LM(p2)*p2
    This is the S-poly provided p1 and p2 are monic
    """
    LM1 = p1.LM
    LM2 = p2.LM
    LCM12 = monomial_lcm(LM1, LM2)
    m1 = monomial_div(LCM12, LM1)
    m2 = monomial_div(LCM12, LM2)
    s1 = p1.mul_monom(m1)
    s2 = p2.mul_monom(m2)
    s = s1 - s2
    return s

def sdm_monomial_lcm(A, B):
    """
    Return the "least common multiple" of ``A`` and ``B``.

    IF `A = M e_j` and `B = N e_j`, where `M` and `N` are polynomial monomials,
    this returns `\lcm(M, N) e_j`. Note that ``A`` and ``B`` involve distinct
    monomials.

    Otherwise the result is undefined.

    >>> from sympy.polys.distributedmodules import sdm_monomial_lcm
    >>> sdm_monomial_lcm((1, 2, 3), (1, 0, 5))
    (1, 2, 5)
    """
    return (A[0],) + monomial_lcm(A[1:], B[1:])

 def prune(P, S, h):
     """
     Prune the pair-set by applying the chain criterion
     [SCA, remark 2.5.11].
     """
     remove = set()
     for (a, b, c) in permutations(S, 3):
         A = sdm_LM(a)
         B = sdm_LM(b)
         C = sdm_LM(c)
         if len(set([A[0], B[0], C[0]])) != 1 or not h in [a, b, c]:
             continue
         if monomial_divides(B[1:], monomial_lcm(A[1:], C[1:])):
             remove.add((tuple(a), tuple(c)))
     return [(f, g) for (f, g) in P if (h not in [f, g]) or \
                 ((tuple(f), tuple(g)) not in remove and \
                  (tuple(g), tuple(f)) not in remove)]

def S_poly(tp1, tp2):
    """expv1,p1 = tp1 with expv1 = p1.leading_expv(), p1 monic; 
    similarly for tp2.
    Compute LCM(LM(p1),LM(p2))/LM(p1)*p1 - LCM(LM(p1),LM(p2))/LM(p2)*p2
    Throw LPolyOverflowError if bits_exp is too small for the result.
    """
    expv1, p1 = tp1
    expv2, p2 = tp2
    lp = p1.lp
    lcm12 = monomial_lcm(expv1, expv2)
    m1 = monomial_div(lcm12, expv1)
    m2 = monomial_div(lcm12, expv2)
    # TODO oprimize
    res = Poly(lp)
    res.iadd_m_mul_q(p1, (m1, 1))
    res.iadd_m_mul_q(p2, (m2, -1))
    return res

def sdm_spoly(f, g, O, K):
    """
    Compute the generalized s-polynomial of ``f`` and ``g``.

    The ground field is assumed to be ``K``, and monomials ordered according to
    ``O``.

    This is invalid if either of ``f`` or ``g`` is zero.

    If the leading terms of `f` and `g` involve different basis elements of
    `F`, their s-poly is defined to be zero. Otherwise it is a certain linear
    combination of `f` and `g` in which the leading terms cancel.
    See [SCA, defn 2.3.6] for details.

    Examples
    ========

    >>> from sympy.polys.distributedmodules import sdm_spoly
    >>> from sympy.polys import QQ, lex
    >>> f = [((2, 1, 1), QQ(1)), ((1, 0, 1), QQ(1))]
    >>> g = [((2, 3, 0), QQ(1))]
    >>> h = [((1, 2, 3), QQ(1))]
    >>> sdm_spoly(f, h, lex, QQ)
    []
    >>> sdm_spoly(f, g, lex, QQ)
    [((1, 2, 1), 1/1)]
    """
    if not f or not g:
        return sdm_zero()
    LM1 = sdm_LM(f)
    LM2 = sdm_LM(g)
    if LM1[0] != LM2[0]:
        return sdm_zero()
    LM1 = LM1[1:]
    LM2 = LM2[1:]
    lcm = monomial_lcm(LM1, LM2)
    return sdm_add(sdm_mul_term(f, (monomial_div(lcm, LM1), K.one), O, K),
                   sdm_mul_term(g, (monomial_div(lcm, LM2),
                                    K.quo(-sdm_LC(f, K), sdm_LC(g, K))), O, K),
                   O, K)

 def select(P):
     # normal selection strategy
     # select the pair with minimum LCM(LM(f), LM(g))
     pr = min(P, key=lambda pair: O(monomial_lcm(sdp_LM(f[pair[0]], u), sdp_LM(f[pair[1]], u))))
     return pr

def test_monomial_lcm():
    assert monomial_lcm((3,4,1), (1,2,0)) == (3,4,1)

    def update(G, B, ih):
        # update G using the set of critical pairs B and h
        # [BW] page 230
        h = f[ih]
        mh = sdp_LM(h, u)

        # filter new pairs (h, g), g in G
        C = G.copy()
        D = set()

        while C:
            # select a pair (h, g) by popping an element from C
            ig = C.pop()
            g = f[ig]
            mg = sdp_LM(g, u)
            LCMhg = monomial_lcm(mh, mg)

            def lcm_divides(ip):
                # LCM(LM(h), LM(p)) divides LCM(LM(h), LM(g))
                m = monomial_lcm(mh, sdp_LM(f[ip], u))
                return monomial_div(LCMhg, m)

            # HT(h) and HT(g) disjoint: mh*mg == LCMhg
            if monomial_mul(mh, mg) == LCMhg or (
                not any(lcm_divides(ipx) for ipx in C) and not any(lcm_divides(pr[1]) for pr in D)
            ):
                D.add((ih, ig))

        E = set()

        while D:
            # select h, g from D (h the same as above)
            ih, ig = D.pop()
            mg = sdp_LM(f[ig], u)
            LCMhg = monomial_lcm(mh, mg)

            if not monomial_mul(mh, mg) == LCMhg:
                E.add((ih, ig))

        # filter old pairs
        B_new = set()

        while B:
            # select g1, g2 from B (-> CP)
            ig1, ig2 = B.pop()
            mg1 = sdp_LM(f[ig1], u)
            mg2 = sdp_LM(f[ig2], u)
            LCM12 = monomial_lcm(mg1, mg2)

            # if HT(h) does not divide lcm(HT(g1), HT(g2))
            if not monomial_div(LCM12, mh) or monomial_lcm(mg1, mh) == LCM12 or monomial_lcm(mg2, mh) == LCM12:
                B_new.add((ig1, ig2))

        B_new |= E

        # filter polynomials
        G_new = set()

        while G:
            ig = G.pop()
            mg = sdp_LM(f[ig], u)

            if not monomial_div(mg, mh):
                G_new.add(ig)

        G_new.add(ih)

        return G_new, B_new

 def lcm_divides(ip):
     # LCM(LM(h), LM(p)) divides LCM(LM(h), LM(g))
     m = monomial_lcm(mh, sdp_LM(f[ip], u))
     return monomial_div(LCMhg, m)

 def lcm_divides(ip):
     # LCM(LM(h), LM(p)) divides LCM(LM(h), LM(g))
     m = monomial_lcm(mh, f[ip].LM)
     return monomial_div(LCMhg, m)

 def select(P):
     # normal selection strategy
     # select the pair with minimum LCM(LM(f), LM(g))
     pr = min(P, key=lambda pair: order(monomial_lcm(f[pair[0]].LM, f[pair[1]].LM)))
     return pr

 def select(P):
     # normal selection strategy
     # select the pair with minimum LCM(LM(f), LM(g))
     pr = minkey(P, key=lambda (i, j): O(monomial_lcm(sdp_LM(f[i], u), sdp_LM(f[j], u))))
     return pr

def sdm_groebner(G, NF, O, K):
    """
    Compute a minimal standard basis of ``G`` with respect to order ``O``.

    The algorithm uses a normal form ``NF``, for example ``sdm_nf_mora``.
    The ground field is assumed to be ``K``, and monomials ordered according
    to ``O``.

    Let `N` denote the submodule generated by elements of `G`. A standard
    basis for `N` is a subset `S` of `N`, such that `in(S) = in(N)`, where for
    any subset `X` of `F`, `in(X)` denotes the submodule generated by the
    initial forms of elements of `X`. [SCA, defn 2.3.2]

    A standard basis is called minimal if no subset of it is a standard basis.

    One may show that standard bases are always generating sets.

    Minimal standard bases are not unique. This algorithm computes a
    deterministic result, depending on the particular order of `G`.

    See [SCA, algorithm 2.3.8, and remark 1.6.3].
    """
    # First compute a standard basis
    S = [f for f in G if f]
    P = list(combinations(S, 2))

    def prune(P, S, h):
        """
        Prune the pair-set by applying the chain criterion
        [SCA, remark 2.5.11].
        """
        remove = set()
        retain = set()
        for (a, b, c) in permutations(S, 3):
            A = sdm_LM(a)
            B = sdm_LM(b)
            C = sdm_LM(c)
            if len(set([A[0], B[0], C[0]])) != 1 or not h in [a, b, c] or \
               any(tuple(x) in retain for x in [a, b, c]):
                continue
            if monomial_divides(B[1:], monomial_lcm(A[1:], C[1:])):
                remove.add((tuple(a), tuple(c)))
                retain.update([tuple(b), tuple(c), tuple(a)])
        return [(f, g) for (f, g) in P if (h not in [f, g]) or \
                    ((tuple(f), tuple(g)) not in remove and \
                     (tuple(g), tuple(f)) not in remove)]

    while P:
        # TODO better data structures!!!
        #print len(P), len(S)
        # Use the "normal selection strategy"
        lcms = [(i, sdm_LM(f)[:1] + monomial_lcm(sdm_LM(f)[1:], sdm_LM(g)[1:])) for \
                i, (f, g) in enumerate(P)]
        i = min(lcms, key=lambda x: O(x[1]))[0]
        f, g = P.pop(i)
        h = NF(sdm_spoly(f, g, O, K), S, O, K)
        if h:
            S.append(h)
            P.extend((h, f) for f in S if sdm_LM(h)[0] == sdm_LM(f)[0])
            P = prune(P, S, h)

    # Now interreduce it. (TODO again, better data structures)
    S = set(tuple(f) for f in S)
    for a, b in permutations(S, 2):
        A = sdm_LM(list(a))
        B = sdm_LM(list(b))
        if sdm_monomial_divides(A, B) and b in S and a in S:
            S.remove(b)

    return sorted((list(f) for f in S), key=lambda f: O(sdm_LM(f)),
                  reverse=True)