Coverage for sympy/solvers/solveset.py : 91%

""" 

This module contains functions to: 

- solve a single equation for a single variable, in any domain either real or complex. 

- solve a system of linear equations with N variables and M equations. 

- solve a system of Non Linear Equations with N variables and M equations 

""" 

from __future__ import print_function, division 

from sympy.core.sympify import sympify 

from sympy.core import S, Pow, Dummy, pi, Expr, Wild, Mul, Equality 

from sympy.core.numbers import I, Number, Rational, oo 

from sympy.core.function import (Lambda, expand, expand_complex) 

from sympy.core.relational import Eq 

from sympy.simplify.simplify import simplify, fraction, trigsimp 

from sympy.core.symbol import Symbol 

from sympy.functions import (log, Abs, tan, cot, sin, cos, sec, csc, exp, 

acos, asin, acsc, asec, arg, 

piecewise_fold) 

from sympy.functions.elementary.trigonometric import (TrigonometricFunction, 

HyperbolicFunction) 

from sympy.functions.elementary.miscellaneous import real_root 

from sympy.sets import (FiniteSet, EmptySet, imageset, Interval, Intersection, 

Union, ConditionSet, ImageSet) 

from sympy.matrices import Matrix 

from sympy.polys import (roots, Poly, degree, together, PolynomialError, 

RootOf) 

from sympy.solvers.solvers import checksol, denoms, unrad, _simple_dens 

from sympy.solvers.polysys import solve_poly_system 

from sympy.solvers.inequalities import solve_univariate_inequality 

from sympy.utilities import filldedent 

from sympy.core.compatibility import ordered, default_sort_key 

def _invert(f_x, y, x, domain=S.Complexes): 

""" 

Reduce the complex valued equation ``f(x) = y`` to a set of equations 

``{g(x) = h_1(y), g(x) = h_2(y), ..., g(x) = h_n(y) }`` where ``g(x)`` is 

a simpler function than ``f(x)``. The return value is a tuple ``(g(x), 

set_h)``, where ``g(x)`` is a function of ``x`` and ``set_h`` is 

the set of function ``{h_1(y), h_2(y), ..., h_n(y)}``. 

Here, ``y`` is not necessarily a symbol. 

The ``set_h`` contains the functions along with the information 

about their domain in which they are valid, through set 

operations. For instance, if ``y = Abs(x) - n``, is inverted 

in the real domain, then, the ``set_h`` doesn't simply return 

`{-n, n}`, as the nature of `n` is unknown; rather it will return: 

`Intersection([0, oo) {n}) U Intersection((-oo, 0], {-n})` 

By default, the complex domain is used but note that inverting even 

seemingly simple functions like ``exp(x)`` can give very different 

result in the complex domain than are obtained in the real domain. 

(In the case of ``exp(x)``, the inversion via ``log`` is multi-valued 

in the complex domain, having infinitely many branches.) 

If you are working with real values only (or you are not sure which 

function to use) you should probably use set the domain to 

``S.Reals`` (or use `invert\_real` which does that automatically). 

Examples 

======== 

>>> from sympy.solvers.solveset import invert_complex, invert_real 

>>> from sympy.abc import x, y 

>>> from sympy import exp, log 

When does exp(x) == y? 

>>> invert_complex(exp(x), y, x) 

(x, ImageSet(Lambda(_n, I*(2*_n*pi + arg(y)) + log(Abs(y))), Integers())) 

>>> invert_real(exp(x), y, x) 

(x, Intersection((-oo, oo), {log(y)})) 

When does exp(x) == 1? 

>>> invert_complex(exp(x), 1, x) 

(x, ImageSet(Lambda(_n, 2*_n*I*pi), Integers())) 

>>> invert_real(exp(x), 1, x) 

(x, {0}) 

See Also 

======== 

invert_real, invert_complex 

""" 

x = sympify(x) 

if not x.is_Symbol: 

raise ValueError("x must be a symbol") 

f_x = sympify(f_x) 

if not f_x.has(x): 

raise ValueError("Inverse of constant function doesn't exist") 

y = sympify(y) 

if y.has(x): 

raise ValueError("y should be independent of x ") 

if domain.is_subset(S.Reals): 

x, s = _invert_real(f_x, FiniteSet(y), x) 

else: 

x, s = _invert_complex(f_x, FiniteSet(y), x) 

return x, s.intersection(domain) if isinstance(s, FiniteSet) else s 

invert_complex = _invert 

def invert_real(f_x, y, x, domain=S.Reals): 

""" 

Inverts a real-valued function. Same as _invert, but sets 

the domain to ``S.Reals`` before inverting. 

""" 

