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), (6, -5, -10, -3)]
    ):
        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)

    # not all symbolic quartics are unresolvable
    eq = Poly(q * x + q / 4 + x ** 4 + x ** 3 + 2 * x ** 2 - Rational(1, 3), x)
    sol = roots_quartic(eq)
    assert all(test_numerically(eq.subs(x, i), 0) for i in sol)
    # but some are (see also iss 1890)
    raises(PolynomialError, lambda: roots_quartic(Poly(y * x ** 4 + x + z, x)))

def _sqrt_symbolic_denest(a, b, r):
    """Given an expression, sqrt(a + b*sqrt(b)), return the denested
    expression or None.

    Algorithm:
    If r = ra + rb*sqrt(rr), try replacing sqrt(rr) in ``a`` with
    (y**2 - ra)/rb, and if the result is a quadratic, ca*y**2 + cb*y + cc, and
    (cb + b)**2 - 4*ca*cc is 0, then sqrt(a + b*sqrt(r)) can be rewritten as
    sqrt(ca*(sqrt(r) + (cb + b)/(2*ca))**2).

    Examples
    ========

    >>> from sympy.simplify.sqrtdenest import _sqrt_symbolic_denest, sqrtdenest
    >>> from sympy import sqrt, Symbol
    >>> from sympy.abc import x

    >>> a, b, r = 16 - 2*sqrt(29), 2, -10*sqrt(29) + 55
    >>> _sqrt_symbolic_denest(a, b, r)
    sqrt(-2*sqrt(29) + 11) + sqrt(5)

    If the expression is numeric, it will be simplified:

    >>> w = sqrt(sqrt(sqrt(3) + 1) + 1) + 1 + sqrt(2)
    >>> sqrtdenest(sqrt((w**2).expand()))
    1 + sqrt(2) + sqrt(1 + sqrt(1 + sqrt(3)))

    Otherwise, it will only be simplified if assumptions allow:

    >>> w = w.subs(sqrt(3), sqrt(x + 3))
    >>> sqrtdenest(sqrt((w**2).expand()))
    sqrt((sqrt(sqrt(sqrt(x + 3) + 1) + 1) + 1 + sqrt(2))**2)

    Notice that the argument of the sqrt is a square. If x is made positive
    then the sqrt of the square is resolved:

    >>> _.subs(x, Symbol('x', positive=True))
    sqrt(sqrt(sqrt(x + 3) + 1) + 1) + 1 + sqrt(2)
    """

    a, b, r = map(sympify, (a, b, r))
    rval = _sqrt_match(r)
    if not rval:
        return None
    ra, rb, rr = rval
    if rb:
        y = Dummy('y', positive=True)
        try:
            newa = Poly(a.subs(sqrt(rr), (y**2 - ra)/rb), y)
        except PolynomialError:
            return None
        if newa.degree() == 2:
            ca, cb, cc = newa.all_coeffs()
            cb += b
            if _mexpand(cb**2 - 4*ca*cc).equals(0):
                z = sqrt(ca*(sqrt(r) + cb/(2*ca))**2)
                if z.is_number:
                    z = _mexpand(Mul._from_args(z.as_content_primitive()))
                return z

def test_issue_8438():
    p = Poly([1, y, -2, -3], x).as_expr()
    roots = roots_cubic(Poly(p, x), x)
    z = -S(3)/2 - 7*I/2  # this will fail in code given in commit msg
    post = [r.subs(y, z) for r in roots]
    assert set(post) == \
    set(roots_cubic(Poly(p.subs(y, z), x)))
    # /!\ if p is not made an expression, this is *very* slow
    assert all(p.subs({y: z, x: i}).n(2, chop=True) == 0 for i in post)

