# Source code for sympy.solvers.solveset

"""
This module contains functions to:

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

- solve a single transcendental equation for a single variable in any domain either real or complex.
(currently supports solving in real domain only)

- 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.containers import Tuple
from sympy.core.facts import InconsistentAssumptions
from sympy.core.numbers import I, Number, Rational, oo
from sympy.core.function import (Lambda, expand_complex, AppliedUndef,
expand_log)
from sympy.core.relational import Eq, Ne
from sympy.core.symbol import Symbol
from sympy.simplify.simplify import simplify, fraction, trigsimp
from sympy.simplify import powdenest, logcombine
from sympy.functions import (log, Abs, tan, cot, sin, cos, sec, csc, exp,
acos, asin, acsc, asec, arg,
piecewise_fold, Piecewise)
from sympy.functions.elementary.trigonometric import (TrigonometricFunction,
HyperbolicFunction)
from sympy.functions.elementary.miscellaneous import real_root
from sympy.logic.boolalg import And
from sympy.sets import (FiniteSet, EmptySet, imageset, Interval, Intersection,
Union, ConditionSet, ImageSet, Complement, Contains)
from sympy.sets.sets import Set
from sympy.matrices import Matrix
from sympy.polys import (roots, Poly, degree, together, PolynomialError,
RootOf, factor)
from sympy.solvers.solvers import (checksol, denoms, unrad,
_simple_dens, recast_to_symbols)
from sympy.solvers.polysys import solve_poly_system
from sympy.solvers.inequalities import solve_univariate_inequality
from sympy.utilities import filldedent
from sympy.utilities.iterables import numbered_symbols
from sympy.calculus.util import periodicity, continuous_domain
from sympy.core.compatibility import ordered, default_sort_key, is_sequence

from types import GeneratorType

"""Return f, with all objects given by atoms replaced with
Dummy symbols, d, and the list of replacements, (d, e),
where e is an object of type given by atoms in which
any other instances of atoms have been recursively replaced with
Dummy symbols, too. The tuples are ordered so that if they are
applied in sequence, the origin f will be restored.

Examples
========

>>> from sympy import cos
>>> from sympy.abc import x

>>> f = cos(cos(x) + 1)
>>> f, reps = _masked(cos(1 + cos(x)), cos)
>>> f
_a1
>>> reps
[(_a1, cos(_a0 + 1)), (_a0, cos(x))]
>>> for d, e in reps:
...     f = f.xreplace({d: e})
>>> f
cos(cos(x) + 1)
"""
sym = numbered_symbols('a', cls=Dummy, real=True)
for a in ordered(f.atoms(*atoms)):
a = a.replace(*i)
for i, (o, n) in enumerate(mask):
f = f.replace(o, n)

def _invert(f_x, y, x, domain=S.Complexes):
r"""
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 the domain in which they are valid, through set
operations. For instance, if y = Abs(x) - n is inverted
in the real domain, then set_h is not simply
{-n, n} as the nature of n is unknown; rather, it is:
Intersection([0, oo) {n}) U Intersection((-oo, 0], {-n})

By default, the complex domain is used which means that inverting even
seemingly simple functions like exp(x) will give very different
results from those 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 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(Reals, {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})

========
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 x not in f_x.free_symbols:
raise ValueError("Inverse of constant function doesn't exist")
y = sympify(y)
if x in y.free_symbols:
raise ValueError("y should be independent of x ")

if domain.is_subset(S.Reals):
x1, s = _invert_real(f_x, FiniteSet(y), x)
else:
x1, s = _invert_complex(f_x, FiniteSet(y), x)

if not isinstance(s, FiniteSet) or x1 != x:
return x1, s

return x1, s.intersection(domain)

invert_complex = _invert

[docs]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):
return _invert_abs(f.args[0], g_ys, symbol)

# 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 denom % 2 == 0:
base_positive = solveset(base >= 0, symbol, S.Reals)
res = imageset(Lambda(n, real_root(n, expo)
), g_ys.intersect(
Interval.Ropen(S.Zero, S.Infinity)))
_inv, _set = _invert_real(base, res, symbol)
return (_inv, _set.intersect(base_positive))

elif 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:
rhs = g_ys.args[0]
if base.is_positive:
return _invert_real(expo,
imageset(Lambda(n, log(n, base, evaluate=False)), g_ys), symbol)
elif base.is_negative:
from sympy.core.power import integer_log
s, b = integer_log(rhs, base)
if b:
return _invert_real(expo, FiniteSet(s), symbol)
else:
return _invert_real(expo, S.EmptySet, symbol)
elif base.is_zero:
one = Eq(rhs, 1)
if one == S.true:
# special case: 0**x - 1
return _invert_real(expo, FiniteSet(0), symbol)
elif one == S.false:
return _invert_real(expo, S.EmptySet, 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')

# 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:
if g in set([S.NegativeInfinity, S.ComplexInfinity, S.Infinity]):
return (h, S.EmptySet)
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, HyperbolicFunction) 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 _invert_abs(f, g_ys, symbol):
"""Helper function for inverting absolute value functions.

Returns the complete result of inverting an absolute value
function along with the conditions which must also be satisfied.

If it is certain that all these conditions are met, a FiniteSet
of all possible solutions is returned. If any condition cannot be
satisfied, an EmptySet is returned. Otherwise, a ConditionSet
of the solutions, with all the required conditions specified, is
returned.

"""
if not g_ys.is_FiniteSet:
# this could be used for FiniteSet, but the
# results are more compact if they aren't, e.g.
# ConditionSet(x, Contains(n, Interval(0, oo)), {-n, n}) vs
# Union(Intersection(Interval(0, oo), {n}), Intersection(Interval(-oo, 0), {-n}))
# for the solution of abs(x) - n
pos = Intersection(g_ys, Interval(0, S.Infinity))
parg = _invert_real(f, pos, symbol)
narg = _invert_real(-f, pos, symbol)
if parg[0] != narg[0]:
raise NotImplementedError
return parg[0], Union(narg[1], parg[1])

# check conditions: all these must be true. If any are unknown
# then return them as conditions which must be satisfied
unknown = []
for a in g_ys.args:
ok = a.is_nonnegative if a.is_Number else a.is_positive
if ok is None:
unknown.append(a)
elif not ok:
return symbol, S.EmptySet
if unknown:
conditions = And(*[Contains(i, Interval(0, oo))
for i in unknown])
else:
conditions = True
n = Dummy('n', real=True)
# this is slightly different than above: instead of solving
# +/-f on positive values, here we solve for f on +/- g_ys
g_x, values = _invert_real(f, Union(
imageset(Lambda(n, n), g_ys),
imageset(Lambda(n, -n), g_ys)), symbol)
return g_x, ConditionSet(g_x, conditions, values)

[docs]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
A.setdefault('finite', A.get('finite', True))
if domain.is_subset(S.Reals):
# if this gets set it will make complex=True, too
A.setdefault('real', True)
else:
# don't change 'real' because being complex implies
A.setdefault('complex', True)
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
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
"""
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):
try:
return _solve_as_poly(g, symbol, domain)
except NotImplementedError:
# The polynomial formed from g could end up having
# coefficients in a ring over which finding roots
# isn't implemented yet, e.g. ZZ[a] for some symbol a
return ConditionSet(symbol, Eq(f, 0), 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):
"""Function to call other helpers to solve trigonometric equations """
sol1 = sol = None
try:
sol1 = _solve_trig1(f, symbol, domain)
except BaseException as error:
pass
if sol1 is None or isinstance(sol1, ConditionSet):
try:
sol = _solve_trig2(f, symbol, domain)
except BaseException as error:
sol = sol1
if isinstance(sol1, ConditionSet) and isinstance(sol, ConditionSet):
if sol1.count_ops() < sol.count_ops():
sol = sol1
else:
sol = sol1
if sol is None:
raise NotImplementedError(filldedent('''
Solution to this kind of trigonometric equations
is yet to be implemented'''))
return sol

def _solve_trig1(f, symbol, domain):
"""Primary 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, ConditionSet):
raise NotImplementedError

if isinstance(solns, FiniteSet):
if any(isinstance(s, RootOf) for s in solns):
raise NotImplementedError
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_trig2(f, symbol, domain):
"""Secondary helper to solve trigonometric equations,
called when first helper fails """
from sympy import ilcm, igcd, expand_trig, degree, simplify
f = trigsimp(f)
f_original = f
trig_functions = f.atoms(sin, cos, tan, sec, cot, csc)
trig_arguments = [e.args[0] for e in trig_functions]
denominators = []
numerators = []

for ar in trig_arguments:
try:
poly_ar = Poly(ar, symbol)

except ValueError:
raise ValueError("give up, we can't solve if this is not a polynomial in x")
if poly_ar.degree() > 1:  # degree >1 still bad
raise ValueError("degree of variable inside polynomial should not exceed one")
if poly_ar.degree() == 0:  # degree 0, don't care
continue
c = poly_ar.all_coeffs()[0]   # got the coefficient of 'symbol'
numerators.append(Rational(c).p)
denominators.append(Rational(c).q)

x = Dummy('x')

# ilcm() and igcd() require more than one argument
if len(numerators) > 1:
mu = Rational(2)*ilcm(*denominators)/igcd(*numerators)
else:
assert len(numerators) == 1
mu = Rational(2)*denominators[0]/numerators[0]

f = f.subs(symbol, mu*x)
f = f.rewrite(tan)
f = expand_trig(f)
f = together(f)

