Coverage for /usr/lib/python3/dist-packages/sympy/solvers/inequalities.py: 6%

472 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

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

2import itertools 

3 

4from sympy.calculus.util import (continuous_domain, periodicity, 

5 function_range) 

6from sympy.core import Symbol, Dummy, sympify 

7from sympy.core.exprtools import factor_terms 

8from sympy.core.relational import Relational, Eq, Ge, Lt 

9from sympy.sets.sets import Interval, FiniteSet, Union, Intersection 

10from sympy.core.singleton import S 

11from sympy.core.function import expand_mul 

12from sympy.functions.elementary.complexes import im, Abs 

13from sympy.logic import And 

14from sympy.polys import Poly, PolynomialError, parallel_poly_from_expr 

15from sympy.polys.polyutils import _nsort 

16from sympy.solvers.solveset import solvify, solveset 

17from sympy.utilities.iterables import sift, iterable 

18from sympy.utilities.misc import filldedent 

19 

20 

21def solve_poly_inequality(poly, rel): 

22 """Solve a polynomial inequality with rational coefficients. 

23 

24 Examples 

25 ======== 

26 

27 >>> from sympy import solve_poly_inequality, Poly 

28 >>> from sympy.abc import x 

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 import Poly 

119 >>> from sympy.solvers.inequalities import solve_poly_inequalities 

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 return Union(*[s for p in polys for s in solve_poly_inequality(*p)]) 

127 

128 

129def solve_rational_inequalities(eqs): 

130 """Solve a system of rational inequalities with rational coefficients. 

131 

132 Examples 

133 ======== 

134 

135 >>> from sympy.abc import x 

136 >>> from sympy import solve_rational_inequalities, Poly 

137 

138 >>> solve_rational_inequalities([[ 

139 ... ((Poly(-x + 1), Poly(1, x)), '>='), 

140 ... ((Poly(-x + 1), Poly(1, x)), '<=')]]) 

141 {1} 

142 

143 >>> solve_rational_inequalities([[ 

144 ... ((Poly(x), Poly(1, x)), '!='), 

145 ... ((Poly(-x + 1), Poly(1, x)), '>=')]]) 

146 Union(Interval.open(-oo, 0), Interval.Lopen(0, 1)) 

147 

148 See Also 

149 ======== 

150 solve_poly_inequality 

151 """ 

152 result = S.EmptySet 

153 

154 for _eqs in eqs: 

155 if not _eqs: 

156 continue 

157 

158 global_intervals = [Interval(S.NegativeInfinity, S.Infinity)] 

159 

160 for (numer, denom), rel in _eqs: 

161 numer_intervals = solve_poly_inequality(numer*denom, rel) 

162 denom_intervals = solve_poly_inequality(denom, '==') 

163 

164 intervals = [] 

165 

166 for numer_interval, global_interval in itertools.product( 

167 numer_intervals, global_intervals): 

168 interval = numer_interval.intersect(global_interval) 

169 

170 if interval is not S.EmptySet: 

171 intervals.append(interval) 

172 

173 global_intervals = intervals 

174 

175 intervals = [] 

176 

177 for global_interval in global_intervals: 

178 for denom_interval in denom_intervals: 

179 global_interval -= denom_interval 

180 

181 if global_interval is not S.EmptySet: 

182 intervals.append(global_interval) 

183 

184 global_intervals = intervals 

185 

186 if not global_intervals: 

187 break 

188 

189 for interval in global_intervals: 

190 result = result.union(interval) 

191 

192 return result 

193 

194 

195def reduce_rational_inequalities(exprs, gen, relational=True): 

196 """Reduce a system of rational inequalities with rational coefficients. 

197 

198 Examples 

199 ======== 

200 

201 >>> from sympy import Symbol 

202 >>> from sympy.solvers.inequalities import reduce_rational_inequalities 

203 

204 >>> x = Symbol('x', real=True) 

205 

206 >>> reduce_rational_inequalities([[x**2 <= 0]], x) 

207 Eq(x, 0) 

208 

209 >>> reduce_rational_inequalities([[x + 2 > 0]], x) 

