def dmp_cancel(f, g, u, K, multout=True):
    """
    Cancel common factors in a rational function ``f/g``.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.euclidtools import dmp_cancel

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

    >>> dmp_cancel(f, g, 1, ZZ)
    ([[2], [2]], [[1], [-1]])

    """
    if dmp_zero_p(f, u) or dmp_zero_p(g, u):
        if multout:
            return f, g
        else:
            return K.one, K.one, f, g

    K0 = None

    if K.has_Field and K.has_assoc_Ring:
        K0, K = K, K.get_ring()

        cq, f = dmp_clear_denoms(f, u, K0, K, convert=True)
        cp, g = dmp_clear_denoms(g, u, K0, K, convert=True)
    else:
        cp, cq = K.one, K.one

    _, p, q = dmp_inner_gcd(f, g, u, K)

    if K0 is not None:
        p = dmp_convert(p, u, K, K0)
        q = dmp_convert(q, u, K, K0)

        K = K0

    p_neg = K.is_negative(dmp_ground_LC(p, u, K))
    q_neg = K.is_negative(dmp_ground_LC(q, u, K))

    if p_neg and q_neg:
        p, q = dmp_neg(p, u, K), dmp_neg(q, u, K)
    elif p_neg:
        cp, p = -cp, dmp_neg(p, u, K)
    elif q_neg:
        cp, q = -cp, dmp_neg(q, u, K)

    if not multout:
        return cp, cq, p, q

    p = dmp_mul_ground(p, cp, u, K)
    q = dmp_mul_ground(q, cq, u, K)

    return p, q

def dmp_cancel(f, g, u, K, include=True):
    """
    Cancel common factors in a rational function `f/g`.

    Examples
    ========

    >>> from sympy.polys import ring, ZZ
    >>> R, x,y = ring("x,y", ZZ)

    >>> R.dmp_cancel(2*x**2 - 2, x**2 - 2*x + 1)
    (2*x + 2, x - 1)

    """
    K0 = None

    if K.has_Field and K.has_assoc_Ring:
        K0, K = K, K.get_ring()

        cq, f = dmp_clear_denoms(f, u, K0, K, convert=True)
        cp, g = dmp_clear_denoms(g, u, K0, K, convert=True)
    else:
        cp, cq = K.one, K.one

    _, p, q = dmp_inner_gcd(f, g, u, K)

    if K0 is not None:
        _, cp, cq = K.cofactors(cp, cq)

        p = dmp_convert(p, u, K, K0)
        q = dmp_convert(q, u, K, K0)

        K = K0

    p_neg = K.is_negative(dmp_ground_LC(p, u, K))
    q_neg = K.is_negative(dmp_ground_LC(q, u, K))

    if p_neg and q_neg:
        p, q = dmp_neg(p, u, K), dmp_neg(q, u, K)
    elif p_neg:
        cp, p = -cp, dmp_neg(p, u, K)
    elif q_neg:
        cp, q = -cp, dmp_neg(q, u, K)

    if not include:
        return cp, cq, p, q

    p = dmp_mul_ground(p, cp, u, K)
    q = dmp_mul_ground(q, cq, u, K)

    return p, q

def dmp_sqf_part(f, u, K):
    """
    Returns square-free part of a polynomial in ``K[X]``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.sqfreetools import dmp_sqf_part

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

    >>> dmp_sqf_part(f, 1, ZZ)
    [[1], [1, 0], []]

    """
    if not u:
        return dup_sqf_part(f, K)

    if not K.has_CharacteristicZero:
        return dmp_gf_sqf_part(f, u, K)

    if dmp_zero_p(f, u):
        return f

    if K.is_negative(dmp_ground_LC(f, u, K)):
        f = dmp_neg(f, u, K)

    gcd = dmp_gcd(f, dmp_diff(f, 1, u, K), u, K)
    sqf = dmp_quo(f, gcd, u, K)

    if K.has_Field or not K.is_Exact:
        return dmp_ground_monic(sqf, u, K)
    else:
        return dmp_ground_primitive(sqf, u, K)[1]

