def test_dup_mul_ground():
    f = dup_normal([], ZZ)

    assert dup_mul_ground(f, ZZ(2), ZZ) == dup_normal([], ZZ)

    f = dup_normal([1,2,3], ZZ)

    assert dup_mul_ground(f, ZZ(0), ZZ) == dup_normal([], ZZ)
    assert dup_mul_ground(f, ZZ(2), ZZ) == dup_normal([2,4,6], ZZ)

def dup_legendre(n, K):
    """Low-level implementation of Legendre polynomials. """
    seq = [[K.one], [K.one, K.zero]]

    for i in range(2, n + 1):
        a = dup_mul_ground(dup_lshift(seq[-1], 1, K), K(2*i - 1, i), K)
        b = dup_mul_ground(seq[-2], K(i - 1, i), K)

        seq.append(dup_sub(a, b, K))

    return seq[n]

def dup_gegenbauer(n, a, K):
    """Low-level implementation of Gegenbauer polynomials. """
    seq = [[K.one], [K(2)*a, K.zero]]

    for i in range(2, n + 1):
        f1 = K(2) * (i + a - K.one) / i
        f2 = (i + K(2)*a - K(2)) / i
        p1 = dup_mul_ground(dup_lshift(seq[-1], 1, K), f1, K)
        p2 = dup_mul_ground(seq[-2], f2, K)
        seq.append(dup_sub(p1, p2, K))

    return seq[n]

def dup_hermite(n, K):
    """Low-level implementation of Hermite polynomials. """
    seq = [[K.one], [K(2), K.zero]]

    for i in range(2, n + 1):
        a = dup_lshift(seq[-1], 1, K)
        b = dup_mul_ground(seq[-2], K(i - 1), K)

        c = dup_mul_ground(dup_sub(a, b, K), K(2), K)

        seq.append(c)

    return seq[n]

def dup_jacobi(n, a, b, K):
    """Low-level implementation of Jacobi polynomials. """
    seq = [[K.one], [(a + b + K(2))/K(2), (a - b)/K(2)]]

    for i in range(2, n + 1):
        den = K(i)*(a + b + i)*(a + b + K(2)*i - K(2))
        f0 = (a + b + K(2)*i - K.one) * (a*a - b*b) / (K(2)*den)
        f1 = (a + b + K(2)*i - K.one) * (a + b + K(2)*i - K(2)) * (a + b + K(2)*i) / (K(2)*den)
        f2 = (a + i - K.one)*(b + i - K.one)*(a + b + K(2)*i) / den
        p0 = dup_mul_ground(seq[-1], f0, K)
        p1 = dup_mul_ground(dup_lshift(seq[-1], 1, K), f1, K)
        p2 = dup_mul_ground(seq[-2], f2, K)
        seq.append(dup_sub(dup_add(p0, p1, K), p2, K))

    return seq[n]