210 -2 < x 

211 >>> reduce_rational_inequalities([[(x + 2, ">")]], x) 

212 -2 < x 

213 >>> reduce_rational_inequalities([[x + 2]], x) 

214 Eq(x, -2) 

215 

216 This function find the non-infinite solution set so if the unknown symbol 

217 is declared as extended real rather than real then the result may include 

218 finiteness conditions: 

219 

220 >>> y = Symbol('y', extended_real=True) 

221 >>> reduce_rational_inequalities([[y + 2 > 0]], y) 

222 (-2 < y) & (y < oo) 

223 """ 

224 exact = True 

225 eqs = [] 

226 solution = S.Reals if exprs else S.EmptySet 

227 for _exprs in exprs: 

228 _eqs = [] 

229 

230 for expr in _exprs: 

231 if isinstance(expr, tuple): 

232 expr, rel = expr 

233 else: 

234 if expr.is_Relational: 

235 expr, rel = expr.lhs - expr.rhs, expr.rel_op 

236 else: 

237 expr, rel = expr, '==' 

238 

239 if expr is S.true: 

240 numer, denom, rel = S.Zero, S.One, '==' 

241 elif expr is S.false: 

242 numer, denom, rel = S.One, S.One, '==' 

243 else: 

244 numer, denom = expr.together().as_numer_denom() 

245 

246 try: 

247 (numer, denom), opt = parallel_poly_from_expr( 

248 (numer, denom), gen) 

249 except PolynomialError: 

250 raise PolynomialError(filldedent(''' 

251 only polynomials and rational functions are 

252 supported in this context. 

253 ''')) 

254 

255 if not opt.domain.is_Exact: 

256 numer, denom, exact = numer.to_exact(), denom.to_exact(), False 

257 

258 domain = opt.domain.get_exact() 

259 

260 if not (domain.is_ZZ or domain.is_QQ): 

261 expr = numer/denom 

262 expr = Relational(expr, 0, rel) 

263 solution &= solve_univariate_inequality(expr, gen, relational=False) 

264 else: 

265 _eqs.append(((numer, denom), rel)) 

266 

267 if _eqs: 

268 eqs.append(_eqs) 

269 

270 if eqs: 

271 solution &= solve_rational_inequalities(eqs) 

272 exclude = solve_rational_inequalities([[((d, d.one), '==') 

273 for i in eqs for ((n, d), _) in i if d.has(gen)]]) 

274 solution -= exclude 

275 

276 if not exact and solution: 

277 solution = solution.evalf() 

278 

279 if relational: 

280 solution = solution.as_relational(gen) 

281 

282 return solution 

283 

284 

285def reduce_abs_inequality(expr, rel, gen): 

286 """Reduce an inequality with nested absolute values. 

287 

288 Examples 

289 ======== 

290 

291 >>> from sympy import reduce_abs_inequality, Abs, Symbol 

292 >>> x = Symbol('x', real=True) 

293 

294 >>> reduce_abs_inequality(Abs(x - 5) - 3, '<', x) 

295 (2 < x) & (x < 8) 

296 

297 >>> reduce_abs_inequality(Abs(x + 2)*3 - 13, '<', x) 

298 (-19/3 < x) & (x < 7/3) 

299 

300 See Also 

301 ======== 

302 

303 reduce_abs_inequalities 