def dmp_content(f, u, K):
    """
    Returns GCD of multivariate coefficients.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.euclidtools import dmp_content

    >>> f = ZZ.map([[2, 6], [4, 12]])

    >>> dmp_content(f, 1, ZZ)
    [2, 6]

    """
    cont, v = dmp_LC(f, K), u-1

    if dmp_zero_p(f, u):
        return cont

    for c in f[1:]:
        cont = dmp_gcd(cont, c, v, K)

        if dmp_one_p(cont, v, K):
            break

    if K.is_negative(dmp_ground_LC(cont, v, K)):
        return dmp_neg(cont, v, K)
    else:
        return cont

def dmp_sqf_part(f, u, K):
    """
    Returns square-free part of a polynomial in ``K[X]``.

    Examples
    ========

    >>> from sympy.polys import ring, ZZ
    >>> R, x,y = ring("x,y", ZZ)

    >>> R.dmp_sqf_part(x**3 + 2*x**2*y + x*y**2)
    x**2 + x*y

    """
    if not u:
        return dup_sqf_part(f, K)

    if K.is_FiniteField:
        return dmp_gf_sqf_part(f, u, K)

    if dmp_zero_p(f, u):
        return f

    if K.is_negative(dmp_ground_LC(f, u, K)):
        f = dmp_neg(f, u, K)

    gcd = dmp_gcd(f, dmp_diff(f, 1, u, K), u, K)
    sqf = dmp_quo(f, gcd, u, K)

    if K.has_Field:
        return dmp_ground_monic(sqf, u, K)
    else:
        return dmp_ground_primitive(sqf, u, K)[1]

def dmp_content(f, u, K):
    """
    Returns GCD of multivariate coefficients.

    Examples
    ========

    >>> from sympy.polys import ring, ZZ
    >>> R, x,y, = ring("x,y", ZZ)

    >>> R.dmp_content(2*x*y + 6*x + 4*y + 12)
    2*y + 6

    """
    cont, v = dmp_LC(f, K), u - 1

    if dmp_zero_p(f, u):
        return cont

    for c in f[1:]:
        cont = dmp_gcd(cont, c, v, K)

        if dmp_one_p(cont, v, K):
            break

    if K.is_negative(dmp_ground_LC(cont, v, K)):
        return dmp_neg(cont, v, K)
    else:
        return cont

def dmp_sqf_list(f, u, K, all=False):
    """
    Return square-free decomposition of a polynomial in ``K[X]``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.sqfreetools import dmp_sqf_list

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

    >>> dmp_sqf_list(f, 1, ZZ)
    (1, [([[1], [1, 0]], 2), ([[1], []], 3)])

    >>> dmp_sqf_list(f, 1, ZZ, all=True)
    (1, [([[1]], 1), ([[1], [1, 0]], 2), ([[1], []], 3)])

    """
    if not u:
        return dup_sqf_list(f, K, all=all)

    if not K.has_CharacteristicZero:
        return dmp_gf_sqf_list(f, u, K, all=all)

    if K.has_Field or not K.is_Exact:
        coeff = dmp_ground_LC(f, u, K)
        f = dmp_ground_monic(f, u, K)
    else:
        coeff, f = dmp_ground_primitive(f, u, K)

        if K.is_negative(dmp_ground_LC(f, u, K)):
            f = dmp_neg(f, u, K)
            coeff = -coeff

    if dmp_degree(f, u) <= 0:
        return coeff, []

    result, i = [], 1

    h = dmp_diff(f, 1, u, K)
    g, p, q = dmp_inner_gcd(f, h, u, K)

    while True:
        d = dmp_diff(p, 1, u, K)
        h = dmp_sub(q, d, u, K)

        if dmp_zero_p(h, u):
            result.append((p, i))
            break

        g, p, q = dmp_inner_gcd(p, h, u, K)

        if all or dmp_degree(g, u) > 0:
            result.append((g, i))

        i += 1

    return coeff, result

