def dup_gff_list(f, K):
    """
    Compute greatest factorial factorization of ``f`` in ``K[x]``.

    Examples
    ========

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

    >>> R.dup_gff_list(x**5 + 2*x**4 - x**3 - 2*x**2)
    [(x, 1), (x + 2, 4)]

    """
    if not f:
        raise ValueError("greatest factorial factorization doesn't exist for a zero polynomial")

    f = dup_monic(f, K)

    if not dup_degree(f):
        return []
    else:
        g = dup_gcd(f, dup_shift(f, K.one, K), K)
        H = dup_gff_list(g, K)

        for i, (h, k) in enumerate(H):
            g = dup_mul(g, dup_shift(h, -K(k), K), K)
            H[i] = (h, k + 1)

        f = dup_quo(f, g, K)

        if not dup_degree(f):
            return H
        else:
            return [(f, 1)] + H

def dup_prem(f, g, K):
    """Polynomial pseudo-remainder in `K[x]`. """
    df = dup_degree(f)
    dg = dup_degree(g)

    r = f

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return r

    N = df - dg + 1
    lc_g = dup_LC(g, K)

    while True:
        dr = dup_degree(r)

        if dr < dg:
            break

        lc_r = dup_LC(r, K)
        j, N = dr-dg, N-1

        R = dup_mul_ground(r, lc_g, K)
        G = dup_mul_term(g, lc_r, j, K)
        r = dup_sub(R, G, K)

    return dup_mul_ground(r, lc_g**N, K)

def dup_sub(f, g, K):
    """
    Subtract dense polynomials in ``K[x]``.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densearith import dup_sub

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

    >>> dup_sub(f, g, ZZ)
    [1, -1, 1]

    """
    if not f:
        return dup_neg(g, K)
    if not g:
        return f

    df = dup_degree(f)
    dg = dup_degree(g)

    if df == dg:
        return dup_strip([ a - b for a, b in zip(f, g) ])
    else:
        k = abs(df - dg)

        if df > dg:
            h, f = f[:k], f[k:]
        else:
            h, g = dup_neg(g[:k], K), g[k:]

        return h + [ a - b for a, b in zip(f, g) ]

def dup_ff_div(f, g, K):
    """Polynomial division with remainder over a field. """
    df = dup_degree(f)
    dg = dup_degree(g)

    q, r = [], f

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return q, r

    lc_g = dup_LC(g, K)

    while True:
        dr = dup_degree(r)

        if dr < dg:
            break

        lc_r = dup_LC(r, K)

        c = K.exquo(lc_r, lc_g)
        j = dr - dg

        q = dup_add_term(q, c, j, K)
        h = dup_mul_term(g, c, j, K)

        r = dup_sub(r, h, K)

        if not K.is_Exact:
            r = dup_normal(r, K)

    return q, r

def dup_rr_div(f, g, K):
    """Univariate division with remainder over a ring. """
    df = dup_degree(f)
    dg = dup_degree(g)

    q, r = [], f

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return q, r

    lc_g = dup_LC(g, K)

    while True:
        dr = dup_degree(r)

        if dr < dg:
            break

        lc_r = dup_LC(r, K)

        if lc_r % lc_g:
            break

        c = K.exquo(lc_r, lc_g)
        j = dr - dg

        q = dup_add_term(q, c, j, K)
        h = dup_mul_term(g, c, j, K)

        r = dup_sub(r, h, K)

    return q, r

def dup_add(f, g, K):
    """
    Add dense polynomials in ``K[x]``.

    Examples
    ========

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

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

    """
    if not f:
        return g
    if not g:
        return f

    df = dup_degree(f)
    dg = dup_degree(g)

    if df == dg:
        return dup_strip([ a + b for a, b in zip(f, g) ])
    else:
        k = abs(df - dg)

        if df > dg:
            h, f = f[:k], f[k:]
        else:
            h, g = g[:k], g[k:]

        return h + [ a + b for a, b in zip(f, g) ]

def dup_sub(f, g, K):
    """
    Subtract dense polynomials in ``K[x]``.

    Examples
    ========

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

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

    """
    if not f:
        return dup_neg(g, K)
    if not g:
        return f

    df = dup_degree(f)
    dg = dup_degree(g)

    if df == dg:
        return dup_strip([ a - b for a, b in zip(f, g) ])
    else:
        k = abs(df - dg)

        if df > dg:
            h, f = f[:k], f[k:]
        else:
            h, g = dup_neg(g[:k], K), g[k:]

        return h + [ a - b for a, b in zip(f, g) ]

def dup_mul(f, g, K):
    """
    Multiply dense polynomials in ``K[x]``.

    Examples
    ========

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

    >>> R.dup_mul(x - 2, x + 2)
    x**2 - 4

    """
    if f == g:
        return dup_sqr(f, K)

    if not (f and g):
        return []

    df = dup_degree(f)
    dg = dup_degree(g)

    h = []

    for i in xrange(0, df + dg + 1):
        coeff = K.zero

        for j in xrange(max(0, i - dg), min(df, i) + 1):
            coeff += f[j]*g[i - j]

        h.append(coeff)

    return dup_strip(h)

