def update(G, CP, h):
        """update G using the set of critical pairs CP and h = (expv,pi)
        see [BW] page 230
        """
        hexpv, hp = f[h]
        # print 'DB10',hp
        # 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
            g = C.pop()
            gexpv = f[g][0]
            LCMhg = lcm_expv(hexpv, gexpv)

            def lcm_divides(p):
                expv = lcm_expv(hexpv, f[p][0])
                # LCM(LM(h), LM(p)) divides LCM(LM(h),LM(g))
                return monomial_div(LCMhg, expv)

            # HT(h) and HT(g) disjoint: hexpv + gexpv == LCMhg
            if monomial_mul(hexpv, gexpv) == LCMhg or (
                not any(lcm_divides(f) for f in C) and not any(lcm_divides(pr[1]) for pr in D)
            ):
                D.add((h, g))

        E = set()
        while D:
            # select h,g from D
            h, g = D.pop()
            gexpv = f[g][0]
            LCMhg = lcm_expv(hexpv, gexpv)
            if not monomial_mul(hexpv, gexpv) == LCMhg:
                E.add((h, g))

        # filter old pairs
        B_new = set()

        while CP:
            # select g1,g2 from CP
            g1, g2 = CP.pop()
            g1expv = f[g1][0]
            g2expv = f[g2][0]
            LCM12 = lcm_expv(g1expv, g2expv)
            # if HT(h) does not divide lcm(HT(g1),HT(g2))
            if not monomial_div(LCM12, hexpv) or lcm_expv(g1expv, hexpv) == LCM12 or lcm_expv(g2expv, hexpv) == LCM12:
                B_new.add((g1, g2))

        B_new |= E

        # filter polynomials
        G_new = set()
        while G:
            g = G.pop()
            if not monomial_div(f[g][0], hexpv):
                G_new.add(g)
        G_new.add(h)

        return G_new, B_new

    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 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_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 _basis(G, ring):
    """
    Computes a list of monomials which are not divisible by the leading
    monomials wrt to ``O`` of ``G``. These monomials are a basis of
    `K[X_1, \ldots, X_n]/(G)`.
    """
    order = ring.order

    leading_monomials = [g.LM for g in G]
    candidates = [ring.zero_monom]
    basis = []

    while candidates:
        t = candidates.pop()
        basis.append(t)

        new_candidates = [_incr_k(t, k) for k in xrange(ring.ngens)
            if all(monomial_div(_incr_k(t, k), lmg) is None
            for lmg in leading_monomials)]
        candidates.extend(new_candidates)
        candidates.sort(key=lambda m: order(m), reverse=True)

    basis = list(set(basis))

    return sorted(basis, key=lambda m: order(m))

def dmp_terms_gcd(f, u, K):
    """
    Remove GCD of terms from ``f`` in ``K[X]``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densebasic import dmp_terms_gcd

    >>> f = ZZ.map([[1, 0], [1, 0, 0], [], []])

    >>> dmp_terms_gcd(f, 1, ZZ)
    ((2, 1), [[1], [1, 0]])

    """
    if dmp_ground_TC(f, u, K) or dmp_zero_p(f, u):
        return (0,)*(u + 1), f

    F = dmp_to_dict(f, u)
    G = monomial_min(*F.keys())

    if all(g == 0 for g in G):
        return G, f

    f = {}

    for monom, coeff in F.iteritems():
        f[monomial_div(monom, G)] = coeff

    return G, dmp_from_dict(f, u, K)

def matrix_fglm(F, u, O_from, O_to, K):
    """
    Converts the reduced Groebner basis ``F`` of a zero-dimensional
    ideal w.r.t. ``O_from`` to a reduced Groebner basis
    w.r.t. ``O_to``.

    **References**
    J.C. Faugere, P. Gianni, D. Lazard, T. Mora (1994). Efficient
    Computation of Zero-dimensional Groebner Bases by Change of
    Ordering

    J.C. Faugere's lecture notes:
    http://www-salsa.lip6.fr/~jcf/Papers/2010_MPRI5e.pdf
    """
    old_basis = _basis(F, u, O_from, K)
    M = _representing_matrices(old_basis, F, u, O_from, K)

    # V contains the normalforms (wrt O_from) of S
    S = [(0,) * (u + 1)]
    V = [[K.one] + [K.zero] * (len(old_basis) - 1)]
    G = []

    L = [(i, 0) for i in xrange(u + 1)]  # (i, j) corresponds to x_i * S[j]
    L.sort(key=lambda (k, l): O_to(_incr_k(S[l], k)), reverse=True)
    t = L.pop()

    P = _identity_matrix(len(old_basis), K)

    while True:
        s = len(S)
        v = _matrix_mul(M[t[0]], V[t[1]], K)
        _lambda = _matrix_mul(P, v, K)

        if all(_lambda[i] == K.zero for i in xrange(s, len(old_basis))):
            # there is a linear combination of v by V

            lt = [(_incr_k(S[t[1]], t[0]), K.one)]
            rest = sdp_strip(sdp_sort([(S[i], _lambda[i]) for i in xrange(s)], O_to))
            g = sdp_sub(lt, rest, u, O_to, K)

            if g != []:
                G.append(g)

        else:
            # v is linearly independant from V
            P = _update(s, _lambda, P, K)
            S.append(_incr_k(S[t[1]], t[0]))
            V.append(v)

            L.extend([(i, s) for i in xrange(u + 1)])
            L = list(set(L))
            L.sort(key=lambda (k, l): O_to(_incr_k(S[l], k)), reverse=True)

        L = [(k, l) for (k, l) in L if all(monomial_div(_incr_k(S[l], k), sdp_LM(g, u)) is None for g in G)]

        if not L:
            G = [sdp_monic(g, K) for g in G]
            return sorted(G, key=lambda g: O_to(sdp_LM(g, u)), reverse=True)

        t = L.pop()