304 """ 

305 if gen.is_extended_real is False: 

306 raise TypeError(filldedent(''' 

307 Cannot solve inequalities with absolute values containing 

308 non-real variables. 

309 ''')) 

310 

311 def _bottom_up_scan(expr): 

312 exprs = [] 

313 

314 if expr.is_Add or expr.is_Mul: 

315 op = expr.func 

316 

317 for arg in expr.args: 

318 _exprs = _bottom_up_scan(arg) 

319 

320 if not exprs: 

321 exprs = _exprs 

322 else: 

323 exprs = [(op(expr, _expr), conds + _conds) for (expr, conds), (_expr, _conds) in 

324 itertools.product(exprs, _exprs)] 

325 elif expr.is_Pow: 

326 n = expr.exp 

327 if not n.is_Integer: 

328 raise ValueError("Only Integer Powers are allowed on Abs.") 

329 

330 exprs.extend((expr**n, conds) for expr, conds in _bottom_up_scan(expr.base)) 

331 elif isinstance(expr, Abs): 

332 _exprs = _bottom_up_scan(expr.args[0]) 

333 

334 for expr, conds in _exprs: 

335 exprs.append(( expr, conds + [Ge(expr, 0)])) 

336 exprs.append((-expr, conds + [Lt(expr, 0)])) 

337 else: 

338 exprs = [(expr, [])] 

339 

340 return exprs 

341 

342 mapping = {'<': '>', '<=': '>='} 

343 inequalities = [] 

344 

345 for expr, conds in _bottom_up_scan(expr): 

346 if rel not in mapping.keys(): 

347 expr = Relational( expr, 0, rel) 

348 else: 

349 expr = Relational(-expr, 0, mapping[rel]) 

350 

351 inequalities.append([expr] + conds) 

352 

353 return reduce_rational_inequalities(inequalities, gen) 

354 

355 

356def reduce_abs_inequalities(exprs, gen): 

357 """Reduce a system of inequalities with nested absolute values. 

358 

359 Examples 

360 ======== 

361 

362 >>> from sympy import reduce_abs_inequalities, Abs, Symbol 

363 >>> x = Symbol('x', extended_real=True) 

364 

365 >>> reduce_abs_inequalities([(Abs(3*x - 5) - 7, '<'), 

366 ... (Abs(x + 25) - 13, '>')], x) 

367 (-2/3 < x) & (x < 4) & (((-oo < x) & (x < -38)) | ((-12 < x) & (x < oo))) 

368 

369 >>> reduce_abs_inequalities([(Abs(x - 4) + Abs(3*x - 5) - 7, '<')], x) 

370 (1/2 < x) & (x < 4) 

371 

372 See Also 

373 ======== 

374 

375 reduce_abs_inequality 

376 """ 

377 return And(*[ reduce_abs_inequality(expr, rel, gen) 

378 for expr, rel in exprs ]) 

379 

380 

381def solve_univariate_inequality(expr, gen, relational=True, domain=S.Reals, continuous=False): 

382 """Solves a real univariate inequality. 

383 

384 Parameters 

385 ========== 

386 

387 expr : Relational 

388 The target inequality 

389 gen : Symbol 

390 The variable for which the inequality is solved 

391 relational : bool 

392 A Relational type output is expected or not 

393 domain : Set 

394 The domain over which the equation is solved 

395 continuous: bool 

396 True if expr is known to be continuous over the given domain 

397 (and so continuous_domain() does not need to be called on it) 

398 

399 Raises 

400 ====== 

401 

402 NotImplementedError 

403 The solution of the inequality cannot be determined due to limitation 

404 in :func:`sympy.solvers.solveset.solvify`. 

405 

406 Notes 

407 ===== 

408 

409 Currently, we cannot solve all the inequalities due to limitations in 

410 :func:`sympy.solvers.solveset.solvify`. Also, the solution returned for trigonometric inequalities 

411 are restricted in its periodic interval. 

412 

413 See Also 

414 ======== 

415 

416 sympy.solvers.solveset.solvify: solver returning solveset solutions with solve's output API 

417 

418 Examples 

419 ======== 

420 

421 >>> from sympy import solve_univariate_inequality, Symbol, sin, Interval, S 

422 >>> x = Symbol('x') 

423 

424 >>> solve_univariate_inequality(x**2 >= 4, x) 

425 ((2 <= x) & (x < oo)) | ((-oo < x) & (x <= -2)) 

426 

427 >>> solve_univariate_inequality(x**2 >= 4, x, relational=False) 

428 Union(Interval(-oo, -2), Interval(2, oo)) 

429 

430 >>> domain = Interval(0, S.Infinity) 

431 >>> solve_univariate_inequality(x**2 >= 4, x, False, domain) 

432 Interval(2, oo) 

433 