def dmp_sqf_list(f, u, K, all=False):
    """
    Return square-free decomposition of a polynomial in ``K[X]``.

    Examples
    ========

    >>> from sympy.polys import ring, ZZ
    >>> R, x,y = ring("x,y", ZZ)

    >>> f = x**5 + 2*x**4*y + x**3*y**2

    >>> R.dmp_sqf_list(f)
    (1, [(x + y, 2), (x, 3)])
    >>> R.dmp_sqf_list(f, all=True)
    (1, [(1, 1), (x + y, 2), (x, 3)])

    """
    if not u:
        return dup_sqf_list(f, K, all=all)

    if K.is_FiniteField:
        return dmp_gf_sqf_list(f, u, K, all=all)

    if K.has_Field:
        coeff = dmp_ground_LC(f, u, K)
        f = dmp_ground_monic(f, u, K)
    else:
        coeff, f = dmp_ground_primitive(f, u, K)

        if K.is_negative(dmp_ground_LC(f, u, K)):
            f = dmp_neg(f, u, K)
            coeff = -coeff

    if dmp_degree(f, u) <= 0:
        return coeff, []

    result, i = [], 1

    h = dmp_diff(f, 1, u, K)
    g, p, q = dmp_inner_gcd(f, h, u, K)

    while True:
        d = dmp_diff(p, 1, u, K)
        h = dmp_sub(q, d, u, K)

        if dmp_zero_p(h, u):
            result.append((p, i))
            break

        g, p, q = dmp_inner_gcd(p, h, u, K)

        if all or dmp_degree(g, u) > 0:
            result.append((g, i))

        i += 1

    return coeff, result

def test_dmp_neg():
    assert dmp_neg([ZZ(-1)], 0, ZZ) == [ZZ(1)]
    assert dmp_neg([QQ(-1,2)], 0, QQ) == [QQ(1,2)]

    assert dmp_neg([[[]]], 2, ZZ) == [[[]]]
    assert dmp_neg([[[ZZ(1)]]], 2, ZZ) == [[[ZZ(-1)]]]
    assert dmp_neg([[[ZZ(-7)]]], 2, ZZ) == [[[ZZ(7)]]]

    assert dmp_neg([[[]]], 2, QQ) == [[[]]]
    assert dmp_neg([[[QQ(1,9)]]], 2, QQ) == [[[QQ(-1,9)]]]
    assert dmp_neg([[[QQ(-7,9)]]], 2, QQ) == [[[QQ(7,9)]]]

    def _parse(cls, rep, dom, lev=None):
        if type(rep) is tuple:
            num, den = rep

            if lev is not None:
                if type(num) is dict:
                    num = dmp_from_dict(num, lev, dom)

                if type(den) is dict:
                    den = dmp_from_dict(den, lev, dom)
            else:
                num, num_lev = dmp_validate(num)
                den, den_lev = dmp_validate(den)

                if num_lev == den_lev:
                    lev = num_lev
                else:
                    raise ValueError('inconsistent number of levels')

            if dmp_zero_p(den, lev):
                raise ZeroDivisionError('fraction denominator')

            if dmp_zero_p(num, lev):
                den = dmp_one(lev, dom)
            else:
                if dmp_negative_p(den, lev, dom):
                    num = dmp_neg(num, lev, dom)
                    den = dmp_neg(den, lev, dom)
        else:
            num = rep

            if lev is not None:
                if type(num) is dict:
                    num = dmp_from_dict(num, lev, dom)
                elif type(num) is not list:
                    num = dmp_ground(dom.convert(num), lev)
            else:
                num, lev = dmp_validate(num)

            den = dmp_one(lev, dom)

        return num, den, lev