return _invert(f_x, y, x, domain) 

def _invert_real(f, g_ys, symbol): 

"""Helper function for _invert.""" 

if f == symbol: 

return (f, g_ys) 

n = Dummy('n', real=True) 

if hasattr(f, 'inverse') and not isinstance(f, ( 

TrigonometricFunction, 

HyperbolicFunction, 

)): 

if len(f.args) > 1: 

raise ValueError("Only functions with one argument are supported.") 

return _invert_real(f.args[0], 

imageset(Lambda(n, f.inverse()(n)), g_ys), 

symbol) 

if isinstance(f, Abs): 

pos = Interval(0, S.Infinity) 

neg = Interval(S.NegativeInfinity, 0) 

return _invert_real(f.args[0], 

Union(imageset(Lambda(n, n), g_ys).intersect(pos), 

imageset(Lambda(n, -n), g_ys).intersect(neg)), symbol) 

if f.is_Add: 

# f = g + h 

g, h = f.as_independent(symbol) 

if g is not S.Zero: 

return _invert_real(h, imageset(Lambda(n, n - g), g_ys), symbol) 

if f.is_Mul: 

# f = g*h 

g, h = f.as_independent(symbol) 

if g is not S.One: 

return _invert_real(h, imageset(Lambda(n, n/g), g_ys), symbol) 

if f.is_Pow: 

base, expo = f.args 

base_has_sym = base.has(symbol) 

expo_has_sym = expo.has(symbol) 

if not expo_has_sym: 

res = imageset(Lambda(n, real_root(n, expo)), g_ys) 

if expo.is_rational: 

numer, denom = expo.as_numer_denom() 

if numer == S.One or numer == - S.One: 

return _invert_real(base, res, symbol) 

else: 

if numer % 2 == 0: 

n = Dummy('n') 

neg_res = imageset(Lambda(n, -n), res) 

return _invert_real(base, res + neg_res, symbol) 

else: 

return _invert_real(base, res, symbol) 

else: 

if not base.is_positive: 

raise ValueError("x**w where w is irrational is not " 

"defined for negative x") 

return _invert_real(base, res, symbol) 

if not base_has_sym: 

return _invert_real(expo, 

imageset(Lambda(n, log(n)/log(base)), g_ys), symbol) 

if isinstance(f, TrigonometricFunction): 

if isinstance(g_ys, FiniteSet): 

def inv(trig): 

if isinstance(f, (sin, csc)): 

F = asin if isinstance(f, sin) else acsc 

return (lambda a: n*pi + (-1)**n*F(a),) 

if isinstance(f, (cos, sec)): 

F = acos if isinstance(f, cos) else asec 

return ( 

lambda a: 2*n*pi + F(a), 

lambda a: 2*n*pi - F(a),) 

if isinstance(f, (tan, cot)): 

return (lambda a: n*pi + f.inverse()(a),) 

n = Dummy('n', integer=True) 

invs = S.EmptySet 

for L in inv(f): 

invs += Union(*[imageset(Lambda(n, L(g)), S.Integers) for g in g_ys]) 

return _invert_real(f.args[0], invs, symbol) 

return (f, g_ys) 

def _invert_complex(f, g_ys, symbol): 

"""Helper function for _invert.""" 

if f == symbol: 

return (f, g_ys) 

n = Dummy('n') 

if f.is_Add: 

# f = g + h 

g, h = f.as_independent(symbol) 

if g is not S.Zero: 

return _invert_complex(h, imageset(Lambda(n, n - g), g_ys), symbol) 

if f.is_Mul: 

# f = g*h 

g, h = f.as_independent(symbol) 

if g is not S.One: 

return _invert_complex(h, imageset(Lambda(n, n/g), g_ys), symbol) 

if hasattr(f, 'inverse') and \ 

not isinstance(f, TrigonometricFunction) and \ 

not isinstance(f, exp): 

if len(f.args) > 1: 

raise ValueError("Only functions with one argument are supported.") 

return _invert_complex(f.args[0], 

imageset(Lambda(n, f.inverse()(n)), g_ys), symbol) 

if isinstance(f, exp): 

if isinstance(g_ys, FiniteSet): 

exp_invs = Union(*[imageset(Lambda(n, I*(2*n*pi + arg(g_y)) + 

log(Abs(g_y))), S.Integers) 

for g_y in g_ys if g_y != 0]) 

return _invert_complex(f.args[0], exp_invs, symbol) 

return (f, g_ys) 

def domain_check(f, symbol, p): 