g, h = fraction(f)
y = Dummy('y')
g, h = g.expand(), h.expand()
g, h = g.subs(tan(x), y), h.subs(tan(x), y)

if g.has(x) or h.has(x):
return ConditionSet(symbol, Eq(f_original, 0), domain)
solns = solveset(g, y, S.Reals) - solveset(h, y, S.Reals)

if isinstance(solns, FiniteSet):
result = Union(*[invert_real(tan(symbol/mu), s, symbol)[1]
for s in solns])
dsol = invert_real(tan(symbol/mu), oo, symbol)[1]
if degree(h) > degree(g):                   # If degree(denom)>degree(num) then there
result = Union(result, dsol)            # would be another sol at Lim(denom-->oo)
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 symbols
# or undefined functions because that makes the solution more complicated.
# For example, expand_complex(a) returns re(a) + I*im(a)
if all([s.atoms(Symbol, AppliedUndef) == 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)

""" Helper function to solve equations with radicals """
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)
for g_y in g_y_s])

if isinstance(result, Complement) or isinstance(result,ConditionSet):
solution_set = result
else:
f_set = []  # solutions for FiniteSet
c_set = []  # solutions for ConditionSet
for s in result:
if checksol(f, symbol, s):
f_set.append(s)
else:
c_set.append(s)
solution_set = FiniteSet(*f_set) + ConditionSet(symbol, Eq(f, 0), FiniteSet(*c_set))

return solution_set

def _solve_abs(f, symbol, domain):
""" Helper function to solve equation involving absolute value function """
if not domain.is_subset(S.Reals):
raise ValueError(filldedent('''
Absolute values cannot be inverted in the
complex domain.'''))
p, q, r = Wild('p'), Wild('q'), Wild('r')
pattern_match = f.match(p*Abs(q) + r) or {}
f_p, f_q, f_r = [pattern_match.get(i, S.Zero) for i in (p, q, r)]

if not (f_p.is_zero or f_q.is_zero):
domain = continuous_domain(f_q, symbol, domain)
q_pos_cond = solve_univariate_inequality(f_q >= 0, symbol,
relational=False, domain=domain, continuous=True)
q_neg_cond = q_pos_cond.complement(domain)

sols_q_pos = solveset_real(f_p*f_q + f_r,
symbol).intersect(q_pos_cond)
sols_q_neg = solveset_real(f_p*(-f_q) + f_r,
symbol).intersect(q_neg_cond)
return Union(sols_q_pos, sols_q_neg)
else:
return ConditionSet(symbol, Eq(f, 0), domain)

def solve_decomposition(f, symbol, domain):
"""
Function to solve equations via the principle of "Decomposition
and Rewriting".

Examples
========
>>> from sympy import exp, sin, Symbol, pprint, S
>>> from sympy.solvers.solveset import solve_decomposition as sd
>>> x = Symbol('x')
>>> f1 = exp(2*x) - 3*exp(x) + 2
>>> sd(f1, x, S.Reals)
{0, log(2)}
>>> f2 = sin(x)**2 + 2*sin(x) + 1
>>> pprint(sd(f2, x, S.Reals), use_unicode=False)
3*pi
{2*n*pi + ---- | n in Integers}
2
>>> f3 = sin(x + 2)
>>> pprint(sd(f3, x, S.Reals), use_unicode=False)
{2*n*pi - 2 | n in Integers} U {2*n*pi - 2 + pi | n in Integers}

"""
from sympy.solvers.decompogen import decompogen
from sympy.calculus.util import function_range
# decompose the given function
g_s = decompogen(f, symbol)
# y_s represents the set of values for which the function g is to be
# solved.
# solutions represent the solutions of the equations g = y_s or
# g = 0 depending on the type of y_s.
# As we are interested in solving the equation: f = 0
y_s = FiniteSet(0)
for g in g_s:
frange = function_range(g, symbol, domain)
y_s = Intersection(frange, y_s)
result = S.EmptySet
if isinstance(y_s, FiniteSet):
for y in y_s:
solutions = solveset(Eq(g, y), symbol, domain)
if not isinstance(solutions, ConditionSet):
result += solutions

else:
if isinstance(y_s, ImageSet):
iter_iset = (y_s,)

elif isinstance(y_s, Union):
iter_iset = y_s.args

elif isinstance(y_s, EmptySet):
# y_s is not in the range of g in g_s, so no solution exists
#in the given domain
return y_s

for iset in iter_iset:
new_solutions = solveset(Eq(iset.lamda.expr, g), symbol, domain)
dummy_var = tuple(iset.lamda.expr.free_symbols)[0]
base_set = iset.base_set
if isinstance(new_solutions, FiniteSet):
new_exprs = new_solutions

elif isinstance(new_solutions, Intersection):
if isinstance(new_solutions.args[1], FiniteSet):
new_exprs = new_solutions.args[1]

for new_expr in new_exprs:
result += ImageSet(Lambda(dummy_var, new_expr), base_set)

if result is S.EmptySet:
return ConditionSet(symbol, Eq(f, 0), domain)

y_s = result

return y_s

def _solveset(f, symbol, domain, _check=False):
"""Helper for solveset to return a result from an expression
that has already been sympify'ed and is known to contain the
given symbol."""
# _check controls whether the answer is checked or not
from sympy.simplify.simplify import signsimp

orig_f = f
if f.is_Mul:
if coeff in set([S.ComplexInfinity, S.NegativeInfinity, S.Infinity]):
f = together(orig_f)
a, h = f.as_independent(symbol)
if m not in set([S.ComplexInfinity, S.Zero, S.Infinity,
S.NegativeInfinity]):
f = a/m + h  # XXX condition m != 0 should be added to soln

# assign the solvers to use
solver = lambda f, x, domain=domain: _solveset(f, x, domain)
inverter = lambda f, rhs, symbol: _invert(f, rhs, symbol, domain)

result = EmptySet()

if f.expand().is_zero:
return domain
elif not f.has(symbol):
return EmptySet()
elif f.is_Mul and all(_is_finite_with_finite_vars(m, domain)
for m in f.args):
# if f(x) and g(x) are both finite we can say that the solution of
# f(x)*g(x) == 0 is same as Union(f(x) == 0, g(x) == 0) is not true in
# general. g(x) can grow to infinitely large for the values where
# f(x) == 0. To be sure that we are not silently allowing any
# wrong solutions we are using this technique only if both f and g are
# finite for a finite input.
result = Union(*[solver(m, symbol) for m in f.args])
elif _is_function_class_equation(TrigonometricFunction, f, symbol) or \
_is_function_class_equation(HyperbolicFunction, f, symbol):
result = _solve_trig(f, symbol, domain)
elif isinstance(f, arg):
a = f.args[0]
result = solveset_real(a > 0, symbol)
elif f.is_Piecewise:
result = EmptySet()
expr_set_pairs = f.as_expr_set_pairs(domain)
for (expr, in_set) in expr_set_pairs:
if in_set.is_Relational:
in_set = in_set.as_set()
solns = solver(expr, symbol, in_set)
result += solns
elif isinstance(f, Eq):
result = solver(Add(f.lhs, - f.rhs, evaluate=False), symbol, domain)

elif f.is_Relational:
if not domain.is_subset(S.Reals):
raise NotImplementedError(filldedent('''
Inequalities in the complex domain are
not supported. Try the real domain by
setting domain=S.Reals'''))
try:
result = solve_univariate_inequality(
f, symbol, domain=domain, relational=False)
except NotImplementedError:
result = ConditionSet(symbol, f, domain)
return result
else:
lhs, rhs_s = inverter(f, 0, symbol)
if lhs == symbol:
# do some very minimal simplification since
# repeated inversion may have left the result
# in a state that other solvers (e.g. poly)
# would have simplified; this is done here
# rather than in the inverter since here it
# is only done once whereas there it would
# be repeated for each step of the inversion
if isinstance(rhs_s, FiniteSet):
rhs_s = FiniteSet(*[Mul(*
signsimp(i).as_content_primitive())
for i in rhs_s])
result = rhs_s

elif isinstance(rhs_s, FiniteSet):
for equation in [lhs - rhs for rhs in rhs_s]:
if equation == f:
if any(_has_rational_power(g, symbol)[0]
for g in equation.args) or _has_rational_power(
equation, symbol)[0]:
symbol,
solver)
elif equation.has(Abs):
result += _solve_abs(f, symbol, domain)
else:
result_rational = _solve_as_rational(equation, symbol, domain)
if isinstance(result_rational, ConditionSet):
# may be a transcendental type equation
result += _transolve(equation, symbol, domain)
else:
result += result_rational
else:
result += solver(equation, symbol)

elif rhs_s is not S.EmptySet:
result = ConditionSet(symbol, Eq(f, 0), domain)

if isinstance(result, ConditionSet):
num, den = f.as_numer_denom()
if den.has(symbol):
_result = _solveset(num, symbol, domain)
if not isinstance(_result, ConditionSet):
singularities = _solveset(den, symbol, domain)
result = _result - singularities

if _check:
if isinstance(result, ConditionSet):
# it wasn't solved or has enumerated all conditions
# -- leave it alone
return result

# whittle away all but the symbol-containing core
# to use this for testing

if isinstance(result, FiniteSet):
# check the result for invalid solutions
result = FiniteSet(*[s for s in result
if isinstance(s, RootOf)
or domain_check(fx, symbol, s)])

return result

