def dup_rr_lcm(f, g, K):
    """
    Computes polynomial LCM over a ring in ``K[x]``.

    **Examples**

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

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

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

    """
    fc, f = dup_primitive(f, K)
    gc, g = dup_primitive(g, K)

    c = K.lcm(fc, gc)

    h = dup_exquo(dup_mul(f, g, K),
                  dup_gcd(f, g, K), K)

    return dup_mul_ground(h, c, K)

def dup_zz_factor_sqf(f, K):
    """Factor square-free (non-primitive) polyomials in `Z[x]`. """
    cont, g = dup_primitive(f, K)

    n = dup_degree(g)

    if dup_LC(g, K) < 0:
        cont, g = -cont, dup_neg(g, K)

    if n <= 0:
        return cont, []
    elif n == 1:
        return cont, [(g, 1)]

    if query('USE_IRREDUCIBLE_IN_FACTOR'):
        if dup_zz_irreducible_p(g, K):
            return cont, [(g, 1)]

    factors = None

    if query('USE_CYCLOTOMIC_FACTOR'):
        factors = dup_zz_cyclotomic_factor(g, K)

    if factors is None:
        factors = dup_zz_zassenhaus(g, K)

    return cont, _sort_factors(factors, multiple=False)

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

    Examples
    ========

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

    >>> R.dup_sqf_part(x**3 - 3*x - 2)
    x**2 - x - 2

    """
    if K.is_FiniteField:
        return dup_gf_sqf_part(f, K)

    if not f:
        return f

    if K.is_negative(dup_LC(f, K)):
        f = dup_neg(f, K)

    gcd = dup_gcd(f, dup_diff(f, 1, K), K)
    sqf = dup_quo(f, gcd, K)

    if K.has_Field:
        return dup_monic(sqf, K)
    else:
        return dup_primitive(sqf, K)[1]

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

    Examples
    ========

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

    >>> dup_sqf_part([ZZ(1), ZZ(0), -ZZ(3), -ZZ(2)], ZZ)
    [1, -1, -2]

    """
    if not K.has_CharacteristicZero:
        return dup_gf_sqf_part(f, K)

    if not f:
        return f

    if K.is_negative(dup_LC(f, K)):
        f = dup_neg(f, K)

    gcd = dup_gcd(f, dup_diff(f, 1, K), K)
    sqf = dup_quo(f, gcd, K)

    if K.has_Field or not K.is_Exact:
        return dup_monic(sqf, K)
    else:
        return dup_primitive(sqf, K)[1]

def dup_factor_list(f, K0):
    """Factor univariate polynomials into irreducibles in `K[x]`. """
    j, f = dup_terms_gcd(f, K0)
    cont, f = dup_primitive(f, K0)

    if K0.is_FiniteField:
        coeff, factors = dup_gf_factor(f, K0)
    elif K0.is_Algebraic:
        coeff, factors = dup_ext_factor(f, K0)
    else:
        if not K0.is_Exact:
            K0_inexact, K0 = K0, K0.get_exact()
            f = dup_convert(f, K0_inexact, K0)
        else:
            K0_inexact = None

        if K0.is_Field:
            K = K0.get_ring()

            denom, f = dup_clear_denoms(f, K0, K)
            f = dup_convert(f, K0, K)
        else:
            K = K0

        if K.is_ZZ:
            coeff, factors = dup_zz_factor(f, K)
        elif K.is_Poly:
            f, u = dmp_inject(f, 0, K)

            coeff, factors = dmp_factor_list(f, u, K.dom)

            for i, (f, k) in enumerate(factors):
                factors[i] = (dmp_eject(f, u, K), k)

            coeff = K.convert(coeff, K.dom)
        else:  # pragma: no cover
            raise DomainError('factorization not supported over %s' % K0)

        if K0.is_Field:
            for i, (f, k) in enumerate(factors):
                factors[i] = (dup_convert(f, K, K0), k)

            coeff = K0.convert(coeff, K)
            coeff = K0.quo(coeff, denom)

            if K0_inexact:
                for i, (f, k) in enumerate(factors):
                    max_norm = dup_max_norm(f, K0)
                    f = dup_quo_ground(f, max_norm, K0)
                    f = dup_convert(f, K0, K0_inexact)
                    factors[i] = (f, k)
                    coeff = K0.mul(coeff, K0.pow(max_norm, k))

                coeff = K0_inexact.convert(coeff, K0)
                K0 = K0_inexact

    if j:
        factors.insert(0, ([K0.one, K0.zero], j))

    return coeff*cont, _sort_factors(factors)

