def gf_berlekamp(f, p, K):
    """Factor a square-free `f` in `GF(p)[x]` for small `p`. """
    Q = gf_Qmatrix(f, p, K)
    V = gf_Qbasis(Q, p, K)

    for i, v in enumerate(V):
        V[i] = gf_strip(list(reversed(v)))

    factors = [f]

    for k in xrange(1, len(V)):
        for f in list(factors):
            s = K.zero

            while s < p:
                g = gf_sub_ground(V[k], s, p, K)
                h = gf_gcd(f, g, p, K)

                if h != [K.one] and h != f:
                    factors.remove(f)

                    f = gf_exquo(f, h, p, K)
                    factors.extend([f, h])

                if len(factors) == len(V):
                    return _sort_factors(factors, multiple=False)

                s += K.one

    return _sort_factors(factors, multiple=False)

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_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 gf_edf_shoup(f, n, p, K):
    """
    Gathen-Shoup: Probabilistic Equal Degree Factorization

    Given a monic square-free polynomial ``f`` in ``GF(p)[x]`` and integer
    ``n`` such that ``n`` divides ``deg(f)``, returns all irreducible factors
    ``f_1,...,f_d`` of ``f``, each of degree ``n``. This is a complete
    factorization over Galois fields.

    This algorithm is an improved version of Zassenhaus algorithm for
    large ``deg(f)`` and modulus ``p`` (especially for ``deg(f) ~ lg(p)``).

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.galoistools import gf_edf_shoup

    >>> gf_edf_shoup(ZZ.map([1, 2837, 2277]), 1, 2917, ZZ)
    [[1, 852], [1, 1985]]

    References
    ==========

    1. [Shoup91]_
    2. [Gathen92]_

    """
    N, q = gf_degree(f), int(p)

    if not N:
        return []
    if N <= n:
        return [f]

    factors, x = [f], [K.one, K.zero]

    r = gf_random(N - 1, p, K)

    h = gf_pow_mod(x, q, f, p, K)
    H = gf_trace_map(r, h, x, n - 1, f, p, K)[1]

    if p == 2:
        h1 = gf_gcd(f, H, p, K)
        h2 = gf_quo(f, h1, p, K)

        factors = gf_edf_shoup(h1, n, p, K) \
            + gf_edf_shoup(h2, n, p, K)
    else:
        h = gf_pow_mod(H, (q - 1)//2, f, p, K)

        h1 = gf_gcd(f, h, p, K)
        h2 = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)
        h3 = gf_quo(f, gf_mul(h1, h2, p, K), p, K)

        factors = gf_edf_shoup(h1, n, p, K) \
            + gf_edf_shoup(h2, n, p, K) \
            + gf_edf_shoup(h3, n, p, K)

    return _sort_factors(factors, multiple=False)

def gf_shoup(f, p, K):
    """Factor a square-free `f` in `GF(p)[x]` for large `p`. """
    factors = []

    for factor, n in gf_ddf_shoup(f, p, K):
        factors += gf_edf_shoup(factor, n, p, K)

    return _sort_factors(factors, multiple=False)

def gf_zassenhaus(f, p, K):
    """Factor a square-free `f` in `GF(p)[x]` for medium `p`. """
    factors = []

    for factor, n in gf_ddf_zassenhaus(f, p, K):
        factors += gf_edf_zassenhaus(factor, n, p, K)

    return _sort_factors(factors, multiple=False)

