solveset_univariate_trig_ineq.py

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

    Examples
    ========

    >>> from sympy.solvers.inequalities import solveset_univariate_trig_inequality
    >>> from sympy.core.symbol import Symbol
    >>> x = Symbol('x')
    >>> solveset((2*cos(x)+1)/(2*cos(x)-1) > 0, x, S.Reals)
 
⎡    ⎧        π        ⎫      ⎧        π        ⎫⎤   ⎛    ⎧        2⋅π        ⎫      ⎧        4⋅π        ⎫⎞
⎢x > ⎨2⋅n⋅π - ─ | n ∊ ℤ⎬, x < ⎨2⋅n⋅π + ─ | n ∊ ℤ⎬⎥ ∪ ⎜x > ⎨2⋅n⋅π + ─── | n ∊ ℤ⎬, x < ⎨2⋅n⋅π + ─── | n ∊ ℤ⎬⎟
⎣    ⎩        3        ⎭      ⎩        3        ⎭⎦   ⎝    ⎩         3         ⎭      ⎩         3         ⎭⎠


    """
    from sympy.solvers.solveset import solveset_real
    from sympy.calculus.singularities import singularities
    from sympy.simplify.simplify import trigsimp
    from sympy.sets import ImageSet
    from sympy.core.function import Lambda

    # This keeps the function independent of the assumptions about `gen`.
    # `solveset` makes sure this function is called only when the domain is
    # real.
    def _extract_expr(imgs, n_new=None):
        expr = []
        expr_extended = {}
        change_dummy = False
        if n_new:
            change_dummy = True
        if isinstance(imgs, ImageSet):
            # value at n = 0
            val = imgs.at(0)
            expr.append(val)
            # there may be a problem when more than 1 Dummy is there
            n_old = (imgs.lamda.expr).atoms(Dummy).pop()
            if change_dummy:
                expr_extended[val] = imgs.subs(n_old, n_new)
            else:
                expr_extended[val] = imgs
        elif isinstance(imgs, Union):
            img_len = len(imgs.args)
            for i in range(0, img_len):
                arg = imgs.args[i]
                val = arg.at(0)
                # list of  expr at n = 0
                expr.append(val)
                # there may be a problem when more than 1 Dummy is there
                n_old = (arg.lamda.expr).atoms(Dummy).pop()
                if change_dummy:
                    expr_extended[val] = imgs.args[i].subs(n_old, n_new)
                else:
                    expr_extended[val] = imgs.args[i]
        if not change_dummy:
            n_new = n_old
        return expr, expr_extended, n_new

    expr = trigsimp(expr)
    if expr is S.true:
        return S.Reals
    elif expr is S.false:
        return S.EmptySet
    else:
        e = expr.lhs - expr.rhs
        n, d = e.as_numer_denom()
        solns_img = solveset_real(n, gen)
        singul_img = solveset_real(d, gen)

        # extract solution expr and singul expr
        # using one dummy `n_new`
        solns, solns_extended, n_new = _extract_expr(solns_img)
        singul, singul_extended,_ = _extract_expr(singul_img, n_new)

        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 = []
        used_reals = None
        try:
            reals = _nsort(set(solns + singul), separated=True)[0]
            used_reals = dict(zip(reals, [ False for i in range(0,len(reals))]))
        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):
                # check  where start and end points are.
                # True means in soln otherwise it is in singul
                start_in_solns = start in solns
                end_in_solns = end in solns

                if start is S.NegativeInfinity:
                    left = start < gen
                elif start_in_solns:
                    left = (solns_extended[start] < gen)
                    used_reals[start] = True
                else:
                    left = (singul_extended[start] < gen)
                    used_reals[start] = True
                if end_in_solns:
                    right = (gen < solns_extended[end])
                    used_reals[end] = True
                else:
                    right = (gen < singul_extended[end])
                    used_reals[end] = True
                sol_sets.append(Interval(left, right, True, True))

            if x in singul:
                singul.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)))
        if valid(pt):
            if start_in_solns:
                left = (solns_extended[start] < gen)
                used_reals[start] = True
            else:
                left = (singul_extended[start] < gen)
                used_reals[start] = True
            if end is S.Infinity:
                right = gen < end
            elif end_in_solns:
                right = (gen < solns_extended[end])
                used_reals[end] = True
            else:
                right = (gen < singul_extended[end])
                used_reals[end] = True
            sol_sets.append(Interval(left, right, True, True))

        # unused values
        unused = []
        for k, v in used_reals.items():
            if v is False:
                unused.append(k)
        # check for -oo and oo and replace with proper Imagset
        first_interval_lhs = sol_sets[0].left.lhs
        if first_interval_lhs is S.NegativeInfinity:
            right = sol_sets[0].right
            for u in unused:
                u_in_solns = u in solns
                img = None
                if u_in_solns:
                    img = solns_extended[u]
                    prev = img.at(-1)
                    if prev < right.rhs.args[0](0):
                        unused.remove(u)
                        exp = img.lamda.expr - u + prev
                        base = img.base_set
                        left = ImageSet(Lambda(n_new, exp), base) < gen
                        sol_sets[0] = Interval(left, right)
                else:
                    img = singul_extended[u]
                    prev = img.at(-1)
                    if prev < right.rhs.args[0](0):
                        unused.remove(u)
                        exp = img.lamda.expr - u + prev
                        base = img.base_set
                        left = ImageSet(Lambda(n_new, exp), base) < gen
                        sol_sets[0] = Interval(left, right)

        # similarly for oo at last Interval
        last_interval_rhs = sol_sets[-1].right.rhs
        if last_interval_rhs is S.Infinity:
            left = sol_sets[-1].lhs
            for u in unused:
                u_in_solns = u in solns
                img = None
                if u_in_solns:
                    img = solns_extended[u]
                    next = img.at(1)
                    if next > left.lhs.args[1](0):
                        unused.remove(u)
                        exp = img.lamda.expr - u + next
                        base = img.base_set
                        right = gen < ImageSet(Lambda(n_new, exp), base)
                        sol_sets[-1] = Interval(left, right)
                else:
                    img = singul_extended[u]
                    prev = img.at(1)
                    if next > right.rhs.args[1](0):
                        unused.remove(u)
                        exp = img.lamda.expr - u + next
                        base = img.base_set
                        right = gen < ImageSet(Lambda(n_new, exp), base)
                        sol_sets[-1] = Interval(left, right)

        rv = Union(*sol_sets)

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