def dup_sqf_list(f, 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 dup_sqf_list

    >>> f = ZZ.map([2, 16, 50, 76, 56, 16])

    >>> dup_sqf_list(f, ZZ)
    (2, [([1, 1], 2), ([1, 2], 3)])

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

    """
    if not K.has_CharacteristicZero:
        return dup_gf_sqf_list(f, K, all=all)

    if K.has_Field or not K.is_Exact:
        coeff = dup_LC(f, K)
        f = dup_monic(f, K)
    else:
        coeff, f = dup_primitive(f, K)

        if K.is_negative(dup_LC(f, K)):
            f = dup_neg(f, K)
            coeff = -coeff

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

    result, i = [], 1

    h = dup_diff(f, 1, K)
    g, p, q = dup_inner_gcd(f, h, K)

    while True:
        d = dup_diff(p, 1, K)
        h = dup_sub(q, d, K)

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

        g, p, q = dup_inner_gcd(p, h, K)

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

        i += 1

    return coeff, result

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

    Examples
    ========

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

    >>> f = 2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16

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

    """
    if K.is_FiniteField:
        return dup_gf_sqf_list(f, K, all=all)

    if K.has_Field:
        coeff = dup_LC(f, K)
        f = dup_monic(f, K)
    else:
        coeff, f = dup_primitive(f, K)

        if K.is_negative(dup_LC(f, K)):
            f = dup_neg(f, K)
            coeff = -coeff

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

    result, i = [], 1

    h = dup_diff(f, 1, K)
    g, p, q = dup_inner_gcd(f, h, K)

    while True:
        d = dup_diff(p, 1, K)
        h = dup_sub(q, d, K)

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

        g, p, q = dup_inner_gcd(p, h, K)

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

        i += 1

    return coeff, result

def dup_rr_prs_gcd(f, g, 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 dup_rr_prs_gcd

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

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

    """
    result = _dup_rr_trivial_gcd(f, g, K)

    if result is not None:
        return result

    fc, F = dup_primitive(f, K)
    gc, G = dup_primitive(g, K)

    c = K.gcd(fc, gc)

    h = dup_subresultants(F, G, K)[-1]
    _, h = dup_primitive(h, K)

    if K.is_negative(dup_LC(h, K)):
        c = -c

    h = dup_mul_ground(h, c, K)

    cff = dup_quo(f, h, K)
    cfg = dup_quo(g, h, K)

    return h, cff, cfg

def dup_rr_prs_gcd(f, g, 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 = ring("x", ZZ)

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

    """
    result = _dup_rr_trivial_gcd(f, g, K)

    if result is not None:
        return result

    fc, F = dup_primitive(f, K)
    gc, G = dup_primitive(g, K)

    c = K.gcd(fc, gc)

    h = dup_subresultants(F, G, K)[-1]
    _, h = dup_primitive(h, K)

    if K.is_negative(dup_LC(h, K)):
        c = -c

    h = dup_mul_ground(h, c, K)

    cff = dup_quo(f, h, K)
    cfg = dup_quo(g, h, K)

    return h, cff, cfg

def dup_rr_lcm(f, g, K):
    """
    Computes polynomial LCM over a ring in `K[x]`.

    Examples
    ========

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

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

    """
    fc, f = dup_primitive(f, K)
    gc, g = dup_primitive(g, K)

    c = K.lcm(fc, gc)

    h = dup_quo(dup_mul(f, g, K), dup_gcd(f, g, K), K)

    return dup_mul_ground(h, c, K)

def dup_primitive_prs(f, g, K):
    """
    Primitive polynomial remainder sequence (PRS) in `K[x]`.

    Examples
    ========

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

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

    >>> prs = dup_primitive_prs(f, g, ZZ)

    >>> prs[0]
    [1, 0, 1, 0, -3, -3, 8, 2, -5]
    >>> prs[1]
    [3, 0, 5, 0, -4, -9, 21]
    >>> prs[2]
    [-5, 0, 1, 0, -3]
    >>> prs[3]
    [13, 25, -49]
    >>> prs[4]
    [4663, -6150]
    >>> prs[5]
    [1]

    """
    prs = [f, g]
    _, h = dup_primitive(dup_prem(f, g, K), K)

    while h:
        prs.append(h)
        f, g = g, h
        _, h = dup_primitive(dup_prem(f, g, K), K)

    return prs

def dup_primitive_prs(f, g, K):
    """
    Primitive polynomial remainder sequence (PRS) in `K[x]`.

    Examples
    ========

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

    >>> f = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    >>> g = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21

    >>> prs = R.dup_primitive_prs(f, g)

    >>> prs[0]
    x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    >>> prs[1]
    3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    >>> prs[2]
    -5*x**4 + x**2 - 3
    >>> prs[3]
    13*x**2 + 25*x - 49
    >>> prs[4]
    4663*x - 6150
    >>> prs[5]
    1

    """
    prs = [f, g]
    _, h = dup_primitive(dup_prem(f, g, K), K)

    while h:
        prs.append(h)
        f, g = g, h
        _, h = dup_primitive(dup_prem(f, g, K), K)

    return prs

def dup_zz_factor_sqf(f, K, **args):
    """Factor square-free (non-primitive) polyomials in `Z[x]`. """
    cont, g = dup_primitive(f, K)

    n = dup_degree(g)

    if dup_LC(g, K) < 0:
        cont, g = -cont, dup_neg(g, K)

    if n <= 0:
        return cont, []

    if n == 1 or dup_zz_irreducible_p(g, K):
        return cont, [(g, 1)]

    factors = []

    if args.get('cyclotomic', True):
        factors = dup_zz_cyclotomic_factor(g, K)

    if factors is None:
        factors = dup_zz_zassenhaus(g, K)

    return cont, _sort_factors(factors, multiple=False)

def dmp_zz_wang_test_points(f, T, ct, A, u, K):
    """Wang/EEZ: Test evaluation points for suitability. """
    if not dmp_eval_tail(dmp_LC(f, K), A, u-1, K):
        raise EvaluationFailed('no luck')

    g = dmp_eval_tail(f, A, u, K)

    if not dup_sqf_p(g, K):
        raise EvaluationFailed('no luck')

    c, h = dup_primitive(g, K)

    if K.is_negative(dup_LC(h, K)):
        c, h = -c, dup_neg(h, K)

    v = u-1

    E = [ dmp_eval_tail(t, A, v, K) for t, _ in T ]
    D = dmp_zz_wang_non_divisors(E, c, ct, K)

    if D is not None:
        return c, h, E
    else:
        raise EvaluationFailed('no luck')

def test_dup_primitive():
    assert dup_primitive([], ZZ) == (ZZ(0), [])
    assert dup_primitive([ZZ(1)], ZZ) == (ZZ(1), [ZZ(1)])
    assert dup_primitive([ZZ(1), ZZ(1)], ZZ) == (ZZ(1), [ZZ(1), ZZ(1)])
    assert dup_primitive([ZZ(2), ZZ(2)], ZZ) == (ZZ(2), [ZZ(1), ZZ(1)])
    assert dup_primitive(
        [ZZ(1), ZZ(2), ZZ(1)], ZZ) == (ZZ(1), [ZZ(1), ZZ(2), ZZ(1)])
    assert dup_primitive(
        [ZZ(2), ZZ(4), ZZ(2)], ZZ) == (ZZ(2), [ZZ(1), ZZ(2), ZZ(1)])

    assert dup_primitive([], QQ) == (QQ(0), [])
    assert dup_primitive([QQ(1)], QQ) == (QQ(1), [QQ(1)])
    assert dup_primitive([QQ(1), QQ(1)], QQ) == (QQ(1), [QQ(1), QQ(1)])
    assert dup_primitive([QQ(2), QQ(2)], QQ) == (QQ(2), [QQ(1), QQ(1)])
    assert dup_primitive(
        [QQ(1), QQ(2), QQ(1)], QQ) == (QQ(1), [QQ(1), QQ(2), QQ(1)])
    assert dup_primitive(
        [QQ(2), QQ(4), QQ(2)], QQ) == (QQ(2), [QQ(1), QQ(2), QQ(1)])

    assert dup_primitive(
        [QQ(2, 3), QQ(4, 9)], QQ) == (QQ(2, 9), [QQ(3), QQ(2)])
    assert dup_primitive(
        [QQ(2, 3), QQ(4, 5)], QQ) == (QQ(2, 15), [QQ(5), QQ(6)])

def test_dmp_zz_wang():
    p = ZZ(nextprime(dmp_zz_mignotte_bound(w_1, 2, ZZ)))

    assert p == ZZ(6291469)

    t_1, k_1, e_1 = dmp_normal([[1],[]], 1, ZZ), 1, ZZ(-14)
    t_2, k_2, e_2 = dmp_normal([[1, 0]], 1, ZZ), 2, ZZ(3)
    t_3, k_3, e_3 = dmp_normal([[1],[ 1, 0]], 1, ZZ), 2, ZZ(-11)
    t_4, k_4, e_4 = dmp_normal([[1],[-1, 0]], 1, ZZ), 1, ZZ(-17)

    T = [t_1, t_2, t_3, t_4]
    K = [k_1, k_2, k_3, k_4]
    E = [e_1, e_2, e_3, e_4]

    T = zip(T, K)

    A = [ZZ(-14), ZZ(3)]

    S = dmp_eval_tail(w_1, A, 2, ZZ)
    cs, s = dup_primitive(S, ZZ)

    assert cs == 1 and s == S == \
        dup_normal([1036728, 915552, 55748, 105621, -17304, -26841, -644], ZZ)

    assert dmp_zz_wang_non_divisors(E, cs, 4, ZZ) == [7, 3, 11, 17]
    assert dup_sqf_p(s, ZZ) and dup_degree(s) == dmp_degree(w_1, 2)

    _, H = dup_zz_factor_sqf(s, ZZ)

    h_1 = dup_normal([44,  42,   1], ZZ)
    h_2 = dup_normal([126, -9,  28], ZZ)
    h_3 = dup_normal([187,  0, -23], ZZ)

    assert H == [h_1, h_2, h_3]

    lc_1 = dmp_normal([[-4], [-4,0]], 1, ZZ)
    lc_2 = dmp_normal([[-1,0,0], []], 1, ZZ)
    lc_3 = dmp_normal([[1], [], [-1,0,0]], 1, ZZ)

    LC = [lc_1, lc_2, lc_3]

    assert dmp_zz_wang_lead_coeffs(w_1, T, cs, E, H, A, 2, ZZ) == (w_1, H, LC)

    H_1 = [ dmp_normal(t, 0, ZZ) for t in [[44L,42L,1L],[126L,-9L,28L],[187L,0L,-23L]] ]
    H_2 = [ dmp_normal(t, 1, ZZ) for t in [[[-4,-12],[-3,0],[1]],[[-9,0],[-9],[-2,0]],[[1,0,-9],[],[1,-9]]] ]
    H_3 = [ dmp_normal(t, 1, ZZ) for t in [[[-4,-12],[-3,0],[1]],[[-9,0],[-9],[-2,0]],[[1,0,-9],[],[1,-9]]] ]

    c_1 = dmp_normal([-70686,-5863,-17826,2009,5031,74], 0, ZZ)
    c_2 = dmp_normal([[9,12,-45,-108,-324],[18,-216,-810,0],[2,9,-252,-288,-945],[-30,-414,0],[2,-54,-3,81],[12,0]], 1, ZZ)
    c_3 = dmp_normal([[-36,-108,0],[-27,-36,-108],[-8,-42,0],[-6,0,9],[2,0]], 1, ZZ)

    T_1 = [ dmp_normal(t, 0, ZZ) for t in [[-3,0],[-2],[1]] ]
    T_2 = [ dmp_normal(t, 1, ZZ) for t in [[[-1,0],[]],[[-3],[]],[[-6]]] ]
    T_3 = [ dmp_normal(t, 1, ZZ) for t in [[[]],[[]],[[-1]]] ]

    assert dmp_zz_diophantine(H_1, c_1,        [], 5, p, 0, ZZ) == T_1
    assert dmp_zz_diophantine(H_2, c_2, [ZZ(-14)], 5, p, 1, ZZ) == T_2
    assert dmp_zz_diophantine(H_3, c_3, [ZZ(-14)], 5, p, 1, ZZ) == T_3

    factors = dmp_zz_wang_hensel_lifting(w_1, H, LC, A, p, 2, ZZ)

    assert dmp_expand(factors, 2, ZZ) == w_1

def dup_zz_zassenhaus(f, K):
    """Factor primitive square-free polynomials in `Z[x]`. """
    n = dup_degree(f)

    if n == 1:
        return [f]

    fc = f[-1]
    A = dup_max_norm(f, K)
    b = dup_LC(f, K)
    B = int(abs(K.sqrt(K(n + 1))*2**n*A*b))
    C = int((n + 1)**(2*n)*A**(2*n - 1))
    gamma = int(_ceil(2*_log(C, 2)))
    bound = int(2*gamma*_log(gamma))
    a = []
    # choose a prime number `p` such that `f` be square free in Z_p
    # if there are many factors in Z_p, choose among a few different `p`
    # the one with fewer factors
    for px in range(3, bound + 1):
        if not isprime(px) or b % px == 0:
            continue

        px = K.convert(px)

        F = gf_from_int_poly(f, px)

        if not gf_sqf_p(F, px, K):
            continue
        fsqfx = gf_factor_sqf(F, px, K)[1]
        a.append((px, fsqfx))
        if len(fsqfx) < 15 or len(a) > 4:
            break
    p, fsqf = min(a, key=lambda x: len(x[1]))

    l = int(_ceil(_log(2*B + 1, p)))

    modular = [gf_to_int_poly(ff, p) for ff in fsqf]

    g = dup_zz_hensel_lift(p, f, modular, l, K)

    sorted_T = range(len(g))
    T = set(sorted_T)
    factors, s = [], 1
    pl = p**l

    while 2*s <= len(T):
        for S in subsets(sorted_T, s):
            # lift the constant coefficient of the product `G` of the factors
            # in the subset `S`; if it is does not divide `fc`, `G` does
            # not divide the input polynomial

            if b == 1:
                q = 1
                for i in S:
                    q = q*g[i][-1]
                q = q % pl
                if not _test_pl(fc, q, pl):
                    continue
            else:
                G = [b]
                for i in S:
                    G = dup_mul(G, g[i], K)
                G = dup_trunc(G, pl, K)
                G = dup_primitive(G, K)[1]
                q = G[-1]
                if q and fc % q != 0:
                    continue

            H = [b]
            S = set(S)
            T_S = T - S

            if b == 1:
                G = [b]
                for i in S:
                    G = dup_mul(G, g[i], K)
                G = dup_trunc(G, pl, K)

            for i in T_S:
                H = dup_mul(H, g[i], K)

            H = dup_trunc(H, pl, K)

            G_norm = dup_l1_norm(G, K)
            H_norm = dup_l1_norm(H, K)

            if G_norm*H_norm <= B:
                T = T_S
                sorted_T = [i for i in sorted_T if i not in S]

                G = dup_primitive(G, K)[1]
                f = dup_primitive(H, K)[1]

                factors.append(G)
                b = dup_LC(f, K)

                break
        else:
            s += 1

#.........这里部分代码省略.........

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

    Given a univariate 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
    Zassenhaus 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**4 - 2`::

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

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

    In result we got the following factorization::

                 f = 2 (x - 1) (x + 1) (x**2 + 1)

    Note that this is a complete factorization over integers,
    however over Gaussian integers we can factor the last term.

    By default, polynomials `x**n - 1` and `x**n + 1` are factored
    using cyclotomic decomposition to speedup computations. To
    disable this behaviour set cyclotomic=False.

    **References**

    1. [Gathen99]_

    """
    cont, g = dup_primitive(f, K)

    n = dup_degree(g)

    if dup_LC(g, K) < 0:
        cont, g = -cont, dup_neg(g, K)

    if n <= 0:
        return cont, []
    elif n == 1:
        return cont, [(g, 1)]

    if query('USE_IRREDUCIBLE_IN_FACTOR'):
        if dup_zz_irreducible_p(g, K):
            return cont, [(g, 1)]

    g = dup_sqf_part(g, K)
    H, factors = None, []

    if query('USE_CYCLOTOMIC_FACTOR'):
        H = dup_zz_cyclotomic_factor(g, K)

    if H is None:
        H = dup_zz_zassenhaus(g, K)

    for h in H:
        k = 0

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

            if not r:
                f, k = q, k+1
            else:
                break

        factors.append((h, k))

    return cont, _sort_factors(factors)

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

    Given a univariate 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
    Zassenhaus 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**4 - 2`::

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

        >>> R.dup_zz_factor(2*x**4 - 2)
        (2, [(x - 1, 1), (x + 1, 1), (x**2 + 1, 1)])

    In result we got the following factorization::

                 f = 2 (x - 1) (x + 1) (x**2 + 1)

    Note that this is a complete factorization over integers,
    however over Gaussian integers we can factor the last term.

    By default, polynomials `x**n - 1` and `x**n + 1` are factored
    using cyclotomic decomposition to speedup computations. To
    disable this behaviour set cyclotomic=False.

    References
    ==========

    1. [Gathen99]_

    """
    cont, g = dup_primitive(f, K)

    n = dup_degree(g)

    if dup_LC(g, K) < 0:
        cont, g = -cont, dup_neg(g, K)

    if n <= 0:
        return cont, []
    elif n == 1:
        return cont, [(g, 1)]

    if query('USE_IRREDUCIBLE_IN_FACTOR'):
        if dup_zz_irreducible_p(g, K):
            return cont, [(g, 1)]

    g = dup_sqf_part(g, K)
    H = None

    if query('USE_CYCLOTOMIC_FACTOR'):
        H = dup_zz_cyclotomic_factor(g, K)

    if H is None:
        H = dup_zz_zassenhaus(g, K)

    factors = dup_trial_division(f, H, K)
    return cont, factors

def dup_zz_zassenhaus(f, K):
    """Factor primitive square-free polynomials in `Z[x]`. """
    n = dup_degree(f)

    if n == 1:
        return [f]

    A = dup_max_norm(f, K)
    b = dup_LC(f, K)
    B = int(abs(K.sqrt(K(n+1))*2**n*A*b))
    C = int((n+1)**(2*n)*A**(2*n-1))
    gamma = int(ceil(2*log(C, 2)))
    bound = int(2*gamma*log(gamma))

    for p in xrange(3, bound+1):
        if not isprime(p) or b % p == 0:
            continue

        p = K.convert(p)

        F = gf_from_int_poly(f, p)

        if gf_sqf_p(F, p, K):
            break

    l = int(ceil(log(2*B + 1, p)))

    modular = []

    for ff in gf_factor_sqf(F, p, K)[1]:
        modular.append(gf_to_int_poly(ff, p))

    g = dup_zz_hensel_lift(p, f, modular, l, K)

    T = set(range(len(g)))
    factors, s = [], 1

    while 2*s <= len(T):
        for S in subsets(T, s):
            G, H = [b], [b]

            S = set(S)

            for i in S:
                G = dup_mul(G, g[i], K)
            for i in T-S:
                H = dup_mul(H, g[i], K)

            G = dup_trunc(G, p**l, K)
            H = dup_trunc(H, p**l, K)

            G_norm = dup_l1_norm(G, K)
            H_norm = dup_l1_norm(H, K)

            if G_norm*H_norm <= B:
                T = T - S

                G = dup_primitive(G, K)[1]
                f = dup_primitive(H, K)[1]

                factors.append(G)
                b = dup_LC(f, K)

                break
        else:
            s += 1

    return factors + [f]