def gf_factor(f, p, K, **args):
    """Factor (non square-free) polynomials in `GF(p)[x]`.

       Given a possibly non square-free polynomial `f` in `GF(p)[x]`, returns
       its complete factorization into irreducibles::

                     f_1(x)**e_1 f_2(x)**e_2 ... f_d(x)**e_d

       where each `f_i` is a monic polynomial and `gcd(f_i, f_j) == 1`, for
       `i != j`.  The result is given as a tuple consisting of the leading
       coefficient of `f` and a list of factors with their multiplicities.

       The algorithm proceeds by first computing square-free decomposition
       of `f` and then iteratively factoring each of the square-free factors.

       Consider a non square-free polynomial `f = (7*x + 1) (x + 2)**2` in
       `GF(11)[x]`. We obtain its factorization into irreducibles as follows::

           >>> from sympy.polys.galoistools import gf_factor
           >>> from sympy.polys.algebratools import ZZ

           >>> gf_factor([5, 2, 7, 2], 11, ZZ)
           (5, [([1, 2], 1), ([1, 8], 2)])

       We arrived with factorization `f = 5 (x + 2) (x + 8)**2`. We didn't
       recover exact form of the input polynomial because we requested to
       get monic factors of `f` and its leading coefficient separately.

       Square-free factors of `f` can be factored into irreducibles over
       `GF(p)` using three very different methods:

       1. Berlekamp - efficient for very small values of `p` (usually `p < 25`)
       2. Cantor-Zassenhaus - efficient on average input and with "typical" `p`
       3. Shoup-Kaltofen-Gathen - efficient with very large inputs and modulus

       If you want to use a specific factorization method, instead of relying,
       on the algorithm to choose one for you, specify `method` keyword and
       set it to one of `berlekamp`, `zassenhaus` or `shoup` values.

       References
       ==========

       .. [Gathen99] J. von zur Gathen, J. Gerhard, Modern Computer Algebra,
           First Edition, Cambridge University Press, 1999, pp. 365

    """
    lc, f = gf_monic(f, p, K)

    if gf_degree(f) < 1:
        return lc, []

    factors = []

    for g, n in gf_sqf_list(f, p, K)[1]:
        for h in gf_factor_sqf(g, p, K, **args)[1]:
            factors.append((h, n))

    return lc, _sort_factors(factors)

def gf_edf_shoup(f, n, p, K):
    """Gathen-Shoup: Probabilistic Equal Degree Factorization

       Given a monic square-free polynomial `f` in `GF(p)[x]` and integer
       `n` such that `n` divides `deg(f)`, returns all irreducible factors
       `f_1 ... f_d` of `f`, each of degree `n`. This is a complete
       factorization over Galois fields.

       This algorithm is an improved version of Zassenhaus algorithm for
       large `deg(f)` and modulus `p` (especially for `deg(f) ~ lg(p)`).

       References
       ==========

       .. [Shoup91] V. Shoup, A Fast Deterministic Algorithm for Factoring
           Polynomials over Finite Fields of Small Characteristic, In
           Proceedings of International Symposium on Symbolic and
           Algebraic Computation, 1991, pp. 14-21

       .. [Gathen92] J. von zur Gathen, V. Shoup, Computing Frobenius Maps
           and Factoring Polynomials, ACM Symposium on Theory of Computing,
           1992, pp. 187-224

    """
    N, q = gf_degree(f), int(p)

    if not N:
        return []
    if N <= n:
        return [f]

    factors, x = [f], [K.one, K.zero]

    r = gf_random(N-1, p, K)

    h = gf_pow_mod(x, q, f, p, K)
    H = gf_trace_map(r, h, x, n-1, f, p, K)[1]

    if p == 2:
        h1 = gf_gcd(f, H, p, K)
        h2 = gf_exquo(f, h1, p, K)

        factors = gf_edf_shoup(h1, n, p, K) \
                + gf_edf_shoup(h2, n, p, K)
    else:
        h = gf_pow_mod(H, (q-1)//2, f, p, K)

        h1 = gf_gcd(f, h, p, K)
        h2 = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)
        h3 = gf_exquo(f, gf_mul(h1, h2, p, K), p, K)

        factors = gf_edf_shoup(h1, n, p, K) \
                + gf_edf_shoup(h2, n, p, K) \
                + gf_edf_shoup(h3, n, p, K)

    return _sort_factors(factors, multiple=False)

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

    if not K0.has_CharacteristicZero:
        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.has_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.has_Field:
            for i, (f, k) in enumerate(factors):
                factors[i] = (dup_convert(f, K, K0), k)

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

            coeff = K0.quo(coeff, denom)

        if K0_inexact is not None:
            for i, (f, k) in enumerate(factors):
                factors[i] = (dup_convert(f, K0, K0_inexact), k)

            coeff = K0_inexact.convert(coeff, K0)

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

    return coeff, _sort_factors(factors)