def reduce_poly_inequalities(exprs, gen, assume=True, relational=True):
    """Reduce a system of polynomial inequalities with rational coefficients. """
    exact = True
    polys = []

    for _exprs in exprs:
        _polys = []

        for expr in _exprs:
            if isinstance(expr, tuple):
                expr, rel = expr
            else:
                if expr.is_Relational:
                    expr, rel = expr.lhs - expr.rhs, expr.rel_op
                else:
                    expr, rel = expr, '=='

            poly = Poly(expr, gen)

            if not poly.get_domain().is_Exact:
                poly, exact = poly.to_exact(), False

            domain = poly.get_domain()

            if not (domain.is_ZZ or domain.is_QQ):
                raise NotImplementedError("inequality solving is not supported over %s" % domain)

            _polys.append((poly, rel))

        polys.append(_polys)

    solution = solve_poly_inequalities(polys)

    if isinstance(solution, Union):
        intervals = list(solution.args)
    elif isinstance(solution, Interval):
        intervals = [solution]
    else:
        intervals = []

    if not exact:
        intervals = map(interval_evalf, intervals)

    if not relational:
        return intervals

    real = ask(gen, 'real', assume)

    def relationalize(gen):
        return Or(*[ i.as_relational(gen) for i in intervals ])

    if not real:
        result = And(relationalize(re(gen)), Eq(im(gen), 0))
    else:
        result = relationalize(gen)

    return result

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),
                                      (6, -5, -10, -3)]):
        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)

    # not all symbolic quartics are unresolvable
    eq = Poly(q*x + q/4 + x**4 + x**3 + 2*x**2 - Rational(1, 3), x)
    sol = roots_quartic(eq)
    assert all(verify_numerically(eq.subs(x, i), 0) for i in sol)
    z = symbols('z', negative=True)
    eq = x**4 + 2*x**3 + 3*x**2 + x*(z + 11) + 5
    zans = roots_quartic(Poly(eq, x))
    assert all([verify_numerically(eq.subs(((x, i), (z, -1))), 0) for i in zans])
    # but some are (see also issue 4989)
    # it's ok if the solution is not Piecewise, but the tests below should pass
    eq = Poly(y*x**4 + x**3 - x + z, x)
    ans = roots_quartic(eq)
    assert all(type(i) == Piecewise for i in ans)
    reps = (
        dict(y=-Rational(1, 3), z=-Rational(1, 4)),  # 4 real
        dict(y=-Rational(1, 3), z=-Rational(1, 2)),  # 2 real
        dict(y=-Rational(1, 3), z=-2))  # 0 real
    for rep in reps:
        sol = roots_quartic(Poly(eq.subs(rep), x))
        assert all([verify_numerically(w.subs(rep) - s, 0) for w, s in zip(ans, sol)])

def _is_negative_or_zero(term):
    if getattr(term, 'is_number', False):
        return term <= 0
    elif isinstance(term, Pow):
        if term.args[1]%2==1:
            return _is_negative_or_zero(term.args[0])

    t = Poly(term).as_dict()
    if (all(c < 0 for c in t.values()) and
        all(i % 2 == 0 for d in t.keys() for i in d)):
        return True
    return ask_is_negative(term)

def _solve_as_poly(f, symbol, solveset_solver, invert_func):
    """
    Solve the equation using polynomial techniques if it already is a
    polynomial equation or, with a change of variables, can be made so.
    """
    result = None
    if f.is_polynomial(symbol):

        solns = roots(f, symbol, cubics=True, quartics=True,
                      quintics=True, domain='EX')
        num_roots = sum(solns.values())
        if degree(f, symbol) <= num_roots:
            result = FiniteSet(*solns.keys())
        else:
            poly = Poly(f, symbol)
            solns = poly.all_roots()
            if poly.degree() <= len(solns):
                result = FiniteSet(*solns)
            else:
                result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
    else:
        poly = Poly(f)
        if poly is None:
            result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
        gens = [g for g in poly.gens if g.has(symbol)]

        if len(gens) == 1:
            poly = Poly(poly, gens[0])
            gen = poly.gen
            deg = poly.degree()
            poly = Poly(poly.as_expr(), poly.gen, composite=True)
            poly_solns = FiniteSet(*roots(poly, cubics=True, quartics=True,
                                          quintics=True).keys())

            if len(poly_solns) < deg:
                result = ConditionSet(symbol, Eq(f, 0), S.Complexes)

            if gen != symbol:
                y = Dummy('y')
                lhs, rhs_s = invert_func(gen, y, symbol)
                if lhs is symbol:
                    result = Union(*[rhs_s.subs(y, s) for s in poly_solns])
                else:
                    result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
        else:
            result = ConditionSet(symbol, Eq(f, 0), S.Complexes)

    if result is not None:
        if isinstance(result, FiniteSet):
            # this is to simplify solutions like -sqrt(-I) to sqrt(2)/2
            # - sqrt(2)*I/2. We are not expanding for solution with free
            # variables because that makes the solution more complicated. For
            # example expand_complex(a) returns re(a) + I*im(a)
            if all([s.free_symbols == set() and not isinstance(s, RootOf)
                    for s in result]):
                s = Dummy('s')
                result = imageset(Lambda(s, expand_complex(s)), result)
        return result
    else:
        return ConditionSet(symbol, Eq(f, 0), S.Complexes)