"""Returns False if point p is infinite or any subexpression of f 

is infinite or becomes so after replacing symbol with p. If none of 

these conditions is met then True will be returned. 

Examples 

======== 

>>> from sympy import Mul, oo 

>>> from sympy.abc import x 

>>> from sympy.solvers.solveset import domain_check 

>>> g = 1/(1 + (1/(x + 1))**2) 

>>> domain_check(g, x, -1) 

False 

>>> domain_check(x**2, x, 0) 

True 

>>> domain_check(1/x, x, oo) 

False 

* The function relies on the assumption that the original form 

of the equation has not been changed by automatic simplification. 

>>> domain_check(x/x, x, 0) # x/x is automatically simplified to 1 

True 

* To deal with automatic evaluations use evaluate=False: 

>>> domain_check(Mul(x, 1/x, evaluate=False), x, 0) 

False 

""" 

f, p = sympify(f), sympify(p) 

if p.is_infinite: 

return False 

return _domain_check(f, symbol, p) 

def _domain_check(f, symbol, p): 

# helper for domain check 

if f.is_Atom and f.is_finite: 

return True 

elif f.subs(symbol, p).is_infinite: 

return False 

else: 

return all([_domain_check(g, symbol, p) 

for g in f.args]) 

def _is_finite_with_finite_vars(f, domain=S.Complexes): 

""" 

Return True if the given expression is finite. For symbols that 

don't assign a value for `complex` and/or `real`, the domain will 

be used to assign a value; symbols that don't assign a value 

for `finite` will be made finite. All other assumptions are 

left unmodified. 

""" 

def assumptions(s): 

A = s.assumptions0 

if A.get('finite', None) is None: 

A['finite'] = True 

A.setdefault('complex', True) 

A.setdefault('real', domain.is_subset(S.Reals)) 

return A 

reps = {s: Dummy(**assumptions(s)) for s in f.free_symbols} 

return f.xreplace(reps).is_finite 

def _is_function_class_equation(func_class, f, symbol): 

""" Tests whether the equation is an equation of the given function class. 

The given equation belongs to the given function class if it is 

comprised of functions of the function class which are multiplied by 

or added to expressions independent of the symbol. In addition, the 

arguments of all such functions must be linear in the symbol as well. 

Examples 

======== 

>>> from sympy.solvers.solveset import _is_function_class_equation 

>>> from sympy import tan, sin, tanh, sinh, exp 

>>> from sympy.abc import x 

>>> from sympy.functions.elementary.trigonometric import (TrigonometricFunction, 

... HyperbolicFunction) 

>>> _is_function_class_equation(TrigonometricFunction, exp(x) + tan(x), x) 

False 

>>> _is_function_class_equation(TrigonometricFunction, tan(x) + sin(x), x) 

True 

>>> _is_function_class_equation(TrigonometricFunction, tan(x**2), x) 

False 

>>> _is_function_class_equation(TrigonometricFunction, tan(x + 2), x) 

True 

>>> _is_function_class_equation(HyperbolicFunction, tanh(x) + sinh(x), x) 

True 

""" 

if f.is_Mul or f.is_Add: 

return all(_is_function_class_equation(func_class, arg, symbol) 

for arg in f.args) 

if f.is_Pow: 

if not f.exp.has(symbol): 

return _is_function_class_equation(func_class, f.base, symbol) 

else: 

return False 

if not f.has(symbol): 

return True 

if isinstance(f, func_class): 

try: 

g = Poly(f.args[0], symbol) 

return g.degree() <= 1 

except PolynomialError: 

return False 

else: 

return False 

def _solve_as_rational(f, symbol, domain): 

""" solve rational functions""" 

f = together(f, deep=True) 

g, h = fraction(f) 

if not h.has(symbol): 

return _solve_as_poly(g, symbol, domain) 

else: 

valid_solns = _solveset(g, symbol, domain) 

invalid_solns = _solveset(h, symbol, domain) 

return valid_solns - invalid_solns 

def _solve_trig(f, symbol, domain): 

""" Helper to solve trigonometric equations """ 

f = trigsimp(f) 

f_original = f 

f = f.rewrite(exp) 

f = together(f) 

g, h = fraction(f) 

y = Dummy('y') 

g, h = g.expand(), h.expand() 

g, h = g.subs(exp(I*symbol), y), h.subs(exp(I*symbol), y) 

if g.has(symbol) or h.has(symbol): 

return ConditionSet(symbol, Eq(f, 0), S.Reals) 

solns = solveset_complex(g, y) - solveset_complex(h, y) 

if isinstance(solns, FiniteSet): 