def gf_edf_zassenhaus(f, n, p, K):
    """Cantor-Zassenhaus: Probabilistic Equal Degree Factorization

       Given a monic square-free polynomial `f` in `GF(p)[x]` and integer
       `n` such that `n` divides `deg(f)`, returns all irreducible factors
       `f_1 ... f_d` of `f`, each of degree `n`. This is a complete
       factorization in Galois fields.

       Consider square-free polynomial `f = x**3 + x**2 + x + 1` in
       `GF(5)[x]`. Lets compute its irreducible factors of degree one::

           >>> from sympy.polys.galoistools import gf_edf_zassenhaus
           >>> from sympy.polys.algebratools import ZZ

           >>> gf_edf_zassenhaus([1,1,1,1], 1, 5, ZZ)
           [[1, 1], [1, 2], [1, 3]]

       References
       ==========

       .. [Gathen99] J. von zur Gathen, J. Gerhard, Modern Computer Algebra,
           First Edition, Cambridge University Press, 1999, pp. 358

       .. [Geddes92] K. Geddes, S. Czapor, G. Labahn, Algorithms for Computer
           Algebra, First Edition, Springer, 1992, pp. 371-373

    """
    factors, q = [f], int(p)

    if gf_degree(f) <= n:
        return factors

    N = gf_degree(f) // n

    while len(factors) < N:
        r = gf_random(2*n-1, p, K)

        if p == 2:
            h = r

            for i in xrange(0, 2**(n*N-1)):
                r = gf_pow_mod(r, 2, f, p, K)
                h = gf_add(h, r, p, K)

            g = gf_gcd(f, h, p, K)
        else:
            h = gf_pow_mod(r, (q**n-1) // 2, f, p, K)
            g = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)

        if g != [K.one] and g != f:
            factors = gf_edf_zassenhaus(g, n, p, K) \
                    + gf_edf_zassenhaus(gf_exquo(f, g, p, K), n, p, K)

    return _sort_factors(factors, multiple=False)