def ratint_ratpart(f, g, x):
    """Horowitz-Ostrogradsky algorithm.

       Given a field K and polynomials f and g in K[x], such that f and g
       are coprime and deg(f) < deg(g), returns fractions A and B in K(x),
       such that f/g = A' + B and B has square-free denominator.

    """
    f = Poly(f, x)
    g = Poly(g, x)

    u, v, _ = g.cofactors(g.diff())

    n = u.degree()
    m = v.degree()

    A_coeffs = [ Dummy('a' + str(n-i)) for i in xrange(0, n) ]
    B_coeffs = [ Dummy('b' + str(m-i)) for i in xrange(0, m) ]

    C_coeffs = A_coeffs + B_coeffs

    A = Poly(A_coeffs, x, domain=ZZ[C_coeffs])
    B = Poly(B_coeffs, x, domain=ZZ[C_coeffs])

    H = f - A.diff()*v + A*(u.diff()*v).quo(u) - B*u

    result = solve(H.coeffs(), C_coeffs)

    A = A.as_expr().subs(result)
    B = B.as_expr().subs(result)

    rat_part = cancel(A/u.as_expr(), x)
    log_part = cancel(B/v.as_expr(), x)

    return rat_part, log_part

def test_nroots2():
    p = Poly(x**5+3*x+1, x)

    roots = p.nroots(n=3)
    # The order of roots matters. The roots are ordered by their real
    # components (if they agree, then by their imaginary components).
    assert [str(r) for r in roots] == \
            ['-0.839 - 0.944*I', '-0.839 + 0.944*I', '-0.332',
                '1.01 - 0.937*I', '1.01 + 0.937*I']

    roots = p.nroots(n=5)
    assert [str(r) for r in roots] == \
            ['-0.83907 - 0.94385*I', '-0.83907 + 0.94385*I',
                '-0.33199', '1.0051 - 0.93726*I', '1.0051 + 0.93726*I']

def max_onepiece(x, f: Poly, g: Poly, l, u):
    roots = sorted(set((f - g).real_roots()))
    new_polynomial_pieces = []
    new_bounds = [l]
    for r in roots:
        if l < r < u:
            m = (r + new_bounds[-1]) / 2
            if f.subs(x, m) >= g.subs(x, m):
                new_polynomial_pieces.append(f)
            else:
                new_polynomial_pieces.append(g)
            new_bounds.append(r)
    new_bounds.append(u)
    return PiecewisePolynomial(new_polynomial_pieces, new_bounds)