def _term_factors(f):
"""
Iterator to get the factors of all terms present
in the given equation.

Parameters
==========
f : Expr
Equation that needs to be addressed

Returns
=======
Factors of all terms present in the equation.

Examples
========

>>> from sympy import symbols
>>> from sympy.solvers.solveset import _term_factors
>>> x = symbols('x')
>>> list(_term_factors(-2 - x**2 + x*(x + 1)))
[-2, -1, x**2, x, x + 1]
"""
yield mul_arg

[docs]def _solve_exponential(lhs, rhs, symbol, domain):
r"""
Helper function for solving (supported) exponential equations.

Exponential equations are the sum of (currently) at most
two terms with one or both of them having a power with a
symbol-dependent exponent.

For example

.. math:: 5^{2x + 3} - 5^{3x - 1}

.. math:: 4^{5 - 9x} - e^{2 - x}

Parameters
==========

lhs, rhs : Expr
The exponential equation to be solved, lhs = rhs

symbol : Symbol
The variable in which the equation is solved

domain : Set
A set over which the equation is solved.

Returns
=======

A set of solutions satisfying the given equation.
A ConditionSet if the equation is unsolvable or
if the assumptions are not properly defined, in that case
a different style of ConditionSet is returned having the
solution(s) of the equation with the desired assumptions.

Examples
========

>>> from sympy.solvers.solveset import _solve_exponential as solve_expo
>>> from sympy import symbols, S
>>> x = symbols('x', real=True)
>>> a, b = symbols('a b')
>>> solve_expo(2**x + 3**x - 5**x, 0, x, S.Reals)  # not solvable
ConditionSet(x, Eq(2**x + 3**x - 5**x, 0), Reals)
>>> solve_expo(a**x - b**x, 0, x, S.Reals)  # solvable but incorrect assumptions
ConditionSet(x, (a > 0) & (b > 0), {0})
>>> solve_expo(3**(2*x) - 2**(x + 3), 0, x, S.Reals)
{-3*log(2)/(-2*log(3) + log(2))}
>>> solve_expo(2**x - 4**x, 0, x, S.Reals)
{0}

* Proof of correctness of the method

The logarithm function is the inverse of the exponential function.
The defining relation between exponentiation and logarithm is:

.. math:: {log_b x} = y \enspace if \enspace b^y = x

Therefore if we are given an equation with exponent terms, we can
convert every term to its corresponding log form. This is achieved by
taking logarithms and expanding the equation using log identities
so that it can easily be handled by solveset.

For example:

.. math:: 3^{2x} = 2^{x + 3}

Taking log both sides will reduce the equation to

.. math:: (2x)*log(3) = (x + 3)*log(2)

This form can be easily handed by solveset.
"""
unsolved_result = ConditionSet(symbol, Eq(lhs - rhs), domain)
newlhs = powdenest(lhs)
if lhs != newlhs:
# it may also be advantageous to factor the new expr
return _solveset(factor(newlhs - rhs), symbol, domain)  # try again with _solveset

if not (isinstance(lhs, Add) and len(lhs.args) == 2):
# solving for the sum of more than two powers is possible
# but not yet implemented
return unsolved_result

if rhs != 0:
return unsolved_result

a, b = list(ordered(lhs.args))
a_term = a.as_independent(symbol)[1]
b_term = b.as_independent(symbol)[1]

a_base, a_exp = a_term.base, a_term.exp
b_base, b_exp = b_term.base, b_term.exp

from sympy.functions.elementary.complexes import im

if domain.is_subset(S.Reals):
conditions = And(
a_base > 0,
b_base > 0,
Eq(im(a_exp), 0),
Eq(im(b_exp), 0))
else:
conditions = And(
Ne(a_base, 0),
Ne(b_base, 0))

L, R = map(lambda i: expand_log(log(i), force=True), (a, -b))
solutions = _solveset(L - R, symbol, domain)

return ConditionSet(symbol, conditions, solutions)

[docs]def _is_exponential(f, symbol):
r"""
Return True if one or more terms contain symbol only in exponents,
else False.

Parameters
==========

f : Expr
The equation to be checked

symbol : Symbol
The variable in which the equation is checked

Examples
========

>>> from sympy import symbols, cos, exp
>>> from sympy.solvers.solveset import _is_exponential as check
>>> x, y = symbols('x y')
>>> check(y, y)
False
>>> check(x**y - 1, y)
True
>>> check(x**y*2**y - 1, y)
True
>>> check(exp(x + 3) + 3**x, x)
True
>>> check(cos(2**x), x)
False

* Philosophy behind the helper

The function extracts each term of the equation and checks if it is
of exponential form w.r.t symbol.
"""
rv = False
for expr_arg in _term_factors(f):
if symbol not in expr_arg.free_symbols:
continue
if (isinstance(expr_arg, Pow) and
symbol not in expr_arg.base.free_symbols or
isinstance(expr_arg, exp)):
rv = True  # symbol in exponent
else:
return False  # dependent on symbol in non-exponential way
return rv

[docs]def _solve_logarithm(lhs, rhs, symbol, domain):
r"""
Helper to solve logarithmic equations which are reducible
to a single instance of log.

Logarithmic equations are (currently) the equations that contains
log terms which can be reduced to a single log term or
a constant using various logarithmic identities.

For example:

.. math:: log(x) + log(x - 4)

can be reduced to:

.. math:: log(x(x - 4))

Parameters
==========

lhs, rhs : Expr
The logarithmic equation to be solved, lhs = rhs

symbol : Symbol
The variable in which the equation is solved

domain : Set
A set over which the equation is solved.

Returns
=======

A set of solutions satisfying the given equation.
A ConditionSet if the equation is unsolvable.

Examples
========

>>> from sympy import symbols, log, S
>>> from sympy.solvers.solveset import _solve_logarithm as solve_log
>>> x = symbols('x')
>>> f = log(x - 3) + log(x + 3)
>>> solve_log(f, 0, x, S.Reals)
{-sqrt(10), sqrt(10)}

* Proof of correctness

A logarithm is another way to write exponent and is defined by

.. math:: {log_b x} = y \enspace if \enspace b^y = x

When one side of the equation contains a single logarithm, the
equation can be solved by rewriting the equation as an equivalent
exponential equation as defined above. But if one side contains
more than one logarithm, we need to use the properties of logarithm
to condense it into a single logarithm.

Take for example

.. math:: log(2*x) - 15 = 0

contains single logarithm, therefore we can directly rewrite it to
exponential form as

.. math:: x = exp(15)/2

But if the equation has more than one logarithm as

.. math:: log(x - 3) + log(x + 3) = 0

we use log identities to convert into a reduced form

Using,

.. math:: log(a) + log(b) = log(a*b)

therefore the equation becomes,

.. math:: log((x - 3)*(x  + 3))

This equation contains one logarithm and can be solved by rewriting
to exponents.
"""
new_lhs = logcombine(lhs, force=True)
new_f = new_lhs - rhs

return _solveset(new_f, symbol, domain)

[docs]def _is_logarithmic(f, symbol):
r"""
Return True if the equation is in the form
a*log(f(x)) + b*log(g(x)) + ... + c else False.

Parameters
==========

f : Expr
The equation to be checked

symbol : Symbol
The variable in which the equation is checked

Returns
=======

True if the equation is logarithmic otherwise False.

Examples
========

>>> from sympy import symbols, tan, log
>>> from sympy.solvers.solveset import _is_logarithmic as check
>>> x, y = symbols('x y')
>>> check(log(x + 2) - log(x + 3), x)
True
>>> check(tan(log(2*x)), x)
False
>>> check(x*log(x), x)
False
>>> check(x + log(x), x)
False
>>> check(y + log(x), x)
True

* Philosophy behind the helper

The function extracts each term and checks whether it is
logarithmic w.r.t symbol.
"""
rv = False
saw_log = False
for term_arg in Mul.make_args(term):
if symbol not in term_arg.free_symbols:
continue
if isinstance(term_arg, log):
if saw_log:
saw_log = True
else:
return False  # dependent on symbol in non-log way
if saw_log:
rv = True
return rv