434 >>> solve_univariate_inequality(sin(x) > 0, x, relational=False) 

435 Interval.open(0, pi) 

436 

437 """ 

438 from sympy.solvers.solvers import denoms 

439 

440 if domain.is_subset(S.Reals) is False: 

441 raise NotImplementedError(filldedent(''' 

442 Inequalities in the complex domain are 

443 not supported. Try the real domain by 

444 setting domain=S.Reals''')) 

445 elif domain is not S.Reals: 

446 rv = solve_univariate_inequality( 

447 expr, gen, relational=False, continuous=continuous).intersection(domain) 

448 if relational: 

449 rv = rv.as_relational(gen) 

450 return rv 

451 else: 

452 pass # continue with attempt to solve in Real domain 

453 

454 # This keeps the function independent of the assumptions about `gen`. 

455 # `solveset` makes sure this function is called only when the domain is 

456 # real. 

457 _gen = gen 

458 _domain = domain 

459 if gen.is_extended_real is False: 

460 rv = S.EmptySet 

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

462 elif gen.is_extended_real is None: 

463 gen = Dummy('gen', extended_real=True) 

464 try: 

465 expr = expr.xreplace({_gen: gen}) 

466 except TypeError: 

467 raise TypeError(filldedent(''' 

468 When gen is real, the relational has a complex part 

469 which leads to an invalid comparison like I < 0. 

470 ''')) 

471 

472 rv = None 

473 

474 if expr is S.true: 

475 rv = domain 

476 

477 elif expr is S.false: 

478 rv = S.EmptySet 

479 

480 else: 

481 e = expr.lhs - expr.rhs 

482 period = periodicity(e, gen) 

483 if period == S.Zero: 

484 e = expand_mul(e) 

485 const = expr.func(e, 0) 

486 if const is S.true: 

487 rv = domain 

488 elif const is S.false: 

489 rv = S.EmptySet 

490 elif period is not None: 

491 frange = function_range(e, gen, domain) 

492 

493 rel = expr.rel_op 

494 if rel in ('<', '<='): 

495 if expr.func(frange.sup, 0): 

496 rv = domain 

497 elif not expr.func(frange.inf, 0): 

498 rv = S.EmptySet 

499 

500 elif rel in ('>', '>='): 

501 if expr.func(frange.inf, 0): 

502 rv = domain 

503 elif not expr.func(frange.sup, 0): 

504 rv = S.EmptySet 

505 

506 inf, sup = domain.inf, domain.sup 

507 if sup - inf is S.Infinity: 

508 domain = Interval(0, period, False, True).intersect(_domain) 

509 _domain = domain 

510 

511 if rv is None: 

512 n, d = e.as_numer_denom() 

513 try: 

514 if gen not in n.free_symbols and len(e.free_symbols) > 1: 

515 raise ValueError 

516 # this might raise ValueError on its own 

517 # or it might give None... 

518 solns = solvify(e, gen, domain) 

519 if solns is None: 

520 # in which case we raise ValueError 

521 raise ValueError 

522 except (ValueError, NotImplementedError): 

523 # replace gen with generic x since it's 

524 # univariate anyway 

525 raise NotImplementedError(filldedent(''' 

526 The inequality, %s, cannot be solved using 

527 solve_univariate_inequality. 

528 ''' % expr.subs(gen, Symbol('x')))) 

529 

530 expanded_e = expand_mul(e) 

531 def valid(x): 

532 # this is used to see if gen=x satisfies the 

533 # relational by substituting it into the 

534 # expanded form and testing against 0, e.g. 

535 # if expr = x*(x + 1) < 2 then e = x*(x + 1) - 2 

536 # and expanded_e = x**2 + x - 2; the test is 

537 # whether a given value of x satisfies 

538 # x**2 + x - 2 < 0 

539 # 

540 # expanded_e, expr and gen used from enclosing scope 

541 v = expanded_e.subs(gen, expand_mul(x)) 

542 try: 

543 r = expr.func(v, 0) 

544 except TypeError: 

545 r = S.false 

546 if r in (S.true, S.false): 

547 return r 

548 if v.is_extended_real is False: 

549 return S.false 

550 else: 

551 v = v.n(2) 

552 if v.is_comparable: 

553 return expr.func(v, 0) 

554 # not comparable or couldn't be evaluated 

555 raise NotImplementedError( 

556 'relationship did not evaluate: %s' % r) 

557 

558 singularities = [] 

559 for d in denoms(expr, gen): 

560 singularities.extend(solvify(d, gen, domain)) 

561 if not continuous: 

562 domain = continuous_domain(expanded_e, gen, domain) 

563 

564 include_x = '=' in expr.rel_op and expr.rel_op != '!=' 

565 

566 try: 

567 discontinuities = set(domain.boundary - 

568 FiniteSet(domain.inf, domain.sup)) 

569 # remove points that are not between inf and sup of domain 

570 critical_points = FiniteSet(*(solns + singularities + list( 

571 discontinuities))).intersection( 

572 Interval(domain.inf, domain.sup, 

573 domain.inf not in domain, domain.sup not in domain)) 

574 if all(r.is_number for r in critical_points): 

575 reals = _nsort(critical_points, separated=True)[0] 

576 else: 

577 sifted = sift(critical_points, lambda x: x.is_extended_real) 

578 if sifted[None]: 

579 # there were some roots that weren't known 

580 # to be real 

581 raise NotImplementedError 

582 try: 

583 reals = sifted[True] 

584 if len(reals) > 1: 

585 reals = sorted(reals) 

586 except TypeError: 

587 raise NotImplementedError 

588 except NotImplementedError: 

589 raise NotImplementedError('sorting of these roots is not supported') 

590 

591 # If expr contains imaginary coefficients, only take real 

592 # values of x for which the imaginary part is 0 

593 make_real = S.Reals 

594 if im(expanded_e) != S.Zero: 

595 check = True 

596 im_sol = FiniteSet() 

597 try: 

598 a = solveset(im(expanded_e), gen, domain) 

599 if not isinstance(a, Interval): 

600 for z in a: 

601 if z not in singularities and valid(z) and z.is_extended_real: 

602 im_sol += FiniteSet(z) 

603 else: 

604 start, end = a.inf, a.sup 

605 for z in _nsort(critical_points + FiniteSet(end)): 

606 valid_start = valid(start) 

607 if start != end: 

608 valid_z = valid(z) 

609 pt = _pt(start, z) 

610 if pt not in singularities and pt.is_extended_real and valid(pt): 

611 if valid_start and valid_z: 

612 im_sol += Interval(start, z) 

613 elif valid_start: 

614 im_sol += Interval.Ropen(start, z) 

615 elif valid_z: 

616 im_sol += Interval.Lopen(start, z) 

617 else: 

618 im_sol += Interval.open(start, z) 

619 start = z 

620 for s in singularities: 

621 im_sol -= FiniteSet(s) 

622 except (TypeError): 

623 im_sol = S.Reals 

624 check = False 

625 

626 if im_sol is S.EmptySet: 

627 raise ValueError(filldedent(''' 

628 %s contains imaginary parts which cannot be 

629 made 0 for any value of %s satisfying the 

630 inequality, leading to relations like I < 0. 

631 ''' % (expr.subs(gen, _gen), _gen))) 

632 

633 make_real = make_real.intersect(im_sol) 

634 

635 sol_sets = [S.EmptySet] 

636 

637 start = domain.inf 

638 if start in domain and valid(start) and start.is_finite: 

639 sol_sets.append(FiniteSet(start)) 

640 

641 for x in reals: 

642 end = x 

643 

644 if valid(_pt(start, end)): 

645 sol_sets.append(Interval(start, end, True, True)) 

646 

647 if x in singularities: 

648 singularities.remove(x) 

649 else: 

650 if x in discontinuities: 

651 discontinuities.remove(x) 

652 _valid = valid(x) 

653 else: # it's a solution 

654 _valid = include_x 

655 if _valid: 

656 sol_sets.append(FiniteSet(x)) 

657 

658 start = end 

659 

660 end = domain.sup 

661 if end in domain and valid(end) and end.is_finite: 

662 sol_sets.append(FiniteSet(end)) 

663 

664 if valid(_pt(start, end)): 

665 sol_sets.append(Interval.open(start, end)) 

666 

667 if im(expanded_e) != S.Zero and check: 

668 rv = (make_real).intersect(_domain) 

669 else: 

670 rv = Intersection( 

671 (Union(*sol_sets)), make_real, _domain).subs(gen, _gen) 

672 

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

674 

675 

676def _pt(start, end): 

677 """Return a point between start and end""" 

678 if not start.is_infinite and not end.is_infinite: 

679 pt = (start + end)/2 

680 elif start.is_infinite and end.is_infinite: 

681 pt = S.Zero 

682 else: 

683 if (start.is_infinite and start.is_extended_positive is None or 

684 end.is_infinite and end.is_extended_positive is None): 

685 raise ValueError('cannot proceed with unsigned infinite values') 

686 if (end.is_infinite and end.is_extended_negative or 

687 start.is_infinite and start.is_extended_positive): 

688 start, end = end, start 

689 # if possible, use a multiple of self which has 

690 # better behavior when checking assumptions than 

691 # an expression obtained by adding or subtracting 1 

692 if end.is_infinite: 

693 if start.is_extended_positive: 

694 pt = start*2 

695 elif start.is_extended_negative: 

696 pt = start*S.Half 

697 else: 

698 pt = start + 1 

699 elif start.is_infinite: 

700 if end.is_extended_positive: 

701 pt = end*S.Half 

702 elif end.is_extended_negative: 

703 pt = end*2 

704 else: 

705 pt = end - 1 

706 return pt 

707 

708 

709def _solve_inequality(ie, s, linear=False): 

710 """Return the inequality with s isolated on the left, if possible. 

