def dmp_sub(f, g, u, K):
    """Subtract dense polynomials in `K[X]`. """
    if not u:
        return dup_sub(f, g, K)

    df = dmp_degree(f, u)

    if df < 0:
        return dmp_neg(g, u, K)

    dg = dmp_degree(g, u)

    if dg < 0:
        return f

    v = u-1

    if df == dg:
        return dmp_strip([ dmp_sub(a, b, v, K) for a, b in zip(f, g) ], u)
    else:
        k = abs(df - dg)

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

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

def dmp_mul(f, g, u, K):
    """Multiply dense polynomials in `K[X]`. """
    if not u:
        return dup_mul(f, g, K)

    if f == g:
        return dmp_sqr(f, u, K)

    df = dmp_degree(f, u)

    if df < 0:
        return f

    dg = dmp_degree(g, u)

    if dg < 0:
        return g

    h, v = [], u-1

    for i in xrange(0, df+dg+1):
        coeff = dmp_zero(v)

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

        h.append(coeff)

    return h

def dmp_prem(f, g, u, K):
    """Polynomial pseudo-remainder in `K[X]`. """
    if not u:
        return dup_prem(f, g, K)

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    r = f

    if df < dg:
        return r

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

    while True:
        dr = dmp_degree(r, u)

        if dr < dg:
            break

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

        R = dmp_mul_term(r, lc_g, 0, u, K)
        G = dmp_mul_term(g, lc_r, j, u, K)
        r = dmp_sub(R, G, u, K)

    c = dmp_pow(lc_g, N, u-1, K)

    return dmp_mul_term(r, c, 0, u, K)

def _dmp_simplify_gcd(f, g, u, K):
    """Try to eliminate ``x_0`` from GCD computation in ``K[X]``. """
    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if df > 0 and dg > 0:
        return None

    if not (df or dg):
        F = dmp_LC(f, K)
        G = dmp_LC(g, K)
    else:
        if not df:
            F = dmp_LC(f, K)
            G = dmp_content(g, u, K)
        else:
            F = dmp_content(f, u, K)
            G = dmp_LC(g, K)

    v = u - 1
    h = dmp_gcd(F, G, v, K)

    cff = [ dmp_exquo(cf, h, v, K) for cf in f ]
    cfg = [ dmp_exquo(cg, h, v, K) for cg in g ]

    return [h], cff, cfg

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 dmp_zz_collins_resultant(f, g, u, K):
    """
    Collins's modular resultant algorithm in `Z[X]`.

    Examples
    ========

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

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

    >>> R.dmp_zz_collins_resultant(f, g)
    -2*y**2 - 5*y + 1

    """

    n = dmp_degree(f, u)
    m = dmp_degree(g, u)

    if n < 0 or m < 0:
        return dmp_zero(u - 1)

    A = dmp_max_norm(f, u, K)
    B = dmp_max_norm(g, u, K)

    a = dmp_ground_LC(f, u, K)
    b = dmp_ground_LC(g, u, K)

    v = u - 1

    B = K(2)*K.factorial(K(n + m))*A**m*B**n
    r, p, P = dmp_zero(v), K.one, K.one

    while P <= B:
        p = K(nextprime(p))

        while not (a % p) or not (b % p):
            p = K(nextprime(p))

        F = dmp_ground_trunc(f, p, u, K)
        G = dmp_ground_trunc(g, p, u, K)

        try:
            R = dmp_zz_modular_resultant(F, G, p, u, K)
        except HomomorphismFailed:
            continue

        if K.is_one(P):
            r = R
        else:
            r = dmp_apply_pairs(r, R, _collins_crt, (P, p, K), v, K)

        P *= p

    return r

def dmp_zz_collins_resultant(f, g, u, K):
    """
    Collins's modular resultant algorithm in `Z[X]`.

    Examples
    ========

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

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

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

    """

    n = dmp_degree(f, u)
    m = dmp_degree(g, u)

    if n < 0 or m < 0:
        return dmp_zero(u-1)

    A = dmp_max_norm(f, u, K)
    B = dmp_max_norm(g, u, K)

    a = dmp_ground_LC(f, u, K)
    b = dmp_ground_LC(g, u, K)

    v = u - 1

    B = K(2)*K.factorial(n+m)*A**m*B**n
    r, p, P = dmp_zero(v), K.one, K.one

    while P <= B:
        p = K(nextprime(p))

        while not (a % p) or not (b % p):
            p = K(nextprime(p))

        F = dmp_ground_trunc(f, p, u, K)
        G = dmp_ground_trunc(g, p, u, K)

        try:
            R = dmp_zz_modular_resultant(F, G, p, u, K)
        except HomomorphismFailed:
            continue

        if K.is_one(P):
            r = R
        else:
            r = dmp_apply_pairs(r, R, _collins_crt, (P, p, K), v, K)

        P *= p

    return r

