1"""Tools for solving inequalities and systems of inequalities. """ 2 3from sympy.core import Symbol, Dummy, sympify 4from sympy.core.compatibility import iterable 5from sympy.core.exprtools import factor_terms 6from sympy.core.relational import Relational, Eq, Ge, Lt 7from sympy.sets import Interval 8from sympy.sets.sets import FiniteSet, Union, EmptySet, Intersection 9from sympy.core.singleton import S 10from sympy.core.function import expand_mul 11 12from sympy.functions import Abs 13from sympy.logic import And 14from sympy.polys import Poly, PolynomialError, parallel_poly_from_expr 15from sympy.polys.polyutils import _nsort 16from sympy.utilities.iterables import sift 17from sympy.utilities.misc import filldedent 18 19 20def solve_poly_inequality(poly, rel): 21 """Solve a polynomial inequality with rational coefficients. 22 23 Examples 24 ======== 25 26 >>> from sympy import Poly 27 >>> from sympy.abc import x 28 >>> from sympy.solvers.inequalities import solve_poly_inequality 29 30 >>> solve_poly_inequality(Poly(x, x, domain='ZZ'), '==') 31 [{0}] 32 33 >>> solve_poly_inequality(Poly(x**2 - 1, x, domain='ZZ'), '!=') 34 [Interval.open(-oo, -1), Interval.open(-1, 1), Interval.open(1, oo)] 35 36 >>> solve_poly_inequality(Poly(x**2 - 1, x, domain='ZZ'), '==') 37 [{-1}, {1}] 38 39 See Also 40 ======== 41 solve_poly_inequalities 42 """ 43 if not isinstance(poly, Poly): 44 raise ValueError( 45 'For efficiency reasons, `poly` should be a Poly instance') 46 if poly.as_expr().is_number: 47 t = Relational(poly.as_expr(), 0, rel) 48 if t is S.true: 49 return [S.Reals] 50 elif t is S.false: 51 return [S.EmptySet] 52 else: 53 raise NotImplementedError( 54 "could not determine truth value of %s" % t) 55 56 reals, intervals = poly.real_roots(multiple=False), [] 57 58 if rel == '==': 59 for root, _ in reals: 60 interval = Interval(root, root) 61 intervals.append(interval) 62 elif rel == '!=': 63 left = S.NegativeInfinity 64 65 for right, _ in reals + [(S.Infinity, 1)]: 66 interval = Interval(left, right, True, True) 67 intervals.append(interval) 68 left = right 69 else: 70 if poly.LC() > 0: 71 sign = +1 72 else: 73 sign = -1 74 75 eq_sign, equal = None, False 76 77 if rel == '>': 78 eq_sign = +1 79 elif rel == '<': 80 eq_sign = -1 81 elif rel == '>=': 82 eq_sign, equal = +1, True 83 elif rel == '<=': 84 eq_sign, equal = -1, True 85 else: 86 raise ValueError("'%s' is not a valid relation" % rel) 87 88 right, right_open = S.Infinity, True 89 90 for left, multiplicity in reversed(reals): 91 if multiplicity % 2: 92 if sign == eq_sign: 93 intervals.insert( 94 0, Interval(left, right, not equal, right_open)) 95 96 sign, right, right_open = -sign, left, not equal 97 else: 98 if sign == eq_sign and not equal: 99 intervals.insert( 100 0, Interval(left, right, True, right_open)) 101 right, right_open = left, True 102 elif sign != eq_sign and equal: 103 intervals.insert(0, Interval(left, left)) 104 105 if sign == eq_sign: 106 intervals.insert( 107 0, Interval(S.NegativeInfinity, right, True, right_open)) 108 109 return intervals 110 111 112def solve_poly_inequalities(polys): 113 """Solve polynomial inequalities with rational coefficients. 114 115 Examples 116 ======== 117 118 >>> from sympy.solvers.inequalities import solve_poly_inequalities 119 >>> from sympy.polys import Poly 120 >>> from sympy.abc import x 121 >>> solve_poly_inequalities((( 122 ... Poly(x**2 - 3), ">"), ( 123 ... Poly(-x**2 + 1), ">"))) 124 Union(Interval.open(-oo, -sqrt(3)), Interval.open(-1, 1), Interval.open(sqrt(3), oo)) 125 """ 126 from sympy import Union 127 return Union(*[s for p in polys for s in solve_poly_inequality(*p)]) 128 129 130def solve_rational_inequalities(eqs): 131 """Solve a system of rational inequalities with rational coefficients. 132 133 Examples 134 ======== 135 136 >>> from sympy.abc import x 137 >>> from sympy import Poly 138 >>> from sympy.solvers.inequalities import solve_rational_inequalities 139 140 >>> solve_rational_inequalities([[ 141 ... ((Poly(-x + 1), Poly(1, x)), '>='), 142 ... ((Poly(-x + 1), Poly(1, x)), '<=')]]) 143 {1} 144 145 >>> solve_rational_inequalities([[ 146 ... ((Poly(x), Poly(1, x)), '!='), 147 ... ((Poly(-x + 1), Poly(1, x)), '>=')]]) 148 Union(Interval.open(-oo, 0), Interval.Lopen(0, 1)) 149 150 See Also 151 ======== 152 solve_poly_inequality 153 """ 154 result = S.EmptySet 155 156 for _eqs in eqs: 157 if not _eqs: 158 continue 159 160 global_intervals = [Interval(S.NegativeInfinity, S.Infinity)] 161 162 for (numer, denom), rel in _eqs: 163 numer_intervals = solve_poly_inequality(numer*denom, rel) 164 denom_intervals = solve_poly_inequality(denom, '==') 165 166 intervals = [] 167 168 for numer_interval in numer_intervals: 169 for global_interval in global_intervals: 170 interval = numer_interval.intersect(global_interval) 171 172 if interval is not S.EmptySet: 173 intervals.append(interval) 174 175 global_intervals = intervals 176 177 intervals = [] 178 179 for global_interval in global_intervals: 180 for denom_interval in denom_intervals: 181 global_interval -= denom_interval 182 183 if global_interval is not S.EmptySet: 184 intervals.append(global_interval) 185 186 global_intervals = intervals 187 188 if not global_intervals: 189 break 190 191 for interval in global_intervals: 192 result = result.union(interval) 193 194 return result 195 196 197def reduce_rational_inequalities(exprs, gen, relational=True): 198 """Reduce a system of rational inequalities with rational coefficients. 199 200 Examples 201 ======== 202 203 >>> from sympy import Symbol 204 >>> from sympy.solvers.inequalities import reduce_rational_inequalities 205 206 >>> x = Symbol('x', real=True) 207 208 >>> reduce_rational_inequalities([[x**2 <= 0]], x) 209 Eq(x, 0) 210 211 >>> reduce_rational_inequalities([[x + 2 > 0]], x) 212 -2 < x 213 >>> reduce_rational_inequalities([[(x + 2, ">")]], x) 214 -2 < x 215 >>> reduce_rational_inequalities([[x + 2]], x) 216 Eq(x, -2) 217 218 This function find the non-infinite solution set so if the unknown symbol 219 is declared as extended real rather than real then the result may include 220 finiteness conditions: 221 222 >>> y = Symbol('y', extended_real=True) 223 >>> reduce_rational_inequalities([[y + 2 > 0]], y) 224 (-2 < y) & (y < oo) 225 """ 226 exact = True 227 eqs = [] 228 solution = S.Reals if exprs else S.EmptySet 229 for _exprs in exprs: 230 _eqs = [] 231 232 for expr in _exprs: 233 if isinstance(expr, tuple): 234 expr, rel = expr 235 else: 236 if expr.is_Relational: 237 expr, rel = expr.lhs - expr.rhs, expr.rel_op 238 else: 239 expr, rel = expr, '==' 240 241 if expr is S.true: 242 numer, denom, rel = S.Zero, S.One, '==' 243 elif expr is S.false: 244 numer, denom, rel = S.One, S.One, '==' 245 else: 246 numer, denom = expr.together().as_numer_denom() 247 248 try: 249 (numer, denom), opt = parallel_poly_from_expr( 250 (numer, denom), gen) 251 except PolynomialError: 252 raise PolynomialError(filldedent(''' 253 only polynomials and rational functions are 254 supported in this context. 255 ''')) 256 257 if not opt.domain.is_Exact: 258 numer, denom, exact = numer.to_exact(), denom.to_exact(), False 259 260 domain = opt.domain.get_exact() 261 262 if not (domain.is_ZZ or domain.is_QQ): 263 expr = numer/denom 264 expr = Relational(expr, 0, rel) 265 solution &= solve_univariate_inequality(expr, gen, relational=False) 266 else: 267 _eqs.append(((numer, denom), rel)) 268 269 if _eqs: 270 eqs.append(_eqs) 271 272 if eqs: 273 solution &= solve_rational_inequalities(eqs) 274 exclude = solve_rational_inequalities([[((d, d.one), '==') 275 for i in eqs for ((n, d), _) in i if d.has(gen)]]) 276 solution -= exclude 277 278 if not exact and solution: 279 solution = solution.evalf() 280 281 if relational: 282 solution = solution.as_relational(gen) 283 284 return solution 285 286 287def reduce_abs_inequality(expr, rel, gen): 288 """Reduce an inequality with nested absolute values. 289 290 Examples 291 ======== 292 293 >>> from sympy import Abs, Symbol 294 >>> from sympy.solvers.inequalities import reduce_abs_inequality 295 >>> x = Symbol('x', real=True) 296 297 >>> reduce_abs_inequality(Abs(x - 5) - 3, '<', x) 298 (2 < x) & (x < 8) 299 300 >>> reduce_abs_inequality(Abs(x + 2)*3 - 13, '<', x) 301 (-19/3 < x) & (x < 7/3) 302 303 See Also 304 ======== 305 306 reduce_abs_inequalities 307 """ 308 if gen.is_extended_real is False: 309 raise TypeError(filldedent(''' 310 can't solve inequalities with absolute values containing 311 non-real variables. 312 ''')) 313 314 def _bottom_up_scan(expr): 315 exprs = [] 316 317 if expr.is_Add or expr.is_Mul: 318 op = expr.func 319 320 for arg in expr.args: 321 _exprs = _bottom_up_scan(arg) 322 323 if not exprs: 324 exprs = _exprs 325 else: 326 args = [] 327 328 for expr, conds in exprs: 329 for _expr, _conds in _exprs: 330 args.append((op(expr, _expr), conds + _conds)) 331 332 exprs = args 333 elif expr.is_Pow: 334 n = expr.exp 335 if not n.is_Integer: 336 raise ValueError("Only Integer Powers are allowed on Abs.") 337 338 _exprs = _bottom_up_scan(expr.base) 339 340 for expr, conds in _exprs: 341 exprs.append((expr**n, conds)) 342 elif isinstance(expr, Abs): 343 _exprs = _bottom_up_scan(expr.args[0]) 344 345 for expr, conds in _exprs: 346 exprs.append(( expr, conds + [Ge(expr, 0)])) 347 exprs.append((-expr, conds + [Lt(expr, 0)])) 348 else: 349 exprs = [(expr, [])] 350 351 return exprs 352 353 exprs = _bottom_up_scan(expr) 354 355 mapping = {'<': '>', '<=': '>='} 356 inequalities = [] 357 358 for expr, conds in exprs: 359 if rel not in mapping.keys(): 360 expr = Relational( expr, 0, rel) 361 else: 362 expr = Relational(-expr, 0, mapping[rel]) 363 364 inequalities.append([expr] + conds) 365 366 return reduce_rational_inequalities(inequalities, gen) 367 368 369def reduce_abs_inequalities(exprs, gen): 370 """Reduce a system of inequalities with nested absolute values. 371 372 Examples 373 ======== 374 375 >>> from sympy import Abs, Symbol 376 >>> from sympy.solvers.inequalities import reduce_abs_inequalities 377 >>> x = Symbol('x', extended_real=True) 378 379 >>> reduce_abs_inequalities([(Abs(3*x - 5) - 7, '<'), 380 ... (Abs(x + 25) - 13, '>')], x) 381 (-2/3 < x) & (x < 4) & (((-oo < x) & (x < -38)) | ((-12 < x) & (x < oo))) 382 383 >>> reduce_abs_inequalities([(Abs(x - 4) + Abs(3*x - 5) - 7, '<')], x) 384 (1/2 < x) & (x < 4) 385 386 See Also 387 ======== 388 389 reduce_abs_inequality 390 """ 391 return And(*[ reduce_abs_inequality(expr, rel, gen) 392 for expr, rel in exprs ]) 393 394 395def solve_univariate_inequality(expr, gen, relational=True, domain=S.Reals, continuous=False): 396 """Solves a real univariate inequality. 397 398 Parameters 399 ========== 400 401 expr : Relational 402 The target inequality 403 gen : Symbol 404 The variable for which the inequality is solved 405 relational : bool 406 A Relational type output is expected or not 407 domain : Set 408 The domain over which the equation is solved 409 continuous: bool 410 True if expr is known to be continuous over the given domain 411 (and so continuous_domain() doesn't need to be called on it) 412 413 Raises 414 ====== 415 416 NotImplementedError 417 The solution of the inequality cannot be determined due to limitation 418 in :func:`sympy.solvers.solveset.solvify`. 419 420 Notes 421 ===== 422 423 Currently, we cannot solve all the inequalities due to limitations in 424 :func:`sympy.solvers.solveset.solvify`. Also, the solution returned for trigonometric inequalities 425 are restricted in its periodic interval. 426 427 See Also 428 ======== 429 430 sympy.solvers.solveset.solvify: solver returning solveset solutions with solve's output API 431 432 Examples 433 ======== 434 435 >>> from sympy.solvers.inequalities import solve_univariate_inequality 436 >>> from sympy import Symbol, sin, Interval, S 437 >>> x = Symbol('x') 438 439 >>> solve_univariate_inequality(x**2 >= 4, x) 440 ((2 <= x) & (x < oo)) | ((x <= -2) & (-oo < x)) 441 442 >>> solve_univariate_inequality(x**2 >= 4, x, relational=False) 443 Union(Interval(-oo, -2), Interval(2, oo)) 444 445 >>> domain = Interval(0, S.Infinity) 446 >>> solve_univariate_inequality(x**2 >= 4, x, False, domain) 447 Interval(2, oo) 448 449 >>> solve_univariate_inequality(sin(x) > 0, x, relational=False) 450 Interval.open(0, pi) 451 452 """ 453 from sympy import im 454 from sympy.calculus.util import (continuous_domain, periodicity, 455 function_range) 456 from sympy.solvers.solvers import denoms 457 from sympy.solvers.solveset import solvify, solveset 458 459 if domain.is_subset(S.Reals) is False: 460 raise NotImplementedError(filldedent(''' 461 Inequalities in the complex domain are 462 not supported. Try the real domain by 463 setting domain=S.Reals''')) 464 elif domain is not S.Reals: 465 rv = solve_univariate_inequality( 466 expr, gen, relational=False, continuous=continuous).intersection(domain) 467 if relational: 468 rv = rv.as_relational(gen) 469 return rv 470 else: 471 pass # continue with attempt to solve in Real domain 472 473 # This keeps the function independent of the assumptions about `gen`. 474 # `solveset` makes sure this function is called only when the domain is 475 # real. 476 _gen = gen 477 _domain = domain 478 if gen.is_extended_real is False: 479 rv = S.EmptySet 480 return rv if not relational else rv.as_relational(_gen) 481 elif gen.is_extended_real is None: 482 gen = Dummy('gen', extended_real=True) 483 try: 484 expr = expr.xreplace({_gen: gen}) 485 except TypeError: 486 raise TypeError(filldedent(''' 487 When gen is real, the relational has a complex part 488 which leads to an invalid comparison like I < 0. 489 ''')) 490 491 rv = None 492 493 if expr is S.true: 494 rv = domain 495 496 elif expr is S.false: 497 rv = S.EmptySet 498 499 else: 500 e = expr.lhs - expr.rhs 501 period = periodicity(e, gen) 502 if period == S.Zero: 503 e = expand_mul(e) 504 const = expr.func(e, 0) 505 if const is S.true: 506 rv = domain 507 elif const is S.false: 508 rv = S.EmptySet 509 elif period is not None: 510 frange = function_range(e, gen, domain) 511 512 rel = expr.rel_op 513 if rel == '<' or rel == '<=': 514 if expr.func(frange.sup, 0): 515 rv = domain 516 elif not expr.func(frange.inf, 0): 517 rv = S.EmptySet 518 519 elif rel == '>' or rel == '>=': 520 if expr.func(frange.inf, 0): 521 rv = domain 522 elif not expr.func(frange.sup, 0): 523 rv = S.EmptySet 524 525 inf, sup = domain.inf, domain.sup 526 if sup - inf is S.Infinity: 527 domain = Interval(0, period, False, True).intersect(_domain) 528 _domain = domain 529 530 if rv is None: 531 n, d = e.as_numer_denom() 532 try: 533 if gen not in n.free_symbols and len(e.free_symbols) > 1: 534 raise ValueError 535 # this might raise ValueError on its own 536 # or it might give None... 537 solns = solvify(e, gen, domain) 538 if solns is None: 539 # in which case we raise ValueError 540 raise ValueError 541 except (ValueError, NotImplementedError): 542 # replace gen with generic x since it's 543 # univariate anyway 544 raise NotImplementedError(filldedent(''' 545 The inequality, %s, cannot be solved using 546 solve_univariate_inequality. 547 ''' % expr.subs(gen, Symbol('x')))) 548 549 expanded_e = expand_mul(e) 550 def valid(x): 551 # this is used to see if gen=x satisfies the 552 # relational by substituting it into the 553 # expanded form and testing against 0, e.g. 554 # if expr = x*(x + 1) < 2 then e = x*(x + 1) - 2 555 # and expanded_e = x**2 + x - 2; the test is 556 # whether a given value of x satisfies 557 # x**2 + x - 2 < 0 558 # 559 # expanded_e, expr and gen used from enclosing scope 560 v = expanded_e.subs(gen, expand_mul(x)) 561 try: 562 r = expr.func(v, 0) 563 except TypeError: 564 r = S.false 565 if r in (S.true, S.false): 566 return r 567 if v.is_extended_real is False: 568 return S.false 569 else: 570 v = v.n(2) 571 if v.is_comparable: 572 return expr.func(v, 0) 573 # not comparable or couldn't be evaluated 574 raise NotImplementedError( 575 'relationship did not evaluate: %s' % r) 576 577 singularities = [] 578 for d in denoms(expr, gen): 579 singularities.extend(solvify(d, gen, domain)) 580 if not continuous: 581 domain = continuous_domain(expanded_e, gen, domain) 582 583 include_x = '=' in expr.rel_op and expr.rel_op != '!=' 584 585 try: 586 discontinuities = set(domain.boundary - 587 FiniteSet(domain.inf, domain.sup)) 588 # remove points that are not between inf and sup of domain 589 critical_points = FiniteSet(*(solns + singularities + list( 590 discontinuities))).intersection( 591 Interval(domain.inf, domain.sup, 592 domain.inf not in domain, domain.sup not in domain)) 593 if all(r.is_number for r in critical_points): 594 reals = _nsort(critical_points, separated=True)[0] 595 else: 596 sifted = sift(critical_points, lambda x: x.is_extended_real) 597 if sifted[None]: 598 # there were some roots that weren't known 599 # to be real 600 raise NotImplementedError 601 try: 602 reals = sifted[True] 603 if len(reals) > 1: 604 reals = list(sorted(reals)) 605 except TypeError: 606 raise NotImplementedError 607 except NotImplementedError: 608 raise NotImplementedError('sorting of these roots is not supported') 609 610 # If expr contains imaginary coefficients, only take real 611 # values of x for which the imaginary part is 0 612 make_real = S.Reals 613 if im(expanded_e) != S.Zero: 614 check = True 615 im_sol = FiniteSet() 616 try: 617 a = solveset(im(expanded_e), gen, domain) 618 if not isinstance(a, Interval): 619 for z in a: 620 if z not in singularities and valid(z) and z.is_extended_real: 621 im_sol += FiniteSet(z) 622 else: 623 start, end = a.inf, a.sup 624 for z in _nsort(critical_points + FiniteSet(end)): 625 valid_start = valid(start) 626 if start != end: 627 valid_z = valid(z) 628 pt = _pt(start, z) 629 if pt not in singularities and pt.is_extended_real and valid(pt): 630 if valid_start and valid_z: 631 im_sol += Interval(start, z) 632 elif valid_start: 633 im_sol += Interval.Ropen(start, z) 634 elif valid_z: 635 im_sol += Interval.Lopen(start, z) 636 else: 637 im_sol += Interval.open(start, z) 638 start = z 639 for s in singularities: 640 im_sol -= FiniteSet(s) 641 except (TypeError): 642 im_sol = S.Reals 643 check = False 644 645 if isinstance(im_sol, EmptySet): 646 raise ValueError(filldedent(''' 647 %s contains imaginary parts which cannot be 648 made 0 for any value of %s satisfying the 649 inequality, leading to relations like I < 0. 650 ''' % (expr.subs(gen, _gen), _gen))) 651 652 make_real = make_real.intersect(im_sol) 653 654 sol_sets = [S.EmptySet] 655 656 start = domain.inf 657 if start in domain and valid(start) and start.is_finite: 658 sol_sets.append(FiniteSet(start)) 659 660 for x in reals: 661 end = x 662 663 if valid(_pt(start, end)): 664 sol_sets.append(Interval(start, end, True, True)) 665 666 if x in singularities: 667 singularities.remove(x) 668 else: 669 if x in discontinuities: 670 discontinuities.remove(x) 671 _valid = valid(x) 672 else: # it's a solution 673 _valid = include_x 674 if _valid: 675 sol_sets.append(FiniteSet(x)) 676 677 start = end 678 679 end = domain.sup 680 if end in domain and valid(end) and end.is_finite: 681 sol_sets.append(FiniteSet(end)) 682 683 if valid(_pt(start, end)): 684 sol_sets.append(Interval.open(start, end)) 685 686 if im(expanded_e) != S.Zero and check: 687 rv = (make_real).intersect(_domain) 688 else: 689 rv = Intersection( 690 (Union(*sol_sets)), make_real, _domain).subs(gen, _gen) 691 692 return rv if not relational else rv.as_relational(_gen) 693 694 695def _pt(start, end): 696 """Return a point between start and end""" 697 if not start.is_infinite and not end.is_infinite: 698 pt = (start + end)/2 699 elif start.is_infinite and end.is_infinite: 700 pt = S.Zero 701 else: 702 if (start.is_infinite and start.is_extended_positive is None or 703 end.is_infinite and end.is_extended_positive is None): 704 raise ValueError('cannot proceed with unsigned infinite values') 705 if (end.is_infinite and end.is_extended_negative or 706 start.is_infinite and start.is_extended_positive): 707 start, end = end, start 708 # if possible, use a multiple of self which has 709 # better behavior when checking assumptions than 710 # an expression obtained by adding or subtracting 1 711 if end.is_infinite: 712 if start.is_extended_positive: 713 pt = start*2 714 elif start.is_extended_negative: 715 pt = start*S.Half 716 else: 717 pt = start + 1 718 elif start.is_infinite: 719 if end.is_extended_positive: 720 pt = end*S.Half 721 elif end.is_extended_negative: 722 pt = end*2 723 else: 724 pt = end - 1 725 return pt 726 727 728def _solve_inequality(ie, s, linear=False): 729 """Return the inequality with s isolated on the left, if possible. 730 If the relationship is non-linear, a solution involving And or Or 731 may be returned. False or True are returned if the relationship 732 is never True or always True, respectively. 733 734 If `linear` is True (default is False) an `s`-dependent expression 735 will be isolated on the left, if possible 736 but it will not be solved for `s` unless the expression is linear 737 in `s`. Furthermore, only "safe" operations which don't change the 738 sense of the relationship are applied: no division by an unsigned 739 value is attempted unless the relationship involves Eq or Ne and 740 no division by a value not known to be nonzero is ever attempted. 741 742 Examples 743 ======== 744 745 >>> from sympy import Eq, Symbol 746 >>> from sympy.solvers.inequalities import _solve_inequality as f 747 >>> from sympy.abc import x, y 748 749 For linear expressions, the symbol can be isolated: 750 751 >>> f(x - 2 < 0, x) 752 x < 2 753 >>> f(-x - 6 < x, x) 754 x > -3 755 756 Sometimes nonlinear relationships will be False 757 758 >>> f(x**2 + 4 < 0, x) 759 False 760 761 Or they may involve more than one region of values: 762 763 >>> f(x**2 - 4 < 0, x) 764 (-2 < x) & (x < 2) 765 766 To restrict the solution to a relational, set linear=True 767 and only the x-dependent portion will be isolated on the left: 768 769 >>> f(x**2 - 4 < 0, x, linear=True) 770 x**2 < 4 771 772 Division of only nonzero quantities is allowed, so x cannot 773 be isolated by dividing by y: 774 775 >>> y.is_nonzero is None # it is unknown whether it is 0 or not 776 True 777 >>> f(x*y < 1, x) 778 x*y < 1 779 780 And while an equality (or inequality) still holds after dividing by a 781 non-zero quantity 782 783 >>> nz = Symbol('nz', nonzero=True) 784 >>> f(Eq(x*nz, 1), x) 785 Eq(x, 1/nz) 786 787 the sign must be known for other inequalities involving > or <: 788 789 >>> f(x*nz <= 1, x) 790 nz*x <= 1 791 >>> p = Symbol('p', positive=True) 792 >>> f(x*p <= 1, x) 793 x <= 1/p 794 795 When there are denominators in the original expression that 796 are removed by expansion, conditions for them will be returned 797 as part of the result: 798 799 >>> f(x < x*(2/x - 1), x) 800 (x < 1) & Ne(x, 0) 801 """ 802 from sympy.solvers.solvers import denoms 803 if s not in ie.free_symbols: 804 return ie 805 if ie.rhs == s: 806 ie = ie.reversed 807 if ie.lhs == s and s not in ie.rhs.free_symbols: 808 return ie 809 810 def classify(ie, s, i): 811 # return True or False if ie evaluates when substituting s with 812 # i else None (if unevaluated) or NaN (when there is an error 813 # in evaluating) 814 try: 815 v = ie.subs(s, i) 816 if v is S.NaN: 817 return v 818 elif v not in (True, False): 819 return 820 return v 821 except TypeError: 822 return S.NaN 823 824 rv = None 825 oo = S.Infinity 826 expr = ie.lhs - ie.rhs 827 try: 828 p = Poly(expr, s) 829 if p.degree() == 0: 830 rv = ie.func(p.as_expr(), 0) 831 elif not linear and p.degree() > 1: 832 # handle in except clause 833 raise NotImplementedError 834 except (PolynomialError, NotImplementedError): 835 if not linear: 836 try: 837 rv = reduce_rational_inequalities([[ie]], s) 838 except PolynomialError: 839 rv = solve_univariate_inequality(ie, s) 840 # remove restrictions wrt +/-oo that may have been 841 # applied when using sets to simplify the relationship 842 okoo = classify(ie, s, oo) 843 if okoo is S.true and classify(rv, s, oo) is S.false: 844 rv = rv.subs(s < oo, True) 845 oknoo = classify(ie, s, -oo) 846 if (oknoo is S.true and 847 classify(rv, s, -oo) is S.false): 848 rv = rv.subs(-oo < s, True) 849 rv = rv.subs(s > -oo, True) 850 if rv is S.true: 851 rv = (s <= oo) if okoo is S.true else (s < oo) 852 if oknoo is not S.true: 853 rv = And(-oo < s, rv) 854 else: 855 p = Poly(expr) 856 857 conds = [] 858 if rv is None: 859 e = p.as_expr() # this is in expanded form 860 # Do a safe inversion of e, moving non-s terms 861 # to the rhs and dividing by a nonzero factor if 862 # the relational is Eq/Ne; for other relationals 863 # the sign must also be positive or negative 864 rhs = 0 865 b, ax = e.as_independent(s, as_Add=True) 866 e -= b 867 rhs -= b 868 ef = factor_terms(e) 869 a, e = ef.as_independent(s, as_Add=False) 870 if (a.is_zero != False or # don't divide by potential 0 871 a.is_negative == 872 a.is_positive is None and # if sign is not known then 873 ie.rel_op not in ('!=', '==')): # reject if not Eq/Ne 874 e = ef 875 a = S.One 876 rhs /= a 877 if a.is_positive: 878 rv = ie.func(e, rhs) 879 else: 880 rv = ie.reversed.func(e, rhs) 881 882 # return conditions under which the value is 883 # valid, too. 884 beginning_denoms = denoms(ie.lhs) | denoms(ie.rhs) 885 current_denoms = denoms(rv) 886 for d in beginning_denoms - current_denoms: 887 c = _solve_inequality(Eq(d, 0), s, linear=linear) 888 if isinstance(c, Eq) and c.lhs == s: 889 if classify(rv, s, c.rhs) is S.true: 890 # rv is permitting this value but it shouldn't 891 conds.append(~c) 892 for i in (-oo, oo): 893 if (classify(rv, s, i) is S.true and 894 classify(ie, s, i) is not S.true): 895 conds.append(s < i if i is oo else i < s) 896 897 conds.append(rv) 898 return And(*conds) 899 900 901def _reduce_inequalities(inequalities, symbols): 902 # helper for reduce_inequalities 903 904 poly_part, abs_part = {}, {} 905 other = [] 906 907 for inequality in inequalities: 908 909 expr, rel = inequality.lhs, inequality.rel_op # rhs is 0 910 911 # check for gens using atoms which is more strict than free_symbols to 912 # guard against EX domain which won't be handled by 913 # reduce_rational_inequalities 914 gens = expr.atoms(Symbol) 915 916 if len(gens) == 1: 917 gen = gens.pop() 918 else: 919 common = expr.free_symbols & symbols 920 if len(common) == 1: 921 gen = common.pop() 922 other.append(_solve_inequality(Relational(expr, 0, rel), gen)) 923 continue 924 else: 925 raise NotImplementedError(filldedent(''' 926 inequality has more than one symbol of interest. 927 ''')) 928 929 if expr.is_polynomial(gen): 930 poly_part.setdefault(gen, []).append((expr, rel)) 931 else: 932 components = expr.find(lambda u: 933 u.has(gen) and ( 934 u.is_Function or u.is_Pow and not u.exp.is_Integer)) 935 if components and all(isinstance(i, Abs) for i in components): 936 abs_part.setdefault(gen, []).append((expr, rel)) 937 else: 938 other.append(_solve_inequality(Relational(expr, 0, rel), gen)) 939 940 poly_reduced = [] 941 abs_reduced = [] 942 943 for gen, exprs in poly_part.items(): 944 poly_reduced.append(reduce_rational_inequalities([exprs], gen)) 945 946 for gen, exprs in abs_part.items(): 947 abs_reduced.append(reduce_abs_inequalities(exprs, gen)) 948 949 return And(*(poly_reduced + abs_reduced + other)) 950 951 952def reduce_inequalities(inequalities, symbols=[]): 953 """Reduce a system of inequalities with rational coefficients. 954 955 Examples 956 ======== 957 958 >>> from sympy.abc import x, y 959 >>> from sympy.solvers.inequalities import reduce_inequalities 960 961 >>> reduce_inequalities(0 <= x + 3, []) 962 (-3 <= x) & (x < oo) 963 964 >>> reduce_inequalities(0 <= x + y*2 - 1, [x]) 965 (x < oo) & (x >= 1 - 2*y) 966 """ 967 if not iterable(inequalities): 968 inequalities = [inequalities] 969 inequalities = [sympify(i) for i in inequalities] 970 971 gens = set().union(*[i.free_symbols for i in inequalities]) 972 973 if not iterable(symbols): 974 symbols = [symbols] 975 symbols = (set(symbols) or gens) & gens 976 if any(i.is_extended_real is False for i in symbols): 977 raise TypeError(filldedent(''' 978 inequalities cannot contain symbols that are not real. 979 ''')) 980 981 # make vanilla symbol real 982 recast = {i: Dummy(i.name, extended_real=True) 983 for i in gens if i.is_extended_real is None} 984 inequalities = [i.xreplace(recast) for i in inequalities] 985 symbols = {i.xreplace(recast) for i in symbols} 986 987 # prefilter 988 keep = [] 989 for i in inequalities: 990 if isinstance(i, Relational): 991 i = i.func(i.lhs.as_expr() - i.rhs.as_expr(), 0) 992 elif i not in (True, False): 993 i = Eq(i, 0) 994 if i == True: 995 continue 996 elif i == False: 997 return S.false 998 if i.lhs.is_number: 999 raise NotImplementedError( 1000 "could not determine truth value of %s" % i) 1001 keep.append(i) 1002 inequalities = keep 1003 del keep 1004 1005 # solve system 1006 rv = _reduce_inequalities(inequalities, symbols) 1007 1008 # restore original symbols and return 1009 return rv.xreplace({v: k for k, v in recast.items()}) 1010