def square(self):
        """square of a polynomial

        Examples
        ========

        >>> from sympy.polys.rings import ring
        >>> from sympy.polys.domains import ZZ
        >>> _, x, y = ring('x, y', ZZ)
        >>> p = x + y**2
        >>> p.square()
        x**2 + 2*x*y**2 + y**4
        """
        ring = self.ring
        p = ring.zero
        get = p.get
        keys = self.keys()
        for i in range(len(keys)):
            k1 = keys[i]
            pk = self[k1]
            for j in range(i):
                k2 = keys[j]
                exp = monomial_mul(k1, k2)
                p[exp] = get(exp, 0) + pk*self[k2]
        p = p.imul_num(2)
        get = p.get
        for k, v in self.iteritems():
            k2 = monomial_mul(k, k)
            p[k2] = get(k2, 0) + v**2
        p.strip_zero()
        return p

    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 sdp_mul_term(f, term, u, O, K):
    """Multiply a distributed polynomial by a term. """
    M, c = term

    if not f or not c:
        return []
    else:
        if K.is_one(c):
            return [ (monomial_mul(f_M, M), f_c) for f_M, f_c in f ]
        else:
            return [ (monomial_mul(f_M, M), f_c * c) for f_M, f_c in f ]

    def _iadd_poly_monom(p1, p2, mc):
        """add to self the product of (p)*(coeff*x0**i0*x1**i1*...)
        unless self is a generator -- then just return the sum of the two.

        mc is a tuple, (monom, coeff), where monomial is (i0, i1, ...)

        Examples
        ========

        >>> from sympy.polys.rings import ring
        >>> from sympy.polys.domains import ZZ

        >>> _, x, y, z = ring('x, y, z', ZZ)
        >>> p1 = x**4 + 2*y
        >>> p2 = y + z
        >>> m = (1, 2, 3)
        >>> p1 = p1._iadd_poly_monom(p2, (m, 3))
        >>> p1
        x**4 + 3*x*y**3*z**3 + 3*x*y**2*z**4 + 2*y

        """
        if p1 in p1.ring.gens:
            p1 = p1.copy()
        (m, c) = mc
        get = p1.get
        zero = p1.ring.domain.zero
        for k, v in p2.iteritems():
            ka = monomial_mul(k, m)
            coeff = get(ka, zero) + v*c
            if coeff:
                p1[ka] = coeff
            else:
                del p1[ka]
        return p1

def is_zero_dimensional(F, u, O, K):
    """
    Checks if the ideal generated by ``F`` is zero-dimensional.

    The algorithm checks if the set of monomials not divisible by a
    leading monomial of any element of ``F`` is bounded. In general
    ``F`` has to be a Groebner basis w.r.t. ``O`` but if ``True``
    is returned, then the ideal is zero-dimensional.

    **References**
    David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties
    and Algorithms, 3rd edition, p. 234
    """
    def single_var(m):
        n = 0
        for e in m:
            if e != 0:
                n += 1
        return n == 1

    exponents = (0,) * (u + 1)

    for f in F:
        if single_var(sdp_LM(f, u)):
            exponents = monomial_mul(exponents, sdp_LM(f, O))  # == sum of exponent vectors

    product = 1
    for e in exponents:
        product *= e

    # If product == 0, then there's a variable for which there's
    # no degree bound.
    return product != 0

def sdp_mul(f, g, u, O, K):
    """Multiply distributed polynomials in `K[X]`. """
    if sdp_term_p(f):
        if not f:
            return f
        else:
            return sdp_mul_term(g, f[0], u, O, K)

    if sdp_term_p(g):
        if not g:
            return g
        else:
            return sdp_mul_term(f, g[0], u, O, K)

    h = {}

    for fm, fc in f:
        for gm, gc in g:
            monom = monomial_mul(fm, gm)
            coeff = fc * gc

            if monom in h:
                coeff += h[monom]

                if not coeff:
                    del h[monom]
                    continue

            h[monom] = coeff

    return sdp_from_dict(h, O)