def _dmp_rr_trivial_gcd(f, g, u, K):
    """Handle trivial cases in GCD algorithm over a ring. """
    zero_f = dmp_zero_p(f, u)
    zero_g = dmp_zero_p(g, u)

    if zero_f and zero_g:
        return tuple(dmp_zeros(3, u, K))
    elif zero_f:
        if K.is_nonnegative(dmp_ground_LC(g, u, K)):
            return g, dmp_zero(u), dmp_one(u, K)
        else:
            return dmp_neg(g, u, K), dmp_zero(u), dmp_ground(-K.one, u)
    elif zero_g:
        if K.is_nonnegative(dmp_ground_LC(f, u, K)):
            return f, dmp_one(u, K), dmp_zero(u)
        else:
            return dmp_neg(f, u, K), dmp_ground(-K.one, u), dmp_zero(u)
    elif query('USE_SIMPLIFY_GCD'):
        return _dmp_simplify_gcd(f, g, u, K)
    else:
        return None

def test_dmp_sqf():
    assert dmp_sqf_part([[]], 1, ZZ) == [[]]
    assert dmp_sqf_p([[]], 1, ZZ) == True

    assert dmp_sqf_part([[7]], 1, ZZ) == [[1]]
    assert dmp_sqf_p([[7]], 1, ZZ) == True

    assert dmp_sqf_p(f_0, 2, ZZ) == True
    assert dmp_sqf_p(dmp_sqr(f_0, 2, ZZ), 2, ZZ) == False
    assert dmp_sqf_p(f_1, 2, ZZ) == True
    assert dmp_sqf_p(dmp_sqr(f_1, 2, ZZ), 2, ZZ) == False
    assert dmp_sqf_p(f_2, 2, ZZ) == True
    assert dmp_sqf_p(dmp_sqr(f_2, 2, ZZ), 2, ZZ) == False
    assert dmp_sqf_p(f_3, 2, ZZ) == True
    assert dmp_sqf_p(dmp_sqr(f_3, 2, ZZ), 2, ZZ) == False
    assert dmp_sqf_p(f_5, 2, ZZ) == False
    assert dmp_sqf_p(dmp_sqr(f_5, 2, ZZ), 2, ZZ) == False

    assert dmp_sqf_p(f_4, 2, ZZ) == True
    assert dmp_sqf_part(f_4, 2, ZZ) == dmp_neg(f_4, 2, ZZ)
    assert dmp_sqf_p(f_6, 3, ZZ) == True
    assert dmp_sqf_part(f_6, 3, ZZ) == f_6

    assert dmp_sqf_part(f_5, 2, ZZ) == [[[1]], [[1], [-1, 0]]]

    assert dup_sqf_list([], ZZ) == (ZZ(0), [])
    assert dup_sqf_list_include([], ZZ) == [([], 1)]

    assert dmp_sqf_list([[ZZ(3)]], 1, ZZ) == (ZZ(3), [])
    assert dmp_sqf_list_include([[ZZ(3)]], 1, ZZ) == [([[ZZ(3)]], 1)]

    f = [-1,1,0,0,1,-1]

    assert dmp_sqf_list(f, 0, ZZ) == \
        (-1, [([1,1,1,1], 1), ([1,-1], 2)])
    assert dmp_sqf_list_include(f, 0, ZZ) == \
        [([-1,-1,-1,-1], 1), ([1,-1], 2)]

    f = [[-1],[1],[],[],[1],[-1]]

    assert dmp_sqf_list(f, 1, ZZ) == \
        (-1, [([[1],[1],[1],[1]], 1), ([[1],[-1]], 2)])
    assert dmp_sqf_list_include(f, 1, ZZ) == \
        [([[-1],[-1],[-1],[-1]], 1), ([[1],[-1]], 2)]

    K = FF(2)

    f = [[-1], [2], [-1]]

    assert dmp_sqf_list_include(f, 1, ZZ) == \
        [([[-1]], 1), ([[1], [-1]], 2)]

    raises(DomainError, "dmp_sqf_list([[K(1), K(0), K(1)]], 1, K)")

