def dmp_deflate(f, u, K):
    """Map `x_i**m_i` to `y_i` in a polynomial in `K[X]`. """
    if dmp_zero_p(f, u):
        return (1,)*(u+1), f

    F = dmp_to_dict(f, u)
    B = [0]*(u+1)

    for M in F.iterkeys():
        for i, m in enumerate(M):
            B[i] = igcd(B[i], m)

    for i, b in enumerate(B):
        if not b:
            B[i] = 1

    B = tuple(B)

    if all([ b == 1 for b in B ]):
        return B, f

    H = {}

    for A, coeff in F.iteritems():
        N = [ a // b for a, b in zip(A, B) ]
        H[tuple(N)] = coeff

    return B, dmp_from_dict(H, u, K)

    def __new__(cls, function, *symbols, **assumptions):

        # Any embedded piecewise functions need to be brought out to the
        # top level so that integration can go into piecewise mode at the
        # earliest possible moment.
        function = piecewise_fold(sympify(function))

        if function is S.NaN:
            return S.NaN

        if symbols:
            limits, sign = _process_limits(*symbols)
        else:
            # no symbols provided -- let's compute full anti-derivative
            limits, sign = [Tuple(s) for s in function.free_symbols], 1

            if len(limits) != 1:
                raise ValueError("specify integration variables to integrate %s" % function)

        while isinstance(function, Integral):
            # denest the integrand
            limits = list(function.limits) + limits
            function = function.function

        obj = Expr.__new__(cls, **assumptions)
        arglist = [sign*function]
        arglist.extend(limits)
        obj._args = tuple(arglist)
        obj.is_commutative = all(s.is_commutative for s in obj.free_symbols)

        return obj

def roots_binomial(f):
    """Returns a list of roots of a binomial polynomial."""
    n = f.degree()

    a, b = f.nth(n), f.nth(0)
    alpha = (-cancel(b/a))**Rational(1, n)

    if alpha.is_number:
        alpha = alpha.expand(complex=True)

    roots, I = [], S.ImaginaryUnit

    for k in xrange(n):
        zeta = exp(2*k*S.Pi*I/n).expand(complex=True)
        roots.append((alpha*zeta).expand(power_base=False))

    if all([ r.is_number for r in roots ]):
        reals, complexes = [], []

        for root in roots:
            if root.is_real:
                reals.append(root)
            else:
                complexes.append(root)

        roots = sorted(reals) + sorted(complexes, key=lambda r: (re(r), -im(r)))

    return roots

def _separatevars_dict(expr, *symbols):
    if symbols:
        assert all((t.is_Atom for t in symbols)), "symbols must be Atoms."
    ret = dict(((i,sympify(1)) for i in symbols))
    ret['coeff'] = sympify(1)
    if expr.is_Mul:
        for i in expr.args:
            expsym = i.atoms(Symbol)
            if len(set(symbols).intersection(expsym)) > 1:
                return None
            if len(set(symbols).intersection(expsym)) == 0:
                # There are no symbols, so it is part of the coefficient
                ret['coeff'] *= i
            else:
                ret[expsym.pop()] *= i
    else:
        expsym = expr.atoms(Symbol)
        if len(set(symbols).intersection(expsym)) > 1:
            return None
        if len(set(symbols).intersection(expsym)) == 0:
            # There are no symbols, so it is part of the coefficient
            ret['coeff'] *= expr
        else:
            ret[expsym.pop()] *= expr

    return ret

def test_roots_quartic():
    assert roots_quartic(Poly(x**4, x)) == [0, 0, 0, 0]
    assert roots_quartic(Poly(x**4 + x**3, x)) in [
        [-1,0,0,0],
        [0,-1,0,0],
        [0,0,-1,0],
        [0,0,0,-1]
    ]
    assert roots_quartic(Poly(x**4 - x**3, x)) in [
        [1,0,0,0],
        [0,1,0,0],
        [0,0,1,0],
        [0,0,0,1]
    ]

    lhs = roots_quartic(Poly(x**4 + x, x))
    rhs = [S.Half + I*sqrt(3)/2, S.Half - I*sqrt(3)/2, S.Zero, -S.One]

    assert sorted(lhs, key=hash) == sorted(rhs, key=hash)

    # test of all branches of roots quartic
    for i, (a, b, c, d) in enumerate([(1, 2, 3, 0),
                                      (3, -7, -9, 9),
                                      (1, 2, 3, 4),
                                      (1, 2, 3, 4),
                                      (-7, -3, 3, -6),
                                      (-3, 5, -6, -4)]):
        if i == 2:
            c = -a*(a**2/S(8) - b/S(2))
        elif i == 3:
            d = a*(a*(3*a**2/S(256) - b/S(16)) + c/S(4))
        eq = x**4 + a*x**3 + b*x**2 + c*x + d
        ans = roots_quartic(Poly(eq, x))
        assert all([eq.subs(x, ai).n(chop=True) == 0 for ai in ans])

def monomial_div(A, B):
    """
    Division of tuples representing monomials.

    Lets divide `x**3*y**4*z` by `x*y**2`::

        >>> from sympy.polys.monomialtools import monomial_div

        >>> monomial_div((3, 4, 1), (1, 2, 0))
        (2, 2, 1)

    which gives `x**2*y**2*z`. However::

        >>> monomial_div((3, 4, 1), (1, 2, 2)) is None
        True

    `x*y**2*z**2` does not divide `x**3*y**4*z`.

    """
    C = [ a - b for a, b in zip(A, B) ]

    if all([ c >= 0 for c in C ]):
        return tuple(C)
    else:
        return None

def dmp_terms_gcd(f, u, K):
    """
    Remove GCD of terms from ``f`` in ``K[X]``.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densebasic import dmp_terms_gcd

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

    >>> dmp_terms_gcd(f, 1, ZZ)
    ((2, 1), [[1], [1, 0]])

    """
    if dmp_ground_TC(f, u, K) or dmp_zero_p(f, u):
        return (0,)*(u+1), f

    F = dmp_to_dict(f, u)
    G = monomial_min(*F.keys())

    if all([ g == 0 for g in G ]):
        return G, f

    f = {}

    for monom, coeff in F.iteritems():
        f[monomial_div(monom, G)] = coeff

    return G, dmp_from_dict(f, u, K)

def dmp_ext_factor(f, u, K):
    """Factor multivariate polynomials over algebraic number fields. """
    if not u:
        return dup_ext_factor(f, K)

    lc = dmp_ground_LC(f, u, K)
    f = dmp_ground_monic(f, u, K)

    if all([ d <= 0 for d in dmp_degree_list(f, u) ]):
        return lc, []

    f, F = dmp_sqf_part(f, u, K), f
    s, g, r = dmp_sqf_norm(f, u, K)

    factors = dmp_factor_list_include(r, u, K.dom)

    if len(factors) == 1:
        coeff, factors = lc, [f]
    else:
        H = dmp_raise([K.one, s*K.unit], u, 0, K)

        for i, (factor, _) in enumerate(factors):
            h = dmp_convert(factor, u, K.dom, K)
            h, _, g = dmp_inner_gcd(h, g, u, K)
            h = dmp_compose(h, H, u, K)
            factors[i] = h

    return lc, dmp_trial_division(F, factors, u, K)

def test_simplify():
    x, y, z, k, n, m, w, f, s, A = symbols('x,y,z,k,n,m,w,f,s,A')

    assert all(simplify(tmp) == tmp for tmp in [I, E, oo, x, -x, -oo, -E, -I])

    e = 1/x + 1/y
    assert e != (x+y)/(x*y)
    assert simplify(e) == (x+y)/(x*y)

    e = A**2*s**4/(4*pi*k*m**3)
    assert simplify(e) == e

    e = (4+4*x-2*(2+2*x))/(2+2*x)
    assert simplify(e) == 0

    e = (-4*x*y**2-2*y**3-2*x**2*y)/(x+y)**2
    assert simplify(e) == -2*y

    e = -x-y-(x+y)**(-1)*y**2+(x+y)**(-1)*x**2
    assert simplify(e) == -2*y

    e = (x+x*y)/x
    assert simplify(e) == 1 + y

    e = (f(x)+y*f(x))/f(x)
    assert simplify(e) == 1 + y

    e = (2 * (1/n - cos(n * pi)/n))/pi
    assert simplify(e) == 2*((1 - 1*cos(pi*n))/(pi*n))

    e = integrate(1/(x**3+1), x).diff(x)
    assert simplify(e) == 1/(x**3+1)

    e = integrate(x/(x**2+3*x+1), x).diff(x)
    assert simplify(e) == x/(x**2+3*x+1)

    A = Matrix([[2*k-m*w**2, -k], [-k, k-m*w**2]]).inv()

    assert simplify((A*Matrix([0,f]))[1]) == \
        (f*(2*k - m*w**2))/(k**2 - 3*k*m*w**2 + m**2*w**4)

    a, b, c, d, e, f, g, h, i = symbols('a,b,c,d,e,f,g,h,i')

    f_1 = x*a + y*b + z*c - 1
    f_2 = x*d + y*e + z*f - 1
    f_3 = x*g + y*h + z*i - 1

    solutions = solve([f_1, f_2, f_3], x, y, z, simplified=False)

    assert simplify(solutions[y]) == \
        (a*i+c*d+f*g-a*f-c*g-d*i)/(a*e*i+b*f*g+c*d*h-a*f*h-b*d*i-c*e*g)

    f = -x + y/(z + t) + z*x/(z + t) + z*a/(z + t) + t*x/(z + t)

    assert simplify(f) == (y + a*z)/(z + t)

    A, B = symbols('A,B', commutative=False)

    assert simplify(A*B - B*A) == A*B - B*A

def denester(nested):
    """
    Denests a list of expressions that contain nested square roots.
    This method should not be called directly - use 'denest' instead.
    This algorithm is based on <http://www.almaden.ibm.com/cs/people/fagin/symb85.pdf>.

    It is assumed that all of the elements of 'nested' share the same
    bottom-level radicand. (This is stated in the paper, on page 177, in
    the paragraph immediately preceding the algorithm.)

    When evaluating all of the arguments in parallel, the bottom-level
    radicand only needs to be denested once. This means that calling
    denester with x arguments results in a recursive invocation with x+1
    arguments; hence denester has polynomial complexity.

    However, if the arguments were evaluated separately, each call would
    result in two recursive invocations, and the algorithm would have
    exponential complexity.

    This is discussed in the paper in the middle paragraph of page 179.
    """
    if all((n ** 2).is_Number for n in nested):  # If none of the arguments are nested
        for f in subsets(len(nested)):  # Test subset 'f' of nested
            p = prod(nested[i] ** 2 for i in range(len(f)) if f[i]).expand()
            if 1 in f and f.count(1) > 1 and f[-1]:
                p = -p
            if sqrt(p).is_Number:
                return sqrt(p), f  # If we got a perfect square, return its square root.
        return nested[-1], [0] * len(nested)  # Otherwise, return the radicand from the previous invocation.
    else:
        a, b, r, R = Wild("a"), Wild("b"), Wild("r"), None
        values = [expr.match(sqrt(a + b * sqrt(r))) for expr in nested]
        for v in values:
            if r in v:  # Since if b=0, r is not defined
                if R is not None:
                    assert R == v[r]  # All the 'r's should be the same.
                else:
                    R = v[r]
        d, f = denester([sqrt((v[a] ** 2).expand() - (R * v[b] ** 2).expand()) for v in values] + [sqrt(R)])
        if not any([f[i] for i in range(len(nested))]):  # If f[i]=0 for all i < len(nested)
            v = values[-1]
            return sqrt(v[a] + v[b] * d), f
        else:
            v = prod(nested[i] ** 2 for i in range(len(nested)) if f[i]).expand().match(a + b * sqrt(r))
            if 1 in f and f.index(1) < len(nested) - 1 and f[len(nested) - 1]:
                v[a] = -1 * v[a]
                v[b] = -1 * v[b]
            if not f[len(nested)]:  # Solution denests with square roots
                return (
                    (
                        sqrt((v[a] + d).expand() / 2) + sign(v[b]) * sqrt((v[b] ** 2 * R / (2 * (v[a] + d))).expand())
                    ).expand(),
                    f,
                )
            else:  # Solution requires a fourth root
                FR, s = (R.expand() ** Rational(1, 4)), sqrt((v[b] * R).expand() + d)
                return (s / (sqrt(2) * FR) + v[a] * FR / (sqrt(2) * s)).expand(), f

def dmp_inflate(f, M, u, K):
    """Map `y_i` to `x_i**k_i` in a polynomial in `K[X]`. """
    if not u:
        return dup_inflate(f, M[0], K)

    if all([ m == 1 for m in M ]):
        return f
    else:
        return _rec_inflate(f, M, u, 0, K)

    def __new__(cls, function, *symbols, **assumptions):
        # Any embedded piecewise functions need to be brought out to the
        # top level so that integration can go into piecewise mode at the
        # earliest possible moment.
        function = piecewise_fold(sympify(function))

        if function is S.NaN:
            return S.NaN

        symbols = list(symbols)
        if not symbols:
            # no symbols provided -- let's compute full anti-derivative
            symbols = sorted(function.free_symbols, Basic.compare)
            if not symbols:
                raise ValueError('An integration variable is required.')

        while isinstance(function, Integral):
            # denest the integrand
            symbols = list(function.limits) + symbols
            function = function.function

        limits = []
        for V in symbols:
            if isinstance(V, Symbol):
                limits.append(Tuple(V))
                continue
            elif isinstance(V, (tuple, list, Tuple)):
                V = sympify(flatten(V))
                if V[0].is_Symbol:
                    newsymbol = V[0]
                    if len(V) == 3:
                        if V[1] is None and V[2] is not None:
                            nlim = [V[2]]
                        elif V[1] is not None and V[2] is None:
                            function = -function
                            nlim = [V[1]]
                        elif V[1] is None and V[2] is None:
                            nlim = []
                        else:
                            nlim = V[1:]
                        limits.append(Tuple(newsymbol, *nlim ))
                        continue
                    elif len(V) == 1 or (len(V) == 2 and V[1] is None):
                        limits.append(Tuple(newsymbol))
                        continue
                    elif len(V) == 2:
                        limits.append(Tuple(newsymbol, V[1]))
                        continue
            raise ValueError("Invalid integration variable or limits: %s" % str(symbols))

        obj = Expr.__new__(cls, **assumptions)
        obj._args = tuple([function] + limits)
        obj.is_commutative = all(s.is_commutative for s in obj.free_symbols)

        return obj

def test_dup_random():
    f = dup_random(0, -10, 10, ZZ)

    assert dup_degree(f) == 0
    assert all([ -10 <= c <= 10 for c in f ])

    f = dup_random(1, -20, 20, ZZ)

    assert dup_degree(f) == 1
    assert all([ -20 <= c <= 20 for c in f ])

    f = dup_random(2, -30, 30, ZZ)

    assert dup_degree(f) == 2
    assert all([ -30 <= c <= 30 for c in f ])

    f = dup_random(3, -40, 40, ZZ)

    assert dup_degree(f) == 3
    assert all([ -40 <= c <= 40 for c in f ])

def dmp_multi_deflate(polys, u, K):
    """
    Map ``x_i**m_i`` to ``y_i`` in a set of polynomials in ``K[X]``.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densebasic import dmp_multi_deflate

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

    >>> dmp_multi_deflate((f, g), 1, ZZ)
    ((2, 1), ([[1, 0, 0, 2], [3, 0, 0, 4]], [[1, 0, 2], [3, 0, 4]]))

    """
    if not u:
        M, H = dup_multi_deflate(polys, K)
        return (M,), H

    F, B = [], [0]*(u+1)

    for p in polys:
        f = dmp_to_dict(p, u)

        if not dmp_zero_p(p, u):
            for M in f.iterkeys():
                for i, m in enumerate(M):
                    B[i] = igcd(B[i], m)

        F.append(f)

    for i, b in enumerate(B):
        if not b:
            B[i] = 1

    B = tuple(B)

    if all([ b == 1 for b in B ]):
        return B, polys

    H = []

    for f in F:
        h = {}

        for A, coeff in f.iteritems():
            N = [ a // b for a, b in zip(A, B) ]
            h[tuple(N)] = coeff

        H.append(dmp_from_dict(h, u, K))

    return B, tuple(H)

def dmp_terms_gcd(f, u, K):
    """Remove GCD of terms from `f` in `K[X]`. """
    if dmp_ground_TC(f, u, K) or dmp_zero_p(f, u):
        return (0,)*(u+1), f

    F = dmp_to_dict(f, u)
    G = monomial_min(*F.keys())

    if all([ g == 0 for g in G ]):
        return G, f

    f = {}

    for monom, coeff in F.iteritems():
        f[monomial_div(monom, G)] = coeff

    return G, dmp_from_dict(f, u, K)

def dmp_deflate(f, u, K):
    """
    Map ``x_i**m_i`` to ``y_i`` in a polynomial in ``K[X]``.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densebasic import dmp_deflate

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

    >>> dmp_deflate(f, 1, ZZ)
    ((2, 3), [[1, 2], [3, 4]])

    """
    if dmp_zero_p(f, u):
        return (1,)*(u+1), f

    F = dmp_to_dict(f, u)
    B = [0]*(u+1)

    for M in F.iterkeys():
        for i, m in enumerate(M):
            B[i] = igcd(B[i], m)

    for i, b in enumerate(B):
        if not b:
            B[i] = 1

    B = tuple(B)

    if all([ b == 1 for b in B ]):
        return B, f

    H = {}

    for A, coeff in F.iteritems():
        N = [ a // b for a, b in zip(A, B) ]
        H[tuple(N)] = coeff

    return B, dmp_from_dict(H, u, K)

def _get_indices_Add(expr):
    """Determine outer indices of an Add object.

    In a sum, each term must have the same set of outer indices.  A valid
    expression could be

        x(i)*y(j) - x(j)*y(i)

    But we do not allow expressions like:

        x(i)*y(j) - z(j)*z(j)

    FIXME: Add support for Numpy broadcasting

    >>> from sympy.tensor.index_methods import _get_indices_Add
    >>> from sympy.tensor.indexed import IndexedBase, Idx
    >>> i, j, k = map(Idx, ['i', 'j', 'k'])
    >>> x = IndexedBase('x')
    >>> y = IndexedBase('y')
    >>> _get_indices_Add(x[i] + x[k]*y[i, k])
    (set([i]), {})

    """

    inds = map(get_indices, expr.args)
    inds, syms = zip(*inds)

    # allow broadcast of scalars
    non_scalars = filter(lambda x: x != set(), inds)
    if not non_scalars:
        return set(), {}

    if not all(map(lambda x: x == non_scalars[0], non_scalars[1:])):
        raise IndexConformanceException("Indices are not consistent: %s"%expr)
    if not reduce(lambda x, y: x!=y or y, syms):
        symmetries = syms[0]
    else:
        # FIXME: search for symmetries
        symmetries = {}

    return non_scalars[0], symmetries

def monomial_div(a, b):
    """Division of tuples representing monomials.

       Lets divide x**3*y**4*z by x*y**2:

       >>> monomial_div((3, 4, 1), (1, 2, 0))
       (2, 2, 1)

       which gives x**2*y**2*z. However

       >>> monomial_div((3, 4, 1), (1, 2, 2))
       None

       x*y**2*z**2 does not divide x**3*y**4*z.

    """
    result = [x - y for x, y in zip(a, b)]

    if all(e >= 0 for e in result):
        return tuple(result)
    else:
        return None

def dmp_multi_deflate(polys, u, K):
    """Map `x_i**m_i` to `y_i` in a set of polynomials in `K[X]`. """
    if not u:
        M, H = dup_multi_deflate(polys, K)
        return (M,), H

    F, B = [], [0]*(u+1)

    for p in polys:
        f = dmp_to_dict(p, u)

        if not dmp_zero_p(p, u):
            for M in f.iterkeys():
                for i, m in enumerate(M):
                    B[i] = igcd(B[i], m)

        F.append(f)

    for i, b in enumerate(B):
        if not b:
            B[i] = 1

    B = tuple(B)

    if all([ b == 1 for b in B ]):
        return B, polys

    H = []

    for f in F:
        h = {}

        for A, coeff in f.iteritems():
            N = [ a // b for a, b in zip(A, B) ]
            h[tuple(N)] = coeff

        H.append(dmp_from_dict(h, u, K))

    return B, tuple(H)

def dmp_inflate(f, M, u, K):
    """
    Map ``y_i`` to ``x_i**k_i`` in a polynomial in ``K[X]``.

    **Examples**

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densebasic import dmp_inflate

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

    >>> dmp_inflate(f, (2, 3), 1, ZZ)
    [[1, 0, 0, 2], [], [3, 0, 0, 4]]

    """
    if not u:
        return dup_inflate(f, M[0], K)

    if all([ m == 1 for m in M ]):
        return f
    else:
        return _rec_inflate(f, M, u, 0, K)