[docs]def _transolve(f, symbol, domain):
r"""
Function to solve transcendental equations. It is a helper to
solveset and should be used internally. \_transolve
currently supports the following class of equations:

- Exponential equations
- Logarithmic equations

Parameters
==========

f : Any transcendental equation that needs to be solved.
This needs to be an expression, which is assumed
to be equal to 0.

symbol : The variable for which the equation is solved.
This needs to be of class Symbol.

domain : A set over which the equation is solved.
This needs to be of class Set.

Returns
=======

Set
A set of values for symbol for which f is equal to
zero. An EmptySet is returned if f does not have solutions
in respective domain. A ConditionSet is returned as unsolved
object if algorithms to evaluate complete solution are not
yet implemented.

How to use \_transolve
========================

\_transolve should not be used as an independent function, because
it assumes that the equation (f) and the symbol comes from
solveset and might have undergone a few modification(s).
To use \_transolve as an independent function the equation (f)
and the symbol should be passed as they would have been by
solveset.

Examples
========

>>> from sympy.solvers.solveset import _transolve as transolve
>>> from sympy.solvers.solvers import _tsolve as tsolve
>>> from sympy import symbols, S, pprint
>>> x = symbols('x', real=True) # assumption added
>>> transolve(5**(x - 3) - 3**(2*x + 1), x, S.Reals)
{-(log(3) + 3*log(5))/(-log(5) + 2*log(3))}

How \_transolve works
=======================

\_transolve uses two types of helper functions to solve equations
of a particular class:

Identifying helpers: To determine whether a given equation
belongs to a certain class of equation or not. Returns either
True or False.

Solving helpers: Once an equation is identified, a corresponding
helper either solves the equation or returns a form of the equation
that solveset might better be able to handle.

* Philosophy behind the module

The purpose of \_transolve is to take equations which are not
already polynomial in their generator(s) and to either recast them
as such through a valid transformation or to solve them outright.
A pair of helper functions for each class of supported
transcendental functions are employed for this purpose. One
identifies the transcendental form of an equation and the other
either solves it or recasts it into a tractable form that can be
solved by  solveset.
For example, an equation in the form a*b**f(x) - c*d**g(x) = 0
can be transformed to
log(a) + f(x)*log(b) - log(c) - g(x)*log(d) = 0
(under certain assumptions) and this can be solved with solveset
if f(x) and g(x) are in polynomial form.

How \_transolve is better than \_tsolve
===========================================

1) Better output

\_transolve provides expressions in a more simplified form.

Consider a simple exponential equation

>>> f = 3**(2*x) - 2**(x + 3)
>>> pprint(transolve(f, x, S.Reals), use_unicode=False)
-3*log(2)
{------------------}
-2*log(3) + log(2)
>>> pprint(tsolve(f, x), use_unicode=False)
/   3     \
| --------|
| log(2/9)|
[-log\2         /]

2) Extensible

The API of \_transolve is designed such that it is easily
extensible, i.e. the code that solves a given class of
equations is encapsulated in a helper and not mixed in with
the code of \_transolve itself.

3) Modular

\_transolve is designed to be modular i.e, for every class of
equation a separate helper for identification and solving is
implemented. This makes it easy to change or modify any of the
method implemented directly in the helpers without interfering
with the actual structure of the API.

4) Faster Computation

Solving equation via \_transolve is much faster as compared to
\_tsolve. In solve, attempts are made computing every possibility
to get the solutions. This series of attempts makes solving a bit
slow. In \_transolve, computation begins only after a particular
type of equation is identified.

How to add new class of equations
=================================

Adding a new class of equation solver is a three-step procedure:

- Identify the type of the equations

Determine the type of the class of equations to which they belong:
it could be of Add, Pow, etc. types. Separate internal functions
are used for each type. Write identification and solving helpers
and use them from within the routine for the given type of equation
(after adding it, if necessary). Something like:

.. code-block:: python

....
if _is_exponential(lhs, x):
new_eq = _solve_exponential(lhs, rhs, x)
....
rhs, lhs = eq.as_independent(x)

- Define the identification helper.

- Define the solving helper.

Apart from this, a few other things needs to be taken care while

- Naming conventions:
Name of the *identification* *helper* should be as
\_is\_class where class will be the name or abbreviation
of the class of equation. The *solving* *helper* will be named as
\_solve\_class.
For example: for exponential equations it becomes
\_is\_exponential and \_solve\_expo.
- The identifying helpersshould take two input parameters,
the equation to be checked and the variable for which a solution
is being sought, while solving helpers would require an additional
domain parameter.
- Be sure to consider corner cases.
- Add tests for each helper.
implemented.
The documentation of the helpers should identify:

- the purpose of the helper,
- the method used to identify and solve the equation,
- a proof of correctness
- the return values of the helpers
"""

"""
Helper for _transolve to handle equations of
Add type, i.e. equations taking the form as
a*f(x) + b*g(x) + .... = c.
For example: 4**x + 8**x = 0
"""
result = ConditionSet(symbol, Eq(lhs - rhs, 0), domain)

# check if it is exponential type equation
if _is_exponential(lhs, symbol):
result = _solve_exponential(lhs, rhs, symbol, domain)
# check if it is logarithmic type equation
elif _is_logarithmic(lhs, symbol):
result = _solve_logarithm(lhs, rhs, symbol, domain)

return result

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

# invert_complex handles the call to the desired inverter based
# on the domain specified.
lhs, rhs_s = invert_complex(f, 0, symbol, domain)

if isinstance(rhs_s, FiniteSet):
assert (len(rhs_s.args)) == 1
rhs = rhs_s.args[0]

result = add_type(lhs, rhs, symbol, domain)
else:
result = rhs_s

return result

[docs]def solveset(f, symbol=None, domain=S.Complexes):
r"""Solves a given inequality or equation with set as output

Parameters
==========

f : Expr or a relational.
The target equation or inequality
symbol : Symbol
The variable for which the equation is solved
domain : Set
The domain over which the equation is solved

Returns
=======

Set
A set of values for symbol for which f is True or is equal to
zero. An EmptySet is returned if f is False or nonzero.
A ConditionSet is returned as unsolved object if algorithms
to evaluate complete solution are not yet implemented.

solveset claims to be complete in the solution set that it returns.

Raises
======

NotImplementedError
The algorithms to solve inequalities in complex domain  are
not yet implemented.
ValueError
The input is not valid.
RuntimeError
It is a bug, please report to the github issue tracker.

Notes
=====

Python interprets 0 and 1 as False and True, respectively, but
in this function they refer to solutions of an expression. So 0 and 1
return the Domain and EmptySet, respectively, while True and False
return the opposite (as they are assumed to be solutions of relational
expressions).

========

solveset_real: solver for real domain
solveset_complex: solver for complex domain

Examples
========

>>> from sympy import exp, sin, Symbol, pprint, S
>>> from sympy.solvers.solveset import solveset, solveset_real

* The default domain is complex. Not specifying a domain will lead
to the solving of the equation in the complex domain (and this
is not affected by the assumptions on the symbol):

>>> x = Symbol('x')
>>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
{2*n*I*pi | n in Integers}

>>> x = Symbol('x', real=True)
>>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
{2*n*I*pi | n in Integers}

* If you want to use solveset to solve the equation in the
real domain, provide a real domain. (Using solveset\_real
does this automatically.)

>>> R = S.Reals
>>> x = Symbol('x')
>>> solveset(exp(x) - 1, x, R)
{0}
>>> solveset_real(exp(x) - 1, x)
{0}

The solution is mostly unaffected by assumptions on the symbol,
but there may be some slight difference:

>>> pprint(solveset(sin(x)/x,x), use_unicode=False)
({2*n*pi | n in Integers} \ {0}) U ({2*n*pi + pi | n in Integers} \ {0})

>>> p = Symbol('p', positive=True)
>>> pprint(solveset(sin(p)/p, p), use_unicode=False)
{2*n*pi | n in Integers} U {2*n*pi + pi | n in Integers}

* Inequalities can be solved over the real domain only. Use of a complex

>>> solveset(exp(x) > 1, x, R)
Interval.open(0, oo)

"""
f = sympify(f)
symbol = sympify(symbol)

if f is S.true:
return domain

if f is S.false:
return S.EmptySet

if not isinstance(f, (Expr, Number)):
raise ValueError("%s is not a valid SymPy expression" % f)

if not isinstance(symbol, Expr) and  symbol is not None:
raise ValueError("%s is not a valid SymPy symbol" % symbol)

if not isinstance(domain, Set):
raise ValueError("%s is not a valid domain" %(domain))

free_symbols = f.free_symbols

if symbol is None and not free_symbols:
b = Eq(f, 0)
if b is S.true:
return domain
elif b is S.false:
return S.EmptySet
else:
raise NotImplementedError(filldedent('''
relationship between value and 0 is unknown: %s''' % b))

if symbol is None:
if len(free_symbols) == 1:
symbol = free_symbols.pop()
elif free_symbols:
raise ValueError(filldedent('''
The independent variable must be specified for a
multivariate equation.'''))
elif not isinstance(symbol, Symbol):
f, s, swap = recast_to_symbols([f], [symbol])
# the xreplace will be needed if a ConditionSet is returned
return solveset(f[0], s[0], domain).xreplace(swap)

if domain.is_subset(S.Reals):
if not symbol.is_real:
assumptions = symbol.assumptions0
assumptions['real'] = True
try:
r = Dummy('r', **assumptions)
return solveset(f.xreplace({symbol: r}), r, domain
).xreplace({r: symbol})
except InconsistentAssumptions:
pass
# Abs has its own handling method which avoids the
# rewriting property that the first piece of abs(x)
# is for x >= 0 and the 2nd piece for x < 0 -- solutions
# can look better if the 2nd condition is x <= 0. Since
# the solution is a set, duplication of results is not
# an issue, e.g. {y, -y} when y is 0 will be {0}
f = f.rewrite(Piecewise) # everything that's not an Abs
# everything *in* an Abs
e = e.func(e.args[0].rewrite(Piecewise))
f = f.xreplace({d: e})
f = piecewise_fold(f)

return _solveset(f, symbol, domain, _check=True)

[docs]def solveset_real(f, symbol):
return solveset(f, symbol, S.Reals)

[docs]def solveset_complex(f, symbol):
return solveset(f, symbol, S.Complexes)

def solvify(f, symbol, domain):
"""Solves an equation using solveset and returns the solution in accordance
with the solve output API.

Returns
=======

We classify the output based on the type of solution returned by solveset.

Solution    |    Output
----------------------------------------
FiniteSet   | list

ImageSet,   | list (if f is periodic)
Union       |

EmptySet    | empty list

Others      | None

Raises
======

NotImplementedError
A ConditionSet is the input.

Examples
========

>>> from sympy.solvers.solveset import solvify, solveset
>>> from sympy.abc import x
>>> from sympy import S, tan, sin, exp
>>> solvify(x**2 - 9, x, S.Reals)
[-3, 3]
>>> solvify(sin(x) - 1, x, S.Reals)
[pi/2]
>>> solvify(tan(x), x, S.Reals)
[0]
>>> solvify(exp(x) - 1, x, S.Complexes)

>>> solvify(exp(x) - 1, x, S.Reals)
[0]

"""
solution_set = solveset(f, symbol, domain)
result = None
if solution_set is S.EmptySet:
result = []