def sig_mult(s, m):
    """
    Multiply a signature by a monomial.

    The product of a signature (m, i) and a monomial n is defined as
    (m * t, i).
    """
    return sig(monomial_mul(s[0], m), s[1])

    def mul_term(f, term):
        monom, coeff = term

        if not f or not coeff:
            return f.ring.zero
        elif monom == f.ring.zero_monom:
            return f.mul_ground(coeff)

        terms = [ (monomial_mul(f_monom, monom), f_coeff*coeff) for f_monom, f_coeff in f.iteritems() ]
        return f.new(terms)

    def representing_matrix(m):
        M = [[K.zero] * len(basis) for _ in xrange(len(basis))]

        for i, v in enumerate(basis):
            r = sdp_rem([(monomial_mul(m, v), K.one)], G, u, O, K)

            for term in r:
                j = basis.index(term[0])
                M[j][i] = term[1]

        return M

    def representing_matrix(m):
        M = [[domain.zero] * len(basis) for _ in xrange(len(basis))]

        for i, v in enumerate(basis):
            r = ring.term_new(monomial_mul(m, v), domain.one).rem(G)

            for monom, coeff in r.terms():
                j = basis.index(monom)
                M[j][i] = coeff

        return M

def sdm_monomial_mul(M, X):
    """
    Multiply tuple ``X`` representing a monomial of `K[X]` into the tuple
    ``M`` representing a monomial of `F`.

    Examples
    ========

    Multiplying `xy^3` into `x f_1` yields `x^2 y^3 f_1`:

    >>> from sympy.polys.distributedmodules import sdm_monomial_mul
    >>> sdm_monomial_mul((1, 1, 0), (1, 3))
    (1, 2, 3)
    """
    return (M[0],) + monomial_mul(X, M[1:])

def sdp_rem(f, G, u, O, K):
    """Returns polynomial remainder in `K[X]`. """
    r = {}

    if K.has_Field:
        term_div = _term_ff_div
    else:
        term_div = _term_rr_div

    ltf = sdp_LT(f, u, K)
    f = dict(f)
    get = f.get
    while f:
        for g in G:
            tq = term_div(ltf, sdp_LT(g, u, K), K)

            if tq is not None:
                m, c = tq
                for mg, cg in g:
                    m1 = monomial_mul(mg, m)
                    c1 = get(m1, 0) - c*cg
                    if not c1:
                        del f[m1]
                    else:
                        f[m1] = c1
                if f:
                    if O == lex:
                        ltm = max(f)
                    else:
                        ltm = max(f, key=lambda mx: O(mx))
                    ltf = ltm, f[ltm]

                break
        else:
            ltm, ltc = ltf
            if ltm in r:
                r[ltm] += ltc
            else:
                r[ltm] = ltc
            del f[ltm]
            if f:
                if O == lex:
                    ltm = max(f)
                else:
                    ltm = max(f, key=lambda mx: O(mx))
                ltf = ltm, f[ltm]
    return sdp_from_dict(r, O)

    def __mul__(p1, p2):
        """Multiply two polynomials.

        Examples
        ========

        >>> from sympy.polys.domains import QQ
        >>> from sympy.polys.rings import ring
        >>> _, x, y = ring('x, y', QQ)
        >>> p1 = x + y
        >>> p2 = x - y
        >>> p1*p2
        x**2 - y**2

        """
        ring = p1.ring
        p = ring.zero
        if not p2:
            return p
        if isinstance(p2, PolyElement):
            if ring == p2.ring:
                get = p.get
                p2it = p2.items()
                for exp1, v1 in p1.iteritems():
                    for exp2, v2 in p2it:
                        exp = monomial_mul(exp1, exp2)
                        p[exp] = get(exp, 0) + v1*v2
                p.strip_zero()
                return p
            elif isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
                pass
            elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
                return p2.__rmul__(p1)
            else:
                return NotImplemented

        try:
            p2 = ring.domain_new(p2)
        except CoercionFailed:
            return NotImplemented
        else:
            for exp1, v1 in p1.iteritems():
                v = v1*p2
                if v:
                    p[exp1] = v
            return p

    def rem(f, G):
        if isinstance(G, PolyElement):
            G = [G]
        domain = f.ring.domain
        order = f.ring.order
        r = f.ring.zero
        term_div = f._term_div()
        ltf = f.LT
        f = f.copy()
        get = f.get
        while f:
            for g in G:
                tq = term_div(ltf, g.LT)
                if tq is not None:
                    m, c = tq
                    for mg, cg in g.terms():
                        m1 = monomial_mul(mg, m)
                        c1 = get(m1, 0) - c*cg
                        if not c1:
                            del f[m1]
                        else:
                            f[m1] = c1
                    if f:
                        if order is lex:
                            ltm = max(f)
                        else:
                            ltm = max(f, key=order)
                        ltf = ltm, f[ltm]

                    break
            else:
                ltm, ltc = ltf
                if ltm in r:
                    r[ltm] += ltc
                else:
                    r[ltm] = ltc
                del f[ltm]
                if f:
                    if order is lex:
                        ltm = max(f)
                    else:
                        ltm = max(f, key=order)
                    ltf = ltm, f[ltm]
        return r

