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),
         -e*(a + c)/(a - c) + sqrt((a*b + c*d - a*d - b*c + 4*a*c*e**2))/(a - c)]

    # 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) == \
        [1,y**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 test_issue_8289():
    roots = (Poly(x**2 + 2)*Poly(x**4 + 2)).all_roots()
    assert roots == _nsort(roots)
    roots = Poly(x**6 + 3*x**3 + 2, x).all_roots()
    assert roots == _nsort(roots)
    roots = Poly(x**6 - x + 1).all_roots()
    assert roots == _nsort(roots)
    # all imaginary roots
    roots = Poly(x**4 + 4*x**2 + 4, x).all_roots()
    assert roots == _nsort(roots)

def test_issue_8285():
    roots = (Poly(4*x**8 - 1, x)*Poly(x**2 + 1)).all_roots()
    assert roots == _nsort(roots)
    f = Poly(x**4 + 5*x**2 + 6, x)
    ro = [RootOf(f, i) for i in range(4)]
    roots = Poly(x**4 + 5*x**2 + 6, x).all_roots()
    assert roots == ro
    assert roots == _nsort(roots)
    # more than 2 complex roots from which to identify the
    # imaginary ones
    roots = Poly(2*x**8 - 1).all_roots()
    assert roots == _nsort(roots)
    assert len(Poly(2*x**10 - 1).all_roots()) == 10  # doesn't fail

def test_roots_binomial():
    assert roots_binomial(Poly(5*x, x)) == [0]
    assert roots_binomial(Poly(5*x**4, x)) == [0, 0, 0, 0]
    assert roots_binomial(Poly(5*x + 2, x)) == [-Rational(2, 5)]

    A = 10**Rational(3, 4)/10

    assert roots_binomial(Poly(5*x**4 + 2, x)) == \
        [-A - A*I, -A + A*I, A - A*I, A + A*I]

    a1 = Symbol('a1', nonnegative=True)
    b1 = Symbol('b1', nonnegative=True)

    r0 = roots_quadratic(Poly(a1*x**2 + b1, x))
    r1 = roots_binomial(Poly(a1*x**2 + b1, x))

    assert powsimp(r0[0]) == powsimp(r1[0])
    assert powsimp(r0[1]) == powsimp(r1[1])
    for a, b, s, n in cartes((1, 2), (1, 2), (-1, 1), (2, 3, 4, 5)):
        if a == b and a != 1:  # a == b == 1 is sufficient
            continue
        p = Poly(a*x**n + s*b)
        ans = roots_binomial(p)
        assert ans == _nsort(ans)

    # issue 8813
    assert roots(Poly(2*x**3 - 16*y**3, x)) == {
        2*y*(-S(1)/2 - sqrt(3)*I/2): 1,
        2*y: 1,
        2*y*(-S(1)/2 + sqrt(3)*I/2): 1}

def test__nsort():
    # issue 6137
    r = S('''[3/2 + sqrt(-14/3 - 2*(-415/216 + 13*I/12)**(1/3) - 4/sqrt(-7/3 +
    61/(18*(-415/216 + 13*I/12)**(1/3)) + 2*(-415/216 + 13*I/12)**(1/3)) -
    61/(18*(-415/216 + 13*I/12)**(1/3)))/2 - sqrt(-7/3 + 61/(18*(-415/216
    + 13*I/12)**(1/3)) + 2*(-415/216 + 13*I/12)**(1/3))/2, 3/2 - sqrt(-7/3
    + 61/(18*(-415/216 + 13*I/12)**(1/3)) + 2*(-415/216 +
    13*I/12)**(1/3))/2 - sqrt(-14/3 - 2*(-415/216 + 13*I/12)**(1/3) -
    4/sqrt(-7/3 + 61/(18*(-415/216 + 13*I/12)**(1/3)) + 2*(-415/216 +
    13*I/12)**(1/3)) - 61/(18*(-415/216 + 13*I/12)**(1/3)))/2, 3/2 +
    sqrt(-14/3 - 2*(-415/216 + 13*I/12)**(1/3) + 4/sqrt(-7/3 +
    61/(18*(-415/216 + 13*I/12)**(1/3)) + 2*(-415/216 + 13*I/12)**(1/3)) -
    61/(18*(-415/216 + 13*I/12)**(1/3)))/2 + sqrt(-7/3 + 61/(18*(-415/216
    + 13*I/12)**(1/3)) + 2*(-415/216 + 13*I/12)**(1/3))/2, 3/2 + sqrt(-7/3
    + 61/(18*(-415/216 + 13*I/12)**(1/3)) + 2*(-415/216 +
    13*I/12)**(1/3))/2 - sqrt(-14/3 - 2*(-415/216 + 13*I/12)**(1/3) +
    4/sqrt(-7/3 + 61/(18*(-415/216 + 13*I/12)**(1/3)) + 2*(-415/216 +
    13*I/12)**(1/3)) - 61/(18*(-415/216 + 13*I/12)**(1/3)))/2]''')
    ans = [r[1], r[0], r[-1], r[-2]]
    assert _nsort(r) == ans
    assert len(_nsort(r, separated=True)[0]) == 0
    b, c, a = exp(-1000), exp(-999), exp(-1001)
    assert _nsort((b, c, a)) == [a, b, c]

    def normal_lines(self, p, prec=None):
        """Normal lines between `p` and the ellipse.

        Parameters
        ==========

        p : Point

        Returns
        =======

        normal_lines : list with 1, 2 or 4 Lines

        Examples
        ========

        >>> from sympy import Line, Point, Ellipse
        >>> e = Ellipse((0, 0), 2, 3)
        >>> c = e.center
        >>> e.normal_lines(c + Point(1, 0))
        [Line2D(Point2D(0, 0), Point2D(1, 0))]
        >>> e.normal_lines(c)
        [Line2D(Point2D(0, 0), Point2D(0, 1)), Line2D(Point2D(0, 0), Point2D(1, 0))]

        Off-axis points require the solution of a quartic equation. This
        often leads to very large expressions that may be of little practical
        use. An approximate solution of `prec` digits can be obtained by
        passing in the desired value:

        >>> e.normal_lines((3, 3), prec=2)
        [Line2D(Point2D(-0.81, -2.7), Point2D(0.19, -1.2)),
        Line2D(Point2D(1.5, -2.0), Point2D(2.5, -2.7))]

        Whereas the above solution has an operation count of 12, the exact
        solution has an operation count of 2020.
        """
        p = Point(p, dim=2)

        # XXX change True to something like self.angle == 0 if the arbitrarily
        # rotated ellipse is introduced.
        # https://github.com/sympy/sympy/issues/2815)
        if True:
            rv = []
            if p.x == self.center.x:
                rv.append(Line(self.center, slope=oo))
            if p.y == self.center.y:
                rv.append(Line(self.center, slope=0))
            if rv:
                # at these special orientations of p either 1 or 2 normals
                # exist and we are done
                return rv

        # find the 4 normal points and construct lines through them with
        # the corresponding slope
        x, y = Dummy('x', real=True), Dummy('y', real=True)
        eq = self.equation(x, y)
        dydx = idiff(eq, y, x)
        norm = -1/dydx
        slope = Line(p, (x, y)).slope
        seq = slope - norm

        # TODO: Replace solve with solveset, when this line is tested
        yis = solve(seq, y)[0]
        xeq = eq.subs(y, yis).as_numer_denom()[0].expand()
        if len(xeq.free_symbols) == 1:
            try:
                # this is so much faster, it's worth a try
                xsol = Poly(xeq, x).real_roots()
            except (DomainError, PolynomialError, NotImplementedError):
                # TODO: Replace solve with solveset, when these lines are tested
                xsol = _nsort(solve(xeq, x), separated=True)[0]
            points = [Point(i, solve(eq.subs(x, i), y)[0]) for i in xsol]
        else:
            raise NotImplementedError(
                'intersections for the general ellipse are not supported')
        slopes = [norm.subs(zip((x, y), pt.args)) for pt in points]
        if prec is not None:
            points = [pt.n(prec) for pt in points]
            slopes = [i if _not_a_coeff(i) else i.n(prec) for i in slopes]
        return [Line(pt, slope=s) for pt,s in zip(points, slopes)]

def solve_univariate_inequality(expr, gen, relational=True):
    """Solves a real univariate inequality.

    Examples
    ========

    >>> from sympy.solvers.inequalities import solve_univariate_inequality
    >>> from sympy.core.symbol import Symbol
    >>> x = Symbol('x')

    >>> solve_univariate_inequality(x**2 >= 4, x)
    Or(And(-oo < x, x <= -2), And(2 <= x, x < oo))

    >>> solve_univariate_inequality(x**2 >= 4, x, relational=False)
    (-oo, -2] U [2, oo)

    """

    from sympy.solvers.solvers import solve, denoms

    # This keeps the function independent of the assumptions about `gen`.
    # `solveset` makes sure this function is called only when the domain is
    # real.
    d = Dummy(real=True)
    expr = expr.subs(gen, d)
    _gen = gen
    gen = d

    if expr is S.true:
        rv = S.Reals
    elif expr is S.false:
        rv = S.EmptySet
    else:
        e = expr.lhs - expr.rhs
        parts = n, d = e.as_numer_denom()
        if all(i.is_polynomial(gen) for i in parts):
            solns = solve(n, gen, check=False)
            singularities = solve(d, gen, check=False)
        else:
            solns = solve(e, gen, check=False)
            singularities = []
            for d in denoms(e):
                singularities.extend(solve(d, gen))

        include_x = expr.func(0, 0)

        def valid(x):
            v = e.subs(gen, x)
            try:
                r = expr.func(v, 0)
            except TypeError:
                r = S.false
            if r in (S.true, S.false):
                return r
            if v.is_real is False:
                return S.false
            else:
                v = v.n(2)
                if v.is_comparable:
                    return expr.func(v, 0)
                return S.false

        start = S.NegativeInfinity
        sol_sets = [S.EmptySet]
        try:
            reals = _nsort(set(solns + singularities), separated=True)[0]
        except NotImplementedError:
            raise NotImplementedError('sorting of these roots is not supported')
        for x in reals:
            end = x

            if end in [S.NegativeInfinity, S.Infinity]:
                if valid(S(0)):
                    sol_sets.append(Interval(start, S.Infinity, True, True))
                    break

            pt = ((start + end)/2 if start is not S.NegativeInfinity else
                (end/2 if end.is_positive else
                (2*end if end.is_negative else
                end - 1)))
            if valid(pt):
                sol_sets.append(Interval(start, end, True, True))

            if x in singularities:
                singularities.remove(x)
            elif include_x:
                sol_sets.append(FiniteSet(x))

            start = end

        end = S.Infinity

        # in case start == -oo then there were no solutions so we just
        # check a point between -oo and oo (e.g. 0) else pick a point
        # past the last solution (which is start after the end of the
        # for-loop above
        pt = (0 if start is S.NegativeInfinity else
            (start/2 if start.is_negative else
            (2*start if start.is_positive else
            start + 1)))
#.........这里部分代码省略.........

def solve_univariate_inequality(expr, gen, relational=True):
    """Solves a real univariate inequality.

    Examples
    ========

    >>> from sympy.solvers.inequalities import solve_univariate_inequality
    >>> from sympy.core.symbol import Symbol
    >>> x = Symbol('x', real=True)

    >>> solve_univariate_inequality(x**2 >= 4, x)
    Or(And(-oo < x, x <= -2), And(2 <= x, x < oo))

    >>> solve_univariate_inequality(x**2 >= 4, x, relational=False)
    (-oo, -2] U [2, oo)

    """

    from sympy.solvers.solvers import solve, denoms

    e = expr.lhs - expr.rhs
    parts = n, d = e.as_numer_denom()
    if all(i.is_polynomial(gen) for i in parts):
        solns = solve(n, gen, check=False)
        singularities = solve(d, gen, check=False)
    else:
        solns = solve(e, gen, check=False)
        singularities = []
        for d in denoms(e):
            singularities.extend(solve(d, gen))

    include_x = expr.func(0, 0)

    def valid(x):
        v = e.subs(gen, x)
        try:
            r = expr.func(v, 0)
        except TypeError:
            r = S.false
        if r in (S.true, S.false):
            return r
        if v.is_real is False:
            return S.false
        else:
            v = v.n(2)
            if v.is_comparable:
                return expr.func(v, 0)
            return S.false

    start = S.NegativeInfinity
    sol_sets = [S.EmptySet]
    try:
        reals = _nsort(set(solns + singularities), separated=True)[0]
    except NotImplementedError:
        raise NotImplementedError('sorting of these roots is not supported')
    for x in reals:
        end = x

        if end in [S.NegativeInfinity, S.Infinity]:
            if valid(S(0)):
                sol_sets.append(Interval(start, S.Infinity, True, True))
                break

        if valid((start + end)/2 if start != S.NegativeInfinity else end - 1):
            sol_sets.append(Interval(start, end, True, True))

        if x in singularities:
            singularities.remove(x)
        elif include_x:
            sol_sets.append(FiniteSet(x))

        start = end

    end = S.Infinity

    if valid(start + 1):
        sol_sets.append(Interval(start, end, True, True))

    rv = Union(*sol_sets)
    return rv if not relational else rv.as_relational(gen)

#.........这里部分代码省略.........
                try:
                    r = expr.func(v, 0)
                except TypeError:
                    r = S.false
                if r in (S.true, S.false):
                    return r
                if v.is_real is False:
                    return S.false
                else:
                    v = v.n(2)
                    if v.is_comparable:
                        return expr.func(v, 0)
                    # not comparable or couldn't be evaluated
                    raise NotImplementedError(
                        'relationship did not evaluate: %s' % r)

            singularities = []
            for d in denoms(expr, gen):
                singularities.extend(solvify(d, gen, domain))
            if not continuous:
                domain = continuous_domain(e, gen, domain)

            include_x = '=' in expr.rel_op and expr.rel_op != '!='

            try:
                discontinuities = set(domain.boundary -
                    FiniteSet(domain.inf, domain.sup))
                # remove points that are not between inf and sup of domain
                critical_points = FiniteSet(*(solns + singularities + list(
                    discontinuities))).intersection(
                    Interval(domain.inf, domain.sup,
                    domain.inf not in domain, domain.sup not in domain))
                if all(r.is_number for r in critical_points):
                    reals = _nsort(critical_points, separated=True)[0]
                else:
                    from sympy.utilities.iterables import sift
                    sifted = sift(critical_points, lambda x: x.is_real)
                    if sifted[None]:
                        # there were some roots that weren't known
                        # to be real
                        raise NotImplementedError
                    try:
                        reals = sifted[True]
                        if len(reals) > 1:
                            reals = list(sorted(reals))
                    except TypeError:
                        raise NotImplementedError
            except NotImplementedError:
                raise NotImplementedError('sorting of these roots is not supported')

            #If expr contains imaginary coefficients
            #Only real values of x for which the imaginary part is 0 are taken
            make_real = S.Reals
            if im(expanded_e) != S.Zero:
                check = True
                im_sol = FiniteSet()
                try:
                    a = solveset(im(expanded_e), gen, domain)
                    if not isinstance(a, Interval):
                        for z in a:
                            if z not in singularities and valid(z) and z.is_real:
                                im_sol += FiniteSet(z)
                    else:
                        start, end = a.inf, a.sup
                        for z in _nsort(critical_points + FiniteSet(end)):
                            valid_start = valid(start)

#.........这里部分代码省略.........
                if expr.func(frange.sup, 0):
                    rv = domain
                elif not expr.func(frange.inf, 0):
                    rv = S.EmptySet

            elif rel == '>' or rel == '>=':
                if expr.func(frange.inf, 0):
                    rv = domain
                elif not expr.func(frange.sup, 0):
                    rv = S.EmptySet

            inf, sup = domain.inf, domain.sup
            if sup - inf is S.Infinity:
                domain = Interval(0, period, False, True)

        if rv is None:
            singularities = []
            for d in denoms(e):
                singularities.extend(solvify(d, gen, domain))
            if not continuous:
                domain = continuous_domain(e, gen, domain)
            solns = solvify(e, gen, domain)

            if solns is None:
                raise NotImplementedError(filldedent('''The inequality cannot be
                    solved using solve_univariate_inequality.'''))

            include_x = expr.func(0, 0)

            def valid(x):
                v = e.subs(gen, x)
                try:
                    r = expr.func(v, 0)
                except TypeError:
                    r = S.false
                if r in (S.true, S.false):
                    return r
                if v.is_real is False:
                    return S.false
                else:
                    v = v.n(2)
                    if v.is_comparable:
                        return expr.func(v, 0)
                    return S.false

            start = domain.inf
            sol_sets = [S.EmptySet]
            try:
                discontinuities = domain.boundary - FiniteSet(domain.inf, domain.sup)
                critical_points = set(solns + singularities + list(discontinuities))
                reals = _nsort(critical_points, separated=True)[0]

            except NotImplementedError:
                raise NotImplementedError('sorting of these roots is not supported')

            if valid(start) and start.is_finite:
                sol_sets.append(FiniteSet(start))

            for x in reals:
                end = x

                if end in [S.NegativeInfinity, S.Infinity]:
                    if valid(S(0)):
                        sol_sets.append(Interval(start, S.Infinity, True, True))
                        break

                pt = ((start + end)/2 if start is not S.NegativeInfinity else
                    (end/2 if end.is_positive else
                    (2*end if end.is_negative else
                    end - 1)))
                if valid(pt):
                    sol_sets.append(Interval(start, end, True, True))

                if x in singularities:
                    singularities.remove(x)
                elif include_x:
                    sol_sets.append(FiniteSet(x))

                start = end

            end = domain.sup

            # in case start == -oo then there were no solutions so we just
            # check a point between -oo and oo (e.g. 0) else pick a point
            # past the last solution (which is start after the end of the
            # for-loop above
            pt = (0 if start is S.NegativeInfinity else
                (start/2 if start.is_negative else
                (2*start if start.is_positive else
                start + 1)))

            if pt >= end:
                pt = (start + end)/2

            if valid(pt):
                sol_sets.append(Interval(start, end, True, True))

            rv = Union(*sol_sets).subs(gen, _gen)

    return rv if not relational else rv.as_relational(_gen)