def dmp_pdiv(f, g, u, K):
    """
    Polynomial pseudo-division in ``K[X]``.

    Examples
    ========

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

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

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

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    q, r = dmp_zero(u), f

    if df < dg:
        return q, r

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

    while True:
        dr = dmp_degree(r, u)

        if dr < dg:
            break

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

        Q = dmp_mul_term(q, lc_g, 0, u, K)
        q = dmp_add_term(Q, lc_r, j, u, K)

        R = dmp_mul_term(r, lc_g, 0, u, K)
        G = dmp_mul_term(g, lc_r, j, u, K)
        r = dmp_sub(R, G, u, K)

    c = dmp_pow(lc_g, N, u-1, K)

    q = dmp_mul_term(q, c, 0, u, K)
    r = dmp_mul_term(r, c, 0, u, K)

    return q, r

def dmp_pdiv(f, g, u, K):
    """
    Polynomial pseudo-division in ``K[X]``.

    Examples
    ========

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

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

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    q, r, dr = dmp_zero(u), f, df

    if df < dg:
        return q, r

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

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

        Q = dmp_mul_term(q, lc_g, 0, u, K)
        q = dmp_add_term(Q, lc_r, j, u, K)

        R = dmp_mul_term(r, lc_g, 0, u, K)
        G = dmp_mul_term(g, lc_r, j, u, K)
        r = dmp_sub(R, G, u, K)

        _dr, dr = dr, dmp_degree(r, u)

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

    c = dmp_pow(lc_g, N, u - 1, K)

    q = dmp_mul_term(q, c, 0, u, K)
    r = dmp_mul_term(r, c, 0, u, K)

    return q, r

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

    Examples
    ========

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

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

    >>> dmp_ff_div(f, g, 1, QQ)
    ([[1/2], [1/2, -1/2]], [[-1/1, 1/1]])

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    q, r = dmp_zero(u), f

    if df < dg:
        return q, r

    lc_g, v = dmp_LC(g, K), u-1

    while True:
        dr = dmp_degree(r, u)

        if dr < dg:
            break

        lc_r = dmp_LC(r, K)

        c, R = dmp_ff_div(lc_r, lc_g, v, K)

        if not dmp_zero_p(R, v):
            break

        j = dr - dg

        q = dmp_add_term(q, c, j, u, K)
        h = dmp_mul_term(g, c, j, u, K)

        r = dmp_sub(r, h, u, K)

    return q, r

def dmp_rr_div(f, g, u, K):
    """
    Multivariate division with remainder over a ring.

    **Examples**

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

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

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

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    q, r = dmp_zero(u), f

    if df < dg:
        return q, r

    lc_g, v = dmp_LC(g, K), u-1

    while True:
        dr = dmp_degree(r, u)

        if dr < dg:
            break

        lc_r = dmp_LC(r, K)

        c, R = dmp_rr_div(lc_r, lc_g, v, K)

        if not dmp_zero_p(R, v):
            break

        j = dr - dg

        q = dmp_add_term(q, c, j, u, K)
        h = dmp_mul_term(g, c, j, u, K)

        r = dmp_sub(r, h, u, K)

    return q, r

def dmp_rr_div(f, g, u, K):
    """
    Multivariate division with remainder over a ring.

    Examples
    ========

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

    >>> R.dmp_rr_div(x**2 + x*y, 2*x + 2)
    (0, x**2 + x*y)

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    q, r = dmp_zero(u), f

    if df < dg:
        return q, r

    lc_g, v = dmp_LC(g, K), u - 1

    while True:
        dr = dmp_degree(r, u)

        if dr < dg:
            break

        lc_r = dmp_LC(r, K)

        c, R = dmp_rr_div(lc_r, lc_g, v, K)

        if not dmp_zero_p(R, v):
            break

        j = dr - dg

        q = dmp_add_term(q, c, j, u, K)
        h = dmp_mul_term(g, c, j, u, K)

        r = dmp_sub(r, h, u, K)

    return q, r

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

    Examples
    ========

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

    >>> R.dmp_ff_div(x**2 + x*y, 2*x + 2)
    (1/2*x + 1/2*y - 1/2, -y + 1)

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    q, r, dr = dmp_zero(u), f, df

    if df < dg:
        return q, r

    lc_g, v = dmp_LC(g, K), u - 1

    while True:
        lc_r = dmp_LC(r, K)
        c, R = dmp_ff_div(lc_r, lc_g, v, K)

        if not dmp_zero_p(R, v):
            break

        j = dr - dg

        q = dmp_add_term(q, c, j, u, K)
        h = dmp_mul_term(g, c, j, u, K)
        r = dmp_sub(r, h, u, K)

        _dr, dr = dr, dmp_degree(r, u)

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

    return q, r