711 If the relationship is non-linear, a solution involving And or Or 

712 may be returned. False or True are returned if the relationship 

713 is never True or always True, respectively. 

714 

715 If `linear` is True (default is False) an `s`-dependent expression 

716 will be isolated on the left, if possible 

717 but it will not be solved for `s` unless the expression is linear 

718 in `s`. Furthermore, only "safe" operations which do not change the 

719 sense of the relationship are applied: no division by an unsigned 

720 value is attempted unless the relationship involves Eq or Ne and 

721 no division by a value not known to be nonzero is ever attempted. 

722 

723 Examples 

724 ======== 

725 

726 >>> from sympy import Eq, Symbol 

727 >>> from sympy.solvers.inequalities import _solve_inequality as f 

728 >>> from sympy.abc import x, y 

729 

730 For linear expressions, the symbol can be isolated: 

731 

732 >>> f(x - 2 < 0, x) 

733 x < 2 

734 >>> f(-x - 6 < x, x) 

735 x > -3 

736 

737 Sometimes nonlinear relationships will be False 

738 

739 >>> f(x**2 + 4 < 0, x) 

740 False 

741 

742 Or they may involve more than one region of values: 

743 

744 >>> f(x**2 - 4 < 0, x) 

745 (-2 < x) & (x < 2) 