elif isinstance(solution_set, ConditionSet):
raise NotImplementedError('solveset is unable to solve this equation.')

elif isinstance(solution_set, FiniteSet):
result = list(solution_set)

else:
period = periodicity(f, symbol)
if period is not None:
solutions = S.EmptySet
iter_solutions = ()
if isinstance(solution_set, ImageSet):
iter_solutions = (solution_set,)
elif isinstance(solution_set, Union):
if all(isinstance(i, ImageSet) for i in solution_set.args):
iter_solutions = solution_set.args

for solution in iter_solutions:
solutions += solution.intersect(Interval(0, period, False, True))

if isinstance(solutions, FiniteSet):
result = list(solutions)

else:
solution = solution_set.intersect(domain)
if isinstance(solution, FiniteSet):
result += solution

return result

###############################################################################
################################ LINSOLVE #####################################
###############################################################################

[docs]def linear_eq_to_matrix(equations, *symbols):
r"""
Converts a given System of Equations into Matrix form.
Here equations must be a linear system of equations in
symbols. The order of symbols in input symbols will
determine the order of coefficients in the returned
Matrix.

The Matrix form corresponds to the augmented matrix form.
For example:

.. math:: 4x + 2y + 3z  = 1
.. math:: 3x +  y +  z  = -6
.. math:: 2x + 4y + 9z  = 2

This system would return A & b as given below:

::

[ 4  2  3 ]          [ 1 ]
A = [ 3  1  1 ]   b  =   [-6 ]
[ 2  4  9 ]          [ 2 ]

Examples
========

>>> from sympy import linear_eq_to_matrix, symbols
>>> x, y, z = symbols('x, y, z')
>>> eqns = [x + 2*y + 3*z - 1, 3*x + y + z + 6, 2*x + 4*y + 9*z - 2]
>>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
>>> A
Matrix([
[1, 2, 3],
[3, 1, 1],
[2, 4, 9]])
>>> b
Matrix([
[ 1],
[-6],
[ 2]])
>>> eqns = [x + z - 1, y + z, x - y]
>>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
>>> A
Matrix([
[1,  0, 1],
[0,  1, 1],
[1, -1, 0]])
>>> b
Matrix([
[1],
[0],
[0]])

* Symbolic coefficients are also supported

>>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
>>> eqns = [a*x + b*y - c, d*x + e*y - f]
>>> A, B = linear_eq_to_matrix(eqns, x, y)
>>> A
Matrix([
[a, b],
[d, e]])
>>> B
Matrix([
[c],
[f]])

"""

if not symbols:
raise ValueError('Symbols must be given, for which coefficients \
are to be found.')

if hasattr(symbols[0], '__iter__'):
symbols = symbols[0]

M = Matrix([symbols])
# initialize Matrix with symbols + 1 columns
M = M.col_insert(len(symbols), Matrix([1]))
row_no = 1

for equation in equations:
f = sympify(equation)
if isinstance(f, Equality):
f = f.lhs - f.rhs

# Extract coeff of symbols
coeff_list = []
for symbol in symbols:
coeff_list.append(f.coeff(symbol))

# append constant term (term free from symbols)

# insert equations coeff's into rows
M = M.row_insert(row_no, Matrix([coeff_list]))
row_no += 1

# delete the initialized (Ist) trivial row
M.row_del(0)
A, b = M[:, :-1], M[:, -1:]
return A, b

[docs]def linsolve(system, *symbols):
r"""
Solve system of N linear equations with M variables, which
means both under - and overdetermined systems are supported.
The possible number of solutions is zero, one or infinite.
Zero solutions throws a ValueError, where as infinite
solutions are represented parametrically in terms of given
symbols. For unique solution a FiniteSet of ordered tuple
is returned.

All Standard input formats are supported:
For the given set of Equations, the respective input types
are given below:

.. math:: 3x + 2y -   z = 1
.. math:: 2x - 2y + 4z = -2
.. math:: 2x -   y + 2z = 0

* Augmented Matrix Form, system given below:

::

[3   2  -1  1]
system = [2  -2   4 -2]
[2  -1   2  0]

* List Of Equations Form

system  =  [3x + 2y - z - 1, 2x - 2y + 4z + 2, 2x - y + 2z]

* Input A & b Matrix Form (from Ax = b) are given as below:

::

[3   2  -1 ]         [  1 ]
A = [2  -2   4 ]    b =  [ -2 ]
[2  -1   2 ]         [  0 ]

system = (A, b)

Symbols can always be passed but are actually only needed
when 1) a system of equations is being passed and 2) the
system is passed as an underdetermined matrix and one wants
to control the name of the free variables in the result.
An error is raised if no symbols are used for case 1, but if
no symbols are provided for case 2, internally generated symbols
will be provided. When providing symbols for case 2, there should
be at least as many symbols are there are columns in matrix A.

The algorithm used here is Gauss-Jordan elimination, which
results, after elimination, in a row echelon form matrix.

Returns
=======

A FiniteSet containing an ordered tuple of values for the
unknowns for which the system has a solution. (Wrapping
the tuple in FiniteSet is used to maintain a consistent
output format throughout solveset.)

Returns EmptySet(), if the linear system is inconsistent.

Raises
======

ValueError
The input is not valid.
The symbols are not given.

Examples
========

>>> from sympy import Matrix, S, linsolve, symbols
>>> x, y, z = symbols("x, y, z")
>>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
>>> b = Matrix([3, 6, 9])
>>> A
Matrix([
[1, 2,  3],
[4, 5,  6],
[7, 8, 10]])
>>> b
Matrix([
[3],
[6],
[9]])
>>> linsolve((A, b), [x, y, z])
{(-1, 2, 0)}

* Parametric Solution: In case the system is underdetermined, the
function will return a parametric solution in terms of the given
symbols. Those that are free will be returned unchanged. e.g. in
the system below, z is returned as the solution for variable z;
it can take on any value.

>>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> b = Matrix([3, 6, 9])
>>> linsolve((A, b), x, y, z)
{(z - 1, -2*z + 2, z)}

If no symbols are given, internally generated symbols will be used.
The tau0 in the 3rd position indicates (as before) that the 3rd
variable -- whatever it's named -- can take on any value:

>>> linsolve((A, b))
{(tau0 - 1, -2*tau0 + 2, tau0)}

* List of Equations as input

>>> Eqns = [3*x + 2*y - z - 1, 2*x - 2*y + 4*z + 2, - x + y/2 - z]
>>> linsolve(Eqns, x, y, z)
{(1, -2, -2)}

* Augmented Matrix as input

>>> aug = Matrix([[2, 1, 3, 1], [2, 6, 8, 3], [6, 8, 18, 5]])
>>> aug
Matrix([
[2, 1,  3, 1],
[2, 6,  8, 3],
[6, 8, 18, 5]])
>>> linsolve(aug, x, y, z)
{(3/10, 2/5, 0)}

* Solve for symbolic coefficients

>>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
>>> eqns = [a*x + b*y - c, d*x + e*y - f]
>>> linsolve(eqns, x, y)
{((-b*f + c*e)/(a*e - b*d), (a*f - c*d)/(a*e - b*d))}

* A degenerate system returns solution as set of given
symbols.

>>> system = Matrix(([0, 0, 0], [0, 0, 0], [0, 0, 0]))
>>> linsolve(system, x, y)
{(x, y)}

* For an empty system linsolve returns empty set

>>> linsolve([], x)
EmptySet()

"""
if not system:
return S.EmptySet

# If second argument is an iterable
if symbols and hasattr(symbols[0], '__iter__'):
symbols = symbols[0]
sym_gen = isinstance(symbols, GeneratorType)

swap = {}
b = None  # if we don't get b the input was bad
syms_needed_msg = None

# unpack system

if hasattr(system, '__iter__'):

# 1). (A, b)
if len(system) == 2 and isinstance(system[0], Matrix):
A, b = system

# 2). (eq1, eq2, ...)
if not isinstance(system[0], Matrix):
if sym_gen or not symbols:
raise ValueError(filldedent('''
When passing a system of equations, the explicit
symbols for which a solution is being sought must
be given as a sequence, too.
'''))
system = list(system)
for i, eq in enumerate(system):
try:
# since we are checking it, we might as well take the
# expanded expr that it will give
system[i] = Poly(eq, symbols).as_expr()
if any (degree(eq, sym) > 1 for sym in symbols):
raise PolynomialError
except PolynomialError:
raise ValueError(filldedent('''
%s contains non-linear terms in the
variables to be evaluated
''') % eq)
system, symbols, swap = recast_to_symbols(system, symbols)
A, b = linear_eq_to_matrix(system, symbols)
syms_needed_msg = 'free symbols in the equations provided'

elif isinstance(system, Matrix) and not (
symbols and not isinstance(symbols, GeneratorType) and
isinstance(symbols[0], Matrix)):
# 3). A augmented with b
A, b = system[:, :-1], system[:, -1:]

if b is None:
raise ValueError("Invalid arguments")

syms_needed_msg  = syms_needed_msg or 'columns of A'

if sym_gen:
symbols = [next(symbols) for i in range(A.cols)]
if any(set(symbols) & (A.free_symbols | b.free_symbols)):
raise ValueError(filldedent('''
At least one of the symbols provided
already appears in the system to be solved.
One way to avoid this is to use Dummy symbols in
the generator, e.g. numbered_symbols('%s', cls=Dummy)
''' % symbols[0].name.rstrip('1234567890')))