def sdp_sqr(f, u, O, K):
    """Square a distributed polynomial in `K[X]`. """
    h = {}

    for fm, fc in f:
        for Fm, Fc in f:
            monom = monomial_mul(fm, Fm)
            coeff = fc * Fc

            if monom in h:
                coeff += h[monom]

                if not coeff:
                    del h[monom]
                    continue

            h[monom] = coeff

    return sdp_from_dict(h, O)

def test_monomial_mul():
    assert monomial_mul((3,4,1), (1,2,0)) == (4,6,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

#.........这里部分代码省略.........
                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)

    while B:
        k, M = B.items()[0]

        for l, N in B.iteritems():
            if O(M, N) == 1:
                k, M = l, N

        del B[k]

        i, j = k[0], k[1]
        p, q = F[i], F[j]

        p_LM, q_LM = sdp_LM(p, u), sdp_LM(q, u)

        if M == monomial_mul(p_LM, q_LM):
            continue

        criterion = False

        for g in G:
            if g == i or g == j:
                continue

            if (min(i, g), max(i, g)) not in B:
                continue

            if (min(j, g), max(j, g)) not in B:
                continue

            if not monomial_div(M, sdp_LM(F[g], u)):
                continue

            criterion = True
            break

        if criterion:
            continue

        p = sdp_mul_term(p, (monomial_div(M, p_LM), K.quo(K.one, sdp_LC(p, K))), u, O, K)
        q = sdp_mul_term(q, (monomial_div(M, q_LM), K.quo(K.one, sdp_LC(q, K))), u, O, K)

        h = normal(sdp_sub(p, q, u, O, K), G)

        if h is not None:
            k, LM = h

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

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

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

            R, P, G, B = generate(R, P, G, B)

    if K.has_Field:
        basis = [ sdp_monic(F[g], K) for g in G ]
    else:
        basis = []

        for g in G:
            _, g = sdp_primitive(F[g], K)

            if K.is_negative(sdp_LC(g, K)):
                basis.append(sdp_neg(g, u, O, K))
            else:
                basis.append(g)

    return list(sorted(basis, O, lambda p: sdp_LM(p, u), True))

                return f[:i] + [(M, coeff)] + f[i+1:]
        else:
            if j > 0:
                hi = i - 1
            else:
                lo = i + 1
    else:
        return f[:i] + [(M, -c)] + f[i+1:]

def sdp_mul_term(f, (M, c), u, O, K):
    """Multiply a distributed polynomial by a term. """
    if not f or not c:
        return []
    else:
        if K.is_one(c):
            return [ (monomial_mul(f_M, M), f_c) for f_M, f_c in f ]
        else:
            return [ (monomial_mul(f_M, M), f_c*c) for f_M, f_c in f ]

def sdp_add(f, g, u, O, K):
    """Add distributed polynomials in `K[X]`. """
    h = dict(f)

    for monom, c in g:
        if h.has_key(monom):
            coeff = h[monom] + c

            if not coeff:
                del h[monom]
            else:
                h[monom] = coeff

 def mul_monom(f, monom):
     terms = [ (monomial_mul(f_monom, monom), f_coeff) for f_monom, f_coeff in f.iteritems() ]
     return f.new(terms)