def _dmp_zz_gcd_interpolate(h, x, v, K):
    """Interpolate polynomial GCD from integer GCD. """
    f = []

    while not dmp_zero_p(h, v):
        g = dmp_ground_trunc(h, x, v, K)
        f.insert(0, g)

        h = dmp_sub(h, g, v, K)
        h = dmp_exquo_ground(h, x, v, K)

    if K.is_negative(dmp_ground_LC(f, v+1, K)):
        return dmp_neg(f, v+1, K)
    else:
        return f

def dmp_fateman_poly_F_3(n, K):
    """Fateman's GCD benchmark: sparse inputs (deg f ~ vars f) """
    u = dup_from_raw_dict({n+1: K.one}, K)

    for i in xrange(0, n-1):
        u = dmp_add_term([u], dmp_one(i, K), n+1, i+1, K)

    v = dmp_add_term(u, dmp_ground(K(2), n-2), 0, n, K)

    f = dmp_sqr(dmp_add_term([dmp_neg(v, n-1, K)], dmp_one(n-1, K), n+1, n, K), n, K)
    g = dmp_sqr(dmp_add_term([v], dmp_one(n-1, K), n+1, n, K), n, K)

    v = dmp_add_term(u, dmp_one(n-2, K), 0, n-1, K)

    h = dmp_sqr(dmp_add_term([v], dmp_one(n-1, K), n+1, n, K), n, K)

    return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h

def dmp_rr_prs_gcd(f, g, u, K):
    """
    Computes polynomial GCD using subresultants over a ring.

    Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``, ``cff = quo(f, h)``,
    and ``cfg = quo(g, h)``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.euclidtools import dmp_rr_prs_gcd

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

    >>> dmp_rr_prs_gcd(f, g, 1, ZZ)
    ([[1], [1, 0]], [[1], [1, 0]], [[1], []])

    """
    if not u:
        return dup_rr_prs_gcd(f, g, K)

    result = _dmp_rr_trivial_gcd(f, g, u, K)

    if result is not None:
        return result

    fc, F = dmp_primitive(f, u, K)
    gc, G = dmp_primitive(g, u, K)

    h = dmp_subresultants(F, G, u, K)[-1]
    c, _, _ = dmp_rr_prs_gcd(fc, gc, u-1, K)

    if K.is_negative(dmp_ground_LC(h, u, K)):
        h = dmp_neg(h, u, K)

    _, h = dmp_primitive(h, u, K)
    h = dmp_mul_term(h, c, 0, u, K)

    cff = dmp_quo(f, h, u, K)
    cfg = dmp_quo(g, h, u, K)

    return h, cff, cfg

def dmp_rr_prs_gcd(f, g, u, K):
    """
    Computes polynomial GCD using subresultants over a ring.

    Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``, ``cff = quo(f, h)``,
    and ``cfg = quo(g, h)``.

    Examples
    ========

    >>> from sympy.polys import ring, ZZ
    >>> R, x,y, = ring("x,y", ZZ)

    >>> f = x**2 + 2*x*y + y**2
    >>> g = x**2 + x*y

    >>> R.dmp_rr_prs_gcd(f, g)
    (x + y, x + y, x)

    """
    if not u:
        return dup_rr_prs_gcd(f, g, K)

    result = _dmp_rr_trivial_gcd(f, g, u, K)

    if result is not None:
        return result

    fc, F = dmp_primitive(f, u, K)
    gc, G = dmp_primitive(g, u, K)

    h = dmp_subresultants(F, G, u, K)[-1]
    c, _, _ = dmp_rr_prs_gcd(fc, gc, u - 1, K)

    if K.is_negative(dmp_ground_LC(h, u, K)):
        h = dmp_neg(h, u, K)

    _, h = dmp_primitive(h, u, K)
    h = dmp_mul_term(h, c, 0, u, K)

    cff = dmp_quo(f, h, u, K)
    cfg = dmp_quo(g, h, u, K)

    return h, cff, cfg