746 

747 To restrict the solution to a relational, set linear=True 

748 and only the x-dependent portion will be isolated on the left: 

749 

750 >>> f(x**2 - 4 < 0, x, linear=True) 

751 x**2 < 4 

752 

753 Division of only nonzero quantities is allowed, so x cannot 

754 be isolated by dividing by y: 

755 

756 >>> y.is_nonzero is None # it is unknown whether it is 0 or not 

757 True 

758 >>> f(x*y < 1, x) 

759 x*y < 1 

760 

761 And while an equality (or inequality) still holds after dividing by a 

762 non-zero quantity 

763 

764 >>> nz = Symbol('nz', nonzero=True) 

765 >>> f(Eq(x*nz, 1), x) 

766 Eq(x, 1/nz) 

767 

768 the sign must be known for other inequalities involving > or <: 

769 

770 >>> f(x*nz <= 1, x) 

771 nz*x <= 1 

772 >>> p = Symbol('p', positive=True) 

773 >>> f(x*p <= 1, x) 

774 x <= 1/p 

775 

776 When there are denominators in the original expression that 

777 are removed by expansion, conditions for them will be returned 

778 as part of the result: 

779 

780 >>> f(x < x*(2/x - 1), x) 

781 (x < 1) & Ne(x, 0) 