def dmp_prem(f, g, u, K):
    """
    Polynomial pseudo-remainder in ``K[X]``.

    **Examples**

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

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

    >>> dmp_prem(f, g, 1, ZZ)
    [[-4, 4]]

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    r = f

    if df < dg:
        return r

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

    while True:
        dr = dmp_degree(r, u)

        if dr < dg:
            break

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

        R = dmp_mul_term(r, lc_g, 0, u, K)
        G = dmp_mul_term(g, lc_r, j, u, K)
        r = dmp_sub(R, G, u, K)

    c = dmp_pow(lc_g, N, u-1, K)

    return dmp_mul_term(r, c, 0, u, K)

def dmp_discriminant(f, u, K):
    """
    Computes discriminant of a polynomial in `K[X]`.

    Examples
    ========

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

    >>> R.dmp_discriminant(x**2*y + x*z + t)
    -4*y*t + z**2

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

    d, v = dmp_degree(f, u), u - 1

    if d <= 0:
        return dmp_zero(v)
    else:
        s = (-1)**((d*(d - 1)) // 2)
        c = dmp_LC(f, K)

        r = dmp_resultant(f, dmp_diff(f, 1, u, K), u, K)
        c = dmp_mul_ground(c, K(s), v, K)

        return dmp_quo(r, c, v, K)

def dmp_discriminant(f, u, K):
    """
    Computes discriminant of a polynomial in ``K[X]``.

    **Examples**

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

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

    >>> dmp_discriminant(f, 3, ZZ)
    [[[-4, 0]], [[1], [], []]]

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

    d, v = dmp_degree(f, u), u-1

    if d <= 0:
        return dmp_zero(v)
    else:
        s = (-1)**((d*(d-1)) // 2)
        c = dmp_LC(f, K)

        r = dmp_resultant(f, dmp_diff(f, 1, u, K), u, K)
        c = dmp_mul_ground(c, K(s), v, K)

        return dmp_exquo(r, c, v, K)

def dmp_sqr(f, u, K):
    """Square dense polynomials in `K[X]`. """
    if not u:
        return dup_sqr(f, K)

    df = dmp_degree(f, u)

    if df < 0:
        return f

    h, v = [], u-1

    for i in xrange(0, 2*df+1):
        c = dmp_zero(v)

        jmin = max(0, i-df)
        jmax = min(i, df)

        n = jmax - jmin + 1

        jmax = jmin + n // 2 - 1

        for j in xrange(jmin, jmax+1):
            c = dmp_add(c, dmp_mul(f[j], f[i-j], v, K), v, K)

        c = dmp_mul_ground(c, 2, v, K)

        if n & 1:
            elem = dmp_sqr(f[jmax+1], v, K)
            c = dmp_add(c, elem, v, K)

        h.append(c)

    return h

def dmp_prem(f, g, u, K):
    """
    Polynomial pseudo-remainder in ``K[X]``.

    Examples
    ========

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

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

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

    df = dmp_degree(f, u)
    dg = dmp_degree(g, u)

    if dg < 0:
        raise ZeroDivisionError("polynomial division")

    r = f

    if df < dg:
        return r

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

    while True:
        dr = dmp_degree(r, u)

        if dr < dg:
            break

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

        R = dmp_mul_term(r, lc_g, 0, u, K)
        G = dmp_mul_term(g, lc_r, j, u, K)
        r = dmp_sub(R, G, u, K)

    c = dmp_pow(lc_g, N, u - 1, K)

    return dmp_mul_term(r, c, 0, u, K)

def dmp_mul(f, g, u, K):
    """
    Multiply dense polynomials in ``K[X]``.

    Examples
    ========

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

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

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

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

    if f == g:
        return dmp_sqr(f, u, K)

    df = dmp_degree(f, u)

    if df < 0:
        return f

    dg = dmp_degree(g, u)

    if dg < 0:
        return g

    h, v = [], u-1

    for i in xrange(0, df+dg+1):
        coeff = dmp_zero(v)

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

        h.append(coeff)

    return dmp_strip(h, u)