result = Union(*[invert_complex(exp(I*symbol), s, symbol)[1] 

for s in solns]) 

return Intersection(result, domain) 

elif solns is S.EmptySet: 

return S.EmptySet 

else: 

return ConditionSet(symbol, Eq(f_original, 0), S.Reals) 

def _solve_as_poly(f, symbol, domain=S.Complexes): 

""" 

Solve the equation using polynomial techniques if it already is a 

polynomial equation or, with a change of variables, can be made so. 

""" 

result = None 

if f.is_polynomial(symbol): 

solns = roots(f, symbol, cubics=True, quartics=True, 

quintics=True, domain='EX') 

num_roots = sum(solns.values()) 

if degree(f, symbol) <= num_roots: 

result = FiniteSet(*solns.keys()) 

else: 

poly = Poly(f, symbol) 

solns = poly.all_roots() 

if poly.degree() <= len(solns): 

result = FiniteSet(*solns) 

else: 

result = ConditionSet(symbol, Eq(f, 0), domain) 

else: 

poly = Poly(f) 

if poly is None: 

result = ConditionSet(symbol, Eq(f, 0), domain) 

gens = [g for g in poly.gens if g.has(symbol)] 

if len(gens) == 1: 

poly = Poly(poly, gens[0]) 

gen = poly.gen 

deg = poly.degree() 

poly = Poly(poly.as_expr(), poly.gen, composite=True) 

poly_solns = FiniteSet(*roots(poly, cubics=True, quartics=True, 

quintics=True).keys()) 

if len(poly_solns) < deg: 

result = ConditionSet(symbol, Eq(f, 0), domain) 

if gen != symbol: 

y = Dummy('y') 

inverter = invert_real if domain.is_subset(S.Reals) else invert_complex 

lhs, rhs_s = inverter(gen, y, symbol) 

if lhs == symbol: 

result = Union(*[rhs_s.subs(y, s) for s in poly_solns]) 

else: 

result = ConditionSet(symbol, Eq(f, 0), domain) 

else: 

result = ConditionSet(symbol, Eq(f, 0), domain) 

if result is not None: 

if isinstance(result, FiniteSet): 

# this is to simplify solutions like -sqrt(-I) to sqrt(2)/2 

# - sqrt(2)*I/2. We are not expanding for solution with free 

# variables because that makes the solution more complicated. For 

# example expand_complex(a) returns re(a) + I*im(a) 

if all([s.free_symbols == set() and not isinstance(s, RootOf) 

for s in result]): 

s = Dummy('s') 

result = imageset(Lambda(s, expand_complex(s)), result) 

if isinstance(result, FiniteSet): 

result = result.intersection(domain) 

return result 

else: 

return ConditionSet(symbol, Eq(f, 0), domain) 

def _has_rational_power(expr, symbol): 

""" 

Returns (bool, den) where bool is True if the term has a 

non-integer rational power and den is the denominator of the 

expression's exponent. 

Examples 

======== 

>>> from sympy.solvers.solveset import _has_rational_power 

>>> from sympy import sqrt 

>>> from sympy.abc import x 

>>> _has_rational_power(sqrt(x), x) 

(True, 2) 

>>> _has_rational_power(x**2, x) 

(False, 1) 

""" 

a, p, q = Wild('a'), Wild('p'), Wild('q') 

pattern_match = expr.match(a*p**q) or {} 

if pattern_match.get(a, S.Zero) is S.Zero: 

return (False, S.One) 

elif p not in pattern_match.keys(): 

return (False, S.One) 

elif isinstance(pattern_match[q], Rational) \ 

and pattern_match[p].has(symbol): 

if not pattern_match[q].q == S.One: 

return (True, pattern_match[q].q) 

if not isinstance(pattern_match[a], Pow) \ 

or isinstance(pattern_match[a], Mul): 

return (False, S.One) 

else: 

return _has_rational_power(pattern_match[a], symbol) 

def _solve_radical(f, symbol, solveset_solver): 

""" Helper function to solve equations with radicals """ 

eq, cov = unrad(f) 

if not cov: 

result = solveset_solver(eq, symbol) - \ 

Union(*[solveset_solver(g, symbol) for g in denoms(f, [symbol])]) 

else: 

y, yeq = cov 

if not solveset_solver(y - I, y): 

yreal = Dummy('yreal', real=True) 

yeq = yeq.xreplace({y: yreal}) 

eq = eq.xreplace({y: yreal}) 

y = yreal 

g_y_s = solveset_solver(yeq, symbol) 

f_y_sols = solveset_solver(eq, y) 

result = Union(*[imageset(Lambda(y, g_y), f_y_sols)