try:
solution, params, free_syms = A.gauss_jordan_solve(b, freevar=True)
except ValueError:
# No solution
return S.EmptySet

# Replace free parameters with free symbols
if params:
if not symbols:
symbols = [_ for _ in params]
# re-use the parameters but put them in order
# params       [x, y, z]
# free_symbols [2, 0, 4]
# idx          [1, 0, 2]
idx = list(zip(*sorted(zip(free_syms, range(len(free_syms))))))[1]
# simultaneous replacements {y: x, x: y, z: z}
replace_dict = dict(zip(symbols, [symbols[i] for i in idx]))
elif len(symbols) >= A.cols:
replace_dict = {v: symbols[free_syms[k]] for k, v in enumerate(params)}
else:
raise IndexError(filldedent('''
the number of symbols passed should have a length
equal to the number of %s.
''' % syms_needed_msg))
solution = [sol.xreplace(replace_dict) for sol in solution]

solution = [simplify(sol).xreplace(swap) for sol in solution]
return FiniteSet(tuple(solution))

##############################################################################
# ------------------------------nonlinsolve ---------------------------------#
##############################################################################

def _return_conditionset(eqs, symbols):
# return conditionset
condition_set = ConditionSet(
Tuple(*symbols),
FiniteSet(*eqs),
S.Complexes)
return condition_set

def substitution(system, symbols, result=[{}], known_symbols=[],
exclude=[], all_symbols=None):
r"""
Solves the system using substitution method. It is used in
nonlinsolve. This will be called from nonlinsolve when any
equation(s) is non polynomial equation.

Parameters
==========

system : list of equations
The target system of equations
symbols : list of symbols to be solved.
The variable(s) for which the system is solved
known_symbols : list of solved symbols
Values are known for these variable(s)
result : An empty list or list of dict
If No symbol values is known then empty list otherwise
symbol as keys and corresponding value in dict.
exclude : Set of expression.
Mostly denominator expression(s) of the equations of the system.
Final solution should not satisfy these expressions.
all_symbols : known_symbols + symbols(unsolved).

Returns
=======

A FiniteSet of ordered tuple of values of all_symbols for which the
system has solution. Order of values in the tuple is same as symbols
present in the parameter all_symbols. If parameter all_symbols is None
then same as symbols present in the parameter symbols.

Please note that general FiniteSet is unordered, the solution returned
here is not simply a FiniteSet of solutions, rather it is a FiniteSet of
ordered tuple, i.e. the first & only argument to FiniteSet is a tuple of
solutions, which is ordered, & hence the returned solution is ordered.

Also note that solution could also have been returned as an ordered tuple,
FiniteSet is just a wrapper {} around the tuple. It has no other
significance except for the fact it is just used to maintain a consistent
output format throughout the solveset.

Raises
======

ValueError
The input is not valid.
The symbols are not given.
AttributeError
The input symbols are not Symbol type.

Examples
========

>>> from sympy.core.symbol import symbols
>>> x, y = symbols('x, y', real=True)
>>> from sympy.solvers.solveset import substitution
>>> substitution([x + y], [x], [{y: 1}], [y], set([]), [x, y])
{(-1, 1)}

* when you want soln should not satisfy eq x + 1 = 0

>>> substitution([x + y], [x], [{y: 1}], [y], set([x + 1]), [y, x])
EmptySet()
>>> substitution([x + y], [x], [{y: 1}], [y], set([x - 1]), [y, x])
{(1, -1)}
>>> substitution([x + y - 1, y - x**2 + 5], [x, y])
{(-3, 4), (2, -1)}

* Returns both real and complex solution

>>> x, y, z = symbols('x, y, z')
>>> from sympy import exp, sin
>>> substitution([exp(x) - sin(y), y**2 - 4], [x, y])
{(log(sin(2)), 2), (ImageSet(Lambda(_n, I*(2*_n*pi + pi) +
log(sin(2))), Integers), -2), (ImageSet(Lambda(_n, 2*_n*I*pi +
Mod(log(sin(2)), 2*I*pi)), Integers), 2)}

>>> eqs = [z**2 + exp(2*x) - sin(y), -3 + exp(-y)]
>>> substitution(eqs, [y, z])
{(-log(3), -sqrt(-exp(2*x) - sin(log(3)))),
(-log(3), sqrt(-exp(2*x) - sin(log(3)))),
(ImageSet(Lambda(_n, 2*_n*I*pi + Mod(-log(3), 2*I*pi)), Integers),
ImageSet(Lambda(_n, -sqrt(-exp(2*x) + sin(2*_n*I*pi +
Mod(-log(3), 2*I*pi)))), Integers)),
(ImageSet(Lambda(_n, 2*_n*I*pi + Mod(-log(3), 2*I*pi)), Integers),
ImageSet(Lambda(_n, sqrt(-exp(2*x) + sin(2*_n*I*pi +
Mod(-log(3), 2*I*pi)))), Integers))}

"""

from sympy import Complement
from sympy.core.compatibility import is_sequence

if not system:
return S.EmptySet

if not symbols:
msg = ('Symbols must be given, for which solution of the '
'system is to be found.')
raise ValueError(filldedent(msg))

if not is_sequence(symbols):
msg = ('symbols should be given as a sequence, e.g. a list.'
'Not type %s: %s')
raise TypeError(filldedent(msg % (type(symbols), symbols)))

try:
sym = symbols[0].is_Symbol
except AttributeError:
sym = False

if not sym:
msg = ('Iterable of symbols must be given as '
'second argument, not type %s: %s')
raise ValueError(filldedent(msg % (type(symbols[0]), symbols[0])))

# By default all_symbols will be same as symbols
if all_symbols is None:
all_symbols = symbols

old_result = result
# storing complements and intersection for particular symbol
complements = {}
intersections = {}

# when total_solveset_call is equals to total_conditionset
# means solvest fail to solve all the eq.
total_conditionset = -1
total_solveset_call = -1

def _unsolved_syms(eq, sort=False):
"""Returns the unsolved symbol present
in the equation eq.
"""
free = eq.free_symbols
unsolved = (free - set(known_symbols)) & set(all_symbols)
if sort:
unsolved = list(unsolved)
unsolved.sort(key=default_sort_key)
return unsolved
# end of _unsolved_syms()

# sort such that equation with the fewest potential symbols is first.
# means eq with less number of variable first in the list.
eqs_in_better_order = list(
ordered(system, lambda _: len(_unsolved_syms(_))))

# If solveset have returned some intersection/complement
# for any symbol. It will be added in final solution.
final_result = []
for res in result:
res_copy = res
for key_res, value_res in res.items():
# Intersection/complement is in Interval or Set.
intersection_true = flags.get('Intersection', True)
complements_true = flags.get('Complement', True)
for key_sym, value_sym in sym_set.items():
if key_sym == key_res:
if intersection_true:
# testcase is not added for this line(intersection)
new_value = \
Intersection(FiniteSet(value_res), value_sym)
if new_value is not S.EmptySet:
res_copy[key_res] = new_value
if complements_true:
new_value = \
Complement(FiniteSet(value_res), value_sym)
if new_value is not S.EmptySet:
res_copy[key_res] = new_value
final_result.append(res_copy)
return final_result

def _extract_main_soln(sol, soln_imageset):
"""separate the Complements, Intersections, ImageSet lambda expr
and it's base_set.
"""
# if there is union, then need to check
# Complement, Intersection, Imageset.
# Order should not be changed.
if isinstance(sol, Complement):
# extract solution and complement
complements[sym] = sol.args[1]
sol = sol.args[0]
# complement will be added at the end
# using add_intersection_complement method
if isinstance(sol, Intersection):
# Interval/Set will be at 0th index always
if sol.args[0] != Interval(-oo, oo):
# sometimes solveset returns soln
# with intersection S.Reals, to confirm that
# soln is in domain=S.Reals or not. We don't consider
# that intersection.
intersections[sym] = sol.args[0]
sol = sol.args[1]
# after intersection and complement Imageset should
# be checked.
if isinstance(sol, ImageSet):
soln_imagest = sol
expr2 = sol.lamda.expr
sol = FiniteSet(expr2)
soln_imageset[expr2] = soln_imagest

# if there is union of Imageset or other in soln.
# no testcase is written for this if block
if isinstance(sol, Union):
sol_args = sol.args
sol = S.EmptySet
# We need in sequence so append finteset elements
# and then imageset or other.
for sol_arg2 in sol_args:
if isinstance(sol_arg2, FiniteSet):
sol += sol_arg2
else:
# ImageSet, Intersection, complement then
# append them directly
sol += FiniteSet(sol_arg2)

if not isinstance(sol, FiniteSet):
sol = FiniteSet(sol)
return sol, soln_imageset
# end of def _extract_main_soln()

# helper function for _append_new_soln
def _check_exclude(rnew, imgset_yes):
rnew_ = rnew
if imgset_yes:
# replace all dummy variables (Imageset lambda variables)
# with zero before checksol. Considering fundamental soln
# for checksol.
rnew_copy = rnew.copy()
dummy_n = imgset_yes[0]
for key_res, value_res in rnew_copy.items():
rnew_copy[key_res] = value_res.subs(dummy_n, 0)
rnew_ = rnew_copy
# satisfy_exclude == true if it satisfies the expr of exclude list.
try:
# something like : Mod(-log(3), 2*I*pi) can't be
# simplified right now, so checksol returns TypeError.
# when this issue is fixed this try block should be
# removed. Mod(-log(3), 2*I*pi) == -log(3)
satisfy_exclude = any(
checksol(d, rnew_) for d in exclude)
except TypeError:
satisfy_exclude = None
return satisfy_exclude
# end of def _check_exclude()