def dmp_fateman_poly_F_2(n, K):
    """Fateman's GCD benchmark: linearly dense quartic inputs """
    u = [K(1), K(0)]

    for i in xrange(0, n - 1):
        u = [dmp_one(i, K), u]

    m = n - 1

    v = dmp_add_term(u, dmp_ground(K(2), m - 1), 0, n, K)

    f = dmp_sqr([dmp_one(m, K), dmp_neg(v, m, K)], n, K)
    g = dmp_sqr([dmp_one(m, K), v], n, K)

    v = dmp_add_term(u, dmp_one(m - 1, K), 0, n, K)

    h = dmp_sqr([dmp_one(m, K), v], n, K)

    return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h

def dmp_zz_factor(f, u, K):
    """
    Factor (non square-free) polynomials in `Z[X]`.

    Given a multivariate polynomial `f` in `Z[x]` computes its complete
    factorization `f_1, ..., f_n` into irreducibles over integers::

                 f = content(f) f_1**k_1 ... f_n**k_n

    The factorization is computed by reducing the input polynomial
    into a primitive square-free polynomial and factoring it using
    Enhanced Extended Zassenhaus (EEZ) algorithm. Trial division
    is used to recover the multiplicities of factors.

    The result is returned as a tuple consisting of::

             (content(f), [(f_1, k_1), ..., (f_n, k_n))

    Consider polynomial `f = 2*(x**2 - y**2)`::

        >>> from sympy.polys.factortools import dmp_zz_factor
        >>> from sympy.polys.domains import ZZ

        >>> dmp_zz_factor([[2], [], [-2, 0, 0]], 1, ZZ)
        (2, [([[1], [-1, 0]], 1), ([[1], [1, 0]], 1)])

    In result we got the following factorization::

                    f = 2 (x - y) (x + y)

    **References**

    1. [Gathen99]_

    """
    if not u:
        return dup_zz_factor(f, K)

    if dmp_zero_p(f, u):
        return K.zero, []

    cont, g = dmp_ground_primitive(f, u, K)

    if dmp_ground_LC(g, u, K) < 0:
        cont, g = -cont, dmp_neg(g, u, K)

    if all([ d <= 0 for d in dmp_degree_list(g, u) ]):
        return cont, []

    G, g = dmp_primitive(g, u, K)

    factors = []

    if dmp_degree(g, u) > 0:
        g = dmp_sqf_part(g, u, K)
        H = dmp_zz_wang(g, u, K)

        for h in H:
            k = 0

            while True:
                q, r = dmp_div(f, h, u, K)

                if dmp_zero_p(r, u):
                    f, k = q, k+1
                else:
                    break

            factors.append((h, k))

    for g, k in dmp_zz_factor(G, u-1, K)[1]:
        factors.insert(0, ([g], k))

    return cont, _sort_factors(factors)

