123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999 |
- """Tools for solving inequalities and systems of inequalities. """
- from sympy.calculus.util import (continuous_domain, periodicity,
- function_range)
- from sympy.core import Symbol, Dummy, sympify
- from sympy.core.exprtools import factor_terms
- from sympy.core.relational import Relational, Eq, Ge, Lt
- from sympy.sets.sets import Interval, FiniteSet, Union, Intersection
- from sympy.core.singleton import S
- from sympy.core.function import expand_mul
- from sympy.functions.elementary.complexes import im, Abs
- from sympy.logic import And
- from sympy.polys import Poly, PolynomialError, parallel_poly_from_expr
- from sympy.polys.polyutils import _nsort
- from sympy.solvers.solveset import solvify, solveset
- from sympy.utilities.iterables import sift, iterable
- from sympy.utilities.misc import filldedent
- def solve_poly_inequality(poly, rel):
- """Solve a polynomial inequality with rational coefficients.
- Examples
- ========
- >>> from sympy import solve_poly_inequality, Poly
- >>> from sympy.abc import x
- >>> 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.as_expr().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 import Poly
- >>> from sympy.solvers.inequalities import solve_poly_inequalities
- >>> 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))
- """
- return Union(*[s for p in polys for s in solve_poly_inequality(*p)])
- def solve_rational_inequalities(eqs):
- """Solve a system of rational inequalities with rational coefficients.
- Examples
- ========
- >>> from sympy.abc import x
- >>> from sympy import solve_rational_inequalities, Poly
- >>> 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
- def reduce_rational_inequalities(exprs, gen, relational=True):
- """Reduce a system of rational inequalities with rational coefficients.
- Examples
- ========
- >>> from sympy import 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
- >>> reduce_rational_inequalities([[(x + 2, ">")]], x)
- -2 < x
- >>> reduce_rational_inequalities([[x + 2]], x)
- Eq(x, -2)
- This function find the non-infinite solution set so if the unknown symbol
- is declared as extended real rather than real then the result may include
- finiteness conditions:
- >>> y = Symbol('y', extended_real=True)
- >>> reduce_rational_inequalities([[y + 2 > 0]], y)
- (-2 < y) & (y < oo)
- """
- 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
- def reduce_abs_inequality(expr, rel, gen):
- """Reduce an inequality with nested absolute values.
- Examples
- ========
- >>> from sympy import reduce_abs_inequality, Abs, Symbol
- >>> 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_extended_real is False:
- raise TypeError(filldedent('''
- Cannot 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)
- def reduce_abs_inequalities(exprs, gen):
- """Reduce a system of inequalities with nested absolute values.
- Examples
- ========
- >>> from sympy import reduce_abs_inequalities, Abs, Symbol
- >>> x = Symbol('x', extended_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 ])
- 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 :func:`sympy.solvers.solveset.solvify`.
- Notes
- =====
- Currently, we cannot solve all the inequalities due to limitations in
- :func:`sympy.solvers.solveset.solvify`. Also, the solution returned for trigonometric inequalities
- are restricted in its periodic interval.
- See Also
- ========
- sympy.solvers.solveset.solvify: solver returning solveset solutions with solve's output API
- Examples
- ========
- >>> from sympy import solve_univariate_inequality, Symbol, sin, Interval, S
- >>> x = Symbol('x')
- >>> solve_univariate_inequality(x**2 >= 4, x)
- ((2 <= x) & (x < oo)) | ((-oo < x) & (x <= -2))
- >>> 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.solvers.solvers import denoms
- if domain.is_subset(S.Reals) is False:
- raise NotImplementedError(filldedent('''
- Inequalities in the complex domain are
- not supported. Try the real domain by
- setting domain=S.Reals'''))
- elif domain is not S.Reals:
- rv = solve_univariate_inequality(
- expr, gen, relational=False, continuous=continuous).intersection(domain)
- if relational:
- rv = rv.as_relational(gen)
- return rv
- else:
- pass # continue with attempt to solve in Real domain
- # 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_extended_real is False:
- rv = S.EmptySet
- return rv if not relational else rv.as_relational(_gen)
- elif gen.is_extended_real is None:
- gen = Dummy('gen', extended_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 == 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 in ('<', '<='):
- if expr.func(frange.sup, 0):
- rv = domain
- elif not expr.func(frange.inf, 0):
- rv = S.EmptySet
- elif rel in ('>', '>='):
- 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).intersect(_domain)
- _domain = domain
- 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_extended_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_extended_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_extended_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_extended_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 im_sol is S.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)
- sol_sets = [S.EmptySet]
- start = domain.inf
- if start in domain and 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 end in domain and 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_extended_positive is None or
- end.is_infinite and end.is_extended_positive is None):
- raise ValueError('cannot proceed with unsigned infinite values')
- if (end.is_infinite and end.is_extended_negative or
- start.is_infinite and start.is_extended_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_extended_positive:
- pt = start*2
- elif start.is_extended_negative:
- pt = start*S.Half
- else:
- pt = start + 1
- elif start.is_infinite:
- if end.is_extended_positive:
- pt = end*S.Half
- elif end.is_extended_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 do not 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))
- def reduce_inequalities(inequalities, symbols=[]):
- """Reduce a system of inequalities with rational coefficients.
- Examples
- ========
- >>> from sympy.abc import x, y
- >>> from sympy 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_extended_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, extended_real=True)
- for i in gens if i.is_extended_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()})
|