def ratint_ratpart(f, g, x):
    """
    Horowitz-Ostrogradsky algorithm.

    Given a field K and polynomials f and g in K[x], such that f and g
    are coprime and deg(f) < deg(g), returns fractions A and B in K(x),
    such that f/g = A' + B and B has square-free denominator.

    Examples
    ========

        >>> from sympy.integrals.rationaltools import ratint_ratpart
        >>> from sympy.abc import x, y
        >>> from sympy import Poly
        >>> ratint_ratpart(Poly(1, x, domain='ZZ'),
        ... Poly(x + 1, x, domain='ZZ'), x)
        (0, 1/(x + 1))
        >>> ratint_ratpart(Poly(1, x, domain='EX'),
        ... Poly(x**2 + y**2, x, domain='EX'), x)
        (0, 1/(x**2 + y**2))
        >>> ratint_ratpart(Poly(36, x, domain='ZZ'),
        ... Poly(x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2, x, domain='ZZ'), x)
        ((12*x + 6)/(x**2 - 1), 12/(x**2 - x - 2))

    See Also
    ========

    ratint, ratint_logpart
    """
    from sympy import solve

    f = Poly(f, x)
    g = Poly(g, x)

    u, v, _ = g.cofactors(g.diff())

    n = u.degree()
    m = v.degree()

    A_coeffs = [ Dummy('a' + str(n - i)) for i in range(0, n) ]
    B_coeffs = [ Dummy('b' + str(m - i)) for i in range(0, m) ]

    C_coeffs = A_coeffs + B_coeffs

    A = Poly(A_coeffs, x, domain=ZZ[C_coeffs])
    B = Poly(B_coeffs, x, domain=ZZ[C_coeffs])

    H = f - A.diff()*v + A*(u.diff()*v).quo(u) - B*u

    result = solve(H.coeffs(), C_coeffs)

    A = A.as_expr().subs(result)
    B = B.as_expr().subs(result)

    rat_part = cancel(A/u.as_expr(), x)
    log_part = cancel(B/v.as_expr(), x)

    return rat_part, log_part

def _is_function_class_equation(func_class, f, symbol):
    """ Tests whether the equation is an equation of the given function class.

    The given equation belongs to the given function class if it is
    comprised of functions of the function class which are multiplied by
    or added to expressions independent of the symbol. In addition, the
    arguments of all such functions must be linear in the symbol as well.

    Examples
    ========

    >>> from sympy.solvers.solveset import _is_function_class_equation
    >>> from sympy import tan, sin, tanh, sinh, exp
    >>> from sympy.abc import x
    >>> from sympy.functions.elementary.trigonometric import (TrigonometricFunction,
    ... HyperbolicFunction)
    >>> _is_function_class_equation(TrigonometricFunction, exp(x) + tan(x), x)
    False
    >>> _is_function_class_equation(TrigonometricFunction, tan(x) + sin(x), x)
    True
    >>> _is_function_class_equation(TrigonometricFunction, tan(x**2), x)
    False
    >>> _is_function_class_equation(TrigonometricFunction, tan(x + 2), x)
    True
    >>> _is_function_class_equation(HyperbolicFunction, tanh(x) + sinh(x), x)
    True
    """
    if f.is_Mul or f.is_Add:
        return all(_is_function_class_equation(func_class, arg, symbol)
                   for arg in f.args)

    if f.is_Pow:
        if not f.exp.has(symbol):
            return _is_function_class_equation(func_class, f.base, symbol)
        else:
            return False

    if not f.has(symbol):
        return True

    if isinstance(f, func_class):
        try:
            g = Poly(f.args[0], symbol)
            return g.degree() <= 1
        except PolynomialError:
            return False
    else:
        return False

def cancel(f, *symbols):
    """Cancel common factors in a given rational function.

       Given a quotient of polynomials, performing only gcd and quo
       operations in polynomial algebra,  return rational function
       with numerator and denominator of minimal total degree in
       an expanded form.

       For all other kinds of expressions the input is returned in
       an unchanged form. Note however, that 'cancel' function can
       thread over sums and relational operators.

       Additionally you can specify a list of variables to perform
       cancelation more efficiently using only those symbols.

       >>> from sympy import *
       >>> x,y = symbols('xy')

       >>> cancel((x**2-1)/(x-1))
       1 + x

       >>> cancel((x**2-y**2)/(x-y), x)
       x + y

       >>> cancel((x**2-2)/(x+sqrt(2)))
       x - 2**(1/2)

    """
    return Poly.cancel(f, *symbols)