def gf_edf_zassenhaus(f, n, p, K):
    """
    Cantor-Zassenhaus: Probabilistic Equal Degree Factorization

    Given a monic square-free polynomial ``f`` in ``GF(p)[x]`` and
    an integer ``n``, such that ``n`` divides ``deg(f)``, returns all
    irreducible factors ``f_1,...,f_d`` of ``f``, each of degree ``n``.
    EDF procedure gives complete factorization over Galois fields.

    Consider the square-free polynomial ``f = x**3 + x**2 + x + 1`` in
    ``GF(5)[x]``. Let's compute its irreducible factors of degree one::

       >>> from sympy.polys.domains import ZZ
       >>> from sympy.polys.galoistools import gf_edf_zassenhaus

       >>> gf_edf_zassenhaus([1,1,1,1], 1, 5, ZZ)
       [[1, 1], [1, 2], [1, 3]]

    References
    ==========

    1. [Gathen99]_
    2. [Geddes92]_

    """
    factors, q = [f], int(p)

    if gf_degree(f) <= n:
        return factors

    N = gf_degree(f) // n

    while len(factors) < N:
        r = gf_random(2*n - 1, p, K)

        if p == 2:
            h = r

            for i in xrange(0, 2**(n*N - 1)):
                r = gf_pow_mod(r, 2, f, p, K)
                h = gf_add(h, r, p, K)

            g = gf_gcd(f, h, p, K)
        else:
            h = gf_pow_mod(r, (q**n - 1) // 2, f, p, K)
            g = gf_gcd(f, gf_sub_ground(h, K.one, p, K), p, K)

        if g != [K.one] and g != f:
            factors = gf_edf_zassenhaus(g, n, p, K) \
                + gf_edf_zassenhaus(gf_quo(f, g, p, K), n, p, K)

    return _sort_factors(factors, multiple=False)

def test__sort_factors():
    assert _sort_factors([], multiple=True) == []
    assert _sort_factors([], multiple=False) == []

    F = [[1, 2, 3], [1, 2], [1]]
    G = [[1], [1, 2], [1, 2, 3]]

    assert _sort_factors(F, multiple=False) == G

    F = [[1, 2], [1, 2, 3], [1, 2], [1]]
    G = [[1], [1, 2], [1, 2], [1, 2, 3]]

    assert _sort_factors(F, multiple=False) == G

    F = [[2, 2], [1, 2, 3], [1, 2], [1]]
    G = [[1], [1, 2], [2, 2], [1, 2, 3]]

    assert _sort_factors(F, multiple=False) == G

    F = [([1, 2, 3], 1), ([1, 2], 1), ([1], 1)]
    G = [([1], 1), ([1, 2], 1), ([1, 2, 3], 1)]

    assert _sort_factors(F, multiple=True) == G

    F = [([1, 2], 1), ([1, 2, 3], 1), ([1, 2], 1), ([1], 1)]
    G = [([1], 1), ([1, 2], 1), ([1, 2], 1), ([1, 2, 3], 1)]

    assert _sort_factors(F, multiple=True) == G

    F = [([2, 2], 1), ([1, 2, 3], 1), ([1, 2], 1), ([1], 1)]
    G = [([1], 1), ([1, 2], 1), ([2, 2], 1), ([1, 2, 3], 1)]

    assert _sort_factors(F, multiple=True) == G

    F = [([2, 2], 1), ([1, 2, 3], 1), ([1, 2], 2), ([1], 1)]
    G = [([1], 1), ([2, 2], 1), ([1, 2], 2), ([1, 2, 3], 1)]

    assert _sort_factors(F, multiple=True) == G

def dup_trial_division(f, factors, K):
    """Determine multiplicities of factors using trial division. """
    result = []

    for factor in factors:
        k = 0

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

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

        result.append((factor, k))

    return _sort_factors(result)

def gf_shoup(f, p, K):
    """
    Factor a square-free ``f`` in ``GF(p)[x]`` for large ``p``.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.galoistools import gf_shoup

    >>> gf_shoup([1, 4, 3], 5, ZZ)
    [[1, 1], [1, 3]]

    """
    factors = []

    for factor, n in gf_ddf_shoup(f, p, K):
        factors += gf_edf_shoup(factor, n, p, K)

    return _sort_factors(factors, multiple=False)

def dmp_trial_division(f, factors, u, K):
    """Determine multiplicities of factors using trial division. """
    result = []

    for factor in factors:
        k = 0

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

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

        result.append((factor, k))

    return _sort_factors(result)
    return sort_factors_if_mult(result)

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_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)

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 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 import ring, ZZ
        >>> R, x,y = ring("x,y", ZZ)

        >>> R.dmp_zz_factor(2*x**2 - 2*y**2)
        (2, [(x - y, 1), (x + y, 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)
        factors = dmp_trial_division(f, H, u, K)

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

    return cont, _sort_factors(factors)

def dmp_factor_list(f, u, K0):
    """Factor polynomials into irreducibles in `K[X]`. """
    if not u:
        return dup_factor_list(f, K0)

    J, f = dmp_terms_gcd(f, u, K0)

    if not K0.has_CharacteristicZero: # pragma: no cover
        coeff, factors = dmp_gf_factor(f, u, K0)
    elif K0.is_Algebraic:
        coeff, factors = dmp_ext_factor(f, u, K0)
    else:
        if not K0.is_Exact:
            K0_inexact, K0 = K0, K0.get_exact()
            f = dmp_convert(f, u, K0_inexact, K0)
        else:
            K0_inexact = None

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

            denom, f = dmp_clear_denoms(f, u, K0, K)
            f = dmp_convert(f, u, K0, K)
        else:
            K = K0

        if K.is_ZZ:
            levels, f, v = dmp_exclude(f, u, K)
            coeff, factors = dmp_zz_factor(f, v, K)

            for i, (f, k) in enumerate(factors):
                factors[i] = (dmp_include(f, levels, v, K), k)
        elif K.is_Poly:
            f, v = dmp_inject(f, u, K)

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

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

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

        if K0.has_Field:
            for i, (f, k) in enumerate(factors):
                factors[i] = (dmp_convert(f, u, K, K0), k)

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

            coeff = K0.quo(coeff, denom)

        if K0_inexact is not None:
            for i, (f, k) in enumerate(factors):
                factors[i] = (dmp_convert(f, u, K0, K0_inexact), k)

            coeff = K0_inexact.convert(coeff, K0)

    for i, j in enumerate(reversed(J)):
        if not j:
            continue

        term = {(0,)*(u-i) + (1,) + (0,)*i: K0.one}
        factors.insert(0, (dmp_from_dict(term, u, K0), j))

    return coeff, _sort_factors(factors)