782 """ 

783 from sympy.solvers.solvers import denoms 

784 if s not in ie.free_symbols: 

785 return ie 

786 if ie.rhs == s: 

787 ie = ie.reversed 

788 if ie.lhs == s and s not in ie.rhs.free_symbols: 

789 return ie 

790 

791 def classify(ie, s, i): 

792 # return True or False if ie evaluates when substituting s with 

793 # i else None (if unevaluated) or NaN (when there is an error 

794 # in evaluating) 

795 try: 

796 v = ie.subs(s, i) 

797 if v is S.NaN: 

798 return v 

799 elif v not in (True, False): 

800 return 

801 return v 

802 except TypeError: 

803 return S.NaN 

804 

805 rv = None 

806 oo = S.Infinity 

807 expr = ie.lhs - ie.rhs 

808 try: 

809 p = Poly(expr, s) 

810 if p.degree() == 0: 

811 rv = ie.func(p.as_expr(), 0) 

812 elif not linear and p.degree() > 1: 

813 # handle in except clause 

814 raise NotImplementedError 

815 except (PolynomialError, NotImplementedError): 

816 if not linear: 

817 try: 

818 rv = reduce_rational_inequalities([[ie]], s) 

819 except PolynomialError: 

820 rv = solve_univariate_inequality(ie, s) 

821 # remove restrictions wrt +/-oo that may have been 

822 # applied when using sets to simplify the relationship 

823 okoo = classify(ie, s, oo) 

824 if okoo is S.true and classify(rv, s, oo) is S.false: 

825 rv = rv.subs(s < oo, True) 

826 oknoo = classify(ie, s, -oo) 

827 if (oknoo is S.true and 

828 classify(rv, s, -oo) is S.false): 

829 rv = rv.subs(-oo < s, True) 

830 rv = rv.subs(s > -oo, True) 

831 if rv is S.true: 

832 rv = (s <= oo) if okoo is S.true else (s < oo) 

833 if oknoo is not S.true: 

834 rv = And(-oo < s, rv) 

835 else: 

836 p = Poly(expr) 

837 

838 conds = [] 

839 if rv is None: 

840 e = p.as_expr() # this is in expanded form 

841 # Do a safe inversion of e, moving non-s terms 

842 # to the rhs and dividing by a nonzero factor if 

843 # the relational is Eq/Ne; for other relationals 

844 # the sign must also be positive or negative 

845 rhs = 0 

846 b, ax = e.as_independent(s, as_Add=True) 

847 e -= b 

848 rhs -= b 

849 ef = factor_terms(e) 

850 a, e = ef.as_independent(s, as_Add=False) 

851 if (a.is_zero != False or # don't divide by potential 0 

852 a.is_negative == 

853 a.is_positive is None and # if sign is not known then 

854 ie.rel_op not in ('!=', '==')): # reject if not Eq/Ne 

855 e = ef 

856 a = S.One 

857 rhs /= a 

858 if a.is_positive: 

859 rv = ie.func(e, rhs) 

860 else: 

861 rv = ie.reversed.func(e, rhs) 

862 

863 # return conditions under which the value is 

864 # valid, too. 

865 beginning_denoms = denoms(ie.lhs) | denoms(ie.rhs) 

866 current_denoms = denoms(rv) 

867 for d in beginning_denoms - current_denoms: 

868 c = _solve_inequality(Eq(d, 0), s, linear=linear) 

869 if isinstance(c, Eq) and c.lhs == s: 

870 if classify(rv, s, c.rhs) is S.true: 

871 # rv is permitting this value but it shouldn't 

872 conds.append(~c) 

873 for i in (-oo, oo): 

874 if (classify(rv, s, i) is S.true and 

875 classify(ie, s, i) is not S.true): 

876 conds.append(s < i if i is oo else i < s) 

877 

878 conds.append(rv) 

879 return And(*conds) 

880 

881 

882def _reduce_inequalities(inequalities, symbols): 

883 # helper for reduce_inequalities 

884 

885 poly_part, abs_part = {}, {} 

886 other = [] 

887 

888 for inequality in inequalities: 

889 

890 expr, rel = inequality.lhs, inequality.rel_op # rhs is 0 

891 

892 # check for gens using atoms which is more strict than free_symbols to 

893 # guard against EX domain which won't be handled by 

894 # reduce_rational_inequalities 

895 gens = expr.atoms(Symbol) 

896 

897 if len(gens) == 1: 

898 gen = gens.pop() 

899 else: 

900 common = expr.free_symbols & symbols 

901 if len(common) == 1: 

902 gen = common.pop() 

903 other.append(_solve_inequality(Relational(expr, 0, rel), gen)) 

904 continue 

905 else: 

906 raise NotImplementedError(filldedent(''' 

907 inequality has more than one symbol of interest. 

908 ''')) 

909 

910 if expr.is_polynomial(gen): 

911 poly_part.setdefault(gen, []).append((expr, rel)) 

912 else: 

913 components = expr.find(lambda u: 

914 u.has(gen) and ( 

915 u.is_Function or u.is_Pow and not u.exp.is_Integer)) 

916 if components and all(isinstance(i, Abs) for i in components): 

917 abs_part.setdefault(gen, []).append((expr, rel)) 

918 else: 

919 other.append(_solve_inequality(Relational(expr, 0, rel), gen)) 

920 

921 poly_reduced = [reduce_rational_inequalities([exprs], gen) for gen, exprs in poly_part.items()] 

922 abs_reduced = [reduce_abs_inequalities(exprs, gen) for gen, exprs in abs_part.items()] 

923 

924 return And(*(poly_reduced + abs_reduced + other)) 

925 

926 

927def reduce_inequalities(inequalities, symbols=[]): 

928 """Reduce a system of inequalities with rational coefficients. 