def _is_negative(term):
    if getattr(term, 'is_number', False):
        return term < 0
    elif isinstance(term, Pow):
        if term.args[1]%2==1:
            return _is_negative(term.args[0])
        else:
            return False

    l = len(term.free_symbols)
    t = Poly(term).as_dict()
    if (all(c < 0 for c in t.values()) and
        all(i % 2 == 0 for d in t.keys() for i in d) and
        (0,)*l in t.keys()):
        return True
    return ask_is_negative(term)

def reduce_poly_inequalities(exprs, gen, assume=True, relational=True):
    """Reduce a system of polynomial inequalities with rational coefficients. """
    exact = True
    polys = []

    for _exprs in exprs:
        _polys = []

        for expr in _exprs:
            if isinstance(expr, tuple):
                expr, rel = expr
            else:
                if expr.is_Relational:
                    expr, rel = expr.lhs - expr.rhs, expr.rel_op
                else:
                    expr, rel = expr, "=="

            poly = Poly(expr, gen)

            if not poly.get_domain().is_Exact:
                poly, exact = poly.to_exact(), False

            domain = poly.get_domain()

            if not (domain.is_ZZ or domain.is_QQ):
                raise NotImplementedError("inequality solving is not supported over %s" % domain)

            _polys.append((poly, rel))

        polys.append(_polys)

    solution = solve_poly_inequalities(polys)

    if not exact:
        solution = solution.evalf()

    if not relational:
        return solution

    real = ask(Q.real(gen), assumptions=assume)

    if not real:
        result = And(solution.as_relational(re(gen)), Eq(im(gen), 0))
    else:
        result = solution.as_relational(gen)

    return result

def _solve_inequality(ie, s):
    """ A hacky replacement for solve, since the latter only works for
        univariate inequalities. """
    from sympy import Poly
    if not ie.rel_op in ('>', '>=', '<', '<='):
        raise NotImplementedError
    expr = ie.lhs - ie.rhs
    p = Poly(expr, s)
    if p.degree() != 1:
        raise NotImplementedError('%s' % ie)
    a, b = p.all_coeffs()
    if a.is_positive:
        return ie.func(s, -b/a)
    elif a.is_negative:
        return ie.func(-b/a, s)
    else:
        raise NotImplementedError

def weak_normalizer(a, d, DE, z=None):
    """
    Weak normalization.

    Given a derivation D on k[t] and f == a/d in k(t), return q in k[t]
    such that f - Dq/q is weakly normalized with respect to t.

    f in k(t) is said to be "weakly normalized" with respect to t if
    residue_p(f) is not a positive integer for any normal irreducible p
    in k[t] such that f is in R_p (Definition 6.1.1).  If f has an
    elementary integral, this is equivalent to no logarithm of
    integral(f) whose argument depends on t has a positive integer
    coefficient, where the arguments of the logarithms not in k(t) are
    in k[t].

    Returns (q, f - Dq/q)
    """
    z = z or Dummy('z')
    dn, ds = splitfactor(d, DE)

    # Compute d1, where dn == d1*d2**2*...*dn**n is a square-free
    # factorization of d.
    g = gcd(dn, dn.diff(DE.t))
    d_sqf_part = dn.quo(g)
    d1 = d_sqf_part.quo(gcd(d_sqf_part, g))

    a1, b = gcdex_diophantine(d.quo(d1).as_poly(DE.t), d1.as_poly(DE.t),
        a.as_poly(DE.t))
    r = (a - Poly(z, DE.t)*derivation(d1, DE)).as_poly(DE.t).resultant(
        d1.as_poly(DE.t))
    r = Poly(r, z)

    if not r.has(z):
        return (Poly(1, DE.t), (a, d))

    N = [i for i in r.real_roots() if i in ZZ and i > 0]

    q = reduce(mul, [gcd(a - Poly(n, DE.t)*derivation(d1, DE), d1) for n in N],
        Poly(1, DE.t))

    dq = derivation(q, DE)
    sn = q*a - d*dq
    sd = q*d
    sn, sd = sn.cancel(sd, include=True)

    return (q, (sn, sd))