def _basis(G, u, O, K):
    """
    Computes a list of monomials which are not divisible by the leading
    monomials wrt to ``O`` of ``G``. These monomials are a basis of
    `K[X_1, \ldots, X_n]/(G)`.
    """
    leading_monomials = [sdp_LM(g, u) for g in G]
    candidates = [(0,) * (u + 1)]
    basis = []

    while candidates:
        t = candidates.pop()
        basis.append(t)

        new_candidates = [
            _incr_k(t, k)
            for k in xrange(u + 1)
            if all(monomial_div(_incr_k(t, k), lmg) is None for lmg in leading_monomials)
        ]
        candidates.extend(new_candidates)
        candidates.sort(key=lambda m: O(m), reverse=True)

    basis = list(set(basis))

    return sorted(basis, key=lambda m: O(m))

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 _term_rr_div(a, b, K):
    """Division of two terms in over a ring. """
    a_lm, a_lc = a
    b_lm, b_lc = b

    monom = monomial_div(a_lm, b_lm)

    if not (monom is None or a_lc % b_lc):
        return monom, K.quo(a_lc, b_lc)
    else:
        return None

def _term_ff_div(a, b, K):
    """Division of two terms in over a field. """
    a_lm, a_lc = a
    b_lm, b_lc = b

    monom = monomial_div(a_lm, b_lm)

    if monom is not None:
        return monom, K.quo(a_lc, b_lc)
    else:
        return None

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 _term_div(self):
        zm = self.ring.zero_monom
        domain = self.ring.domain
        domain_quo = domain.quo

        if domain.has_Field or not domain.is_Exact:
            def term_div((a_lm, a_lc), (b_lm, b_lc)):
                if b_lm == zm: # apparently this is a very common case
                    monom = a_lm
                else:
                    monom = monomial_div(a_lm, b_lm)
                if monom is not None:
                    return monom, domain_quo(a_lc, b_lc)
                else:
                    return None

def dmp_terms_gcd(f, u, K):
    """Remove GCD of terms from `f` in `K[X]`. """
    if dmp_ground_TC(f, u, K) or dmp_zero_p(f, u):
        return (0,)*(u+1), f

    F = dmp_to_dict(f, u)
    G = monomial_min(*F.keys())

    if all([ g == 0 for g in G ]):
        return G, f

    f = {}

    for monom, coeff in F.iteritems():
        f[monomial_div(monom, G)] = coeff

    return G, dmp_from_dict(f, u, K)

def test_monomial_div():
    assert monomial_div((3,4,1), (1,2,0)) == (2,2,1)

 def lcm_divides(p):
     expv = lcm_expv(hexpv, f[p][0])
     # LCM(LM(h), LM(p)) divides LCM(LM(h),LM(g))
     return monomial_div(LCMhg, expv)

def sdp_groebner(F, u, O, K):
    """Computes Groebner basis for a set of polynomials in `K[X]`.

       Given a set of multivariate polynomials `F`, finds another
       set `G`, such that Ideal `F = Ideal G` and `G` is a reduced
       Groebner basis.

       The resulting basis is unique and has monic generators if the
       ground domains is a field. Otherwise the result is non-unique
       but Groebner bases over e.g. integers can be computed (if the
       input polynomials are monic).

       Groebner bases can be used to choose specific generators for a
       polynomial ideal. Because these bases are unique you can check
       for ideal equality by comparing the Groebner bases.  To see if
       one polynomial lies in an ideal, divide by the elements in the
       base and see if the remainder vanishes.

       They can also be used to  solve systems of polynomial equations
       as,  by choosing lexicographic ordering,  you can eliminate one
       variable at a time, provided that the ideal is zero-dimensional
       (finite number of solutions).

       References
       ==========

       .. [Bose03] N.K. Bose, B. Buchberger, J.P. Guiver, Multidimensional
           Systems Theory and Applications, Springer, 2003, pp. 98+

       .. [Giovini91] A. Giovini, T. Mora, "One sugar cube, please" or
           Selection strategies in Buchberger algorithm, ISSAC '91, ACM

       .. [Ajwa95] I.A. Ajwa, Z. Liu, P.S. Wang, Groebner Bases Algorithm,
           http://citeseer.ist.psu.edu/ajwa95grbner.html, 1995

       .. [Cox97] D. Cox, J. Little, D. O'Shea, Ideals, Varieties and
           Algorithms, Springer, Second Edition, 1997, pp. 62

    """
    F = [ f for f in F if f ]

    if not F:
        return [[]]

    R, P, G, B, I = set(), set(), set(), {}, {}

    for i, f in enumerate(F):
        I[tuple(f)] = i
        R.add(i)

    def normal(g, J):
        h = sdp_rem(g, [ F[j] for j in J ], u, O, K)

        if not h:
            return None
        else:
            H = tuple(h)

            if not H in I:
                I[H] = len(F)
                F.append(h)

            return I[H], sdp_LM(h, u)

    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

    R, P, G, B = generate(R, P, G, B)
#.........这里部分代码省略.........

 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 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, f[ip].LM)
     return monomial_div(LCMhg, m)