def dup_sqf_list_include(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_include

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

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

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

    """
    coeff, factors = dup_sqf_list(f, K, all=all)

    if factors and factors[0][1] == 1:
        g = dup_mul_ground(factors[0][0], coeff, K)
        return [(g, 1)] + factors[1:]
    else:
        g = dup_strip([coeff])
        return [(g, 1)] + factors

def dup_transform(f, p, q, K):
    """
    Evaluate functional transformation ``q**n * f(p/q)`` in ``K[x]``.

    Examples
    ========

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

    >>> R.dup_transform(x**2 - 2*x + 1, x**2 + 1, x - 1)
    x**4 - 2*x**3 + 5*x**2 - 4*x + 4

    """
    if not f:
        return []

    n = len(f) - 1
    h, Q = [f[0]], [[K.one]]

    for i in range(0, n):
        Q.append(dup_mul(Q[-1], q, K))

    for c, q in zip(f[1:], Q[1:]):
        h = dup_mul(h, p, K)
        q = dup_mul_ground(q, c, K)
        h = dup_add(h, q, K)

    return h

def dup_transform(f, p, q, K):
    """
    Evaluate functional transformation ``q**n * f(p/q)`` in ``K[x]``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densetools import dup_transform

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

    >>> dup_transform(f, p, q, ZZ)
    [1, -2, 5, -4, 4]

    """
    if not f:
        return []

    n = dup_degree(f)
    h, Q = [f[0]], [[K.one]]

    for i in xrange(0, n):
        Q.append(dup_mul(Q[-1], q, K))

    for c, q in zip(f[1:], Q[1:]):
        h = dup_mul(h, p, K)
        q = dup_mul_ground(q, c, K)
        h = dup_add(h, q, K)

    return h

def dup_sqf_list_include(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_include(f)
    [(2, 1), (x + 1, 2), (x + 2, 3)]
    >>> R.dup_sqf_list_include(f, all=True)
    [(2, 1), (x + 1, 2), (x + 2, 3)]

    """
    coeff, factors = dup_sqf_list(f, K, all=all)

    if factors and factors[0][1] == 1:
        g = dup_mul_ground(factors[0][0], coeff, K)
        return [(g, 1)] + factors[1:]
    else:
        g = dup_strip([coeff])
        return [(g, 1)] + factors

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.domains import QQ
    >>> from sympy.polys.densetools import dup_revert

    >>> f = [-QQ(1,720), QQ(0), QQ(1,24), QQ(0), -QQ(1,2), QQ(0), QQ(1)]

    >>> dup_revert(f, 8, QQ)
    [61/720, 0/1, 5/24, 0/1, 1/2, 0/1, 1/1]

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

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

    for i in xrange(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

def dup_clear_denoms(f, K0, K1=None, convert=False):
    """
    Clear denominators, i.e. transform ``K_0`` to ``K_1``.

    Examples
    ========

    >>> from sympy.polys.domains import QQ, ZZ
    >>> from sympy.polys.densetools import dup_clear_denoms

    >>> f = [QQ(1,2), QQ(1,3)]
    >>> dup_clear_denoms(f, QQ, convert=False)
    (6, [3/1, 2/1])

    >>> f = [QQ(1,2), QQ(1,3)]
    >>> dup_clear_denoms(f, QQ, convert=True)
    (6, [3, 2])

    """
    if K1 is None:
        K1 = K0.get_ring()

    common = K1.one

    for c in f:
        common = K1.lcm(common, K0.denom(c))

    if not K1.is_one(common):
        f = dup_mul_ground(f, common, K0)

    if not convert:
        return common, f
    else:
        return common, dup_convert(f, K0, K1)

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

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_qq_heu_gcd(f, g, K0):
    """
    Heuristic polynomial GCD in `Q[x]`.

    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 QQ
    >>> from sympy.polys.euclidtools import dup_qq_heu_gcd

    >>> f = [QQ(1,2), QQ(7,4), QQ(3,2)]
    >>> g = [QQ(1,2), QQ(1), QQ(0)]

    >>> dup_qq_heu_gcd(f, g, QQ)
    ([1/1, 2/1], [1/2, 3/4], [1/2, 0/1])

    """
    result = _dup_ff_trivial_gcd(f, g, K0)

    if result is not None:
        return result

    K1 = K0.get_ring()

    cf, f = dup_clear_denoms(f, K0, K1)
    cg, g = dup_clear_denoms(g, K0, K1)

    f = dup_convert(f, K0, K1)
    g = dup_convert(g, K0, K1)

    h, cff, cfg = dup_zz_heu_gcd(f, g, K1)

    h = dup_convert(h, K1, K0)

    c = dup_LC(h, K0)
    h = dup_monic(h, K0)

    cff = dup_convert(cff, K1, K0)
    cfg = dup_convert(cfg, K1, K0)

    cff = dup_mul_ground(cff, K0.quo(c, cf), K0)
    cfg = dup_mul_ground(cfg, K0.quo(c, cg), K0)

    return h, cff, cfg

def dup_qq_heu_gcd(f, g, K0):
    """
    Heuristic polynomial GCD in `Q[x]`.

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

    >>> f = QQ(1,2)*x**2 + QQ(7,4)*x + QQ(3,2)
    >>> g = QQ(1,2)*x**2 + x

    >>> R.dup_qq_heu_gcd(f, g)
    (x + 2, 1/2*x + 3/4, 1/2*x)

    """
    result = _dup_ff_trivial_gcd(f, g, K0)

    if result is not None:
        return result

    K1 = K0.get_ring()

    cf, f = dup_clear_denoms(f, K0, K1)
    cg, g = dup_clear_denoms(g, K0, K1)

    f = dup_convert(f, K0, K1)
    g = dup_convert(g, K0, K1)

    h, cff, cfg = dup_zz_heu_gcd(f, g, K1)

    h = dup_convert(h, K1, K0)

    c = dup_LC(h, K0)
    h = dup_monic(h, K0)

    cff = dup_convert(cff, K1, K0)
    cfg = dup_convert(cfg, K1, K0)

    cff = dup_mul_ground(cff, K0.quo(c, cf), K0)
    cfg = dup_mul_ground(cfg, K0.quo(c, cg), K0)

    return h, cff, cfg

def dup_factor_list(f, K0, **args):
    """Factor polynomials into irreducibles in `K[x]`. """
    if not K0.has_CharacteristicZero: # pragma: no cover
        raise DomainError('only characteristic zero allowed')

    if 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_ground_to_ring(f, K0, K)
            f = dup_convert(f, K0, K)
        else:
            K = K0

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

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

            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 not args.get('include', False):
        return coeff, factors
    else:
        if not factors:
            return [(dup_strip([coeff]), 1)]
        else:
            g = dup_mul_ground(factors[0][0], coeff, K)
            return [(g, factors[0][1])] + factors[1:]

def dup_factor_list_include(f, K):
    """Factor polynomials into irreducibles in `K[x]`. """
    coeff, factors = dup_factor_list(f, K)

    if not factors:
        return [(dup_strip([coeff]), 1)]
    else:
        g = dup_mul_ground(factors[0][0], coeff, K)
        return [(g, factors[0][1])] + factors[1:]

def dup_chebyshevt(n, K):
    """Low-level implementation of Chebyshev polynomials of the 1st kind. """
    seq = [[K.one], [K.one, K.zero]]

    for i in range(2, n + 1):
        a = dup_mul_ground(dup_lshift(seq[-1], 1, K), K(2), K)
        seq.append(dup_sub(a, seq[-2], K))

    return seq[n]

def dup_spherical_bessel_fn_minus(n, K):
    """ Low-level implementation of fn(-n, x) """
    seq = [[K.one, K.zero], [K.zero]]

    for i in range(2, n + 1):
        a = dup_mul_ground(dup_lshift(seq[-1], 1, K), K(3 - 2*i), K)
        seq.append(dup_sub(a, seq[-2], K))

    return seq[n]

def dup_zz_hensel_lift(p, f, f_list, l, K):
    """
    Multifactor Hensel lifting in `Z[x]`.

    Given a prime `p`, polynomial `f` over `Z[x]` such that `lc(f)`
    is a unit modulo `p`, monic pair-wise coprime polynomials `f_i`
    over `Z[x]` satisfying::

        f = lc(f) f_1 ... f_r (mod p)

    and a positive integer `l`, returns a list of monic polynomials
    `F_1`, `F_2`, ..., `F_r` satisfying::

       f = lc(f) F_1 ... F_r (mod p**l)

       F_i = f_i (mod p), i = 1..r

    References
    ==========

    1. [Gathen99]_

    """
    r = len(f_list)
    lc = dup_LC(f, K)

    if r == 1:
        F = dup_mul_ground(f, K.gcdex(lc, p**l)[0], K)
        return [ dup_trunc(F, p**l, K) ]

    m = p
    k = r // 2
    d = int(_ceil(_log(l, 2)))

    g = gf_from_int_poly([lc], p)

    for f_i in f_list[:k]:
        g = gf_mul(g, gf_from_int_poly(f_i, p), p, K)

    h = gf_from_int_poly(f_list[k], p)

    for f_i in f_list[k + 1:]:
        h = gf_mul(h, gf_from_int_poly(f_i, p), p, K)

    s, t, _ = gf_gcdex(g, h, p, K)

    g = gf_to_int_poly(g, p)
    h = gf_to_int_poly(h, p)
    s = gf_to_int_poly(s, p)
    t = gf_to_int_poly(t, p)

    for _ in range(1, d + 1):
        (g, h, s, t), m = dup_zz_hensel_step(m, f, g, h, s, t, K), m**2

    return dup_zz_hensel_lift(p, g, f_list[:k], l, K) \
        + dup_zz_hensel_lift(p, h, f_list[k:], l, K)