inequalities.py 33 KB

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