Source code for sympy.solvers.inequalities

"""Tools for solving inequalities and systems of inequalities. """

from __future__ import print_function, division

from sympy.core import Symbol, Dummy, sympify
from sympy.core.compatibility import iterable
from sympy.core.exprtools import factor_terms
from sympy.core.relational import Relational, Eq, Ge, Lt, Ne
from sympy.sets import Interval
from sympy.sets.sets import FiniteSet, Union, EmptySet, Intersection
from sympy.sets.fancysets import ImageSet
from sympy.core.singleton import S
from sympy.core.function import expand_mul

from sympy.functions import Abs
from sympy.logic import And
from sympy.polys import Poly, PolynomialError, parallel_poly_from_expr
from sympy.polys.polyutils import _nsort
from sympy.utilities.iterables import sift
from sympy.utilities.misc import filldedent

[docs]def solve_poly_inequality(poly, rel): """Solve a polynomial inequality with rational coefficients. Examples ======== >>> from sympy import Poly >>> from sympy.abc import x >>> from sympy.solvers.inequalities import solve_poly_inequality >>> solve_poly_inequality(Poly(x, x, domain='ZZ'), '==') [{0}] >>> solve_poly_inequality(Poly(x**2 - 1, x, domain='ZZ'), '!=') [Interval.open(-oo, -1), Interval.open(-1, 1), Interval.open(1, oo)] >>> solve_poly_inequality(Poly(x**2 - 1, x, domain='ZZ'), '==') [{-1}, {1}] See Also ======== solve_poly_inequalities """ if not isinstance(poly, Poly): raise ValueError( 'For efficiency reasons, `poly` should be a Poly instance') if poly.is_number: t = Relational(poly.as_expr(), 0, rel) if t is S.true: return [S.Reals] elif t is S.false: return [S.EmptySet] else: raise NotImplementedError( "could not determine truth value of %s" % t) reals, intervals = poly.real_roots(multiple=False), [] if rel == '==': for root, _ in reals: interval = Interval(root, root) intervals.append(interval) elif rel == '!=': left = S.NegativeInfinity for right, _ in reals + [(S.Infinity, 1)]: interval = Interval(left, right, True, True) intervals.append(interval) left = right else: if poly.LC() > 0: sign = +1 else: sign = -1 eq_sign, equal = None, False if rel == '>': eq_sign = +1 elif rel == '<': eq_sign = -1 elif rel == '>=': eq_sign, equal = +1, True elif rel == '<=': eq_sign, equal = -1, True else: raise ValueError("'%s' is not a valid relation" % rel) right, right_open = S.Infinity, True for left, multiplicity in reversed(reals): if multiplicity % 2: if sign == eq_sign: intervals.insert( 0, Interval(left, right, not equal, right_open)) sign, right, right_open = -sign, left, not equal else: if sign == eq_sign and not equal: intervals.insert( 0, Interval(left, right, True, right_open)) right, right_open = left, True elif sign != eq_sign and equal: intervals.insert(0, Interval(left, left)) if sign == eq_sign: intervals.insert( 0, Interval(S.NegativeInfinity, right, True, right_open)) return intervals
def solve_poly_inequalities(polys): """Solve polynomial inequalities with rational coefficients. Examples ======== >>> from sympy.solvers.inequalities import solve_poly_inequalities >>> from sympy.polys import Poly >>> from sympy.abc import x >>> solve_poly_inequalities((( ... Poly(x**2 - 3), ">"), ( ... Poly(-x**2 + 1), ">"))) Union(Interval.open(-oo, -sqrt(3)), Interval.open(-1, 1), Interval.open(sqrt(3), oo)) """ from sympy import Union return Union(*[s for p in polys for s in solve_poly_inequality(*p)])
[docs]def solve_rational_inequalities(eqs): """Solve a system of rational inequalities with rational coefficients. Examples ======== >>> from sympy.abc import x >>> from sympy import Poly >>> from sympy.solvers.inequalities import solve_rational_inequalities >>> solve_rational_inequalities([[ ... ((Poly(-x + 1), Poly(1, x)), '>='), ... ((Poly(-x + 1), Poly(1, x)), '<=')]]) {1} >>> solve_rational_inequalities([[ ... ((Poly(x), Poly(1, x)), '!='), ... ((Poly(-x + 1), Poly(1, x)), '>=')]]) Union(Interval.open(-oo, 0), Interval.Lopen(0, 1)) See Also ======== solve_poly_inequality """ result = S.EmptySet for _eqs in eqs: if not _eqs: continue global_intervals = [Interval(S.NegativeInfinity, S.Infinity)] for (numer, denom), rel in _eqs: numer_intervals = solve_poly_inequality(numer*denom, rel) denom_intervals = solve_poly_inequality(denom, '==') intervals = [] for numer_interval in numer_intervals: for global_interval in global_intervals: interval = numer_interval.intersect(global_interval) if interval is not S.EmptySet: intervals.append(interval) global_intervals = intervals intervals = [] for global_interval in global_intervals: for denom_interval in denom_intervals: global_interval -= denom_interval if global_interval is not S.EmptySet: intervals.append(global_interval) global_intervals = intervals if not global_intervals: break for interval in global_intervals: result = result.union(interval) return result
[docs]def reduce_rational_inequalities(exprs, gen, relational=True): """Reduce a system of rational inequalities with rational coefficients. Examples ======== >>> from sympy import Poly, Symbol >>> from sympy.solvers.inequalities import reduce_rational_inequalities >>> x = Symbol('x', real=True) >>> reduce_rational_inequalities([[x**2 <= 0]], x) Eq(x, 0) >>> reduce_rational_inequalities([[x + 2 > 0]], x) (-2 < x) & (x < oo) >>> reduce_rational_inequalities([[(x + 2, ">")]], x) (-2 < x) & (x < oo) >>> reduce_rational_inequalities([[x + 2]], x) Eq(x, -2) """ exact = True eqs = [] solution = S.Reals if exprs else S.EmptySet for _exprs in exprs: _eqs = [] 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, '==' if expr is S.true: numer, denom, rel = S.Zero, S.One, '==' elif expr is S.false: numer, denom, rel = S.One, S.One, '==' else: numer, denom = expr.together().as_numer_denom() try: (numer, denom), opt = parallel_poly_from_expr( (numer, denom), gen) except PolynomialError: raise PolynomialError(filldedent(''' only polynomials and rational functions are supported in this context. ''')) if not opt.domain.is_Exact: numer, denom, exact = numer.to_exact(), denom.to_exact(), False domain = opt.domain.get_exact() if not (domain.is_ZZ or domain.is_QQ): expr = numer/denom expr = Relational(expr, 0, rel) solution &= solve_univariate_inequality(expr, gen, relational=False) else: _eqs.append(((numer, denom), rel)) if _eqs: eqs.append(_eqs) if eqs: solution &= solve_rational_inequalities(eqs) exclude = solve_rational_inequalities([[((d, d.one), '==') for i in eqs for ((n, d), _) in i if d.has(gen)]]) solution -= exclude if not exact and solution: solution = solution.evalf() if relational: solution = solution.as_relational(gen) return solution
[docs]def reduce_abs_inequality(expr, rel, gen): """Reduce an inequality with nested absolute values. Examples ======== >>> from sympy import Abs, Symbol >>> from sympy.solvers.inequalities import reduce_abs_inequality >>> x = Symbol('x', real=True) >>> reduce_abs_inequality(Abs(x - 5) - 3, '<', x) (2 < x) & (x < 8) >>> reduce_abs_inequality(Abs(x + 2)*3 - 13, '<', x) (-19/3 < x) & (x < 7/3) See Also ======== reduce_abs_inequalities """ if gen.is_real is False: raise TypeError(filldedent(''' can't solve inequalities with absolute values containing non-real variables. ''')) def _bottom_up_scan(expr): exprs = [] if expr.is_Add or expr.is_Mul: op = expr.func for arg in expr.args: _exprs = _bottom_up_scan(arg) if not exprs: exprs = _exprs else: args = [] for expr, conds in exprs: for _expr, _conds in _exprs: args.append((op(expr, _expr), conds + _conds)) exprs = args elif expr.is_Pow: n = expr.exp if not n.is_Integer: raise ValueError("Only Integer Powers are allowed on Abs.") _exprs = _bottom_up_scan(expr.base) for expr, conds in _exprs: exprs.append((expr**n, conds)) elif isinstance(expr, Abs): _exprs = _bottom_up_scan(expr.args[0]) for expr, conds in _exprs: exprs.append(( expr, conds + [Ge(expr, 0)])) exprs.append((-expr, conds + [Lt(expr, 0)])) else: exprs = [(expr, [])] return exprs exprs = _bottom_up_scan(expr) mapping = {'<': '>', '<=': '>='} inequalities = [] for expr, conds in exprs: if rel not in mapping.keys(): expr = Relational( expr, 0, rel) else: expr = Relational(-expr, 0, mapping[rel]) inequalities.append([expr] + conds) return reduce_rational_inequalities(inequalities, gen)
[docs]def reduce_abs_inequalities(exprs, gen): """Reduce a system of inequalities with nested absolute values. Examples ======== >>> from sympy import Abs, Symbol >>> from sympy.abc import x >>> from sympy.solvers.inequalities import reduce_abs_inequalities >>> x = Symbol('x', real=True) >>> reduce_abs_inequalities([(Abs(3*x - 5) - 7, '<'), ... (Abs(x + 25) - 13, '>')], x) (-2/3 < x) & (x < 4) & (((-oo < x) & (x < -38)) | ((-12 < x) & (x < oo))) >>> reduce_abs_inequalities([(Abs(x - 4) + Abs(3*x - 5) - 7, '<')], x) (1/2 < x) & (x < 4) See Also ======== reduce_abs_inequality """ return And(*[ reduce_abs_inequality(expr, rel, gen) for expr, rel in exprs ])
[docs]def solve_univariate_inequality(expr, gen, relational=True, domain=S.Reals, continuous=False): """Solves a real univariate inequality. Parameters ========== expr : Relational The target inequality gen : Symbol The variable for which the inequality is solved relational : bool A Relational type output is expected or not domain : Set The domain over which the equation is solved continuous: bool True if expr is known to be continuous over the given domain (and so continuous_domain() doesn't need to be called on it) Raises ====== NotImplementedError The solution of the inequality cannot be determined due to limitation in `solvify`. Notes ===== Currently, we cannot solve all the inequalities due to limitations in `solvify`. Also, the solution returned for trigonometric inequalities are restricted in its periodic interval. See Also ======== solvify: solver returning solveset solutions with solve's output API Examples ======== >>> from sympy.solvers.inequalities import solve_univariate_inequality >>> from sympy import Symbol, sin, Interval, S >>> x = Symbol('x') >>> solve_univariate_inequality(x**2 >= 4, x) ((2 <= x) & (x < oo)) | ((x <= -2) & (-oo < x)) >>> solve_univariate_inequality(x**2 >= 4, x, relational=False) Union(Interval(-oo, -2), Interval(2, oo)) >>> domain = Interval(0, S.Infinity) >>> solve_univariate_inequality(x**2 >= 4, x, False, domain) Interval(2, oo) >>> solve_univariate_inequality(sin(x) > 0, x, relational=False) Interval.open(0, pi) """ from sympy import im from sympy.calculus.util import (continuous_domain, periodicity, function_range) from sympy.solvers.solvers import denoms from sympy.solvers.solveset import solveset_real, solvify, solveset from sympy.solvers.solvers import solve # This keeps the function independent of the assumptions about `gen`. # `solveset` makes sure this function is called only when the domain is # real. _gen = gen _domain = domain if gen.is_real is False: rv = S.EmptySet return rv if not relational else rv.as_relational(_gen) elif gen.is_real is None: gen = Dummy('gen', real=True) try: expr = expr.xreplace({_gen: gen}) except TypeError: raise TypeError(filldedent(''' When gen is real, the relational has a complex part which leads to an invalid comparison like I < 0. ''')) rv = None if expr is S.true: rv = domain elif expr is S.false: rv = S.EmptySet else: e = expr.lhs - expr.rhs period = periodicity(e, gen) if period is S.Zero: e = expand_mul(e) const = expr.func(e, 0) if const is S.true: rv = domain elif const is S.false: rv = S.EmptySet elif period is not None: frange = function_range(e, gen, domain) rel = expr.rel_op if rel == '<' or rel == '<=': 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: n, d = e.as_numer_denom() try: if gen not in n.free_symbols and len(e.free_symbols) > 1: raise ValueError # this might raise ValueError on its own # or it might give None... solns = solvify(e, gen, domain) if solns is None: # in which case we raise ValueError raise ValueError except (ValueError, NotImplementedError): # replace gen with generic x since it's # univariate anyway raise NotImplementedError(filldedent(''' The inequality, %s, cannot be solved using solve_univariate_inequality. ''' % expr.subs(gen, Symbol('x')))) expanded_e = expand_mul(e) def valid(x): # this is used to see if gen=x satisfies the # relational by substituting it into the # expanded form and testing against 0, e.g. # if expr = x*(x + 1) < 2 then e = x*(x + 1) - 2 # and expanded_e = x**2 + x - 2; the test is # whether a given value of x satisfies # x**2 + x - 2 < 0 # # expanded_e, expr and gen used from enclosing scope v = expanded_e.subs(gen, expand_mul(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) # 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(expanded_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: 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 take real # values of x for which the imaginary part is 0 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 start != end: valid_z = valid(z) pt = _pt(start, z) if pt not in singularities and pt.is_real and valid(pt): if valid_start and valid_z: im_sol += Interval(start, z) elif valid_start: im_sol += Interval.Ropen(start, z) elif valid_z: im_sol += Interval.Lopen(start, z) else: im_sol += Interval.open(start, z) start = z for s in singularities: im_sol -= FiniteSet(s) except (TypeError): im_sol = S.Reals check = False if isinstance(im_sol, EmptySet): raise ValueError(filldedent(''' %s contains imaginary parts which cannot be made 0 for any value of %s satisfying the inequality, leading to relations like I < 0. ''' % (expr.subs(gen, _gen), _gen))) make_real = make_real.intersect(im_sol) empty = sol_sets = [S.EmptySet] start = domain.inf if valid(start) and start.is_finite: sol_sets.append(FiniteSet(start)) for x in reals: end = x if valid(_pt(start, end)): sol_sets.append(Interval(start, end, True, True)) if x in singularities: singularities.remove(x) else: if x in discontinuities: discontinuities.remove(x) _valid = valid(x) else: # it's a solution _valid = include_x if _valid: sol_sets.append(FiniteSet(x)) start = end end = domain.sup if valid(end) and end.is_finite: sol_sets.append(FiniteSet(end)) if valid(_pt(start, end)): sol_sets.append(Interval.open(start, end)) if im(expanded_e) != S.Zero and check: rv = (make_real).intersect(_domain) else: rv = Intersection( (Union(*sol_sets)), make_real, _domain).subs(gen, _gen) return rv if not relational else rv.as_relational(_gen)
def _pt(start, end): """Return a point between start and end""" if not start.is_infinite and not end.is_infinite: pt = (start + end)/2 elif start.is_infinite and end.is_infinite: pt = S.Zero else: if (start.is_infinite and start.is_positive is None or end.is_infinite and end.is_positive is None): raise ValueError('cannot proceed with unsigned infinite values') if (end.is_infinite and end.is_negative or start.is_infinite and start.is_positive): start, end = end, start # if possible, use a multiple of self which has # better behavior when checking assumptions than # an expression obtained by adding or subtracting 1 if end.is_infinite: if start.is_positive: pt = start*2 elif start.is_negative: pt = start*S.Half else: pt = start + 1 elif start.is_infinite: if end.is_positive: pt = end*S.Half elif end.is_negative: pt = end*2 else: pt = end - 1 return pt def _solve_inequality(ie, s, linear=False): """Return the inequality with s isolated on the left, if possible. If the relationship is non-linear, a solution involving And or Or may be returned. False or True are returned if the relationship is never True or always True, respectively. If `linear` is True (default is False) an `s`-dependent expression will be isolated on the left, if possible but it will not be solved for `s` unless the expression is linear in `s`. Furthermore, only "safe" operations which don't change the sense of the relationship are applied: no division by an unsigned value is attempted unless the relationship involves Eq or Ne and no division by a value not known to be nonzero is ever attempted. Examples ======== >>> from sympy import Eq, Symbol >>> from sympy.solvers.inequalities import _solve_inequality as f >>> from sympy.abc import x, y For linear expressions, the symbol can be isolated: >>> f(x - 2 < 0, x) x < 2 >>> f(-x - 6 < x, x) x > -3 Sometimes nonlinear relationships will be False >>> f(x**2 + 4 < 0, x) False Or they may involve more than one region of values: >>> f(x**2 - 4 < 0, x) (-2 < x) & (x < 2) To restrict the solution to a relational, set linear=True and only the x-dependent portion will be isolated on the left: >>> f(x**2 - 4 < 0, x, linear=True) x**2 < 4 Division of only nonzero quantities is allowed, so x cannot be isolated by dividing by y: >>> y.is_nonzero is None # it is unknown whether it is 0 or not True >>> f(x*y < 1, x) x*y < 1 And while an equality (or inequality) still holds after dividing by a non-zero quantity >>> nz = Symbol('nz', nonzero=True) >>> f(Eq(x*nz, 1), x) Eq(x, 1/nz) the sign must be known for other inequalities involving > or <: >>> f(x*nz <= 1, x) nz*x <= 1 >>> p = Symbol('p', positive=True) >>> f(x*p <= 1, x) x <= 1/p When there are denominators in the original expression that are removed by expansion, conditions for them will be returned as part of the result: >>> f(x < x*(2/x - 1), x) (x < 1) & Ne(x, 0) """ from sympy.solvers.solvers import denoms if s not in ie.free_symbols: return ie if ie.rhs == s: ie = ie.reversed if ie.lhs == s and s not in ie.rhs.free_symbols: return ie def classify(ie, s, i): # return True or False if ie evaluates when substituting s with # i else None (if unevaluated) or NaN (when there is an error # in evaluating) try: v = ie.subs(s, i) if v is S.NaN: return v elif v not in (True, False): return return v except TypeError: return S.NaN rv = None oo = S.Infinity expr = ie.lhs - ie.rhs try: p = Poly(expr, s) if p.degree() == 0: rv = ie.func(p.as_expr(), 0) elif not linear and p.degree() > 1: # handle in except clause raise NotImplementedError except (PolynomialError, NotImplementedError): if not linear: try: rv = reduce_rational_inequalities([[ie]], s) except PolynomialError: rv = solve_univariate_inequality(ie, s) # remove restrictions wrt +/-oo that may have been # applied when using sets to simplify the relationship okoo = classify(ie, s, oo) if okoo is S.true and classify(rv, s, oo) is S.false: rv = rv.subs(s < oo, True) oknoo = classify(ie, s, -oo) if (oknoo is S.true and classify(rv, s, -oo) is S.false): rv = rv.subs(-oo < s, True) rv = rv.subs(s > -oo, True) if rv is S.true: rv = (s <= oo) if okoo is S.true else (s < oo) if oknoo is not S.true: rv = And(-oo < s, rv) else: p = Poly(expr) conds = [] if rv is None: e = p.as_expr() # this is in expanded form # Do a safe inversion of e, moving non-s terms # to the rhs and dividing by a nonzero factor if # the relational is Eq/Ne; for other relationals # the sign must also be positive or negative rhs = 0 b, ax = e.as_independent(s, as_Add=True) e -= b rhs -= b ef = factor_terms(e) a, e = ef.as_independent(s, as_Add=False) if (a.is_zero != False or # don't divide by potential 0 a.is_negative == a.is_positive is None and # if sign is not known then ie.rel_op not in ('!=', '==')): # reject if not Eq/Ne e = ef a = S.One rhs /= a if a.is_positive: rv = ie.func(e, rhs) else: rv = ie.reversed.func(e, rhs) # return conditions under which the value is # valid, too. beginning_denoms = denoms(ie.lhs) | denoms(ie.rhs) current_denoms = denoms(rv) for d in beginning_denoms - current_denoms: c = _solve_inequality(Eq(d, 0), s, linear=linear) if isinstance(c, Eq) and c.lhs == s: if classify(rv, s, c.rhs) is S.true: # rv is permitting this value but it shouldn't conds.append(~c) for i in (-oo, oo): if (classify(rv, s, i) is S.true and classify(ie, s, i) is not S.true): conds.append(s < i if i is oo else i < s) conds.append(rv) return And(*conds) def _reduce_inequalities(inequalities, symbols): # helper for reduce_inequalities poly_part, abs_part = {}, {} other = [] for inequality in inequalities: expr, rel = inequality.lhs, inequality.rel_op # rhs is 0 # check for gens using atoms which is more strict than free_symbols to # guard against EX domain which won't be handled by # reduce_rational_inequalities gens = expr.atoms(Symbol) if len(gens) == 1: gen = gens.pop() else: common = expr.free_symbols & symbols if len(common) == 1: gen = common.pop() other.append(_solve_inequality(Relational(expr, 0, rel), gen)) continue else: raise NotImplementedError(filldedent(''' inequality has more than one symbol of interest. ''')) if expr.is_polynomial(gen): poly_part.setdefault(gen, []).append((expr, rel)) else: components = expr.find(lambda u: u.has(gen) and ( u.is_Function or u.is_Pow and not u.exp.is_Integer)) if components and all(isinstance(i, Abs) for i in components): abs_part.setdefault(gen, []).append((expr, rel)) else: other.append(_solve_inequality(Relational(expr, 0, rel), gen)) poly_reduced = [] abs_reduced = [] for gen, exprs in poly_part.items(): poly_reduced.append(reduce_rational_inequalities([exprs], gen)) for gen, exprs in abs_part.items(): abs_reduced.append(reduce_abs_inequalities(exprs, gen)) return And(*(poly_reduced + abs_reduced + other))
[docs]def reduce_inequalities(inequalities, symbols=[]): """Reduce a system of inequalities with rational coefficients. Examples ======== >>> from sympy import sympify as S, Symbol >>> from sympy.abc import x, y >>> from sympy.solvers.inequalities import reduce_inequalities >>> reduce_inequalities(0 <= x + 3, []) (-3 <= x) & (x < oo) >>> reduce_inequalities(0 <= x + y*2 - 1, [x]) (x < oo) & (x >= 1 - 2*y) """ if not iterable(inequalities): inequalities = [inequalities] inequalities = [sympify(i) for i in inequalities] gens = set().union(*[i.free_symbols for i in inequalities]) if not iterable(symbols): symbols = [symbols] symbols = (set(symbols) or gens) & gens if any(i.is_real is False for i in symbols): raise TypeError(filldedent(''' inequalities cannot contain symbols that are not real. ''')) # make vanilla symbol real recast = {i: Dummy(i.name, real=True) for i in gens if i.is_real is None} inequalities = [i.xreplace(recast) for i in inequalities] symbols = {i.xreplace(recast) for i in symbols} # prefilter keep = [] for i in inequalities: if isinstance(i, Relational): i = i.func(i.lhs.as_expr() - i.rhs.as_expr(), 0) elif i not in (True, False): i = Eq(i, 0) if i == True: continue elif i == False: return S.false if i.lhs.is_number: raise NotImplementedError( "could not determine truth value of %s" % i) keep.append(i) inequalities = keep del keep # solve system rv = _reduce_inequalities(inequalities, symbols) # restore original symbols and return return rv.xreplace({v: k for k, v in recast.items()})