# helper function for _append_new_soln
def _restore_imgset(rnew, original_imageset, newresult):
restore_sym = set(rnew.keys()) & \
set(original_imageset.keys())
for key_sym in restore_sym:
img = original_imageset[key_sym]
rnew[key_sym] = img
if rnew not in newresult:
newresult.append(rnew)
# end of def _restore_imgset()

def _append_eq(eq, result, res, delete_soln, n=None):
u = Dummy('u')
if n:
eq = eq.subs(n, 0)
satisfy = checksol(u, u, eq, minimal=True)
if satisfy is False:
delete_soln = True
res = {}
else:
result.append(res)
return result, res, delete_soln

def _append_new_soln(rnew, sym, sol, imgset_yes, soln_imageset,
original_imageset, newresult, eq=None):
"""If rnew (A dict <symbol: soln>) contains valid soln
append it to newresult list.
imgset_yes is (base, dummy_var) if there was imageset in previously
calculated result(otherwise empty tuple). original_imageset is dict
of imageset expr and imageset from this result.
soln_imageset dict of imageset expr and imageset of new soln.
"""
satisfy_exclude = _check_exclude(rnew, imgset_yes)
delete_soln = False
# soln should not satisfy expr present in exclude list.
if not satisfy_exclude:
local_n = None
# if it is imageset
if imgset_yes:
local_n = imgset_yes[0]
base = imgset_yes[1]
if sym and sol:
# when sym and sol is None means no new
# soln. In that case we will append rnew directly after
# substituting original imagesets in rnew values if present
# (second last line of this function using _restore_imgset)
dummy_list = list(sol.atoms(Dummy))
# use one dummy n which is in
# previous imageset
local_n_list = [
local_n for i in range(
0, len(dummy_list))]

dummy_zip = zip(dummy_list, local_n_list)
lam = Lambda(local_n, sol.subs(dummy_zip))
rnew[sym] = ImageSet(lam, base)
if eq is not None:
newresult, rnew, delete_soln = _append_eq(
eq, newresult, rnew, delete_soln, local_n)
elif eq is not None:
newresult, rnew, delete_soln = _append_eq(
eq, newresult, rnew, delete_soln)
elif soln_imageset:
rnew[sym] = soln_imageset[sol]
# restore original imageset
_restore_imgset(rnew, original_imageset, newresult)
else:
newresult.append(rnew)
elif satisfy_exclude:
delete_soln = True
rnew = {}
_restore_imgset(rnew, original_imageset, newresult)
return newresult, delete_soln
# end of def _append_new_soln()

def _new_order_result(result, eq):
# separate first, second priority. res that makes eq value equals
# to zero, should be used first then other result(second priority).
# If it is not done then we may miss some soln.
first_priority = []
second_priority = []
for res in result:
if not any(isinstance(val, ImageSet) for val in res.values()):
if eq.subs(res) == 0:
first_priority.append(res)
else:
second_priority.append(res)
if first_priority or second_priority:
return first_priority + second_priority
return result

def _solve_using_known_values(result, solver):
"""Solves the system using already known solution
(result contains the dict <symbol: value>).
solver is solveset_complex or solveset_real.
"""
# stores imageset <expr: imageset(Lambda(n, expr), base)>.
soln_imageset = {}
total_solvest_call = 0
total_conditionst = 0

# sort such that equation with the fewest potential symbols is first.
# means eq with less variable first
for index, eq in enumerate(eqs_in_better_order):
newresult = []
original_imageset = {}
# if imageset expr is used to solve other symbol
imgset_yes = False
result = _new_order_result(result, eq)
for res in result:
got_symbol = set()  # symbols solved in one iteration
if soln_imageset:
# find the imageset and use its expr.
for key_res, value_res in res.items():
if isinstance(value_res, ImageSet):
res[key_res] = value_res.lamda.expr
original_imageset[key_res] = value_res
dummy_n = value_res.lamda.expr.atoms(Dummy).pop()
base = value_res.base_set
imgset_yes = (dummy_n, base)
# update eq with everything that is known so far
eq2 = eq.subs(res)
unsolved_syms = _unsolved_syms(eq2, sort=True)
if not unsolved_syms:
if res:
newresult, delete_res = _append_new_soln(
res, None, None, imgset_yes, soln_imageset,
original_imageset, newresult, eq2)
if delete_res:
# delete_res is true, means substituting res in
# eq2 doesn't return zero or deleting the res
# (a soln) since it staisfies expr of exclude
# list.
result.remove(res)
continue  # skip as it's independent of desired symbols
depen = eq2.as_independent(unsolved_syms)[0]
if depen.has(Abs) and solver == solveset_complex:
# Absolute values cannot be inverted in the
# complex domain
continue
soln_imageset = {}
for sym in unsolved_syms:
not_solvable = False
try:
soln = solver(eq2, sym)
total_solvest_call += 1
soln_new = S.EmptySet
if isinstance(soln, Complement):
# separate solution and complement
complements[sym] = soln.args[1]
soln = soln.args[0]
# complement will be added at the end
if isinstance(soln, Intersection):
# Interval will be at 0th index always
if soln.args[0] != Interval(-oo, oo):
# sometimes solveset returns soln
# with intersection S.Reals, to confirm that
# soln is in domain=S.Reals
intersections[sym] = soln.args[0]
soln_new += soln.args[1]
soln = soln_new if soln_new else soln
if index > 0 and solver == solveset_real:
# one symbol's real soln , another symbol may have
# corresponding complex soln.
if not isinstance(soln, (ImageSet, ConditionSet)):
soln += solveset_complex(eq2, sym)
except NotImplementedError:
# If sovleset is not able to solve equation eq2. Next
# time we may get soln using next equation eq2
continue
if isinstance(soln, ConditionSet):
soln = S.EmptySet
# don't do continue we may get soln
# in terms of other symbol(s)
not_solvable = True
total_conditionst += 1

if soln is not S.EmptySet:
soln, soln_imageset = _extract_main_soln(
soln, soln_imageset)

for sol in soln:
# sol is not a Union since we checked it
# before this loop
sol, soln_imageset = _extract_main_soln(
sol, soln_imageset)
sol = set(sol).pop()
free = sol.free_symbols
if got_symbol and any([
ss in free for ss in got_symbol
]):
# sol depends on previously solved symbols
# then continue
continue
rnew = res.copy()
# put each solution in res and append the new  result
# in the new result list (solution for symbol s)
# along with old results.
for k, v in res.items():
if isinstance(v, Expr):
# if any unsolved symbol is present
# Then subs known value
rnew[k] = v.subs(sym, sol)
# and add this new solution
if soln_imageset:
# replace all lambda variables with 0.
imgst = soln_imageset[sol]
rnew[sym] = imgst.lamda(
*[0 for i in range(0, len(
imgst.lamda.variables))])
else:
rnew[sym] = sol
newresult, delete_res = _append_new_soln(
rnew, sym, sol, imgset_yes, soln_imageset,
original_imageset, newresult)
if delete_res:
# deleting the res (a soln) since it staisfies
# eq of exclude list
result.remove(res)
# solution got for sym
if not not_solvable:
# next time use this new soln
if newresult:
result = newresult
return result, total_solvest_call, total_conditionst
# end def _solve_using_know_values()

new_result_real, solve_call1, cnd_call1 = _solve_using_known_values(
old_result, solveset_real)
new_result_complex, solve_call2, cnd_call2 = _solve_using_known_values(
old_result, solveset_complex)

# when total_solveset_call is equals to total_conditionset
# means solvest fails to solve all the eq.
# return conditionset in this case
total_conditionset += (cnd_call1 + cnd_call2)
total_solveset_call += (solve_call1 + solve_call2)

if total_conditionset == total_solveset_call and total_solveset_call != -1:
return _return_conditionset(eqs_in_better_order, all_symbols)

# overall result
result = new_result_real + new_result_complex

result_all_variables = []
result_infinite = []
for res in result:
if not res:
# means {None : None}
continue
# If length < len(all_symbols) means infinite soln.
# Some or all the soln is dependent on 1 symbol.
# eg. {x: y+2} then final soln {x: y+2, y: y}
if len(res) < len(all_symbols):
solved_symbols = res.keys()
unsolved = list(filter(
lambda x: x not in solved_symbols, all_symbols))
for unsolved_sym in unsolved:
res[unsolved_sym] = unsolved_sym
result_infinite.append(res)
if res not in result_all_variables:
result_all_variables.append(res)

if result_infinite:
# we have general soln
# eg : [{x: -1, y : 1}, {x : -y , y: y}] then
# return [{x : -y, y : y}]
result_all_variables = result_infinite
if intersections and complements:
# no testcase is added for this block
result_all_variables, intersections,
Intersection=True, Complement=True)
elif intersections:
result_all_variables, intersections, Intersection=True)
elif complements:
result_all_variables, complements, Complement=True)

# convert to ordered tuple
result = S.EmptySet
for r in result_all_variables:
temp = [r[symb] for symb in all_symbols]
result += FiniteSet(tuple(temp))
return result
# end of def substitution()

def _solveset_work(system, symbols):
soln = solveset(system[0], symbols[0])
if isinstance(soln, FiniteSet):
_soln = FiniteSet(*[tuple((s,)) for s in soln])
return _soln
else:
return FiniteSet(tuple(FiniteSet(soln)))