def test_roots_quadratic():
    assert roots_quadratic(Poly(2*x**2, x)) == [0, 0]
    assert roots_quadratic(Poly(2*x**2 + 3*x, x)) == [-Rational(3, 2), 0]
    assert roots_quadratic(Poly(2*x**2 + 3, x)) == [-I*sqrt(6)/2, I*sqrt(6)/2]
    assert roots_quadratic(Poly(2*x**2 + 4*x + 3, x)) == [-1 - I*sqrt(2)/2, -1 + I*sqrt(2)/2]

    f = x**2 + (2*a*e + 2*c*e)/(a - c)*x + (d - b + a*e**2 - c*e**2)/(a - c)

    assert roots_quadratic(Poly(f, x)) == \
        [-e*(a + c)/(a - c) - sqrt((a*b + c*d - a*d - b*c + 4*a*c*e**2)/(a - c)**2),
         -e*(a + c)/(a - c) + sqrt((a*b + c*d - a*d - b*c + 4*a*c*e**2)/(a - c)**2)]

    # check for simplification
    f = Poly(y*x**2 - 2*x - 2*y, x)
    assert roots_quadratic(f) == \
        [-sqrt(2*y**2 + 1)/y + 1/y, sqrt(2*y**2 + 1)/y + 1/y]
    f = Poly(x**2 + (-y**2 - 2)*x + y**2 + 1, x)
    assert roots_quadratic(f) == \
        [y**2/2 - sqrt(y**4)/2 + 1, y**2/2 + sqrt(y**4)/2 + 1]

    f = Poly(sqrt(2)*x**2 - 1, x)
    r = roots_quadratic(f)
    assert r == _nsort(r)

    # issue 8255
    f = Poly(-24*x**2 - 180*x + 264)
    assert [w.n(2) for w in f.all_roots(radicals=True)] == \
           [w.n(2) for w in f.all_roots(radicals=False)]
    for _a, _b, _c in cartes((-2, 2), (-2, 2), (0, -1)):
        f = Poly(_a*x**2 + _b*x + _c)
        roots = roots_quadratic(f)
        assert roots == _nsort(roots)

def root_factors(f, *gens, **args):
    """Returns all factors of a univariate polynomial.

       Examples
       ========

           >>> from sympy.abc import x, y
           >>> from sympy.polys.polyroots import root_factors

           >>> root_factors(x**2-y, x)
           [x - y**(1/2), x + y**(1/2)]

    """
    F = Poly(f, *gens, **args)

    if not F.is_Poly:
        return [f]

    if F.is_multivariate:
        raise ValueError('multivariate polynomials not supported')

    x = F.gens[0]

    if 'multiple' in args:
        del args['multiple']

    zeros = roots(F, **args)

    if not zeros:
        factors = [F]
    else:
        factors, N = [], 0

        for r, n in zeros.iteritems():
            factors, N = factors + [Poly(x-r, x)]*n, N + n

        if N < F.degree():
            g = reduce(lambda p,q: p*q, factors)
            factors.append(f.exquo(g))

    if not isinstance(f, Poly):
        return [ f.as_basic() for f in factors ]
    else:
        return factors

def _solve_inequality(ie, s):
    """ A hacky replacement for solve, since the latter only works for
        univariate inequalities. """
    if not ie.rel_op in (">", ">=", "<", "<="):
        raise NotImplementedError
    expr = ie.lhs - ie.rhs
    try:
        p = Poly(expr, s)
    except PolynomialError:
        raise NotImplementedError
    if p.degree() != 1:
        raise NotImplementedError("%s" % ie)
    a, b = p.all_coeffs()
    if a.is_positive:
        return ie.func(s, -b / a)
    elif a.is_negative:
        return ie.func(-b / a, s)
    else:
        raise NotImplementedError