def dup_ext_factor(f, K):
    """Factor univariate polynomials over algebraic number fields. """
    n, lc = dup_degree(f), dup_LC(f, K)

    f = dup_monic(f, K)

    if n <= 0:
        return lc, []
    if n == 1:
        return lc, [(f, 1)]

    f, F = dup_sqf_part(f, K), f
    s, g, r = dup_sqf_norm(f, K)

    factors = dup_factor_list_include(r, K.dom)

    if len(factors) == 1:
        return lc, [(f, n//dup_degree(f))]

    H = s*K.unit

    for i, (factor, _) in enumerate(factors):
        h = dup_convert(factor, K.dom, K)
        h, _, g = dup_inner_gcd(h, g, K)
        h = dup_shift(h, H, K)
        factors[i] = h

    factors = dup_trial_division(F, factors, K)

    return lc, 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_mul(f, g, K):
    """
    Multiply dense polynomials in ``K[x]``.

    Examples
    ========

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

    >>> R.dup_mul(x - 2, x + 2)
    x**2 - 4

    """
    if f == g:
        return dup_sqr(f, K)

    if not (f and g):
        return []

    df = dup_degree(f)
    dg = dup_degree(g)

    n = max(df, dg) + 1

    if n < 100:
        h = []

        for i in xrange(0, df + dg + 1):
            coeff = K.zero

            for j in xrange(max(0, i - dg), min(df, i) + 1):
                coeff += f[j]*g[i - j]

            h.append(coeff)

        return dup_strip(h)
    else:
        # Use Karatsuba's algorithm (divide and conquer), see e.g.:
        # Joris van der Hoeven, Relax But Don't Be Too Lazy,
        # J. Symbolic Computation, 11 (2002), section 3.1.1.
        n2 = n//2

        fl, gl = dup_slice(f, 0, n2, K), dup_slice(g, 0, n2, K)

        fh = dup_rshift(dup_slice(f, n2, n, K), n2, K)
        gh = dup_rshift(dup_slice(g, n2, n, K), n2, K)

        lo, hi = dup_mul(fl, gl, K), dup_mul(fh, gh, K)

        mid = dup_mul(dup_add(fl, fh, K), dup_add(gl, gh, K), K)
        mid = dup_sub(mid, dup_add(lo, hi, K), K)

        return dup_add(dup_add(lo, dup_lshift(mid, n2, K), K),
                       dup_lshift(hi, 2*n2, K), K)

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_prs_resultant(f, g, K):
    """
    Resultant algorithm in `K[x]` using subresultant PRS.

    Examples
    ========

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

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

    >>> dup_prs_resultant(f, g, ZZ)
    (4, [[1, 0, 1], [1, 0, -1], [-2]])

    """
    if not f or not g:
        return (K.zero, [])

    R, B, D = dup_inner_subresultants(f, g, K)

    if dup_degree(R[-1]) > 0:
        return (K.zero, R)
    if R[-2] == [K.one]:
        return (dup_LC(R[-1], K), R)

    s, i = 1, 1
    p, q = K.one, K.one

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

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

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

        p *= b**dv * lc**(du - dw)
        q *= lc**(dv*(1 + d))

    if s < 0:
        p = -p

    i = dup_degree(R[-2])

    res = dup_LC(R[-1], K)**i

    res = K.quo(res*p, q)

    return res, R

def dup_pdiv(f, g, K):
    """
    Polynomial pseudo-division in ``K[x]``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densearith import dup_pdiv

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

    >>> dup_pdiv(f, g, ZZ)
    ([2, 4], [20])

    """
    df = dup_degree(f)
    dg = dup_degree(g)

    q, r = [], f

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return q, r

    N = df - dg + 1
    lc_g = dup_LC(g, K)

    while True:
        dr = dup_degree(r)

        if dr < dg:
            break

        lc_r = dup_LC(r, K)
        j, N = dr-dg, N-1

        Q = dup_mul_ground(q, lc_g, K)
        q = dup_add_term(Q, lc_r, j, K)

        R = dup_mul_ground(r, lc_g, K)
        G = dup_mul_term(g, lc_r, j, K)
        r = dup_sub(R, G, K)

    c = lc_g**N

    q = dup_mul_ground(q, c, K)
    r = dup_mul_ground(r, c, K)

    return q, r

def dup_pdiv(f, g, K):
    """
    Polynomial pseudo-division in ``K[x]``.

    Examples
    ========

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

    >>> R.dup_pdiv(x**2 + 1, 2*x - 4)
    (2*x + 4, 20)

    """
    df = dup_degree(f)
    dg = dup_degree(g)

    q, r, dr = [], f, df

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return q, r

    N = df - dg + 1
    lc_g = dup_LC(g, K)

    while True:
        lc_r = dup_LC(r, K)
        j, N = dr - dg, N - 1

        Q = dup_mul_ground(q, lc_g, K)
        q = dup_add_term(Q, lc_r, j, K)

        R = dup_mul_ground(r, lc_g, K)
        G = dup_mul_term(g, lc_r, j, K)
        r = dup_sub(R, G, K)

        _dr, dr = dr, dup_degree(r)

        if dr < dg:
            break
        elif not (dr < _dr):
            raise PolynomialDivisionFailed(f, g, K)

    c = lc_g**N

    q = dup_mul_ground(q, c, K)
    r = dup_mul_ground(r, c, K)

    return q, r

def dup_ff_div(f, g, K):
    """
    Polynomial division with remainder over a field.

    Examples
    ========

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

    >>> R.dup_ff_div(x**2 + 1, 2*x - 4)
    (1/2*x + 1, 5)

    """
    df = dup_degree(f)
    dg = dup_degree(g)

    q, r, dr = [], f, df

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return q, r

    lc_g = dup_LC(g, K)

    while True:
        lc_r = dup_LC(r, K)

        c = K.exquo(lc_r, lc_g)
        j = dr - dg

        q = dup_add_term(q, c, j, K)
        h = dup_mul_term(g, c, j, K)
        r = dup_sub(r, h, K)

        _dr, dr = dr, dup_degree(r)

        if dr < dg:
            break
        elif dr == _dr and not K.is_Exact:
            # remove leading term created by rounding error
            r = dup_strip(r[1:])
            dr = dup_degree(r)
            if dr < dg:
                break
        elif not (dr < _dr):
            raise PolynomialDivisionFailed(f, g, K)

    return q, r

def dup_prs_resultant(f, g, K):
    """
    Resultant algorithm in `K[x]` using subresultant PRS.

    Examples
    ========

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

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

    """
    if not f or not g:
        return (K.zero, [])

    R, B, D = dup_inner_subresultants(f, g, K)

    if dup_degree(R[-1]) > 0:
        return (K.zero, R)
    if R[-2] == [K.one]:
        return (dup_LC(R[-1], K), R)

    s, i = 1, 1
    p, q = K.one, K.one

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

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

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

        p *= b**dv * lc**(du - dw)
        q *= lc**(dv*(1 + d))

    if s < 0:
        p = -p

    i = dup_degree(R[-2])

    res = dup_LC(R[-1], K)**i

    res = K.quo(res*p, q)

    return res, R

def dup_rr_div(f, g, K):
    """
    Univariate division with remainder over a ring.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densearith import dup_rr_div

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

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

    """
    df = dup_degree(f)
    dg = dup_degree(g)

    q, r = [], f

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return q, r

    lc_g = dup_LC(g, K)

    while True:
        dr = dup_degree(r)

        if dr < dg:
            break

        lc_r = dup_LC(r, K)

        if lc_r % lc_g:
            break

        c = K.exquo(lc_r, lc_g)
        j = dr - dg

        q = dup_add_term(q, c, j, K)
        h = dup_mul_term(g, c, j, K)

        r = dup_sub(r, h, K)

    return q, r

def dup_ff_div(f, g, K):
    """
    Polynomial division with remainder over a field.

    Examples
    ========

    >>> from sympy.polys.domains import QQ
    >>> from sympy.polys.densearith import dup_ff_div

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

    >>> dup_ff_div(f, g, QQ)
    ([1/2, 1/1], [5/1])

    """
    df = dup_degree(f)
    dg = dup_degree(g)

    q, r = [], f

    if not g:
        raise ZeroDivisionError("polynomial division")
    elif df < dg:
        return q, r

    lc_g = dup_LC(g, K)

    while True:
        dr = dup_degree(r)

        if dr < dg:
            break

        lc_r = dup_LC(r, K)

        c = K.exquo(lc_r, lc_g)
        j = dr - dg

        q = dup_add_term(q, c, j, K)
        h = dup_mul_term(g, c, j, K)

        r = dup_sub(r, h, K)

        if not K.is_Exact:
            r = dup_normal(r, K)

    return q, r

def dup_revert(f, n, K):
    """
    Compute ``f**(-1)`` mod ``x**n`` using Newton iteration.

    This function computes first ``2**n`` terms of a polynomial that
    is a result of inversion of a polynomial modulo ``x**n``. This is
    useful to efficiently compute series expansion of ``1/f``.

    Examples
    ========

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

    >>> f = -QQ(1,720)*x**6 + QQ(1,24)*x**4 - QQ(1,2)*x**2 + 1

    >>> R.dup_revert(f, 8)
    61/720*x**6 + 5/24*x**4 + 1/2*x**2 + 1

    """
    g = [K.revert(dup_TC(f, K))]
    h = [K.one, K.zero, K.zero]

    N = int(_ceil(_log(n, 2)))

    for i in range(1, N + 1):
        a = dup_mul_ground(g, K(2), K)
        b = dup_mul(f, dup_sqr(g, K), K)
        g = dup_rem(dup_sub(a, b, K), h, K)
        h = dup_lshift(h, dup_degree(h), K)

    return g