929 

930 Examples 

931 ======== 

932 

933 >>> from sympy.abc import x, y 

934 >>> from sympy import reduce_inequalities 

935 

936 >>> reduce_inequalities(0 <= x + 3, []) 

937 (-3 <= x) & (x < oo) 

938 

939 >>> reduce_inequalities(0 <= x + y*2 - 1, [x]) 

940 (x < oo) & (x >= 1 - 2*y) 

941 """ 

942 if not iterable(inequalities): 

943 inequalities = [inequalities] 

944 inequalities = [sympify(i) for i in inequalities] 

945 

946 gens = set().union(*[i.free_symbols for i in inequalities]) 

947 

948 if not iterable(symbols): 

949 symbols = [symbols] 

950 symbols = (set(symbols) or gens) & gens 

951 if any(i.is_extended_real is False for i in symbols): 

952 raise TypeError(filldedent(''' 

953 inequalities cannot contain symbols that are not real. 

954 ''')) 

955 

956 # make vanilla symbol real 

957 recast = {i: Dummy(i.name, extended_real=True) 

958 for i in gens if i.is_extended_real is None} 

959 inequalities = [i.xreplace(recast) for i in inequalities] 

960 symbols = {i.xreplace(recast) for i in symbols} 

961 

962 # prefilter 

963 keep = [] 

964 for i in inequalities: 

965 if isinstance(i, Relational): 

966 i = i.func(i.lhs.as_expr() - i.rhs.as_expr(), 0) 

967 elif i not in (True, False): 

968 i = Eq(i, 0) 

969 if i == True: 

970 continue 

971 elif i == False: 

972 return S.false 

973 if i.lhs.is_number: 

974 raise NotImplementedError( 

975 "could not determine truth value of %s" % i) 

976 keep.append(i) 

977 inequalities = keep 

978 del keep 

979 

980 # solve system 

981 rv = _reduce_inequalities(inequalities, symbols) 

982 

983 # restore original symbols and return 

984 return rv.xreplace({v: k for k, v in recast.items()})