def _handle_positive_dimensional(polys, symbols, denominators):
from sympy.polys.polytools import groebner
# substitution method where new system is groebner basis of the system
_symbols = list(symbols)
_symbols.sort(key=default_sort_key)
basis = groebner(polys, _symbols, polys=True)
new_system = []
for poly_eq in basis:
new_system.append(poly_eq.as_expr())
result = [{}]
result = substitution(
new_system, symbols, result, [],
denominators)
return result
# end of def _handle_positive_dimensional()

def _handle_zero_dimensional(polys, symbols, system):
# solve 0 dimensional poly system using solve_poly_system
result = solve_poly_system(polys, *symbols)
# May be some extra soln is added because
# we used unrad in _separate_poly_nonpoly, so
# need to check and remove if it is not a soln.
result_update = S.EmptySet
for res in result:
dict_sym_value = dict(list(zip(symbols, res)))
if all(checksol(eq, dict_sym_value) for eq in system):
result_update += FiniteSet(res)
return result_update
# end of def _handle_zero_dimensional()

def _separate_poly_nonpoly(system, symbols):
polys = []
polys_expr = []
nonpolys = []
denominators = set()
poly = None
for eq in system:
# Store denom expression if it contains symbol
denominators.update(_simple_dens(eq, symbols))
# try to remove sqrt and rational power
if not cov:
if isinstance(eq, Expr):
eq = eq.as_numer_denom()[0]
poly = eq.as_poly(*symbols, extension=True)
elif simplify(eq).is_number:
continue
if poly is not None:
polys.append(poly)
polys_expr.append(poly.as_expr())
else:
nonpolys.append(eq)
return polys, polys_expr, nonpolys, denominators
# end of def _separate_poly_nonpoly()

[docs]def nonlinsolve(system, *symbols):
r"""
Solve system of N non linear equations with M variables, which means both
under and overdetermined systems are supported. Positive dimensional
system is also supported (A system with infinitely many solutions is said
to be positive-dimensional). In Positive dimensional system solution will
be dependent on at least one symbol. Returns both real solution
and complex solution(If system have). The possible number of solutions
is zero, one or infinite.

Parameters
==========

system : list of equations
The target system of equations
symbols : list of Symbols
symbols should be given as a sequence eg. list

Returns
=======

A FiniteSet of ordered tuple of values of symbols for which the system
has solution. Order of values in the tuple is same as symbols present in
the parameter symbols.

Please note that general FiniteSet is unordered, the solution returned
here is not simply a FiniteSet of solutions, rather it is a FiniteSet of
ordered tuple, i.e. the first & only argument to FiniteSet is a tuple of
solutions, which is ordered, & hence the returned solution is ordered.

Also note that solution could also have been returned as an ordered tuple,
FiniteSet is just a wrapper {} around the tuple. It has no other
significance except for the fact it is just used to maintain a consistent
output format throughout the solveset.

For the given set of Equations, the respective input types
are given below:

.. math:: x*y - 1 = 0
.. math:: 4*x**2 + y**2 - 5 = 0

system  = [x*y - 1, 4*x**2 + y**2 - 5]
symbols = [x, y]

Raises
======

ValueError
The input is not valid.
The symbols are not given.
AttributeError
The input symbols are not Symbol type.

Examples
========

>>> from sympy.core.symbol import symbols
>>> from sympy.solvers.solveset import nonlinsolve
>>> x, y, z = symbols('x, y, z', real=True)
>>> nonlinsolve([x*y - 1, 4*x**2 + y**2 - 5], [x, y])
{(-1, -1), (-1/2, -2), (1/2, 2), (1, 1)}

1. Positive dimensional system and complements:

>>> from sympy import pprint
>>> from sympy.polys.polytools import is_zero_dimensional
>>> a, b, c, d = symbols('a, b, c, d', real=True)
>>> eq1 =  a + b + c + d
>>> eq2 = a*b + b*c + c*d + d*a
>>> eq3 = a*b*c + b*c*d + c*d*a + d*a*b
>>> eq4 = a*b*c*d - 1
>>> system = [eq1, eq2, eq3, eq4]
>>> is_zero_dimensional(system)
False
>>> pprint(nonlinsolve(system, [a, b, c, d]), use_unicode=False)
-1       1               1      -1
{(---, -d, -, {d} \ {0}), (-, -d, ---, {d} \ {0})}
d       d               d       d
>>> nonlinsolve([(x+y)**2 - 4, x + y - 2], [x, y])
{(-y + 2, y)}

2. If some of the equations are non polynomial equation then nonlinsolve
will call substitution function and returns real and complex solutions,
if present.

>>> from sympy import exp, sin
>>> nonlinsolve([exp(x) - sin(y), y**2 - 4], [x, y])
{(log(sin(2)), 2), (ImageSet(Lambda(_n, I*(2*_n*pi + pi) +
log(sin(2))), Integers), -2), (ImageSet(Lambda(_n, 2*_n*I*pi +
Mod(log(sin(2)), 2*I*pi)), Integers), 2)}

3. If system is Non linear polynomial zero dimensional then it returns
both solution (real and complex solutions, if present using
solve_poly_system):

>>> from sympy import sqrt
>>> nonlinsolve([x**2 - 2*y**2 -2, x*y - 2], [x, y])
{(-2, -1), (2, 1), (-sqrt(2)*I, sqrt(2)*I), (sqrt(2)*I, -sqrt(2)*I)}

4. nonlinsolve can solve some linear(zero or positive dimensional)
system (because it is using groebner function to get the
groebner basis and then substitution function basis as the new system).
But it is not recommended to solve linear system using nonlinsolve,
because linsolve is better for all kind of linear system.

>>> nonlinsolve([x + 2*y -z - 3, x - y - 4*z + 9 , y + z - 4], [x, y, z])
{(3*z - 5, -z + 4, z)}

5. System having polynomial equations and only real solution is present
(will be solved using solve_poly_system):

>>> e1 = sqrt(x**2 + y**2) - 10
>>> e2 = sqrt(y**2 + (-x + 10)**2) - 3
>>> nonlinsolve((e1, e2), (x, y))
{(191/20, -3*sqrt(391)/20), (191/20, 3*sqrt(391)/20)}
>>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [x, y])
{(1, 2), (1 + sqrt(5), -sqrt(5) + 2), (-sqrt(5) + 1, 2 + sqrt(5))}
>>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [y, x])
{(2, 1), (2 + sqrt(5), -sqrt(5) + 1), (-sqrt(5) + 2, 1 + sqrt(5))}

6. It is better to use symbols instead of Trigonometric Function or
Function (e.g. replace sin(x) with symbol, replace f(x) with symbol
and so on. Get soln from nonlinsolve and then using solveset get
the value of x)

How nonlinsolve is better than old solver _solve_system :
===========================================================

1. A positive dimensional system solver : nonlinsolve can return
solution for positive dimensional system. It finds the
Groebner Basis of the positive dimensional system(calling it as
basis) then we can start solving equation(having least number of
variable first in the basis) using solveset and substituting that
solved solutions into other equation(of basis) to get solution in
terms of minimum variables. Here the important thing is how we
are substituting the known values and in which equations.

2. Real and Complex both solutions : nonlinsolve returns both real
and complex solution. If all the equations in the system are polynomial
then using solve_poly_system both real and complex solution is returned.
If all the equations in the system are not polynomial equation then goes to
substitution method with this polynomial and non polynomial equation(s),
to solve for unsolved variables. Here to solve for particular variable
solveset_real and solveset_complex is used. For both real and complex
solution function _solve_using_know_values is used inside substitution
function.(substitution function will be called when there is any non
polynomial equation(s) is present). When solution is valid then add its
general solution in the final result.

3. Complement and Intersection will be added if any : nonlinsolve maintains
dict for complements and Intersections. If solveset find complements or/and
Intersection with any Interval or set during the execution of
substitution function ,then complement or/and Intersection for that
variable is added before returning final solution.

"""
from sympy.polys.polytools import is_zero_dimensional

if not system:
return S.EmptySet

if not symbols:
msg = ('Symbols must be given, for which solution of the '
'system is to be found.')
raise ValueError(filldedent(msg))

if hasattr(symbols[0], '__iter__'):
symbols = symbols[0]

if not is_sequence(symbols) or not symbols:
msg = ('Symbols must be given, for which solution of the '
'system is to be found.')
raise IndexError(filldedent(msg))

system, symbols, swap = recast_to_symbols(system, symbols)
if swap:
soln = nonlinsolve(system, symbols)
return FiniteSet(*[tuple(i.xreplace(swap) for i in s) for s in soln])

if len(system) == 1 and len(symbols) == 1:
return _solveset_work(system, symbols)

# main code of def nonlinsolve() starts from here
polys, polys_expr, nonpolys, denominators = _separate_poly_nonpoly(
system, symbols)

if len(symbols) == len(polys):
# If all the equations in the system are poly
if is_zero_dimensional(polys, symbols):
# finite number of soln (Zero dimensional system)
try:
return _handle_zero_dimensional(polys, symbols, system)
except NotImplementedError:
# Right now it doesn't fail for any polynomial system of
# equation. If solve_poly_system fails then substitution
# method will handle it.
result = substitution(
polys_expr, symbols, exclude=denominators)
return result

# positive dimensional system
return _handle_positive_dimensional(polys, symbols, denominators)

else:
# If all the equations are not polynomial.
# Use substitution method for the system
result = substitution(
polys_expr + nonpolys, symbols, exclude=denominators)
return result