#.........这里部分代码省略.........
    if mod is None:
        if u == 1:
            mod = 2
        else:
            mod = 1

    history, configs, A, r = set([]), [], [K.zero]*u, None

    try:
        cs, s, E = dmp_zz_wang_test_points(f, T, ct, A, u, K)

        _, H = dup_zz_factor_sqf(s, K)

        r = len(H)

        if r == 1:
            return [f]

        bad_points = set([tuple(A)])
        configs = [(s, cs, E, H, A)]
    except EvaluationFailed:
        pass

    eez_num_configs = query('EEZ_NUMBER_OF_CONFIGS')
    eez_num_tries = query('EEZ_NUMBER_OF_TRIES')
    eez_mod_step = query('EEZ_MODULUS_STEP')

    while len(configs) < eez_num_configs:
        for _ in xrange(eez_num_tries):
            A = [ K(randint(-mod, mod)) for _ in xrange(u) ]

            if tuple(A) not in history:
                history.add(tuple(A))
            else:
                continue

            try:
                cs, s, E = dmp_zz_wang_test_points(f, T, ct, A, u, K)
            except EvaluationFailed:
                continue

            _, H = dup_zz_factor_sqf(s, K)

            rr = len(H)

            if r is not None:
                if rr != r: # pragma: no cover
                    if rr < r:
                        configs, r = [], rr
                    else:
                        continue
            else:
                r = rr

            if r == 1:
                return [f]

            configs.append((s, cs, E, H, A))

            if len(configs) == eez_num_configs:
                break
        else:
            mod += eez_mod_step

    s_norm, s_arg, i = None, 0, 0

    for s, _, _, _, _ in configs:
        _s_norm = dup_max_norm(s, K)

        if s_norm is not None:
            if _s_norm < s_norm:
                s_norm = _s_norm
                s_arg = i
        else:
            s_norm = _s_norm

        i += 1

    _, cs, E, H, A = configs[s_arg]

    try:
        f, H, LC = dmp_zz_wang_lead_coeffs(f, T, cs, E, H, A, u, K)
        factors = dmp_zz_wang_hensel_lifting(f, H, LC, A, p, u, K)
    except ExtraneousFactors: # pragma: no cover
        if query('EEZ_RESTART_IF_NEEDED'):
            return dmp_zz_wang(f, u, K, mod+1)
        else:
            raise ExtraneousFactors("we need to restart algorithm with better parameters")

    negative, result = 0, []

    for f in factors:
        _, f = dmp_ground_primitive(f, u, K)

        if K.is_negative(dmp_ground_LC(f, u, K)):
            f = dmp_neg(f, u, K)

        result.append(f)

    return result

def dmp_prs_resultant(f, g, u, K):
    """
    Resultant algorithm in ``K[X]`` using subresultant PRS.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.euclidtools import dmp_prs_resultant

    >>> f = ZZ.map([[3, 0], [], [-1, 0, 0, -4]])
    >>> g = ZZ.map([[1], [1, 0, 0, 0], [-9]])

    >>> a = ZZ.map([[3, 0, 0, 0, 0], [1, 0, -27, 4]])
    >>> b = ZZ.map([[-3, 0, 0, -12, 1, 0, -54, 8, 729, -216, 16]])

    >>> dmp_prs_resultant(f, g, 1, ZZ) == (b[0], [f, g, a, b])
    True

    """
    if not u:
        return dup_prs_resultant(f, g, K)

    if dmp_zero_p(f, u) or dmp_zero_p(g, u):
        return (dmp_zero(u-1), [])

    R, B, D = dmp_inner_subresultants(f, g, u, K)

    if dmp_degree(R[-1], u) > 0:
        return (dmp_zero(u-1), R)
    if dmp_one_p(R[-2], u, K):
        return (dmp_LC(R[-1], K), R)

    s, i, v = 1, 1, u-1

    p = dmp_one(v, K)
    q = dmp_one(v, K)

    for b, d in zip(B, D)[:-1]:
        du = dmp_degree(R[i-1], u)
        dv = dmp_degree(R[i  ], u)
        dw = dmp_degree(R[i+1], u)

        if du % 2 and dv % 2:
            s = -s

        lc, i = dmp_LC(R[i], K), i+1

        p = dmp_mul(dmp_mul(p, dmp_pow(b, dv, v, K), v, K),
                               dmp_pow(lc, du-dw, v, K), v, K)
        q = dmp_mul(q, dmp_pow(lc, dv*(1+d), v, K), v, K)

        _, p, q = dmp_inner_gcd(p, q, v, K)

    if s < 0:
        p = dmp_neg(p, v, K)

    i = dmp_degree(R[-2], u)

    res = dmp_pow(dmp_LC(R[-1], K), i, v, K)
    res = dmp_exquo(dmp_mul(res, p, v, K), q, v, K)

    return res, R