Coverage for /usr/lib/python3/dist-packages/sympy/solvers/solveset.py: 6%
1384 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2This module contains functions to:
4 - solve a single equation for a single variable, in any domain either real or complex.
6 - solve a single transcendental equation for a single variable in any domain either real or complex.
7 (currently supports solving in real domain only)
9 - solve a system of linear equations with N variables and M equations.
11 - solve a system of Non Linear Equations with N variables and M equations
12"""
13from sympy.core.sympify import sympify
14from sympy.core import (S, Pow, Dummy, pi, Expr, Wild, Mul, Equality,
15 Add, Basic)
16from sympy.core.containers import Tuple
17from sympy.core.function import (Lambda, expand_complex, AppliedUndef,
18 expand_log, _mexpand, expand_trig, nfloat)
19from sympy.core.mod import Mod
20from sympy.core.numbers import igcd, I, Number, Rational, oo, ilcm
21from sympy.core.power import integer_log
22from sympy.core.relational import Eq, Ne, Relational
23from sympy.core.sorting import default_sort_key, ordered
24from sympy.core.symbol import Symbol, _uniquely_named_symbol
25from sympy.core.sympify import _sympify
26from sympy.polys.matrices.linsolve import _linear_eq_to_dict
27from sympy.polys.polyroots import UnsolvableFactorError
28from sympy.simplify.simplify import simplify, fraction, trigsimp, nsimplify
29from sympy.simplify import powdenest, logcombine
30from sympy.functions import (log, tan, cot, sin, cos, sec, csc, exp,
31 acos, asin, acsc, asec,
32 piecewise_fold, Piecewise)
33from sympy.functions.elementary.complexes import Abs, arg, re, im
34from sympy.functions.elementary.hyperbolic import HyperbolicFunction
35from sympy.functions.elementary.miscellaneous import real_root
36from sympy.functions.elementary.trigonometric import TrigonometricFunction
37from sympy.logic.boolalg import And, BooleanTrue
38from sympy.sets import (FiniteSet, imageset, Interval, Intersection,
39 Union, ConditionSet, ImageSet, Complement, Contains)
40from sympy.sets.sets import Set, ProductSet
41from sympy.matrices import zeros, Matrix, MatrixBase
42from sympy.ntheory import totient
43from sympy.ntheory.factor_ import divisors
44from sympy.ntheory.residue_ntheory import discrete_log, nthroot_mod
45from sympy.polys import (roots, Poly, degree, together, PolynomialError,
46 RootOf, factor, lcm, gcd)
47from sympy.polys.polyerrors import CoercionFailed
48from sympy.polys.polytools import invert, groebner, poly
49from sympy.polys.solvers import (sympy_eqs_to_ring, solve_lin_sys,
50 PolyNonlinearError)
51from sympy.polys.matrices.linsolve import _linsolve
52from sympy.solvers.solvers import (checksol, denoms, unrad,
53 _simple_dens, recast_to_symbols)
54from sympy.solvers.polysys import solve_poly_system
55from sympy.utilities import filldedent
56from sympy.utilities.iterables import (numbered_symbols, has_dups,
57 is_sequence, iterable)
58from sympy.calculus.util import periodicity, continuous_domain, function_range
60from types import GeneratorType
63class NonlinearError(ValueError):
64 """Raised when unexpectedly encountering nonlinear equations"""
65 pass
68_rc = Dummy("R", real=True), Dummy("C", complex=True)
71def _masked(f, *atoms):
72 """Return ``f``, with all objects given by ``atoms`` replaced with
73 Dummy symbols, ``d``, and the list of replacements, ``(d, e)``,
74 where ``e`` is an object of type given by ``atoms`` in which
75 any other instances of atoms have been recursively replaced with
76 Dummy symbols, too. The tuples are ordered so that if they are
77 applied in sequence, the origin ``f`` will be restored.
79 Examples
80 ========
82 >>> from sympy import cos
83 >>> from sympy.abc import x
84 >>> from sympy.solvers.solveset import _masked
86 >>> f = cos(cos(x) + 1)
87 >>> f, reps = _masked(cos(1 + cos(x)), cos)
88 >>> f
89 _a1
90 >>> reps
91 [(_a1, cos(_a0 + 1)), (_a0, cos(x))]
92 >>> for d, e in reps:
93 ... f = f.xreplace({d: e})
94 >>> f
95 cos(cos(x) + 1)
96 """
97 sym = numbered_symbols('a', cls=Dummy, real=True)
98 mask = []
99 for a in ordered(f.atoms(*atoms)):
100 for i in mask:
101 a = a.replace(*i)
102 mask.append((a, next(sym)))
103 for i, (o, n) in enumerate(mask):
104 f = f.replace(o, n)
105 mask[i] = (n, o)
106 mask = list(reversed(mask))
107 return f, mask
110def _invert(f_x, y, x, domain=S.Complexes):
111 r"""
112 Reduce the complex valued equation $f(x) = y$ to a set of equations
114 $$\left\{g(x) = h_1(y),\ g(x) = h_2(y),\ \dots,\ g(x) = h_n(y) \right\}$$
116 where $g(x)$ is a simpler function than $f(x)$. The return value is a tuple
117 $(g(x), \mathrm{set}_h)$, where $g(x)$ is a function of $x$ and $\mathrm{set}_h$ is
118 the set of function $\left\{h_1(y), h_2(y), \dots, h_n(y)\right\}$.
119 Here, $y$ is not necessarily a symbol.
121 $\mathrm{set}_h$ contains the functions, along with the information
122 about the domain in which they are valid, through set
123 operations. For instance, if :math:`y = |x| - n` is inverted
124 in the real domain, then $\mathrm{set}_h$ is not simply
125 $\{-n, n\}$ as the nature of `n` is unknown; rather, it is:
127 $$ \left(\left[0, \infty\right) \cap \left\{n\right\}\right) \cup
128 \left(\left(-\infty, 0\right] \cap \left\{- n\right\}\right)$$
130 By default, the complex domain is used which means that inverting even
131 seemingly simple functions like $\exp(x)$ will give very different
132 results from those obtained in the real domain.
133 (In the case of $\exp(x)$, the inversion via $\log$ is multi-valued
134 in the complex domain, having infinitely many branches.)
136 If you are working with real values only (or you are not sure which
137 function to use) you should probably set the domain to
138 ``S.Reals`` (or use ``invert_real`` which does that automatically).
141 Examples
142 ========
144 >>> from sympy.solvers.solveset import invert_complex, invert_real
145 >>> from sympy.abc import x, y
146 >>> from sympy import exp
148 When does exp(x) == y?
150 >>> invert_complex(exp(x), y, x)
151 (x, ImageSet(Lambda(_n, I*(2*_n*pi + arg(y)) + log(Abs(y))), Integers))
152 >>> invert_real(exp(x), y, x)
153 (x, Intersection({log(y)}, Reals))
155 When does exp(x) == 1?
157 >>> invert_complex(exp(x), 1, x)
158 (x, ImageSet(Lambda(_n, 2*_n*I*pi), Integers))
159 >>> invert_real(exp(x), 1, x)
160 (x, {0})
162 See Also
163 ========
164 invert_real, invert_complex
165 """
166 x = sympify(x)
167 if not x.is_Symbol:
168 raise ValueError("x must be a symbol")
169 f_x = sympify(f_x)
170 if x not in f_x.free_symbols:
171 raise ValueError("Inverse of constant function doesn't exist")
172 y = sympify(y)
173 if x in y.free_symbols:
174 raise ValueError("y should be independent of x ")
176 if domain.is_subset(S.Reals):
177 x1, s = _invert_real(f_x, FiniteSet(y), x)
178 else:
179 x1, s = _invert_complex(f_x, FiniteSet(y), x)
181 if not isinstance(s, FiniteSet) or x1 != x:
182 return x1, s
184 # Avoid adding gratuitous intersections with S.Complexes. Actual
185 # conditions should be handled by the respective inverters.
186 if domain is S.Complexes:
187 return x1, s
188 else:
189 return x1, s.intersection(domain)
192invert_complex = _invert
195def invert_real(f_x, y, x):
196 """
197 Inverts a real-valued function. Same as :func:`invert_complex`, but sets
198 the domain to ``S.Reals`` before inverting.
199 """
200 return _invert(f_x, y, x, S.Reals)
203def _invert_real(f, g_ys, symbol):
204 """Helper function for _invert."""
206 if f == symbol or g_ys is S.EmptySet:
207 return (f, g_ys)
209 n = Dummy('n', real=True)
211 if isinstance(f, exp) or (f.is_Pow and f.base == S.Exp1):
212 return _invert_real(f.exp,
213 imageset(Lambda(n, log(n)), g_ys),
214 symbol)
216 if hasattr(f, 'inverse') and f.inverse() is not None and not isinstance(f, (
217 TrigonometricFunction,
218 HyperbolicFunction,
219 )):
220 if len(f.args) > 1:
221 raise ValueError("Only functions with one argument are supported.")
222 return _invert_real(f.args[0],
223 imageset(Lambda(n, f.inverse()(n)), g_ys),
224 symbol)
226 if isinstance(f, Abs):
227 return _invert_abs(f.args[0], g_ys, symbol)
229 if f.is_Add:
230 # f = g + h
231 g, h = f.as_independent(symbol)
232 if g is not S.Zero:
233 return _invert_real(h, imageset(Lambda(n, n - g), g_ys), symbol)
235 if f.is_Mul:
236 # f = g*h
237 g, h = f.as_independent(symbol)
239 if g is not S.One:
240 return _invert_real(h, imageset(Lambda(n, n/g), g_ys), symbol)
242 if f.is_Pow:
243 base, expo = f.args
244 base_has_sym = base.has(symbol)
245 expo_has_sym = expo.has(symbol)
247 if not expo_has_sym:
249 if expo.is_rational:
250 num, den = expo.as_numer_denom()
252 if den % 2 == 0 and num % 2 == 1 and den.is_zero is False:
253 # Here we have f(x)**(num/den) = y
254 # where den is nonzero and even and y is an element
255 # of the set g_ys.
256 # den is even, so we are only interested in the cases
257 # where both f(x) and y are positive.
258 # Restricting y to be positive (using the set g_ys_pos)
259 # means that y**(den/num) is always positive.
260 # Therefore it isn't necessary to also constrain f(x)
261 # to be positive because we are only going to
262 # find solutions of f(x) = y**(d/n)
263 # where the rhs is already required to be positive.
264 root = Lambda(n, real_root(n, expo))
265 g_ys_pos = g_ys & Interval(0, oo)
266 res = imageset(root, g_ys_pos)
267 _inv, _set = _invert_real(base, res, symbol)
268 return (_inv, _set)
270 if den % 2 == 1:
271 root = Lambda(n, real_root(n, expo))
272 res = imageset(root, g_ys)
273 if num % 2 == 0:
274 neg_res = imageset(Lambda(n, -n), res)
275 return _invert_real(base, res + neg_res, symbol)
276 if num % 2 == 1:
277 return _invert_real(base, res, symbol)
279 elif expo.is_irrational:
280 root = Lambda(n, real_root(n, expo))
281 g_ys_pos = g_ys & Interval(0, oo)
282 res = imageset(root, g_ys_pos)
283 return _invert_real(base, res, symbol)
285 else:
286 # indeterminate exponent, e.g. Float or parity of
287 # num, den of rational could not be determined
288 pass # use default return
290 if not base_has_sym:
291 rhs = g_ys.args[0]
292 if base.is_positive:
293 return _invert_real(expo,
294 imageset(Lambda(n, log(n, base, evaluate=False)), g_ys), symbol)
295 elif base.is_negative:
296 s, b = integer_log(rhs, base)
297 if b:
298 return _invert_real(expo, FiniteSet(s), symbol)
299 else:
300 return (expo, S.EmptySet)
301 elif base.is_zero:
302 one = Eq(rhs, 1)
303 if one == S.true:
304 # special case: 0**x - 1
305 return _invert_real(expo, FiniteSet(0), symbol)
306 elif one == S.false:
307 return (expo, S.EmptySet)
310 if isinstance(f, TrigonometricFunction):
311 if isinstance(g_ys, FiniteSet):
312 def inv(trig):
313 if isinstance(trig, (sin, csc)):
314 F = asin if isinstance(trig, sin) else acsc
315 return (lambda a: n*pi + S.NegativeOne**n*F(a),)
316 if isinstance(trig, (cos, sec)):
317 F = acos if isinstance(trig, cos) else asec
318 return (
319 lambda a: 2*n*pi + F(a),
320 lambda a: 2*n*pi - F(a),)
321 if isinstance(trig, (tan, cot)):
322 return (lambda a: n*pi + trig.inverse()(a),)
324 n = Dummy('n', integer=True)
325 invs = S.EmptySet
326 for L in inv(f):
327 invs += Union(*[imageset(Lambda(n, L(g)), S.Integers) for g in g_ys])
328 return _invert_real(f.args[0], invs, symbol)
330 return (f, g_ys)
333def _invert_complex(f, g_ys, symbol):
334 """Helper function for _invert."""
336 if f == symbol or g_ys is S.EmptySet:
337 return (f, g_ys)
339 n = Dummy('n')
341 if f.is_Add:
342 # f = g + h
343 g, h = f.as_independent(symbol)
344 if g is not S.Zero:
345 return _invert_complex(h, imageset(Lambda(n, n - g), g_ys), symbol)
347 if f.is_Mul:
348 # f = g*h
349 g, h = f.as_independent(symbol)
351 if g is not S.One:
352 if g in {S.NegativeInfinity, S.ComplexInfinity, S.Infinity}:
353 return (h, S.EmptySet)
354 return _invert_complex(h, imageset(Lambda(n, n/g), g_ys), symbol)
356 if f.is_Pow:
357 base, expo = f.args
358 # special case: g**r = 0
359 # Could be improved like `_invert_real` to handle more general cases.
360 if expo.is_Rational and g_ys == FiniteSet(0):
361 if expo.is_positive:
362 return _invert_complex(base, g_ys, symbol)
364 if hasattr(f, 'inverse') and f.inverse() is not None and \
365 not isinstance(f, TrigonometricFunction) and \
366 not isinstance(f, HyperbolicFunction) and \
367 not isinstance(f, exp):
368 if len(f.args) > 1:
369 raise ValueError("Only functions with one argument are supported.")
370 return _invert_complex(f.args[0],
371 imageset(Lambda(n, f.inverse()(n)), g_ys), symbol)
373 if isinstance(f, exp) or (f.is_Pow and f.base == S.Exp1):
374 if isinstance(g_ys, ImageSet):
375 # can solve upto `(d*exp(exp(...(exp(a*x + b))...) + c)` format.
376 # Further can be improved to `(d*exp(exp(...(exp(a*x**n + b*x**(n-1) + ... + f))...) + c)`.
377 g_ys_expr = g_ys.lamda.expr
378 g_ys_vars = g_ys.lamda.variables
379 k = Dummy('k{}'.format(len(g_ys_vars)))
380 g_ys_vars_1 = (k,) + g_ys_vars
381 exp_invs = Union(*[imageset(Lambda((g_ys_vars_1,), (I*(2*k*pi + arg(g_ys_expr))
382 + log(Abs(g_ys_expr)))), S.Integers**(len(g_ys_vars_1)))])
383 return _invert_complex(f.exp, exp_invs, symbol)
385 elif isinstance(g_ys, FiniteSet):
386 exp_invs = Union(*[imageset(Lambda(n, I*(2*n*pi + arg(g_y)) +
387 log(Abs(g_y))), S.Integers)
388 for g_y in g_ys if g_y != 0])
389 return _invert_complex(f.exp, exp_invs, symbol)
391 return (f, g_ys)
394def _invert_abs(f, g_ys, symbol):
395 """Helper function for inverting absolute value functions.
397 Returns the complete result of inverting an absolute value
398 function along with the conditions which must also be satisfied.
400 If it is certain that all these conditions are met, a :class:`~.FiniteSet`
401 of all possible solutions is returned. If any condition cannot be
402 satisfied, an :class:`~.EmptySet` is returned. Otherwise, a
403 :class:`~.ConditionSet` of the solutions, with all the required conditions
404 specified, is returned.
406 """
407 if not g_ys.is_FiniteSet:
408 # this could be used for FiniteSet, but the
409 # results are more compact if they aren't, e.g.
410 # ConditionSet(x, Contains(n, Interval(0, oo)), {-n, n}) vs
411 # Union(Intersection(Interval(0, oo), {n}), Intersection(Interval(-oo, 0), {-n}))
412 # for the solution of abs(x) - n
413 pos = Intersection(g_ys, Interval(0, S.Infinity))
414 parg = _invert_real(f, pos, symbol)
415 narg = _invert_real(-f, pos, symbol)
416 if parg[0] != narg[0]:
417 raise NotImplementedError
418 return parg[0], Union(narg[1], parg[1])
420 # check conditions: all these must be true. If any are unknown
421 # then return them as conditions which must be satisfied
422 unknown = []
423 for a in g_ys.args:
424 ok = a.is_nonnegative if a.is_Number else a.is_positive
425 if ok is None:
426 unknown.append(a)
427 elif not ok:
428 return symbol, S.EmptySet
429 if unknown:
430 conditions = And(*[Contains(i, Interval(0, oo))
431 for i in unknown])
432 else:
433 conditions = True
434 n = Dummy('n', real=True)
435 # this is slightly different than above: instead of solving
436 # +/-f on positive values, here we solve for f on +/- g_ys
437 g_x, values = _invert_real(f, Union(
438 imageset(Lambda(n, n), g_ys),
439 imageset(Lambda(n, -n), g_ys)), symbol)
440 return g_x, ConditionSet(g_x, conditions, values)
443def domain_check(f, symbol, p):
444 """Returns False if point p is infinite or any subexpression of f
445 is infinite or becomes so after replacing symbol with p. If none of
446 these conditions is met then True will be returned.
448 Examples
449 ========
451 >>> from sympy import Mul, oo
452 >>> from sympy.abc import x
453 >>> from sympy.solvers.solveset import domain_check
454 >>> g = 1/(1 + (1/(x + 1))**2)
455 >>> domain_check(g, x, -1)
456 False
457 >>> domain_check(x**2, x, 0)
458 True
459 >>> domain_check(1/x, x, oo)
460 False
462 * The function relies on the assumption that the original form
463 of the equation has not been changed by automatic simplification.
465 >>> domain_check(x/x, x, 0) # x/x is automatically simplified to 1
466 True
468 * To deal with automatic evaluations use evaluate=False:
470 >>> domain_check(Mul(x, 1/x, evaluate=False), x, 0)
471 False
472 """
473 f, p = sympify(f), sympify(p)
474 if p.is_infinite:
475 return False
476 return _domain_check(f, symbol, p)
479def _domain_check(f, symbol, p):
480 # helper for domain check
481 if f.is_Atom and f.is_finite:
482 return True
483 elif f.subs(symbol, p).is_infinite:
484 return False
485 elif isinstance(f, Piecewise):
486 # Check the cases of the Piecewise in turn. There might be invalid
487 # expressions in later cases that don't apply e.g.
488 # solveset(Piecewise((0, Eq(x, 0)), (1/x, True)), x)
489 for expr, cond in f.args:
490 condsubs = cond.subs(symbol, p)
491 if condsubs is S.false:
492 continue
493 elif condsubs is S.true:
494 return _domain_check(expr, symbol, p)
495 else:
496 # We don't know which case of the Piecewise holds. On this
497 # basis we cannot decide whether any solution is in or out of
498 # the domain. Ideally this function would allow returning a
499 # symbolic condition for the validity of the solution that
500 # could be handled in the calling code. In the mean time we'll
501 # give this particular solution the benefit of the doubt and
502 # let it pass.
503 return True
504 else:
505 # TODO : We should not blindly recurse through all args of arbitrary expressions like this
506 return all(_domain_check(g, symbol, p)
507 for g in f.args)
510def _is_finite_with_finite_vars(f, domain=S.Complexes):
511 """
512 Return True if the given expression is finite. For symbols that
513 do not assign a value for `complex` and/or `real`, the domain will
514 be used to assign a value; symbols that do not assign a value
515 for `finite` will be made finite. All other assumptions are
516 left unmodified.
517 """
518 def assumptions(s):
519 A = s.assumptions0
520 A.setdefault('finite', A.get('finite', True))
521 if domain.is_subset(S.Reals):
522 # if this gets set it will make complex=True, too
523 A.setdefault('real', True)
524 else:
525 # don't change 'real' because being complex implies
526 # nothing about being real
527 A.setdefault('complex', True)
528 return A
530 reps = {s: Dummy(**assumptions(s)) for s in f.free_symbols}
531 return f.xreplace(reps).is_finite
534def _is_function_class_equation(func_class, f, symbol):
535 """ Tests whether the equation is an equation of the given function class.
537 The given equation belongs to the given function class if it is
538 comprised of functions of the function class which are multiplied by
539 or added to expressions independent of the symbol. In addition, the
540 arguments of all such functions must be linear in the symbol as well.
542 Examples
543 ========
545 >>> from sympy.solvers.solveset import _is_function_class_equation
546 >>> from sympy import tan, sin, tanh, sinh, exp
547 >>> from sympy.abc import x
548 >>> from sympy.functions.elementary.trigonometric import TrigonometricFunction
549 >>> from sympy.functions.elementary.hyperbolic import HyperbolicFunction
550 >>> _is_function_class_equation(TrigonometricFunction, exp(x) + tan(x), x)
551 False
552 >>> _is_function_class_equation(TrigonometricFunction, tan(x) + sin(x), x)
553 True
554 >>> _is_function_class_equation(TrigonometricFunction, tan(x**2), x)
555 False
556 >>> _is_function_class_equation(TrigonometricFunction, tan(x + 2), x)
557 True
558 >>> _is_function_class_equation(HyperbolicFunction, tanh(x) + sinh(x), x)
559 True
560 """
561 if f.is_Mul or f.is_Add:
562 return all(_is_function_class_equation(func_class, arg, symbol)
563 for arg in f.args)
565 if f.is_Pow:
566 if not f.exp.has(symbol):
567 return _is_function_class_equation(func_class, f.base, symbol)
568 else:
569 return False
571 if not f.has(symbol):
572 return True
574 if isinstance(f, func_class):
575 try:
576 g = Poly(f.args[0], symbol)
577 return g.degree() <= 1
578 except PolynomialError:
579 return False
580 else:
581 return False
584def _solve_as_rational(f, symbol, domain):
585 """ solve rational functions"""
586 f = together(_mexpand(f, recursive=True), deep=True)
587 g, h = fraction(f)
588 if not h.has(symbol):
589 try:
590 return _solve_as_poly(g, symbol, domain)
591 except NotImplementedError:
592 # The polynomial formed from g could end up having
593 # coefficients in a ring over which finding roots
594 # isn't implemented yet, e.g. ZZ[a] for some symbol a
595 return ConditionSet(symbol, Eq(f, 0), domain)
596 except CoercionFailed:
597 # contained oo, zoo or nan
598 return S.EmptySet
599 else:
600 valid_solns = _solveset(g, symbol, domain)
601 invalid_solns = _solveset(h, symbol, domain)
602 return valid_solns - invalid_solns
605class _SolveTrig1Error(Exception):
606 """Raised when _solve_trig1 heuristics do not apply"""
608def _solve_trig(f, symbol, domain):
609 """Function to call other helpers to solve trigonometric equations """
610 sol = None
611 try:
612 sol = _solve_trig1(f, symbol, domain)
613 except _SolveTrig1Error:
614 try:
615 sol = _solve_trig2(f, symbol, domain)
616 except ValueError:
617 raise NotImplementedError(filldedent('''
618 Solution to this kind of trigonometric equations
619 is yet to be implemented'''))
620 return sol
623def _solve_trig1(f, symbol, domain):
624 """Primary solver for trigonometric and hyperbolic equations
626 Returns either the solution set as a ConditionSet (auto-evaluated to a
627 union of ImageSets if no variables besides 'symbol' are involved) or
628 raises _SolveTrig1Error if f == 0 cannot be solved.
630 Notes
631 =====
632 Algorithm:
633 1. Do a change of variable x -> mu*x in arguments to trigonometric and
634 hyperbolic functions, in order to reduce them to small integers. (This
635 step is crucial to keep the degrees of the polynomials of step 4 low.)
636 2. Rewrite trigonometric/hyperbolic functions as exponentials.
637 3. Proceed to a 2nd change of variable, replacing exp(I*x) or exp(x) by y.
638 4. Solve the resulting rational equation.
639 5. Use invert_complex or invert_real to return to the original variable.
640 6. If the coefficients of 'symbol' were symbolic in nature, add the
641 necessary consistency conditions in a ConditionSet.
643 """
644 # Prepare change of variable
645 x = Dummy('x')
646 if _is_function_class_equation(HyperbolicFunction, f, symbol):
647 cov = exp(x)
648 inverter = invert_real if domain.is_subset(S.Reals) else invert_complex
649 else:
650 cov = exp(I*x)
651 inverter = invert_complex
653 f = trigsimp(f)
654 f_original = f
655 trig_functions = f.atoms(TrigonometricFunction, HyperbolicFunction)
656 trig_arguments = [e.args[0] for e in trig_functions]
657 # trigsimp may have reduced the equation to an expression
658 # that is independent of 'symbol' (e.g. cos**2+sin**2)
659 if not any(a.has(symbol) for a in trig_arguments):
660 return solveset(f_original, symbol, domain)
662 denominators = []
663 numerators = []
664 for ar in trig_arguments:
665 try:
666 poly_ar = Poly(ar, symbol)
667 except PolynomialError:
668 raise _SolveTrig1Error("trig argument is not a polynomial")
669 if poly_ar.degree() > 1: # degree >1 still bad
670 raise _SolveTrig1Error("degree of variable must not exceed one")
671 if poly_ar.degree() == 0: # degree 0, don't care
672 continue
673 c = poly_ar.all_coeffs()[0] # got the coefficient of 'symbol'
674 numerators.append(fraction(c)[0])
675 denominators.append(fraction(c)[1])
677 mu = lcm(denominators)/gcd(numerators)
678 f = f.subs(symbol, mu*x)
679 f = f.rewrite(exp)
680 f = together(f)
681 g, h = fraction(f)
682 y = Dummy('y')
683 g, h = g.expand(), h.expand()
684 g, h = g.subs(cov, y), h.subs(cov, y)
685 if g.has(x) or h.has(x):
686 raise _SolveTrig1Error("change of variable not possible")
688 solns = solveset_complex(g, y) - solveset_complex(h, y)
689 if isinstance(solns, ConditionSet):
690 raise _SolveTrig1Error("polynomial has ConditionSet solution")
692 if isinstance(solns, FiniteSet):
693 if any(isinstance(s, RootOf) for s in solns):
694 raise _SolveTrig1Error("polynomial results in RootOf object")
695 # revert the change of variable
696 cov = cov.subs(x, symbol/mu)
697 result = Union(*[inverter(cov, s, symbol)[1] for s in solns])
698 # In case of symbolic coefficients, the solution set is only valid
699 # if numerator and denominator of mu are non-zero.
700 if mu.has(Symbol):
701 syms = (mu).atoms(Symbol)
702 munum, muden = fraction(mu)
703 condnum = munum.as_independent(*syms, as_Add=False)[1]
704 condden = muden.as_independent(*syms, as_Add=False)[1]
705 cond = And(Ne(condnum, 0), Ne(condden, 0))
706 else:
707 cond = True
708 # Actual conditions are returned as part of the ConditionSet. Adding an
709 # intersection with C would only complicate some solution sets due to
710 # current limitations of intersection code. (e.g. #19154)
711 if domain is S.Complexes:
712 # This is a slight abuse of ConditionSet. Ideally this should
713 # be some kind of "PiecewiseSet". (See #19507 discussion)
714 return ConditionSet(symbol, cond, result)
715 else:
716 return ConditionSet(symbol, cond, Intersection(result, domain))
717 elif solns is S.EmptySet:
718 return S.EmptySet
719 else:
720 raise _SolveTrig1Error("polynomial solutions must form FiniteSet")
723def _solve_trig2(f, symbol, domain):
724 """Secondary helper to solve trigonometric equations,
725 called when first helper fails """
726 f = trigsimp(f)
727 f_original = f
728 trig_functions = f.atoms(sin, cos, tan, sec, cot, csc)
729 trig_arguments = [e.args[0] for e in trig_functions]
730 denominators = []
731 numerators = []
733 # todo: This solver can be extended to hyperbolics if the
734 # analogous change of variable to tanh (instead of tan)
735 # is used.
736 if not trig_functions:
737 return ConditionSet(symbol, Eq(f_original, 0), domain)
739 # todo: The pre-processing below (extraction of numerators, denominators,
740 # gcd, lcm, mu, etc.) should be updated to the enhanced version in
741 # _solve_trig1. (See #19507)
742 for ar in trig_arguments:
743 try:
744 poly_ar = Poly(ar, symbol)
745 except PolynomialError:
746 raise ValueError("give up, we cannot solve if this is not a polynomial in x")
747 if poly_ar.degree() > 1: # degree >1 still bad
748 raise ValueError("degree of variable inside polynomial should not exceed one")
749 if poly_ar.degree() == 0: # degree 0, don't care
750 continue
751 c = poly_ar.all_coeffs()[0] # got the coefficient of 'symbol'
752 try:
753 numerators.append(Rational(c).p)
754 denominators.append(Rational(c).q)
755 except TypeError:
756 return ConditionSet(symbol, Eq(f_original, 0), domain)
758 x = Dummy('x')
760 # ilcm() and igcd() require more than one argument
761 if len(numerators) > 1:
762 mu = Rational(2)*ilcm(*denominators)/igcd(*numerators)
763 else:
764 assert len(numerators) == 1
765 mu = Rational(2)*denominators[0]/numerators[0]
767 f = f.subs(symbol, mu*x)
768 f = f.rewrite(tan)
769 f = expand_trig(f)
770 f = together(f)
772 g, h = fraction(f)
773 y = Dummy('y')
774 g, h = g.expand(), h.expand()
775 g, h = g.subs(tan(x), y), h.subs(tan(x), y)
777 if g.has(x) or h.has(x):
778 return ConditionSet(symbol, Eq(f_original, 0), domain)
779 solns = solveset(g, y, S.Reals) - solveset(h, y, S.Reals)
781 if isinstance(solns, FiniteSet):
782 result = Union(*[invert_real(tan(symbol/mu), s, symbol)[1]
783 for s in solns])
784 dsol = invert_real(tan(symbol/mu), oo, symbol)[1]
785 if degree(h) > degree(g): # If degree(denom)>degree(num) then there
786 result = Union(result, dsol) # would be another sol at Lim(denom-->oo)
787 return Intersection(result, domain)
788 elif solns is S.EmptySet:
789 return S.EmptySet
790 else:
791 return ConditionSet(symbol, Eq(f_original, 0), S.Reals)
794def _solve_as_poly(f, symbol, domain=S.Complexes):
795 """
796 Solve the equation using polynomial techniques if it already is a
797 polynomial equation or, with a change of variables, can be made so.
798 """
799 result = None
800 if f.is_polynomial(symbol):
801 solns = roots(f, symbol, cubics=True, quartics=True,
802 quintics=True, domain='EX')
803 num_roots = sum(solns.values())
804 if degree(f, symbol) <= num_roots:
805 result = FiniteSet(*solns.keys())
806 else:
807 poly = Poly(f, symbol)
808 solns = poly.all_roots()
809 if poly.degree() <= len(solns):
810 result = FiniteSet(*solns)
811 else:
812 result = ConditionSet(symbol, Eq(f, 0), domain)
813 else:
814 poly = Poly(f)
815 if poly is None:
816 result = ConditionSet(symbol, Eq(f, 0), domain)
817 gens = [g for g in poly.gens if g.has(symbol)]
819 if len(gens) == 1:
820 poly = Poly(poly, gens[0])
821 gen = poly.gen
822 deg = poly.degree()
823 poly = Poly(poly.as_expr(), poly.gen, composite=True)
824 poly_solns = FiniteSet(*roots(poly, cubics=True, quartics=True,
825 quintics=True).keys())
827 if len(poly_solns) < deg:
828 result = ConditionSet(symbol, Eq(f, 0), domain)
830 if gen != symbol:
831 y = Dummy('y')
832 inverter = invert_real if domain.is_subset(S.Reals) else invert_complex
833 lhs, rhs_s = inverter(gen, y, symbol)
834 if lhs == symbol:
835 result = Union(*[rhs_s.subs(y, s) for s in poly_solns])
836 if isinstance(result, FiniteSet) and isinstance(gen, Pow
837 ) and gen.base.is_Rational:
838 result = FiniteSet(*[expand_log(i) for i in result])
839 else:
840 result = ConditionSet(symbol, Eq(f, 0), domain)
841 else:
842 result = ConditionSet(symbol, Eq(f, 0), domain)
844 if result is not None:
845 if isinstance(result, FiniteSet):
846 # this is to simplify solutions like -sqrt(-I) to sqrt(2)/2
847 # - sqrt(2)*I/2. We are not expanding for solution with symbols
848 # or undefined functions because that makes the solution more complicated.
849 # For example, expand_complex(a) returns re(a) + I*im(a)
850 if all(s.atoms(Symbol, AppliedUndef) == set() and not isinstance(s, RootOf)
851 for s in result):
852 s = Dummy('s')
853 result = imageset(Lambda(s, expand_complex(s)), result)
854 if isinstance(result, FiniteSet) and domain != S.Complexes:
855 # Avoid adding gratuitous intersections with S.Complexes. Actual
856 # conditions should be handled elsewhere.
857 result = result.intersection(domain)
858 return result
859 else:
860 return ConditionSet(symbol, Eq(f, 0), domain)
863def _solve_radical(f, unradf, symbol, solveset_solver):
864 """ Helper function to solve equations with radicals """
865 res = unradf
866 eq, cov = res if res else (f, [])
867 if not cov:
868 result = solveset_solver(eq, symbol) - \
869 Union(*[solveset_solver(g, symbol) for g in denoms(f, symbol)])
870 else:
871 y, yeq = cov
872 if not solveset_solver(y - I, y):
873 yreal = Dummy('yreal', real=True)
874 yeq = yeq.xreplace({y: yreal})
875 eq = eq.xreplace({y: yreal})
876 y = yreal
877 g_y_s = solveset_solver(yeq, symbol)
878 f_y_sols = solveset_solver(eq, y)
879 result = Union(*[imageset(Lambda(y, g_y), f_y_sols)
880 for g_y in g_y_s])
882 def check_finiteset(solutions):
883 f_set = [] # solutions for FiniteSet
884 c_set = [] # solutions for ConditionSet
885 for s in solutions:
886 if checksol(f, symbol, s):
887 f_set.append(s)
888 else:
889 c_set.append(s)
890 return FiniteSet(*f_set) + ConditionSet(symbol, Eq(f, 0), FiniteSet(*c_set))
892 def check_set(solutions):
893 if solutions is S.EmptySet:
894 return solutions
895 elif isinstance(solutions, ConditionSet):
896 # XXX: Maybe the base set should be checked?
897 return solutions
898 elif isinstance(solutions, FiniteSet):
899 return check_finiteset(solutions)
900 elif isinstance(solutions, Complement):
901 A, B = solutions.args
902 return Complement(check_set(A), B)
903 elif isinstance(solutions, Union):
904 return Union(*[check_set(s) for s in solutions.args])
905 else:
906 # XXX: There should be more cases checked here. The cases above
907 # are all those that come up in the test suite for now.
908 return solutions
910 solution_set = check_set(result)
912 return solution_set
915def _solve_abs(f, symbol, domain):
916 """ Helper function to solve equation involving absolute value function """
917 if not domain.is_subset(S.Reals):
918 raise ValueError(filldedent('''
919 Absolute values cannot be inverted in the
920 complex domain.'''))
921 p, q, r = Wild('p'), Wild('q'), Wild('r')
922 pattern_match = f.match(p*Abs(q) + r) or {}
923 f_p, f_q, f_r = [pattern_match.get(i, S.Zero) for i in (p, q, r)]
925 if not (f_p.is_zero or f_q.is_zero):
926 domain = continuous_domain(f_q, symbol, domain)
927 from .inequalities import solve_univariate_inequality
928 q_pos_cond = solve_univariate_inequality(f_q >= 0, symbol,
929 relational=False, domain=domain, continuous=True)
930 q_neg_cond = q_pos_cond.complement(domain)
932 sols_q_pos = solveset_real(f_p*f_q + f_r,
933 symbol).intersect(q_pos_cond)
934 sols_q_neg = solveset_real(f_p*(-f_q) + f_r,
935 symbol).intersect(q_neg_cond)
936 return Union(sols_q_pos, sols_q_neg)
937 else:
938 return ConditionSet(symbol, Eq(f, 0), domain)
941def solve_decomposition(f, symbol, domain):
942 """
943 Function to solve equations via the principle of "Decomposition
944 and Rewriting".
946 Examples
947 ========
948 >>> from sympy import exp, sin, Symbol, pprint, S
949 >>> from sympy.solvers.solveset import solve_decomposition as sd
950 >>> x = Symbol('x')
951 >>> f1 = exp(2*x) - 3*exp(x) + 2
952 >>> sd(f1, x, S.Reals)
953 {0, log(2)}
954 >>> f2 = sin(x)**2 + 2*sin(x) + 1
955 >>> pprint(sd(f2, x, S.Reals), use_unicode=False)
956 3*pi
957 {2*n*pi + ---- | n in Integers}
958 2
959 >>> f3 = sin(x + 2)
960 >>> pprint(sd(f3, x, S.Reals), use_unicode=False)
961 {2*n*pi - 2 | n in Integers} U {2*n*pi - 2 + pi | n in Integers}
963 """
964 from sympy.solvers.decompogen import decompogen
965 # decompose the given function
966 g_s = decompogen(f, symbol)
967 # `y_s` represents the set of values for which the function `g` is to be
968 # solved.
969 # `solutions` represent the solutions of the equations `g = y_s` or
970 # `g = 0` depending on the type of `y_s`.
971 # As we are interested in solving the equation: f = 0
972 y_s = FiniteSet(0)
973 for g in g_s:
974 frange = function_range(g, symbol, domain)
975 y_s = Intersection(frange, y_s)
976 result = S.EmptySet
977 if isinstance(y_s, FiniteSet):
978 for y in y_s:
979 solutions = solveset(Eq(g, y), symbol, domain)
980 if not isinstance(solutions, ConditionSet):
981 result += solutions
983 else:
984 if isinstance(y_s, ImageSet):
985 iter_iset = (y_s,)
987 elif isinstance(y_s, Union):
988 iter_iset = y_s.args
990 elif y_s is S.EmptySet:
991 # y_s is not in the range of g in g_s, so no solution exists
992 #in the given domain
993 return S.EmptySet
995 for iset in iter_iset:
996 new_solutions = solveset(Eq(iset.lamda.expr, g), symbol, domain)
997 dummy_var = tuple(iset.lamda.expr.free_symbols)[0]
998 (base_set,) = iset.base_sets
999 if isinstance(new_solutions, FiniteSet):
1000 new_exprs = new_solutions
1002 elif isinstance(new_solutions, Intersection):
1003 if isinstance(new_solutions.args[1], FiniteSet):
1004 new_exprs = new_solutions.args[1]
1006 for new_expr in new_exprs:
1007 result += ImageSet(Lambda(dummy_var, new_expr), base_set)
1009 if result is S.EmptySet:
1010 return ConditionSet(symbol, Eq(f, 0), domain)
1012 y_s = result
1014 return y_s
1017def _solveset(f, symbol, domain, _check=False):
1018 """Helper for solveset to return a result from an expression
1019 that has already been sympify'ed and is known to contain the
1020 given symbol."""
1021 # _check controls whether the answer is checked or not
1022 from sympy.simplify.simplify import signsimp
1024 if isinstance(f, BooleanTrue):
1025 return domain
1027 orig_f = f
1028 if f.is_Mul:
1029 coeff, f = f.as_independent(symbol, as_Add=False)
1030 if coeff in {S.ComplexInfinity, S.NegativeInfinity, S.Infinity}:
1031 f = together(orig_f)
1032 elif f.is_Add:
1033 a, h = f.as_independent(symbol)
1034 m, h = h.as_independent(symbol, as_Add=False)
1035 if m not in {S.ComplexInfinity, S.Zero, S.Infinity,
1036 S.NegativeInfinity}:
1037 f = a/m + h # XXX condition `m != 0` should be added to soln
1039 # assign the solvers to use
1040 solver = lambda f, x, domain=domain: _solveset(f, x, domain)
1041 inverter = lambda f, rhs, symbol: _invert(f, rhs, symbol, domain)
1043 result = S.EmptySet
1045 if f.expand().is_zero:
1046 return domain
1047 elif not f.has(symbol):
1048 return S.EmptySet
1049 elif f.is_Mul and all(_is_finite_with_finite_vars(m, domain)
1050 for m in f.args):
1051 # if f(x) and g(x) are both finite we can say that the solution of
1052 # f(x)*g(x) == 0 is same as Union(f(x) == 0, g(x) == 0) is not true in
1053 # general. g(x) can grow to infinitely large for the values where
1054 # f(x) == 0. To be sure that we are not silently allowing any
1055 # wrong solutions we are using this technique only if both f and g are
1056 # finite for a finite input.
1057 result = Union(*[solver(m, symbol) for m in f.args])
1058 elif _is_function_class_equation(TrigonometricFunction, f, symbol) or \
1059 _is_function_class_equation(HyperbolicFunction, f, symbol):
1060 result = _solve_trig(f, symbol, domain)
1061 elif isinstance(f, arg):
1062 a = f.args[0]
1063 result = Intersection(_solveset(re(a) > 0, symbol, domain),
1064 _solveset(im(a), symbol, domain))
1065 elif f.is_Piecewise:
1066 expr_set_pairs = f.as_expr_set_pairs(domain)
1067 for (expr, in_set) in expr_set_pairs:
1068 if in_set.is_Relational:
1069 in_set = in_set.as_set()
1070 solns = solver(expr, symbol, in_set)
1071 result += solns
1072 elif isinstance(f, Eq):
1073 result = solver(Add(f.lhs, - f.rhs, evaluate=False), symbol, domain)
1075 elif f.is_Relational:
1076 from .inequalities import solve_univariate_inequality
1077 try:
1078 result = solve_univariate_inequality(
1079 f, symbol, domain=domain, relational=False)
1080 except NotImplementedError:
1081 result = ConditionSet(symbol, f, domain)
1082 return result
1083 elif _is_modular(f, symbol):
1084 result = _solve_modular(f, symbol, domain)
1085 else:
1086 lhs, rhs_s = inverter(f, 0, symbol)
1087 if lhs == symbol:
1088 # do some very minimal simplification since
1089 # repeated inversion may have left the result
1090 # in a state that other solvers (e.g. poly)
1091 # would have simplified; this is done here
1092 # rather than in the inverter since here it
1093 # is only done once whereas there it would
1094 # be repeated for each step of the inversion
1095 if isinstance(rhs_s, FiniteSet):
1096 rhs_s = FiniteSet(*[Mul(*
1097 signsimp(i).as_content_primitive())
1098 for i in rhs_s])
1099 result = rhs_s
1101 elif isinstance(rhs_s, FiniteSet):
1102 for equation in [lhs - rhs for rhs in rhs_s]:
1103 if equation == f:
1104 u = unrad(f, symbol)
1105 if u:
1106 result += _solve_radical(equation, u,
1107 symbol,
1108 solver)
1109 elif equation.has(Abs):
1110 result += _solve_abs(f, symbol, domain)
1111 else:
1112 result_rational = _solve_as_rational(equation, symbol, domain)
1113 if not isinstance(result_rational, ConditionSet):
1114 result += result_rational
1115 else:
1116 # may be a transcendental type equation
1117 t_result = _transolve(equation, symbol, domain)
1118 if isinstance(t_result, ConditionSet):
1119 # might need factoring; this is expensive so we
1120 # have delayed until now. To avoid recursion
1121 # errors look for a non-trivial factoring into
1122 # a product of symbol dependent terms; I think
1123 # that something that factors as a Pow would
1124 # have already been recognized by now.
1125 factored = equation.factor()
1126 if factored.is_Mul and equation != factored:
1127 _, dep = factored.as_independent(symbol)
1128 if not dep.is_Add:
1129 # non-trivial factoring of equation
1130 # but use form with constants
1131 # in case they need special handling
1132 t_results = []
1133 for fac in Mul.make_args(factored):
1134 if fac.has(symbol):
1135 t_results.append(solver(fac, symbol))
1136 t_result = Union(*t_results)
1137 result += t_result
1138 else:
1139 result += solver(equation, symbol)
1141 elif rhs_s is not S.EmptySet:
1142 result = ConditionSet(symbol, Eq(f, 0), domain)
1144 if isinstance(result, ConditionSet):
1145 if isinstance(f, Expr):
1146 num, den = f.as_numer_denom()
1147 if den.has(symbol):
1148 _result = _solveset(num, symbol, domain)
1149 if not isinstance(_result, ConditionSet):
1150 singularities = _solveset(den, symbol, domain)
1151 result = _result - singularities
1153 if _check:
1154 if isinstance(result, ConditionSet):
1155 # it wasn't solved or has enumerated all conditions
1156 # -- leave it alone
1157 return result
1159 # whittle away all but the symbol-containing core
1160 # to use this for testing
1161 if isinstance(orig_f, Expr):
1162 fx = orig_f.as_independent(symbol, as_Add=True)[1]
1163 fx = fx.as_independent(symbol, as_Add=False)[1]
1164 else:
1165 fx = orig_f
1167 if isinstance(result, FiniteSet):
1168 # check the result for invalid solutions
1169 result = FiniteSet(*[s for s in result
1170 if isinstance(s, RootOf)
1171 or domain_check(fx, symbol, s)])
1173 return result
1176def _is_modular(f, symbol):
1177 """
1178 Helper function to check below mentioned types of modular equations.
1179 ``A - Mod(B, C) = 0``
1181 A -> This can or cannot be a function of symbol.
1182 B -> This is surely a function of symbol.
1183 C -> It is an integer.
1185 Parameters
1186 ==========
1188 f : Expr
1189 The equation to be checked.
1191 symbol : Symbol
1192 The concerned variable for which the equation is to be checked.
1194 Examples
1195 ========
1197 >>> from sympy import symbols, exp, Mod
1198 >>> from sympy.solvers.solveset import _is_modular as check
1199 >>> x, y = symbols('x y')
1200 >>> check(Mod(x, 3) - 1, x)
1201 True
1202 >>> check(Mod(x, 3) - 1, y)
1203 False
1204 >>> check(Mod(x, 3)**2 - 5, x)
1205 False
1206 >>> check(Mod(x, 3)**2 - y, x)
1207 False
1208 >>> check(exp(Mod(x, 3)) - 1, x)
1209 False
1210 >>> check(Mod(3, y) - 1, y)
1211 False
1212 """
1214 if not f.has(Mod):
1215 return False
1217 # extract modterms from f.
1218 modterms = list(f.atoms(Mod))
1220 return (len(modterms) == 1 and # only one Mod should be present
1221 modterms[0].args[0].has(symbol) and # B-> function of symbol
1222 modterms[0].args[1].is_integer and # C-> to be an integer.
1223 any(isinstance(term, Mod)
1224 for term in list(_term_factors(f))) # free from other funcs
1225 )
1228def _invert_modular(modterm, rhs, n, symbol):
1229 """
1230 Helper function to invert modular equation.
1231 ``Mod(a, m) - rhs = 0``
1233 Generally it is inverted as (a, ImageSet(Lambda(n, m*n + rhs), S.Integers)).
1234 More simplified form will be returned if possible.
1236 If it is not invertible then (modterm, rhs) is returned.
1238 The following cases arise while inverting equation ``Mod(a, m) - rhs = 0``:
1240 1. If a is symbol then m*n + rhs is the required solution.
1242 2. If a is an instance of ``Add`` then we try to find two symbol independent
1243 parts of a and the symbol independent part gets transferred to the other
1244 side and again the ``_invert_modular`` is called on the symbol
1245 dependent part.
1247 3. If a is an instance of ``Mul`` then same as we done in ``Add`` we separate
1248 out the symbol dependent and symbol independent parts and transfer the
1249 symbol independent part to the rhs with the help of invert and again the
1250 ``_invert_modular`` is called on the symbol dependent part.
1252 4. If a is an instance of ``Pow`` then two cases arise as following:
1254 - If a is of type (symbol_indep)**(symbol_dep) then the remainder is
1255 evaluated with the help of discrete_log function and then the least
1256 period is being found out with the help of totient function.
1257 period*n + remainder is the required solution in this case.
1258 For reference: (https://en.wikipedia.org/wiki/Euler's_theorem)
1260 - If a is of type (symbol_dep)**(symbol_indep) then we try to find all
1261 primitive solutions list with the help of nthroot_mod function.
1262 m*n + rem is the general solution where rem belongs to solutions list
1263 from nthroot_mod function.
1265 Parameters
1266 ==========
1268 modterm, rhs : Expr
1269 The modular equation to be inverted, ``modterm - rhs = 0``
1271 symbol : Symbol
1272 The variable in the equation to be inverted.
1274 n : Dummy
1275 Dummy variable for output g_n.
1277 Returns
1278 =======
1280 A tuple (f_x, g_n) is being returned where f_x is modular independent function
1281 of symbol and g_n being set of values f_x can have.
1283 Examples
1284 ========
1286 >>> from sympy import symbols, exp, Mod, Dummy, S
1287 >>> from sympy.solvers.solveset import _invert_modular as invert_modular
1288 >>> x, y = symbols('x y')
1289 >>> n = Dummy('n')
1290 >>> invert_modular(Mod(exp(x), 7), S(5), n, x)
1291 (Mod(exp(x), 7), 5)
1292 >>> invert_modular(Mod(x, 7), S(5), n, x)
1293 (x, ImageSet(Lambda(_n, 7*_n + 5), Integers))
1294 >>> invert_modular(Mod(3*x + 8, 7), S(5), n, x)
1295 (x, ImageSet(Lambda(_n, 7*_n + 6), Integers))
1296 >>> invert_modular(Mod(x**4, 7), S(5), n, x)
1297 (x, EmptySet)
1298 >>> invert_modular(Mod(2**(x**2 + x + 1), 7), S(2), n, x)
1299 (x**2 + x + 1, ImageSet(Lambda(_n, 3*_n + 1), Naturals0))
1301 """
1302 a, m = modterm.args
1304 if rhs.is_real is False or any(term.is_real is False
1305 for term in list(_term_factors(a))):
1306 # Check for complex arguments
1307 return modterm, rhs
1309 if abs(rhs) >= abs(m):
1310 # if rhs has value greater than value of m.
1311 return symbol, S.EmptySet
1313 if a == symbol:
1314 return symbol, ImageSet(Lambda(n, m*n + rhs), S.Integers)
1316 if a.is_Add:
1317 # g + h = a
1318 g, h = a.as_independent(symbol)
1319 if g is not S.Zero:
1320 x_indep_term = rhs - Mod(g, m)
1321 return _invert_modular(Mod(h, m), Mod(x_indep_term, m), n, symbol)
1323 if a.is_Mul:
1324 # g*h = a
1325 g, h = a.as_independent(symbol)
1326 if g is not S.One:
1327 x_indep_term = rhs*invert(g, m)
1328 return _invert_modular(Mod(h, m), Mod(x_indep_term, m), n, symbol)
1330 if a.is_Pow:
1331 # base**expo = a
1332 base, expo = a.args
1333 if expo.has(symbol) and not base.has(symbol):
1334 # remainder -> solution independent of n of equation.
1335 # m, rhs are made coprime by dividing igcd(m, rhs)
1336 try:
1337 remainder = discrete_log(m / igcd(m, rhs), rhs, a.base)
1338 except ValueError: # log does not exist
1339 return modterm, rhs
1340 # period -> coefficient of n in the solution and also referred as
1341 # the least period of expo in which it is repeats itself.
1342 # (a**(totient(m)) - 1) divides m. Here is link of theorem:
1343 # (https://en.wikipedia.org/wiki/Euler's_theorem)
1344 period = totient(m)
1345 for p in divisors(period):
1346 # there might a lesser period exist than totient(m).
1347 if pow(a.base, p, m / igcd(m, a.base)) == 1:
1348 period = p
1349 break
1350 # recursion is not applied here since _invert_modular is currently
1351 # not smart enough to handle infinite rhs as here expo has infinite
1352 # rhs = ImageSet(Lambda(n, period*n + remainder), S.Naturals0).
1353 return expo, ImageSet(Lambda(n, period*n + remainder), S.Naturals0)
1354 elif base.has(symbol) and not expo.has(symbol):
1355 try:
1356 remainder_list = nthroot_mod(rhs, expo, m, all_roots=True)
1357 if remainder_list == []:
1358 return symbol, S.EmptySet
1359 except (ValueError, NotImplementedError):
1360 return modterm, rhs
1361 g_n = S.EmptySet
1362 for rem in remainder_list:
1363 g_n += ImageSet(Lambda(n, m*n + rem), S.Integers)
1364 return base, g_n
1366 return modterm, rhs
1369def _solve_modular(f, symbol, domain):
1370 r"""
1371 Helper function for solving modular equations of type ``A - Mod(B, C) = 0``,
1372 where A can or cannot be a function of symbol, B is surely a function of
1373 symbol and C is an integer.
1375 Currently ``_solve_modular`` is only able to solve cases
1376 where A is not a function of symbol.
1378 Parameters
1379 ==========
1381 f : Expr
1382 The modular equation to be solved, ``f = 0``
1384 symbol : Symbol
1385 The variable in the equation to be solved.
1387 domain : Set
1388 A set over which the equation is solved. It has to be a subset of
1389 Integers.
1391 Returns
1392 =======
1394 A set of integer solutions satisfying the given modular equation.
1395 A ``ConditionSet`` if the equation is unsolvable.
1397 Examples
1398 ========
1400 >>> from sympy.solvers.solveset import _solve_modular as solve_modulo
1401 >>> from sympy import S, Symbol, sin, Intersection, Interval, Mod
1402 >>> x = Symbol('x')
1403 >>> solve_modulo(Mod(5*x - 8, 7) - 3, x, S.Integers)
1404 ImageSet(Lambda(_n, 7*_n + 5), Integers)
1405 >>> solve_modulo(Mod(5*x - 8, 7) - 3, x, S.Reals) # domain should be subset of integers.
1406 ConditionSet(x, Eq(Mod(5*x + 6, 7) - 3, 0), Reals)
1407 >>> solve_modulo(-7 + Mod(x, 5), x, S.Integers)
1408 EmptySet
1409 >>> solve_modulo(Mod(12**x, 21) - 18, x, S.Integers)
1410 ImageSet(Lambda(_n, 6*_n + 2), Naturals0)
1411 >>> solve_modulo(Mod(sin(x), 7) - 3, x, S.Integers) # not solvable
1412 ConditionSet(x, Eq(Mod(sin(x), 7) - 3, 0), Integers)
1413 >>> solve_modulo(3 - Mod(x, 5), x, Intersection(S.Integers, Interval(0, 100)))
1414 Intersection(ImageSet(Lambda(_n, 5*_n + 3), Integers), Range(0, 101, 1))
1415 """
1416 # extract modterm and g_y from f
1417 unsolved_result = ConditionSet(symbol, Eq(f, 0), domain)
1418 modterm = list(f.atoms(Mod))[0]
1419 rhs = -S.One*(f.subs(modterm, S.Zero))
1420 if f.as_coefficients_dict()[modterm].is_negative:
1421 # checks if coefficient of modterm is negative in main equation.
1422 rhs *= -S.One
1424 if not domain.is_subset(S.Integers):
1425 return unsolved_result
1427 if rhs.has(symbol):
1428 # TODO Case: A-> function of symbol, can be extended here
1429 # in future.
1430 return unsolved_result
1432 n = Dummy('n', integer=True)
1433 f_x, g_n = _invert_modular(modterm, rhs, n, symbol)
1435 if f_x == modterm and g_n == rhs:
1436 return unsolved_result
1438 if f_x == symbol:
1439 if domain is not S.Integers:
1440 return domain.intersect(g_n)
1441 return g_n
1443 if isinstance(g_n, ImageSet):
1444 lamda_expr = g_n.lamda.expr
1445 lamda_vars = g_n.lamda.variables
1446 base_sets = g_n.base_sets
1447 sol_set = _solveset(f_x - lamda_expr, symbol, S.Integers)
1448 if isinstance(sol_set, FiniteSet):
1449 tmp_sol = S.EmptySet
1450 for sol in sol_set:
1451 tmp_sol += ImageSet(Lambda(lamda_vars, sol), *base_sets)
1452 sol_set = tmp_sol
1453 else:
1454 sol_set = ImageSet(Lambda(lamda_vars, sol_set), *base_sets)
1455 return domain.intersect(sol_set)
1457 return unsolved_result
1460def _term_factors(f):
1461 """
1462 Iterator to get the factors of all terms present
1463 in the given equation.
1465 Parameters
1466 ==========
1467 f : Expr
1468 Equation that needs to be addressed
1470 Returns
1471 =======
1472 Factors of all terms present in the equation.
1474 Examples
1475 ========
1477 >>> from sympy import symbols
1478 >>> from sympy.solvers.solveset import _term_factors
1479 >>> x = symbols('x')
1480 >>> list(_term_factors(-2 - x**2 + x*(x + 1)))
1481 [-2, -1, x**2, x, x + 1]
1482 """
1483 for add_arg in Add.make_args(f):
1484 yield from Mul.make_args(add_arg)
1487def _solve_exponential(lhs, rhs, symbol, domain):
1488 r"""
1489 Helper function for solving (supported) exponential equations.
1491 Exponential equations are the sum of (currently) at most
1492 two terms with one or both of them having a power with a
1493 symbol-dependent exponent.
1495 For example
1497 .. math:: 5^{2x + 3} - 5^{3x - 1}
1499 .. math:: 4^{5 - 9x} - e^{2 - x}
1501 Parameters
1502 ==========
1504 lhs, rhs : Expr
1505 The exponential equation to be solved, `lhs = rhs`
1507 symbol : Symbol
1508 The variable in which the equation is solved
1510 domain : Set
1511 A set over which the equation is solved.
1513 Returns
1514 =======
1516 A set of solutions satisfying the given equation.
1517 A ``ConditionSet`` if the equation is unsolvable or
1518 if the assumptions are not properly defined, in that case
1519 a different style of ``ConditionSet`` is returned having the
1520 solution(s) of the equation with the desired assumptions.
1522 Examples
1523 ========
1525 >>> from sympy.solvers.solveset import _solve_exponential as solve_expo
1526 >>> from sympy import symbols, S
1527 >>> x = symbols('x', real=True)
1528 >>> a, b = symbols('a b')
1529 >>> solve_expo(2**x + 3**x - 5**x, 0, x, S.Reals) # not solvable
1530 ConditionSet(x, Eq(2**x + 3**x - 5**x, 0), Reals)
1531 >>> solve_expo(a**x - b**x, 0, x, S.Reals) # solvable but incorrect assumptions
1532 ConditionSet(x, (a > 0) & (b > 0), {0})
1533 >>> solve_expo(3**(2*x) - 2**(x + 3), 0, x, S.Reals)
1534 {-3*log(2)/(-2*log(3) + log(2))}
1535 >>> solve_expo(2**x - 4**x, 0, x, S.Reals)
1536 {0}
1538 * Proof of correctness of the method
1540 The logarithm function is the inverse of the exponential function.
1541 The defining relation between exponentiation and logarithm is:
1543 .. math:: {\log_b x} = y \enspace if \enspace b^y = x
1545 Therefore if we are given an equation with exponent terms, we can
1546 convert every term to its corresponding logarithmic form. This is
1547 achieved by taking logarithms and expanding the equation using
1548 logarithmic identities so that it can easily be handled by ``solveset``.
1550 For example:
1552 .. math:: 3^{2x} = 2^{x + 3}
1554 Taking log both sides will reduce the equation to
1556 .. math:: (2x)\log(3) = (x + 3)\log(2)
1558 This form can be easily handed by ``solveset``.
1559 """
1560 unsolved_result = ConditionSet(symbol, Eq(lhs - rhs, 0), domain)
1561 newlhs = powdenest(lhs)
1562 if lhs != newlhs:
1563 # it may also be advantageous to factor the new expr
1564 neweq = factor(newlhs - rhs)
1565 if neweq != (lhs - rhs):
1566 return _solveset(neweq, symbol, domain) # try again with _solveset
1568 if not (isinstance(lhs, Add) and len(lhs.args) == 2):
1569 # solving for the sum of more than two powers is possible
1570 # but not yet implemented
1571 return unsolved_result
1573 if rhs != 0:
1574 return unsolved_result
1576 a, b = list(ordered(lhs.args))
1577 a_term = a.as_independent(symbol)[1]
1578 b_term = b.as_independent(symbol)[1]
1580 a_base, a_exp = a_term.as_base_exp()
1581 b_base, b_exp = b_term.as_base_exp()
1583 if domain.is_subset(S.Reals):
1584 conditions = And(
1585 a_base > 0,
1586 b_base > 0,
1587 Eq(im(a_exp), 0),
1588 Eq(im(b_exp), 0))
1589 else:
1590 conditions = And(
1591 Ne(a_base, 0),
1592 Ne(b_base, 0))
1594 L, R = (expand_log(log(i), force=True) for i in (a, -b))
1595 solutions = _solveset(L - R, symbol, domain)
1597 return ConditionSet(symbol, conditions, solutions)
1600def _is_exponential(f, symbol):
1601 r"""
1602 Return ``True`` if one or more terms contain ``symbol`` only in
1603 exponents, else ``False``.
1605 Parameters
1606 ==========
1608 f : Expr
1609 The equation to be checked
1611 symbol : Symbol
1612 The variable in which the equation is checked
1614 Examples
1615 ========
1617 >>> from sympy import symbols, cos, exp
1618 >>> from sympy.solvers.solveset import _is_exponential as check
1619 >>> x, y = symbols('x y')
1620 >>> check(y, y)
1621 False
1622 >>> check(x**y - 1, y)
1623 True
1624 >>> check(x**y*2**y - 1, y)
1625 True
1626 >>> check(exp(x + 3) + 3**x, x)
1627 True
1628 >>> check(cos(2**x), x)
1629 False
1631 * Philosophy behind the helper
1633 The function extracts each term of the equation and checks if it is
1634 of exponential form w.r.t ``symbol``.
1635 """
1636 rv = False
1637 for expr_arg in _term_factors(f):
1638 if symbol not in expr_arg.free_symbols:
1639 continue
1640 if (isinstance(expr_arg, Pow) and
1641 symbol not in expr_arg.base.free_symbols or
1642 isinstance(expr_arg, exp)):
1643 rv = True # symbol in exponent
1644 else:
1645 return False # dependent on symbol in non-exponential way
1646 return rv
1649def _solve_logarithm(lhs, rhs, symbol, domain):
1650 r"""
1651 Helper to solve logarithmic equations which are reducible
1652 to a single instance of `\log`.
1654 Logarithmic equations are (currently) the equations that contains
1655 `\log` terms which can be reduced to a single `\log` term or
1656 a constant using various logarithmic identities.
1658 For example:
1660 .. math:: \log(x) + \log(x - 4)
1662 can be reduced to:
1664 .. math:: \log(x(x - 4))
1666 Parameters
1667 ==========
1669 lhs, rhs : Expr
1670 The logarithmic equation to be solved, `lhs = rhs`
1672 symbol : Symbol
1673 The variable in which the equation is solved
1675 domain : Set
1676 A set over which the equation is solved.
1678 Returns
1679 =======
1681 A set of solutions satisfying the given equation.
1682 A ``ConditionSet`` if the equation is unsolvable.
1684 Examples
1685 ========
1687 >>> from sympy import symbols, log, S
1688 >>> from sympy.solvers.solveset import _solve_logarithm as solve_log
1689 >>> x = symbols('x')
1690 >>> f = log(x - 3) + log(x + 3)
1691 >>> solve_log(f, 0, x, S.Reals)
1692 {-sqrt(10), sqrt(10)}
1694 * Proof of correctness
1696 A logarithm is another way to write exponent and is defined by
1698 .. math:: {\log_b x} = y \enspace if \enspace b^y = x
1700 When one side of the equation contains a single logarithm, the
1701 equation can be solved by rewriting the equation as an equivalent
1702 exponential equation as defined above. But if one side contains
1703 more than one logarithm, we need to use the properties of logarithm
1704 to condense it into a single logarithm.
1706 Take for example
1708 .. math:: \log(2x) - 15 = 0
1710 contains single logarithm, therefore we can directly rewrite it to
1711 exponential form as
1713 .. math:: x = \frac{e^{15}}{2}
1715 But if the equation has more than one logarithm as
1717 .. math:: \log(x - 3) + \log(x + 3) = 0
1719 we use logarithmic identities to convert it into a reduced form
1721 Using,
1723 .. math:: \log(a) + \log(b) = \log(ab)
1725 the equation becomes,
1727 .. math:: \log((x - 3)(x + 3))
1729 This equation contains one logarithm and can be solved by rewriting
1730 to exponents.
1731 """
1732 new_lhs = logcombine(lhs, force=True)
1733 new_f = new_lhs - rhs
1735 return _solveset(new_f, symbol, domain)
1738def _is_logarithmic(f, symbol):
1739 r"""
1740 Return ``True`` if the equation is in the form
1741 `a\log(f(x)) + b\log(g(x)) + ... + c` else ``False``.
1743 Parameters
1744 ==========
1746 f : Expr
1747 The equation to be checked
1749 symbol : Symbol
1750 The variable in which the equation is checked
1752 Returns
1753 =======
1755 ``True`` if the equation is logarithmic otherwise ``False``.
1757 Examples
1758 ========
1760 >>> from sympy import symbols, tan, log
1761 >>> from sympy.solvers.solveset import _is_logarithmic as check
1762 >>> x, y = symbols('x y')
1763 >>> check(log(x + 2) - log(x + 3), x)
1764 True
1765 >>> check(tan(log(2*x)), x)
1766 False
1767 >>> check(x*log(x), x)
1768 False
1769 >>> check(x + log(x), x)
1770 False
1771 >>> check(y + log(x), x)
1772 True
1774 * Philosophy behind the helper
1776 The function extracts each term and checks whether it is
1777 logarithmic w.r.t ``symbol``.
1778 """
1779 rv = False
1780 for term in Add.make_args(f):
1781 saw_log = False
1782 for term_arg in Mul.make_args(term):
1783 if symbol not in term_arg.free_symbols:
1784 continue
1785 if isinstance(term_arg, log):
1786 if saw_log:
1787 return False # more than one log in term
1788 saw_log = True
1789 else:
1790 return False # dependent on symbol in non-log way
1791 if saw_log:
1792 rv = True
1793 return rv
1796def _is_lambert(f, symbol):
1797 r"""
1798 If this returns ``False`` then the Lambert solver (``_solve_lambert``) will not be called.
1800 Explanation
1801 ===========
1803 Quick check for cases that the Lambert solver might be able to handle.
1805 1. Equations containing more than two operands and `symbol`s involving any of
1806 `Pow`, `exp`, `HyperbolicFunction`,`TrigonometricFunction`, `log` terms.
1808 2. In `Pow`, `exp` the exponent should have `symbol` whereas for
1809 `HyperbolicFunction`,`TrigonometricFunction`, `log` should contain `symbol`.
1811 3. For `HyperbolicFunction`,`TrigonometricFunction` the number of trigonometric functions in
1812 equation should be less than number of symbols. (since `A*cos(x) + B*sin(x) - c`
1813 is not the Lambert type).
1815 Some forms of lambert equations are:
1816 1. X**X = C
1817 2. X*(B*log(X) + D)**A = C
1818 3. A*log(B*X + A) + d*X = C
1819 4. (B*X + A)*exp(d*X + g) = C
1820 5. g*exp(B*X + h) - B*X = C
1821 6. A*D**(E*X + g) - B*X = C
1822 7. A*cos(X) + B*sin(X) - D*X = C
1823 8. A*cosh(X) + B*sinh(X) - D*X = C
1825 Where X is any variable,
1826 A, B, C, D, E are any constants,
1827 g, h are linear functions or log terms.
1829 Parameters
1830 ==========
1832 f : Expr
1833 The equation to be checked
1835 symbol : Symbol
1836 The variable in which the equation is checked
1838 Returns
1839 =======
1841 If this returns ``False`` then the Lambert solver (``_solve_lambert``) will not be called.
1843 Examples
1844 ========
1846 >>> from sympy.solvers.solveset import _is_lambert
1847 >>> from sympy import symbols, cosh, sinh, log
1848 >>> x = symbols('x')
1850 >>> _is_lambert(3*log(x) - x*log(3), x)
1851 True
1852 >>> _is_lambert(log(log(x - 3)) + log(x-3), x)
1853 True
1854 >>> _is_lambert(cosh(x) - sinh(x), x)
1855 False
1856 >>> _is_lambert((x**2 - 2*x + 1).subs(x, (log(x) + 3*x)**2 - 1), x)
1857 True
1859 See Also
1860 ========
1862 _solve_lambert
1864 """
1865 term_factors = list(_term_factors(f.expand()))
1867 # total number of symbols in equation
1868 no_of_symbols = len([arg for arg in term_factors if arg.has(symbol)])
1869 # total number of trigonometric terms in equation
1870 no_of_trig = len([arg for arg in term_factors \
1871 if arg.has(HyperbolicFunction, TrigonometricFunction)])
1873 if f.is_Add and no_of_symbols >= 2:
1874 # `log`, `HyperbolicFunction`, `TrigonometricFunction` should have symbols
1875 # and no_of_trig < no_of_symbols
1876 lambert_funcs = (log, HyperbolicFunction, TrigonometricFunction)
1877 if any(isinstance(arg, lambert_funcs)\
1878 for arg in term_factors if arg.has(symbol)):
1879 if no_of_trig < no_of_symbols:
1880 return True
1881 # here, `Pow`, `exp` exponent should have symbols
1882 elif any(isinstance(arg, (Pow, exp)) \
1883 for arg in term_factors if (arg.as_base_exp()[1]).has(symbol)):
1884 return True
1885 return False
1888def _transolve(f, symbol, domain):
1889 r"""
1890 Function to solve transcendental equations. It is a helper to
1891 ``solveset`` and should be used internally. ``_transolve``
1892 currently supports the following class of equations:
1894 - Exponential equations
1895 - Logarithmic equations
1897 Parameters
1898 ==========
1900 f : Any transcendental equation that needs to be solved.
1901 This needs to be an expression, which is assumed
1902 to be equal to ``0``.
1904 symbol : The variable for which the equation is solved.
1905 This needs to be of class ``Symbol``.
1907 domain : A set over which the equation is solved.
1908 This needs to be of class ``Set``.
1910 Returns
1911 =======
1913 Set
1914 A set of values for ``symbol`` for which ``f`` is equal to
1915 zero. An ``EmptySet`` is returned if ``f`` does not have solutions
1916 in respective domain. A ``ConditionSet`` is returned as unsolved
1917 object if algorithms to evaluate complete solution are not
1918 yet implemented.
1920 How to use ``_transolve``
1921 =========================
1923 ``_transolve`` should not be used as an independent function, because
1924 it assumes that the equation (``f``) and the ``symbol`` comes from
1925 ``solveset`` and might have undergone a few modification(s).
1926 To use ``_transolve`` as an independent function the equation (``f``)
1927 and the ``symbol`` should be passed as they would have been by
1928 ``solveset``.
1930 Examples
1931 ========
1933 >>> from sympy.solvers.solveset import _transolve as transolve
1934 >>> from sympy.solvers.solvers import _tsolve as tsolve
1935 >>> from sympy import symbols, S, pprint
1936 >>> x = symbols('x', real=True) # assumption added
1937 >>> transolve(5**(x - 3) - 3**(2*x + 1), x, S.Reals)
1938 {-(log(3) + 3*log(5))/(-log(5) + 2*log(3))}
1940 How ``_transolve`` works
1941 ========================
1943 ``_transolve`` uses two types of helper functions to solve equations
1944 of a particular class:
1946 Identifying helpers: To determine whether a given equation
1947 belongs to a certain class of equation or not. Returns either
1948 ``True`` or ``False``.
1950 Solving helpers: Once an equation is identified, a corresponding
1951 helper either solves the equation or returns a form of the equation
1952 that ``solveset`` might better be able to handle.
1954 * Philosophy behind the module
1956 The purpose of ``_transolve`` is to take equations which are not
1957 already polynomial in their generator(s) and to either recast them
1958 as such through a valid transformation or to solve them outright.
1959 A pair of helper functions for each class of supported
1960 transcendental functions are employed for this purpose. One
1961 identifies the transcendental form of an equation and the other
1962 either solves it or recasts it into a tractable form that can be
1963 solved by ``solveset``.
1964 For example, an equation in the form `ab^{f(x)} - cd^{g(x)} = 0`
1965 can be transformed to
1966 `\log(a) + f(x)\log(b) - \log(c) - g(x)\log(d) = 0`
1967 (under certain assumptions) and this can be solved with ``solveset``
1968 if `f(x)` and `g(x)` are in polynomial form.
1970 How ``_transolve`` is better than ``_tsolve``
1971 =============================================
1973 1) Better output
1975 ``_transolve`` provides expressions in a more simplified form.
1977 Consider a simple exponential equation
1979 >>> f = 3**(2*x) - 2**(x + 3)
1980 >>> pprint(transolve(f, x, S.Reals), use_unicode=False)
1981 -3*log(2)
1982 {------------------}
1983 -2*log(3) + log(2)
1984 >>> pprint(tsolve(f, x), use_unicode=False)
1985 / 3 \
1986 | --------|
1987 | log(2/9)|
1988 [-log\2 /]
1990 2) Extensible
1992 The API of ``_transolve`` is designed such that it is easily
1993 extensible, i.e. the code that solves a given class of
1994 equations is encapsulated in a helper and not mixed in with
1995 the code of ``_transolve`` itself.
1997 3) Modular
1999 ``_transolve`` is designed to be modular i.e, for every class of
2000 equation a separate helper for identification and solving is
2001 implemented. This makes it easy to change or modify any of the
2002 method implemented directly in the helpers without interfering
2003 with the actual structure of the API.
2005 4) Faster Computation
2007 Solving equation via ``_transolve`` is much faster as compared to
2008 ``_tsolve``. In ``solve``, attempts are made computing every possibility
2009 to get the solutions. This series of attempts makes solving a bit
2010 slow. In ``_transolve``, computation begins only after a particular
2011 type of equation is identified.
2013 How to add new class of equations
2014 =================================
2016 Adding a new class of equation solver is a three-step procedure:
2018 - Identify the type of the equations
2020 Determine the type of the class of equations to which they belong:
2021 it could be of ``Add``, ``Pow``, etc. types. Separate internal functions
2022 are used for each type. Write identification and solving helpers
2023 and use them from within the routine for the given type of equation
2024 (after adding it, if necessary). Something like:
2026 .. code-block:: python
2028 def add_type(lhs, rhs, x):
2029 ....
2030 if _is_exponential(lhs, x):
2031 new_eq = _solve_exponential(lhs, rhs, x)
2032 ....
2033 rhs, lhs = eq.as_independent(x)
2034 if lhs.is_Add:
2035 result = add_type(lhs, rhs, x)
2037 - Define the identification helper.
2039 - Define the solving helper.
2041 Apart from this, a few other things needs to be taken care while
2042 adding an equation solver:
2044 - Naming conventions:
2045 Name of the identification helper should be as
2046 ``_is_class`` where class will be the name or abbreviation
2047 of the class of equation. The solving helper will be named as
2048 ``_solve_class``.
2049 For example: for exponential equations it becomes
2050 ``_is_exponential`` and ``_solve_expo``.
2051 - The identifying helpers should take two input parameters,
2052 the equation to be checked and the variable for which a solution
2053 is being sought, while solving helpers would require an additional
2054 domain parameter.
2055 - Be sure to consider corner cases.
2056 - Add tests for each helper.
2057 - Add a docstring to your helper that describes the method
2058 implemented.
2059 The documentation of the helpers should identify:
2061 - the purpose of the helper,
2062 - the method used to identify and solve the equation,
2063 - a proof of correctness
2064 - the return values of the helpers
2065 """
2067 def add_type(lhs, rhs, symbol, domain):
2068 """
2069 Helper for ``_transolve`` to handle equations of
2070 ``Add`` type, i.e. equations taking the form as
2071 ``a*f(x) + b*g(x) + .... = c``.
2072 For example: 4**x + 8**x = 0
2073 """
2074 result = ConditionSet(symbol, Eq(lhs - rhs, 0), domain)
2076 # check if it is exponential type equation
2077 if _is_exponential(lhs, symbol):
2078 result = _solve_exponential(lhs, rhs, symbol, domain)
2079 # check if it is logarithmic type equation
2080 elif _is_logarithmic(lhs, symbol):
2081 result = _solve_logarithm(lhs, rhs, symbol, domain)
2083 return result
2085 result = ConditionSet(symbol, Eq(f, 0), domain)
2087 # invert_complex handles the call to the desired inverter based
2088 # on the domain specified.
2089 lhs, rhs_s = invert_complex(f, 0, symbol, domain)
2091 if isinstance(rhs_s, FiniteSet):
2092 assert (len(rhs_s.args)) == 1
2093 rhs = rhs_s.args[0]
2095 if lhs.is_Add:
2096 result = add_type(lhs, rhs, symbol, domain)
2097 else:
2098 result = rhs_s
2100 return result
2103def solveset(f, symbol=None, domain=S.Complexes):
2104 r"""Solves a given inequality or equation with set as output
2106 Parameters
2107 ==========
2109 f : Expr or a relational.
2110 The target equation or inequality
2111 symbol : Symbol
2112 The variable for which the equation is solved
2113 domain : Set
2114 The domain over which the equation is solved
2116 Returns
2117 =======
2119 Set
2120 A set of values for `symbol` for which `f` is True or is equal to
2121 zero. An :class:`~.EmptySet` is returned if `f` is False or nonzero.
2122 A :class:`~.ConditionSet` is returned as unsolved object if algorithms
2123 to evaluate complete solution are not yet implemented.
2125 ``solveset`` claims to be complete in the solution set that it returns.
2127 Raises
2128 ======
2130 NotImplementedError
2131 The algorithms to solve inequalities in complex domain are
2132 not yet implemented.
2133 ValueError
2134 The input is not valid.
2135 RuntimeError
2136 It is a bug, please report to the github issue tracker.
2139 Notes
2140 =====
2142 Python interprets 0 and 1 as False and True, respectively, but
2143 in this function they refer to solutions of an expression. So 0 and 1
2144 return the domain and EmptySet, respectively, while True and False
2145 return the opposite (as they are assumed to be solutions of relational
2146 expressions).
2149 See Also
2150 ========
2152 solveset_real: solver for real domain
2153 solveset_complex: solver for complex domain
2155 Examples
2156 ========
2158 >>> from sympy import exp, sin, Symbol, pprint, S, Eq
2159 >>> from sympy.solvers.solveset import solveset, solveset_real
2161 * The default domain is complex. Not specifying a domain will lead
2162 to the solving of the equation in the complex domain (and this
2163 is not affected by the assumptions on the symbol):
2165 >>> x = Symbol('x')
2166 >>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
2167 {2*n*I*pi | n in Integers}
2169 >>> x = Symbol('x', real=True)
2170 >>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
2171 {2*n*I*pi | n in Integers}
2173 * If you want to use ``solveset`` to solve the equation in the
2174 real domain, provide a real domain. (Using ``solveset_real``
2175 does this automatically.)
2177 >>> R = S.Reals
2178 >>> x = Symbol('x')
2179 >>> solveset(exp(x) - 1, x, R)
2180 {0}
2181 >>> solveset_real(exp(x) - 1, x)
2182 {0}
2184 The solution is unaffected by assumptions on the symbol:
2186 >>> p = Symbol('p', positive=True)
2187 >>> pprint(solveset(p**2 - 4))
2188 {-2, 2}
2190 When a :class:`~.ConditionSet` is returned, symbols with assumptions that
2191 would alter the set are replaced with more generic symbols:
2193 >>> i = Symbol('i', imaginary=True)
2194 >>> solveset(Eq(i**2 + i*sin(i), 1), i, domain=S.Reals)
2195 ConditionSet(_R, Eq(_R**2 + _R*sin(_R) - 1, 0), Reals)
2197 * Inequalities can be solved over the real domain only. Use of a complex
2198 domain leads to a NotImplementedError.
2200 >>> solveset(exp(x) > 1, x, R)
2201 Interval.open(0, oo)
2203 """
2204 f = sympify(f)
2205 symbol = sympify(symbol)
2207 if f is S.true:
2208 return domain
2210 if f is S.false:
2211 return S.EmptySet
2213 if not isinstance(f, (Expr, Relational, Number)):
2214 raise ValueError("%s is not a valid SymPy expression" % f)
2216 if not isinstance(symbol, (Expr, Relational)) and symbol is not None:
2217 raise ValueError("%s is not a valid SymPy symbol" % (symbol,))
2219 if not isinstance(domain, Set):
2220 raise ValueError("%s is not a valid domain" %(domain))
2222 free_symbols = f.free_symbols
2224 if f.has(Piecewise):
2225 f = piecewise_fold(f)
2227 if symbol is None and not free_symbols:
2228 b = Eq(f, 0)
2229 if b is S.true:
2230 return domain
2231 elif b is S.false:
2232 return S.EmptySet
2233 else:
2234 raise NotImplementedError(filldedent('''
2235 relationship between value and 0 is unknown: %s''' % b))
2237 if symbol is None:
2238 if len(free_symbols) == 1:
2239 symbol = free_symbols.pop()
2240 elif free_symbols:
2241 raise ValueError(filldedent('''
2242 The independent variable must be specified for a
2243 multivariate equation.'''))
2244 elif not isinstance(symbol, Symbol):
2245 f, s, swap = recast_to_symbols([f], [symbol])
2246 # the xreplace will be needed if a ConditionSet is returned
2247 return solveset(f[0], s[0], domain).xreplace(swap)
2249 # solveset should ignore assumptions on symbols
2250 if symbol not in _rc:
2251 x = _rc[0] if domain.is_subset(S.Reals) else _rc[1]
2252 rv = solveset(f.xreplace({symbol: x}), x, domain)
2253 # try to use the original symbol if possible
2254 try:
2255 _rv = rv.xreplace({x: symbol})
2256 except TypeError:
2257 _rv = rv
2258 if rv.dummy_eq(_rv):
2259 rv = _rv
2260 return rv
2262 # Abs has its own handling method which avoids the
2263 # rewriting property that the first piece of abs(x)
2264 # is for x >= 0 and the 2nd piece for x < 0 -- solutions
2265 # can look better if the 2nd condition is x <= 0. Since
2266 # the solution is a set, duplication of results is not
2267 # an issue, e.g. {y, -y} when y is 0 will be {0}
2268 f, mask = _masked(f, Abs)
2269 f = f.rewrite(Piecewise) # everything that's not an Abs
2270 for d, e in mask:
2271 # everything *in* an Abs
2272 e = e.func(e.args[0].rewrite(Piecewise))
2273 f = f.xreplace({d: e})
2274 f = piecewise_fold(f)
2276 return _solveset(f, symbol, domain, _check=True)
2279def solveset_real(f, symbol):
2280 return solveset(f, symbol, S.Reals)
2283def solveset_complex(f, symbol):
2284 return solveset(f, symbol, S.Complexes)
2287def _solveset_multi(eqs, syms, domains):
2288 '''Basic implementation of a multivariate solveset.
2290 For internal use (not ready for public consumption)'''
2292 rep = {}
2293 for sym, dom in zip(syms, domains):
2294 if dom is S.Reals:
2295 rep[sym] = Symbol(sym.name, real=True)
2296 eqs = [eq.subs(rep) for eq in eqs]
2297 syms = [sym.subs(rep) for sym in syms]
2299 syms = tuple(syms)
2301 if len(eqs) == 0:
2302 return ProductSet(*domains)
2304 if len(syms) == 1:
2305 sym = syms[0]
2306 domain = domains[0]
2307 solsets = [solveset(eq, sym, domain) for eq in eqs]
2308 solset = Intersection(*solsets)
2309 return ImageSet(Lambda((sym,), (sym,)), solset).doit()
2311 eqs = sorted(eqs, key=lambda eq: len(eq.free_symbols & set(syms)))
2313 for n, eq in enumerate(eqs):
2314 sols = []
2315 all_handled = True
2316 for sym in syms:
2317 if sym not in eq.free_symbols:
2318 continue
2319 sol = solveset(eq, sym, domains[syms.index(sym)])
2321 if isinstance(sol, FiniteSet):
2322 i = syms.index(sym)
2323 symsp = syms[:i] + syms[i+1:]
2324 domainsp = domains[:i] + domains[i+1:]
2325 eqsp = eqs[:n] + eqs[n+1:]
2326 for s in sol:
2327 eqsp_sub = [eq.subs(sym, s) for eq in eqsp]
2328 sol_others = _solveset_multi(eqsp_sub, symsp, domainsp)
2329 fun = Lambda((symsp,), symsp[:i] + (s,) + symsp[i:])
2330 sols.append(ImageSet(fun, sol_others).doit())
2331 else:
2332 all_handled = False
2333 if all_handled:
2334 return Union(*sols)
2337def solvify(f, symbol, domain):
2338 """Solves an equation using solveset and returns the solution in accordance
2339 with the `solve` output API.
2341 Returns
2342 =======
2344 We classify the output based on the type of solution returned by `solveset`.
2346 Solution | Output
2347 ----------------------------------------
2348 FiniteSet | list
2350 ImageSet, | list (if `f` is periodic)
2351 Union |
2353 Union | list (with FiniteSet)
2355 EmptySet | empty list
2357 Others | None
2360 Raises
2361 ======
2363 NotImplementedError
2364 A ConditionSet is the input.
2366 Examples
2367 ========
2369 >>> from sympy.solvers.solveset import solvify
2370 >>> from sympy.abc import x
2371 >>> from sympy import S, tan, sin, exp
2372 >>> solvify(x**2 - 9, x, S.Reals)
2373 [-3, 3]
2374 >>> solvify(sin(x) - 1, x, S.Reals)
2375 [pi/2]
2376 >>> solvify(tan(x), x, S.Reals)
2377 [0]
2378 >>> solvify(exp(x) - 1, x, S.Complexes)
2380 >>> solvify(exp(x) - 1, x, S.Reals)
2381 [0]
2383 """
2384 solution_set = solveset(f, symbol, domain)
2385 result = None
2386 if solution_set is S.EmptySet:
2387 result = []
2389 elif isinstance(solution_set, ConditionSet):
2390 raise NotImplementedError('solveset is unable to solve this equation.')
2392 elif isinstance(solution_set, FiniteSet):
2393 result = list(solution_set)
2395 else:
2396 period = periodicity(f, symbol)
2397 if period is not None:
2398 solutions = S.EmptySet
2399 iter_solutions = ()
2400 if isinstance(solution_set, ImageSet):
2401 iter_solutions = (solution_set,)
2402 elif isinstance(solution_set, Union):
2403 if all(isinstance(i, ImageSet) for i in solution_set.args):
2404 iter_solutions = solution_set.args
2406 for solution in iter_solutions:
2407 solutions += solution.intersect(Interval(0, period, False, True))
2409 if isinstance(solutions, FiniteSet):
2410 result = list(solutions)
2412 else:
2413 solution = solution_set.intersect(domain)
2414 if isinstance(solution, Union):
2415 # concerned about only FiniteSet with Union but not about ImageSet
2416 # if required could be extend
2417 if any(isinstance(i, FiniteSet) for i in solution.args):
2418 result = [sol for soln in solution.args \
2419 for sol in soln.args if isinstance(soln,FiniteSet)]
2420 else:
2421 return None
2423 elif isinstance(solution, FiniteSet):
2424 result += solution
2426 return result
2429###############################################################################
2430################################ LINSOLVE #####################################
2431###############################################################################
2434def linear_coeffs(eq, *syms, dict=False):
2435 """Return a list whose elements are the coefficients of the
2436 corresponding symbols in the sum of terms in ``eq``.
2437 The additive constant is returned as the last element of the
2438 list.
2440 Raises
2441 ======
2443 NonlinearError
2444 The equation contains a nonlinear term
2445 ValueError
2446 duplicate or unordered symbols are passed
2448 Parameters
2449 ==========
2451 dict - (default False) when True, return coefficients as a
2452 dictionary with coefficients keyed to syms that were present;
2453 key 1 gives the constant term
2455 Examples
2456 ========
2458 >>> from sympy.solvers.solveset import linear_coeffs
2459 >>> from sympy.abc import x, y, z
2460 >>> linear_coeffs(3*x + 2*y - 1, x, y)
2461 [3, 2, -1]
2463 It is not necessary to expand the expression:
2465 >>> linear_coeffs(x + y*(z*(x*3 + 2) + 3), x)
2466 [3*y*z + 1, y*(2*z + 3)]
2468 When nonlinear is detected, an error will be raised:
2470 * even if they would cancel after expansion (so the
2471 situation does not pass silently past the caller's
2472 attention)
2474 >>> eq = 1/x*(x - 1) + 1/x
2475 >>> linear_coeffs(eq.expand(), x)
2476 [0, 1]
2477 >>> linear_coeffs(eq, x)
2478 Traceback (most recent call last):
2479 ...
2480 NonlinearError:
2481 nonlinear in given generators
2483 * when there are cross terms
2485 >>> linear_coeffs(x*(y + 1), x, y)
2486 Traceback (most recent call last):
2487 ...
2488 NonlinearError:
2489 symbol-dependent cross-terms encountered
2491 * when there are terms that contain an expression
2492 dependent on the symbols that is not linear
2494 >>> linear_coeffs(x**2, x)
2495 Traceback (most recent call last):
2496 ...
2497 NonlinearError:
2498 nonlinear in given generators
2499 """
2500 eq = _sympify(eq)
2501 if len(syms) == 1 and iterable(syms[0]) and not isinstance(syms[0], Basic):
2502 raise ValueError('expecting unpacked symbols, *syms')
2503 symset = set(syms)
2504 if len(symset) != len(syms):
2505 raise ValueError('duplicate symbols given')
2506 try:
2507 d, c = _linear_eq_to_dict([eq], symset)
2508 d = d[0]
2509 c = c[0]
2510 except PolyNonlinearError as err:
2511 raise NonlinearError(str(err))
2512 if dict:
2513 if c:
2514 d[S.One] = c
2515 return d
2516 rv = [S.Zero]*(len(syms) + 1)
2517 rv[-1] = c
2518 for i, k in enumerate(syms):
2519 if k not in d:
2520 continue
2521 rv[i] = d[k]
2522 return rv
2525def linear_eq_to_matrix(equations, *symbols):
2526 r"""
2527 Converts a given System of Equations into Matrix form.
2528 Here `equations` must be a linear system of equations in
2529 `symbols`. Element ``M[i, j]`` corresponds to the coefficient
2530 of the jth symbol in the ith equation.
2532 The Matrix form corresponds to the augmented matrix form.
2533 For example:
2535 .. math:: 4x + 2y + 3z = 1
2536 .. math:: 3x + y + z = -6
2537 .. math:: 2x + 4y + 9z = 2
2539 This system will return $A$ and $b$ as:
2541 $$ A = \left[\begin{array}{ccc}
2542 4 & 2 & 3 \\
2543 3 & 1 & 1 \\
2544 2 & 4 & 9
2545 \end{array}\right] \ \ b = \left[\begin{array}{c}
2546 1 \\ -6 \\ 2
2547 \end{array}\right] $$
2549 The only simplification performed is to convert
2550 ``Eq(a, b)`` $\Rightarrow a - b$.
2552 Raises
2553 ======
2555 NonlinearError
2556 The equations contain a nonlinear term.
2557 ValueError
2558 The symbols are not given or are not unique.
2560 Examples
2561 ========
2563 >>> from sympy import linear_eq_to_matrix, symbols
2564 >>> c, x, y, z = symbols('c, x, y, z')
2566 The coefficients (numerical or symbolic) of the symbols will
2567 be returned as matrices:
2569 >>> eqns = [c*x + z - 1 - c, y + z, x - y]
2570 >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
2571 >>> A
2572 Matrix([
2573 [c, 0, 1],
2574 [0, 1, 1],
2575 [1, -1, 0]])
2576 >>> b
2577 Matrix([
2578 [c + 1],
2579 [ 0],
2580 [ 0]])
2582 This routine does not simplify expressions and will raise an error
2583 if nonlinearity is encountered:
2585 >>> eqns = [
2586 ... (x**2 - 3*x)/(x - 3) - 3,
2587 ... y**2 - 3*y - y*(y - 4) + x - 4]
2588 >>> linear_eq_to_matrix(eqns, [x, y])
2589 Traceback (most recent call last):
2590 ...
2591 NonlinearError:
2592 symbol-dependent term can be ignored using `strict=False`
2594 Simplifying these equations will discard the removable singularity
2595 in the first and reveal the linear structure of the second:
2597 >>> [e.simplify() for e in eqns]
2598 [x - 3, x + y - 4]
2600 Any such simplification needed to eliminate nonlinear terms must
2601 be done *before* calling this routine.
2602 """
2603 if not symbols:
2604 raise ValueError(filldedent('''
2605 Symbols must be given, for which coefficients
2606 are to be found.
2607 '''))
2609 if hasattr(symbols[0], '__iter__'):
2610 symbols = symbols[0]
2612 if has_dups(symbols):
2613 raise ValueError('Symbols must be unique')
2615 equations = sympify(equations)
2616 if isinstance(equations, MatrixBase):
2617 equations = list(equations)
2618 elif isinstance(equations, (Expr, Eq)):
2619 equations = [equations]
2620 elif not is_sequence(equations):
2621 raise ValueError(filldedent('''
2622 Equation(s) must be given as a sequence, Expr,
2623 Eq or Matrix.
2624 '''))
2626 # construct the dictionaries
2627 try:
2628 eq, c = _linear_eq_to_dict(equations, symbols)
2629 except PolyNonlinearError as err:
2630 raise NonlinearError(str(err))
2631 # prepare output matrices
2632 n, m = shape = len(eq), len(symbols)
2633 ix = dict(zip(symbols, range(m)))
2634 A = zeros(*shape)
2635 for row, d in enumerate(eq):
2636 for k in d:
2637 col = ix[k]
2638 A[row, col] = d[k]
2639 b = Matrix(n, 1, [-i for i in c])
2640 return A, b
2643def linsolve(system, *symbols):
2644 r"""
2645 Solve system of $N$ linear equations with $M$ variables; both
2646 underdetermined and overdetermined systems are supported.
2647 The possible number of solutions is zero, one or infinite.
2648 Zero solutions throws a ValueError, whereas infinite
2649 solutions are represented parametrically in terms of the given
2650 symbols. For unique solution a :class:`~.FiniteSet` of ordered tuples
2651 is returned.
2653 All standard input formats are supported:
2654 For the given set of equations, the respective input types
2655 are given below:
2657 .. math:: 3x + 2y - z = 1
2658 .. math:: 2x - 2y + 4z = -2
2659 .. math:: 2x - y + 2z = 0
2661 * Augmented matrix form, ``system`` given below:
2663 $$ \text{system} = \left[{array}{cccc}
2664 3 & 2 & -1 & 1\\
2665 2 & -2 & 4 & -2\\
2666 2 & -1 & 2 & 0
2667 \end{array}\right] $$
2669 ::
2671 system = Matrix([[3, 2, -1, 1], [2, -2, 4, -2], [2, -1, 2, 0]])
2673 * List of equations form
2675 ::
2677 system = [3x + 2y - z - 1, 2x - 2y + 4z + 2, 2x - y + 2z]
2679 * Input $A$ and $b$ in matrix form (from $Ax = b$) are given as:
2681 $$ A = \left[\begin{array}{ccc}
2682 3 & 2 & -1 \\
2683 2 & -2 & 4 \\
2684 2 & -1 & 2
2685 \end{array}\right] \ \ b = \left[\begin{array}{c}
2686 1 \\ -2 \\ 0
2687 \end{array}\right] $$
2689 ::
2691 A = Matrix([[3, 2, -1], [2, -2, 4], [2, -1, 2]])
2692 b = Matrix([[1], [-2], [0]])
2693 system = (A, b)
2695 Symbols can always be passed but are actually only needed
2696 when 1) a system of equations is being passed and 2) the
2697 system is passed as an underdetermined matrix and one wants
2698 to control the name of the free variables in the result.
2699 An error is raised if no symbols are used for case 1, but if
2700 no symbols are provided for case 2, internally generated symbols
2701 will be provided. When providing symbols for case 2, there should
2702 be at least as many symbols are there are columns in matrix A.
2704 The algorithm used here is Gauss-Jordan elimination, which
2705 results, after elimination, in a row echelon form matrix.
2707 Returns
2708 =======
2710 A FiniteSet containing an ordered tuple of values for the
2711 unknowns for which the `system` has a solution. (Wrapping
2712 the tuple in FiniteSet is used to maintain a consistent
2713 output format throughout solveset.)
2715 Returns EmptySet, if the linear system is inconsistent.
2717 Raises
2718 ======
2720 ValueError
2721 The input is not valid.
2722 The symbols are not given.
2724 Examples
2725 ========
2727 >>> from sympy import Matrix, linsolve, symbols
2728 >>> x, y, z = symbols("x, y, z")
2729 >>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
2730 >>> b = Matrix([3, 6, 9])
2731 >>> A
2732 Matrix([
2733 [1, 2, 3],
2734 [4, 5, 6],
2735 [7, 8, 10]])
2736 >>> b
2737 Matrix([
2738 [3],
2739 [6],
2740 [9]])
2741 >>> linsolve((A, b), [x, y, z])
2742 {(-1, 2, 0)}
2744 * Parametric Solution: In case the system is underdetermined, the
2745 function will return a parametric solution in terms of the given
2746 symbols. Those that are free will be returned unchanged. e.g. in
2747 the system below, `z` is returned as the solution for variable z;
2748 it can take on any value.
2750 >>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
2751 >>> b = Matrix([3, 6, 9])
2752 >>> linsolve((A, b), x, y, z)
2753 {(z - 1, 2 - 2*z, z)}
2755 If no symbols are given, internally generated symbols will be used.
2756 The ``tau0`` in the third position indicates (as before) that the third
2757 variable -- whatever it is named -- can take on any value:
2759 >>> linsolve((A, b))
2760 {(tau0 - 1, 2 - 2*tau0, tau0)}
2762 * List of equations as input
2764 >>> Eqns = [3*x + 2*y - z - 1, 2*x - 2*y + 4*z + 2, - x + y/2 - z]
2765 >>> linsolve(Eqns, x, y, z)
2766 {(1, -2, -2)}
2768 * Augmented matrix as input
2770 >>> aug = Matrix([[2, 1, 3, 1], [2, 6, 8, 3], [6, 8, 18, 5]])
2771 >>> aug
2772 Matrix([
2773 [2, 1, 3, 1],
2774 [2, 6, 8, 3],
2775 [6, 8, 18, 5]])
2776 >>> linsolve(aug, x, y, z)
2777 {(3/10, 2/5, 0)}
2779 * Solve for symbolic coefficients
2781 >>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
2782 >>> eqns = [a*x + b*y - c, d*x + e*y - f]
2783 >>> linsolve(eqns, x, y)
2784 {((-b*f + c*e)/(a*e - b*d), (a*f - c*d)/(a*e - b*d))}
2786 * A degenerate system returns solution as set of given
2787 symbols.
2789 >>> system = Matrix(([0, 0, 0], [0, 0, 0], [0, 0, 0]))
2790 >>> linsolve(system, x, y)
2791 {(x, y)}
2793 * For an empty system linsolve returns empty set
2795 >>> linsolve([], x)
2796 EmptySet
2798 * An error is raised if any nonlinearity is detected, even
2799 if it could be removed with expansion
2801 >>> linsolve([x*(1/x - 1)], x)
2802 Traceback (most recent call last):
2803 ...
2804 NonlinearError: nonlinear term: 1/x
2806 >>> linsolve([x*(y + 1)], x, y)
2807 Traceback (most recent call last):
2808 ...
2809 NonlinearError: nonlinear cross-term: x*(y + 1)
2811 >>> linsolve([x**2 - 1], x)
2812 Traceback (most recent call last):
2813 ...
2814 NonlinearError: nonlinear term: x**2
2815 """
2816 if not system:
2817 return S.EmptySet
2819 # If second argument is an iterable
2820 if symbols and hasattr(symbols[0], '__iter__'):
2821 symbols = symbols[0]
2822 sym_gen = isinstance(symbols, GeneratorType)
2823 dup_msg = 'duplicate symbols given'
2826 b = None # if we don't get b the input was bad
2827 # unpack system
2829 if hasattr(system, '__iter__'):
2831 # 1). (A, b)
2832 if len(system) == 2 and isinstance(system[0], MatrixBase):
2833 A, b = system
2835 # 2). (eq1, eq2, ...)
2836 if not isinstance(system[0], MatrixBase):
2837 if sym_gen or not symbols:
2838 raise ValueError(filldedent('''
2839 When passing a system of equations, the explicit
2840 symbols for which a solution is being sought must
2841 be given as a sequence, too.
2842 '''))
2843 if len(set(symbols)) != len(symbols):
2844 raise ValueError(dup_msg)
2846 #
2847 # Pass to the sparse solver implemented in polys. It is important
2848 # that we do not attempt to convert the equations to a matrix
2849 # because that would be very inefficient for large sparse systems
2850 # of equations.
2851 #
2852 eqs = system
2853 eqs = [sympify(eq) for eq in eqs]
2854 try:
2855 sol = _linsolve(eqs, symbols)
2856 except PolyNonlinearError as exc:
2857 # e.g. cos(x) contains an element of the set of generators
2858 raise NonlinearError(str(exc))
2860 if sol is None:
2861 return S.EmptySet
2863 sol = FiniteSet(Tuple(*(sol.get(sym, sym) for sym in symbols)))
2864 return sol
2866 elif isinstance(system, MatrixBase) and not (
2867 symbols and not isinstance(symbols, GeneratorType) and
2868 isinstance(symbols[0], MatrixBase)):
2869 # 3). A augmented with b
2870 A, b = system[:, :-1], system[:, -1:]
2872 if b is None:
2873 raise ValueError("Invalid arguments")
2874 if sym_gen:
2875 symbols = [next(symbols) for i in range(A.cols)]
2876 symset = set(symbols)
2877 if any(symset & (A.free_symbols | b.free_symbols)):
2878 raise ValueError(filldedent('''
2879 At least one of the symbols provided
2880 already appears in the system to be solved.
2881 One way to avoid this is to use Dummy symbols in
2882 the generator, e.g. numbered_symbols('%s', cls=Dummy)
2883 ''' % symbols[0].name.rstrip('1234567890')))
2884 elif len(symset) != len(symbols):
2885 raise ValueError(dup_msg)
2887 if not symbols:
2888 symbols = [Dummy() for _ in range(A.cols)]
2889 name = _uniquely_named_symbol('tau', (A, b),
2890 compare=lambda i: str(i).rstrip('1234567890')).name
2891 gen = numbered_symbols(name)
2892 else:
2893 gen = None
2895 # This is just a wrapper for solve_lin_sys
2896 eqs = []
2897 rows = A.tolist()
2898 for rowi, bi in zip(rows, b):
2899 terms = [elem * sym for elem, sym in zip(rowi, symbols) if elem]
2900 terms.append(-bi)
2901 eqs.append(Add(*terms))
2903 eqs, ring = sympy_eqs_to_ring(eqs, symbols)
2904 sol = solve_lin_sys(eqs, ring, _raw=False)
2905 if sol is None:
2906 return S.EmptySet
2907 #sol = {sym:val for sym, val in sol.items() if sym != val}
2908 sol = FiniteSet(Tuple(*(sol.get(sym, sym) for sym in symbols)))
2910 if gen is not None:
2911 solsym = sol.free_symbols
2912 rep = {sym: next(gen) for sym in symbols if sym in solsym}
2913 sol = sol.subs(rep)
2915 return sol
2918##############################################################################
2919# ------------------------------nonlinsolve ---------------------------------#
2920##############################################################################
2923def _return_conditionset(eqs, symbols):
2924 # return conditionset
2925 eqs = (Eq(lhs, 0) for lhs in eqs)
2926 condition_set = ConditionSet(
2927 Tuple(*symbols), And(*eqs), S.Complexes**len(symbols))
2928 return condition_set
2931def substitution(system, symbols, result=[{}], known_symbols=[],
2932 exclude=[], all_symbols=None):
2933 r"""
2934 Solves the `system` using substitution method. It is used in
2935 :func:`~.nonlinsolve`. This will be called from :func:`~.nonlinsolve` when any
2936 equation(s) is non polynomial equation.
2938 Parameters
2939 ==========
2941 system : list of equations
2942 The target system of equations
2943 symbols : list of symbols to be solved.
2944 The variable(s) for which the system is solved
2945 known_symbols : list of solved symbols
2946 Values are known for these variable(s)
2947 result : An empty list or list of dict
2948 If No symbol values is known then empty list otherwise
2949 symbol as keys and corresponding value in dict.
2950 exclude : Set of expression.
2951 Mostly denominator expression(s) of the equations of the system.
2952 Final solution should not satisfy these expressions.
2953 all_symbols : known_symbols + symbols(unsolved).
2955 Returns
2956 =======
2958 A FiniteSet of ordered tuple of values of `all_symbols` for which the
2959 `system` has solution. Order of values in the tuple is same as symbols
2960 present in the parameter `all_symbols`. If parameter `all_symbols` is None
2961 then same as symbols present in the parameter `symbols`.
2963 Please note that general FiniteSet is unordered, the solution returned
2964 here is not simply a FiniteSet of solutions, rather it is a FiniteSet of
2965 ordered tuple, i.e. the first & only argument to FiniteSet is a tuple of
2966 solutions, which is ordered, & hence the returned solution is ordered.
2968 Also note that solution could also have been returned as an ordered tuple,
2969 FiniteSet is just a wrapper `{}` around the tuple. It has no other
2970 significance except for the fact it is just used to maintain a consistent
2971 output format throughout the solveset.
2973 Raises
2974 ======
2976 ValueError
2977 The input is not valid.
2978 The symbols are not given.
2979 AttributeError
2980 The input symbols are not :class:`~.Symbol` type.
2982 Examples
2983 ========
2985 >>> from sympy import symbols, substitution
2986 >>> x, y = symbols('x, y', real=True)
2987 >>> substitution([x + y], [x], [{y: 1}], [y], set([]), [x, y])
2988 {(-1, 1)}
2990 * When you want a soln not satisfying $x + 1 = 0$
2992 >>> substitution([x + y], [x], [{y: 1}], [y], set([x + 1]), [y, x])
2993 EmptySet
2994 >>> substitution([x + y], [x], [{y: 1}], [y], set([x - 1]), [y, x])
2995 {(1, -1)}
2996 >>> substitution([x + y - 1, y - x**2 + 5], [x, y])
2997 {(-3, 4), (2, -1)}
2999 * Returns both real and complex solution
3001 >>> x, y, z = symbols('x, y, z')
3002 >>> from sympy import exp, sin
3003 >>> substitution([exp(x) - sin(y), y**2 - 4], [x, y])
3004 {(ImageSet(Lambda(_n, I*(2*_n*pi + pi) + log(sin(2))), Integers), -2),
3005 (ImageSet(Lambda(_n, 2*_n*I*pi + log(sin(2))), Integers), 2)}
3007 >>> eqs = [z**2 + exp(2*x) - sin(y), -3 + exp(-y)]
3008 >>> substitution(eqs, [y, z])
3009 {(-log(3), -sqrt(-exp(2*x) - sin(log(3)))),
3010 (-log(3), sqrt(-exp(2*x) - sin(log(3)))),
3011 (ImageSet(Lambda(_n, 2*_n*I*pi - log(3)), Integers),
3012 ImageSet(Lambda(_n, -sqrt(-exp(2*x) + sin(2*_n*I*pi - log(3)))), Integers)),
3013 (ImageSet(Lambda(_n, 2*_n*I*pi - log(3)), Integers),
3014 ImageSet(Lambda(_n, sqrt(-exp(2*x) + sin(2*_n*I*pi - log(3)))), Integers))}
3016 """
3018 if not system:
3019 return S.EmptySet
3021 if not symbols:
3022 msg = ('Symbols must be given, for which solution of the '
3023 'system is to be found.')
3024 raise ValueError(filldedent(msg))
3026 if not is_sequence(symbols):
3027 msg = ('symbols should be given as a sequence, e.g. a list.'
3028 'Not type %s: %s')
3029 raise TypeError(filldedent(msg % (type(symbols), symbols)))
3031 if not getattr(symbols[0], 'is_Symbol', False):
3032 msg = ('Iterable of symbols must be given as '
3033 'second argument, not type %s: %s')
3034 raise ValueError(filldedent(msg % (type(symbols[0]), symbols[0])))
3036 # By default `all_symbols` will be same as `symbols`
3037 if all_symbols is None:
3038 all_symbols = symbols
3040 old_result = result
3041 # storing complements and intersection for particular symbol
3042 complements = {}
3043 intersections = {}
3045 # when total_solveset_call equals total_conditionset
3046 # it means that solveset failed to solve all eqs.
3047 total_conditionset = -1
3048 total_solveset_call = -1
3050 def _unsolved_syms(eq, sort=False):
3051 """Returns the unsolved symbol present
3052 in the equation `eq`.
3053 """
3054 free = eq.free_symbols
3055 unsolved = (free - set(known_symbols)) & set(all_symbols)
3056 if sort:
3057 unsolved = list(unsolved)
3058 unsolved.sort(key=default_sort_key)
3059 return unsolved
3060 # end of _unsolved_syms()
3062 # sort such that equation with the fewest potential symbols is first.
3063 # means eq with less number of variable first in the list.
3064 eqs_in_better_order = list(
3065 ordered(system, lambda _: len(_unsolved_syms(_))))
3067 def add_intersection_complement(result, intersection_dict, complement_dict):
3068 # If solveset has returned some intersection/complement
3069 # for any symbol, it will be added in the final solution.
3070 final_result = []
3071 for res in result:
3072 res_copy = res
3073 for key_res, value_res in res.items():
3074 intersect_set, complement_set = None, None
3075 for key_sym, value_sym in intersection_dict.items():
3076 if key_sym == key_res:
3077 intersect_set = value_sym
3078 for key_sym, value_sym in complement_dict.items():
3079 if key_sym == key_res:
3080 complement_set = value_sym
3081 if intersect_set or complement_set:
3082 new_value = FiniteSet(value_res)
3083 if intersect_set and intersect_set != S.Complexes:
3084 new_value = Intersection(new_value, intersect_set)
3085 if complement_set:
3086 new_value = Complement(new_value, complement_set)
3087 if new_value is S.EmptySet:
3088 res_copy = None
3089 break
3090 elif new_value.is_FiniteSet and len(new_value) == 1:
3091 res_copy[key_res] = set(new_value).pop()
3092 else:
3093 res_copy[key_res] = new_value
3095 if res_copy is not None:
3096 final_result.append(res_copy)
3097 return final_result
3098 # end of def add_intersection_complement()
3100 def _extract_main_soln(sym, sol, soln_imageset):
3101 """Separate the Complements, Intersections, ImageSet lambda expr and
3102 its base_set. This function returns the unmasks sol from different classes
3103 of sets and also returns the appended ImageSet elements in a
3104 soln_imageset (dict: where key as unmasked element and value as ImageSet).
3105 """
3106 # if there is union, then need to check
3107 # Complement, Intersection, Imageset.
3108 # Order should not be changed.
3109 if isinstance(sol, ConditionSet):
3110 # extracts any solution in ConditionSet
3111 sol = sol.base_set
3113 if isinstance(sol, Complement):
3114 # extract solution and complement
3115 complements[sym] = sol.args[1]
3116 sol = sol.args[0]
3117 # complement will be added at the end
3118 # using `add_intersection_complement` method
3120 # if there is union of Imageset or other in soln.
3121 # no testcase is written for this if block
3122 if isinstance(sol, Union):
3123 sol_args = sol.args
3124 sol = S.EmptySet
3125 # We need in sequence so append finteset elements
3126 # and then imageset or other.
3127 for sol_arg2 in sol_args:
3128 if isinstance(sol_arg2, FiniteSet):
3129 sol += sol_arg2
3130 else:
3131 # ImageSet, Intersection, complement then
3132 # append them directly
3133 sol += FiniteSet(sol_arg2)
3135 if isinstance(sol, Intersection):
3136 # Interval/Set will be at 0th index always
3137 if sol.args[0] not in (S.Reals, S.Complexes):
3138 # Sometimes solveset returns soln with intersection
3139 # S.Reals or S.Complexes. We don't consider that
3140 # intersection.
3141 intersections[sym] = sol.args[0]
3142 sol = sol.args[1]
3143 # after intersection and complement Imageset should
3144 # be checked.
3145 if isinstance(sol, ImageSet):
3146 soln_imagest = sol
3147 expr2 = sol.lamda.expr
3148 sol = FiniteSet(expr2)
3149 soln_imageset[expr2] = soln_imagest
3151 if not isinstance(sol, FiniteSet):
3152 sol = FiniteSet(sol)
3153 return sol, soln_imageset
3154 # end of def _extract_main_soln()
3156 # helper function for _append_new_soln
3157 def _check_exclude(rnew, imgset_yes):
3158 rnew_ = rnew
3159 if imgset_yes:
3160 # replace all dummy variables (Imageset lambda variables)
3161 # with zero before `checksol`. Considering fundamental soln
3162 # for `checksol`.
3163 rnew_copy = rnew.copy()
3164 dummy_n = imgset_yes[0]
3165 for key_res, value_res in rnew_copy.items():
3166 rnew_copy[key_res] = value_res.subs(dummy_n, 0)
3167 rnew_ = rnew_copy
3168 # satisfy_exclude == true if it satisfies the expr of `exclude` list.
3169 try:
3170 # something like : `Mod(-log(3), 2*I*pi)` can't be
3171 # simplified right now, so `checksol` returns `TypeError`.
3172 # when this issue is fixed this try block should be
3173 # removed. Mod(-log(3), 2*I*pi) == -log(3)
3174 satisfy_exclude = any(
3175 checksol(d, rnew_) for d in exclude)
3176 except TypeError:
3177 satisfy_exclude = None
3178 return satisfy_exclude
3179 # end of def _check_exclude()
3181 # helper function for _append_new_soln
3182 def _restore_imgset(rnew, original_imageset, newresult):
3183 restore_sym = set(rnew.keys()) & \
3184 set(original_imageset.keys())
3185 for key_sym in restore_sym:
3186 img = original_imageset[key_sym]
3187 rnew[key_sym] = img
3188 if rnew not in newresult:
3189 newresult.append(rnew)
3190 # end of def _restore_imgset()
3192 def _append_eq(eq, result, res, delete_soln, n=None):
3193 u = Dummy('u')
3194 if n:
3195 eq = eq.subs(n, 0)
3196 satisfy = eq if eq in (True, False) else checksol(u, u, eq, minimal=True)
3197 if satisfy is False:
3198 delete_soln = True
3199 res = {}
3200 else:
3201 result.append(res)
3202 return result, res, delete_soln
3204 def _append_new_soln(rnew, sym, sol, imgset_yes, soln_imageset,
3205 original_imageset, newresult, eq=None):
3206 """If `rnew` (A dict <symbol: soln>) contains valid soln
3207 append it to `newresult` list.
3208 `imgset_yes` is (base, dummy_var) if there was imageset in previously
3209 calculated result(otherwise empty tuple). `original_imageset` is dict
3210 of imageset expr and imageset from this result.
3211 `soln_imageset` dict of imageset expr and imageset of new soln.
3212 """
3213 satisfy_exclude = _check_exclude(rnew, imgset_yes)
3214 delete_soln = False
3215 # soln should not satisfy expr present in `exclude` list.
3216 if not satisfy_exclude:
3217 local_n = None
3218 # if it is imageset
3219 if imgset_yes:
3220 local_n = imgset_yes[0]
3221 base = imgset_yes[1]
3222 if sym and sol:
3223 # when `sym` and `sol` is `None` means no new
3224 # soln. In that case we will append rnew directly after
3225 # substituting original imagesets in rnew values if present
3226 # (second last line of this function using _restore_imgset)
3227 dummy_list = list(sol.atoms(Dummy))
3228 # use one dummy `n` which is in
3229 # previous imageset
3230 local_n_list = [
3231 local_n for i in range(
3232 0, len(dummy_list))]
3234 dummy_zip = zip(dummy_list, local_n_list)
3235 lam = Lambda(local_n, sol.subs(dummy_zip))
3236 rnew[sym] = ImageSet(lam, base)
3237 if eq is not None:
3238 newresult, rnew, delete_soln = _append_eq(
3239 eq, newresult, rnew, delete_soln, local_n)
3240 elif eq is not None:
3241 newresult, rnew, delete_soln = _append_eq(
3242 eq, newresult, rnew, delete_soln)
3243 elif sol in soln_imageset.keys():
3244 rnew[sym] = soln_imageset[sol]
3245 # restore original imageset
3246 _restore_imgset(rnew, original_imageset, newresult)
3247 else:
3248 newresult.append(rnew)
3249 elif satisfy_exclude:
3250 delete_soln = True
3251 rnew = {}
3252 _restore_imgset(rnew, original_imageset, newresult)
3253 return newresult, delete_soln
3254 # end of def _append_new_soln()
3256 def _new_order_result(result, eq):
3257 # separate first, second priority. `res` that makes `eq` value equals
3258 # to zero, should be used first then other result(second priority).
3259 # If it is not done then we may miss some soln.
3260 first_priority = []
3261 second_priority = []
3262 for res in result:
3263 if not any(isinstance(val, ImageSet) for val in res.values()):
3264 if eq.subs(res) == 0:
3265 first_priority.append(res)
3266 else:
3267 second_priority.append(res)
3268 if first_priority or second_priority:
3269 return first_priority + second_priority
3270 return result
3272 def _solve_using_known_values(result, solver):
3273 """Solves the system using already known solution
3274 (result contains the dict <symbol: value>).
3275 solver is :func:`~.solveset_complex` or :func:`~.solveset_real`.
3276 """
3277 # stores imageset <expr: imageset(Lambda(n, expr), base)>.
3278 soln_imageset = {}
3279 total_solvest_call = 0
3280 total_conditionst = 0
3282 # sort such that equation with the fewest potential symbols is first.
3283 # means eq with less variable first
3284 for index, eq in enumerate(eqs_in_better_order):
3285 newresult = []
3286 original_imageset = {}
3287 # if imageset expr is used to solve other symbol
3288 imgset_yes = False
3289 result = _new_order_result(result, eq)
3290 for res in result:
3291 got_symbol = set() # symbols solved in one iteration
3292 # find the imageset and use its expr.
3293 for key_res, value_res in res.items():
3294 if isinstance(value_res, ImageSet):
3295 res[key_res] = value_res.lamda.expr
3296 original_imageset[key_res] = value_res
3297 dummy_n = value_res.lamda.expr.atoms(Dummy).pop()
3298 (base,) = value_res.base_sets
3299 imgset_yes = (dummy_n, base)
3300 # update eq with everything that is known so far
3301 eq2 = eq.subs(res).expand()
3302 unsolved_syms = _unsolved_syms(eq2, sort=True)
3303 if not unsolved_syms:
3304 if res:
3305 newresult, delete_res = _append_new_soln(
3306 res, None, None, imgset_yes, soln_imageset,
3307 original_imageset, newresult, eq2)
3308 if delete_res:
3309 # `delete_res` is true, means substituting `res` in
3310 # eq2 doesn't return `zero` or deleting the `res`
3311 # (a soln) since it staisfies expr of `exclude`
3312 # list.
3313 result.remove(res)
3314 continue # skip as it's independent of desired symbols
3315 depen1, depen2 = (eq2.rewrite(Add)).as_independent(*unsolved_syms)
3316 if (depen1.has(Abs) or depen2.has(Abs)) and solver == solveset_complex:
3317 # Absolute values cannot be inverted in the
3318 # complex domain
3319 continue
3320 soln_imageset = {}
3321 for sym in unsolved_syms:
3322 not_solvable = False
3323 try:
3324 soln = solver(eq2, sym)
3325 total_solvest_call += 1
3326 soln_new = S.EmptySet
3327 if isinstance(soln, Complement):
3328 # separate solution and complement
3329 complements[sym] = soln.args[1]
3330 soln = soln.args[0]
3331 # complement will be added at the end
3332 if isinstance(soln, Intersection):
3333 # Interval will be at 0th index always
3334 if soln.args[0] != Interval(-oo, oo):
3335 # sometimes solveset returns soln
3336 # with intersection S.Reals, to confirm that
3337 # soln is in domain=S.Reals
3338 intersections[sym] = soln.args[0]
3339 soln_new += soln.args[1]
3340 soln = soln_new if soln_new else soln
3341 if index > 0 and solver == solveset_real:
3342 # one symbol's real soln, another symbol may have
3343 # corresponding complex soln.
3344 if not isinstance(soln, (ImageSet, ConditionSet)):
3345 soln += solveset_complex(eq2, sym) # might give ValueError with Abs
3346 except (NotImplementedError, ValueError):
3347 # If solveset is not able to solve equation `eq2`. Next
3348 # time we may get soln using next equation `eq2`
3349 continue
3350 if isinstance(soln, ConditionSet):
3351 if soln.base_set in (S.Reals, S.Complexes):
3352 soln = S.EmptySet
3353 # don't do `continue` we may get soln
3354 # in terms of other symbol(s)
3355 not_solvable = True
3356 total_conditionst += 1
3357 else:
3358 soln = soln.base_set
3360 if soln is not S.EmptySet:
3361 soln, soln_imageset = _extract_main_soln(
3362 sym, soln, soln_imageset)
3364 for sol in soln:
3365 # sol is not a `Union` since we checked it
3366 # before this loop
3367 sol, soln_imageset = _extract_main_soln(
3368 sym, sol, soln_imageset)
3369 sol = set(sol).pop()
3370 free = sol.free_symbols
3371 if got_symbol and any(
3372 ss in free for ss in got_symbol
3373 ):
3374 # sol depends on previously solved symbols
3375 # then continue
3376 continue
3377 rnew = res.copy()
3378 # put each solution in res and append the new result
3379 # in the new result list (solution for symbol `s`)
3380 # along with old results.
3381 for k, v in res.items():
3382 if isinstance(v, Expr) and isinstance(sol, Expr):
3383 # if any unsolved symbol is present
3384 # Then subs known value
3385 rnew[k] = v.subs(sym, sol)
3386 # and add this new solution
3387 if sol in soln_imageset.keys():
3388 # replace all lambda variables with 0.
3389 imgst = soln_imageset[sol]
3390 rnew[sym] = imgst.lamda(
3391 *[0 for i in range(0, len(
3392 imgst.lamda.variables))])
3393 else:
3394 rnew[sym] = sol
3395 newresult, delete_res = _append_new_soln(
3396 rnew, sym, sol, imgset_yes, soln_imageset,
3397 original_imageset, newresult)
3398 if delete_res:
3399 # deleting the `res` (a soln) since it staisfies
3400 # eq of `exclude` list
3401 result.remove(res)
3402 # solution got for sym
3403 if not not_solvable:
3404 got_symbol.add(sym)
3405 # next time use this new soln
3406 if newresult:
3407 result = newresult
3408 return result, total_solvest_call, total_conditionst
3409 # end def _solve_using_know_values()
3411 new_result_real, solve_call1, cnd_call1 = _solve_using_known_values(
3412 old_result, solveset_real)
3413 new_result_complex, solve_call2, cnd_call2 = _solve_using_known_values(
3414 old_result, solveset_complex)
3416 # If total_solveset_call is equal to total_conditionset
3417 # then solveset failed to solve all of the equations.
3418 # In this case we return a ConditionSet here.
3419 total_conditionset += (cnd_call1 + cnd_call2)
3420 total_solveset_call += (solve_call1 + solve_call2)
3422 if total_conditionset == total_solveset_call and total_solveset_call != -1:
3423 return _return_conditionset(eqs_in_better_order, all_symbols)
3425 # don't keep duplicate solutions
3426 filtered_complex = []
3427 for i in list(new_result_complex):
3428 for j in list(new_result_real):
3429 if i.keys() != j.keys():
3430 continue
3431 if all(a.dummy_eq(b) for a, b in zip(i.values(), j.values()) \
3432 if not (isinstance(a, int) and isinstance(b, int))):
3433 break
3434 else:
3435 filtered_complex.append(i)
3436 # overall result
3437 result = new_result_real + filtered_complex
3439 result_all_variables = []
3440 result_infinite = []
3441 for res in result:
3442 if not res:
3443 # means {None : None}
3444 continue
3445 # If length < len(all_symbols) means infinite soln.
3446 # Some or all the soln is dependent on 1 symbol.
3447 # eg. {x: y+2} then final soln {x: y+2, y: y}
3448 if len(res) < len(all_symbols):
3449 solved_symbols = res.keys()
3450 unsolved = list(filter(
3451 lambda x: x not in solved_symbols, all_symbols))
3452 for unsolved_sym in unsolved:
3453 res[unsolved_sym] = unsolved_sym
3454 result_infinite.append(res)
3455 if res not in result_all_variables:
3456 result_all_variables.append(res)
3458 if result_infinite:
3459 # we have general soln
3460 # eg : [{x: -1, y : 1}, {x : -y, y: y}] then
3461 # return [{x : -y, y : y}]
3462 result_all_variables = result_infinite
3463 if intersections or complements:
3464 result_all_variables = add_intersection_complement(
3465 result_all_variables, intersections, complements)
3467 # convert to ordered tuple
3468 result = S.EmptySet
3469 for r in result_all_variables:
3470 temp = [r[symb] for symb in all_symbols]
3471 result += FiniteSet(tuple(temp))
3472 return result
3473# end of def substitution()
3476def _solveset_work(system, symbols):
3477 soln = solveset(system[0], symbols[0])
3478 if isinstance(soln, FiniteSet):
3479 _soln = FiniteSet(*[(s,) for s in soln])
3480 return _soln
3481 else:
3482 return FiniteSet(tuple(FiniteSet(soln)))
3485def _handle_positive_dimensional(polys, symbols, denominators):
3486 from sympy.polys.polytools import groebner
3487 # substitution method where new system is groebner basis of the system
3488 _symbols = list(symbols)
3489 _symbols.sort(key=default_sort_key)
3490 basis = groebner(polys, _symbols, polys=True)
3491 new_system = []
3492 for poly_eq in basis:
3493 new_system.append(poly_eq.as_expr())
3494 result = [{}]
3495 result = substitution(
3496 new_system, symbols, result, [],
3497 denominators)
3498 return result
3499# end of def _handle_positive_dimensional()
3502def _handle_zero_dimensional(polys, symbols, system):
3503 # solve 0 dimensional poly system using `solve_poly_system`
3504 result = solve_poly_system(polys, *symbols)
3505 # May be some extra soln is added because
3506 # we used `unrad` in `_separate_poly_nonpoly`, so
3507 # need to check and remove if it is not a soln.
3508 result_update = S.EmptySet
3509 for res in result:
3510 dict_sym_value = dict(list(zip(symbols, res)))
3511 if all(checksol(eq, dict_sym_value) for eq in system):
3512 result_update += FiniteSet(res)
3513 return result_update
3514# end of def _handle_zero_dimensional()
3517def _separate_poly_nonpoly(system, symbols):
3518 polys = []
3519 polys_expr = []
3520 nonpolys = []
3521 # unrad_changed stores a list of expressions containing
3522 # radicals that were processed using unrad
3523 # this is useful if solutions need to be checked later.
3524 unrad_changed = []
3525 denominators = set()
3526 poly = None
3527 for eq in system:
3528 # Store denom expressions that contain symbols
3529 denominators.update(_simple_dens(eq, symbols))
3530 # Convert equality to expression
3531 if isinstance(eq, Equality):
3532 eq = eq.rewrite(Add)
3533 # try to remove sqrt and rational power
3534 without_radicals = unrad(simplify(eq), *symbols)
3535 if without_radicals:
3536 unrad_changed.append(eq)
3537 eq_unrad, cov = without_radicals
3538 if not cov:
3539 eq = eq_unrad
3540 if isinstance(eq, Expr):
3541 eq = eq.as_numer_denom()[0]
3542 poly = eq.as_poly(*symbols, extension=True)
3543 elif simplify(eq).is_number:
3544 continue
3545 if poly is not None:
3546 polys.append(poly)
3547 polys_expr.append(poly.as_expr())
3548 else:
3549 nonpolys.append(eq)
3550 return polys, polys_expr, nonpolys, denominators, unrad_changed
3551# end of def _separate_poly_nonpoly()
3554def _handle_poly(polys, symbols):
3555 # _handle_poly(polys, symbols) -> (poly_sol, poly_eqs)
3556 #
3557 # We will return possible solution information to nonlinsolve as well as a
3558 # new system of polynomial equations to be solved if we cannot solve
3559 # everything directly here. The new system of polynomial equations will be
3560 # a lex-order Groebner basis for the original system. The lex basis
3561 # hopefully separate some of the variables and equations and give something
3562 # easier for substitution to work with.
3564 # The format for representing solution sets in nonlinsolve and substitution
3565 # is a list of dicts. These are the special cases:
3566 no_information = [{}] # No equations solved yet
3567 no_solutions = [] # The system is inconsistent and has no solutions.
3569 # If there is no need to attempt further solution of these equations then
3570 # we return no equations:
3571 no_equations = []
3573 inexact = any(not p.domain.is_Exact for p in polys)
3574 if inexact:
3575 # The use of Groebner over RR is likely to result incorrectly in an
3576 # inconsistent Groebner basis. So, convert any float coefficients to
3577 # Rational before computing the Groebner basis.
3578 polys = [poly(nsimplify(p, rational=True)) for p in polys]
3580 # Compute a Groebner basis in grevlex order wrt the ordering given. We will
3581 # try to convert this to lex order later. Usually it seems to be more
3582 # efficient to compute a lex order basis by computing a grevlex basis and
3583 # converting to lex with fglm.
3584 basis = groebner(polys, symbols, order='grevlex', polys=False)
3586 #
3587 # No solutions (inconsistent equations)?
3588 #
3589 if 1 in basis:
3591 # No solutions:
3592 poly_sol = no_solutions
3593 poly_eqs = no_equations
3595 #
3596 # Finite number of solutions (zero-dimensional case)
3597 #
3598 elif basis.is_zero_dimensional:
3600 # Convert Groebner basis to lex ordering
3601 basis = basis.fglm('lex')
3603 # Convert polynomial coefficients back to float before calling
3604 # solve_poly_system
3605 if inexact:
3606 basis = [nfloat(p) for p in basis]
3608 # Solve the zero-dimensional case using solve_poly_system if possible.
3609 # If some polynomials have factors that cannot be solved in radicals
3610 # then this will fail. Using solve_poly_system(..., strict=True)
3611 # ensures that we either get a complete solution set in radicals or
3612 # UnsolvableFactorError will be raised.
3613 try:
3614 result = solve_poly_system(basis, *symbols, strict=True)
3615 except UnsolvableFactorError:
3616 # Failure... not fully solvable in radicals. Return the lex-order
3617 # basis for substitution to handle.
3618 poly_sol = no_information
3619 poly_eqs = list(basis)
3620 else:
3621 # Success! We have a finite solution set and solve_poly_system has
3622 # succeeded in finding all solutions. Return the solutions and also
3623 # an empty list of remaining equations to be solved.
3624 poly_sol = [dict(zip(symbols, res)) for res in result]
3625 poly_eqs = no_equations
3627 #
3628 # Infinite families of solutions (positive-dimensional case)
3629 #
3630 else:
3631 # In this case the grevlex basis cannot be converted to lex using the
3632 # fglm method and also solve_poly_system cannot solve the equations. We
3633 # would like to return a lex basis but since we can't use fglm we
3634 # compute the lex basis directly here. The time required to recompute
3635 # the basis is generally significantly less than the time required by
3636 # substitution to solve the new system.
3637 poly_sol = no_information
3638 poly_eqs = list(groebner(polys, symbols, order='lex', polys=False))
3640 if inexact:
3641 poly_eqs = [nfloat(p) for p in poly_eqs]
3643 return poly_sol, poly_eqs
3646def nonlinsolve(system, *symbols):
3647 r"""
3648 Solve system of $N$ nonlinear equations with $M$ variables, which means both
3649 under and overdetermined systems are supported. Positive dimensional
3650 system is also supported (A system with infinitely many solutions is said
3651 to be positive-dimensional). In a positive dimensional system the solution will
3652 be dependent on at least one symbol. Returns both real solution
3653 and complex solution (if they exist).
3655 Parameters
3656 ==========
3658 system : list of equations
3659 The target system of equations
3660 symbols : list of Symbols
3661 symbols should be given as a sequence eg. list
3663 Returns
3664 =======
3666 A :class:`~.FiniteSet` of ordered tuple of values of `symbols` for which the `system`
3667 has solution. Order of values in the tuple is same as symbols present in
3668 the parameter `symbols`.
3670 Please note that general :class:`~.FiniteSet` is unordered, the solution
3671 returned here is not simply a :class:`~.FiniteSet` of solutions, rather it
3672 is a :class:`~.FiniteSet` of ordered tuple, i.e. the first and only
3673 argument to :class:`~.FiniteSet` is a tuple of solutions, which is
3674 ordered, and, hence ,the returned solution is ordered.
3676 Also note that solution could also have been returned as an ordered tuple,
3677 FiniteSet is just a wrapper ``{}`` around the tuple. It has no other
3678 significance except for the fact it is just used to maintain a consistent
3679 output format throughout the solveset.
3681 For the given set of equations, the respective input types
3682 are given below:
3684 .. math:: xy - 1 = 0
3685 .. math:: 4x^2 + y^2 - 5 = 0
3687 ::
3689 system = [x*y - 1, 4*x**2 + y**2 - 5]
3690 symbols = [x, y]
3692 Raises
3693 ======
3695 ValueError
3696 The input is not valid.
3697 The symbols are not given.
3698 AttributeError
3699 The input symbols are not `Symbol` type.
3701 Examples
3702 ========
3704 >>> from sympy import symbols, nonlinsolve
3705 >>> x, y, z = symbols('x, y, z', real=True)
3706 >>> nonlinsolve([x*y - 1, 4*x**2 + y**2 - 5], [x, y])
3707 {(-1, -1), (-1/2, -2), (1/2, 2), (1, 1)}
3709 1. Positive dimensional system and complements:
3711 >>> from sympy import pprint
3712 >>> from sympy.polys.polytools import is_zero_dimensional
3713 >>> a, b, c, d = symbols('a, b, c, d', extended_real=True)
3714 >>> eq1 = a + b + c + d
3715 >>> eq2 = a*b + b*c + c*d + d*a
3716 >>> eq3 = a*b*c + b*c*d + c*d*a + d*a*b
3717 >>> eq4 = a*b*c*d - 1
3718 >>> system = [eq1, eq2, eq3, eq4]
3719 >>> is_zero_dimensional(system)
3720 False
3721 >>> pprint(nonlinsolve(system, [a, b, c, d]), use_unicode=False)
3722 -1 1 1 -1
3723 {(---, -d, -, {d} \ {0}), (-, -d, ---, {d} \ {0})}
3724 d d d d
3725 >>> nonlinsolve([(x+y)**2 - 4, x + y - 2], [x, y])
3726 {(2 - y, y)}
3728 2. If some of the equations are non-polynomial then `nonlinsolve`
3729 will call the ``substitution`` function and return real and complex solutions,
3730 if present.
3732 >>> from sympy import exp, sin
3733 >>> nonlinsolve([exp(x) - sin(y), y**2 - 4], [x, y])
3734 {(ImageSet(Lambda(_n, I*(2*_n*pi + pi) + log(sin(2))), Integers), -2),
3735 (ImageSet(Lambda(_n, 2*_n*I*pi + log(sin(2))), Integers), 2)}
3737 3. If system is non-linear polynomial and zero-dimensional then it
3738 returns both solution (real and complex solutions, if present) using
3739 :func:`~.solve_poly_system`:
3741 >>> from sympy import sqrt
3742 >>> nonlinsolve([x**2 - 2*y**2 -2, x*y - 2], [x, y])
3743 {(-2, -1), (2, 1), (-sqrt(2)*I, sqrt(2)*I), (sqrt(2)*I, -sqrt(2)*I)}
3745 4. ``nonlinsolve`` can solve some linear (zero or positive dimensional)
3746 system (because it uses the :func:`sympy.polys.polytools.groebner` function to get the
3747 groebner basis and then uses the ``substitution`` function basis as the
3748 new `system`). But it is not recommended to solve linear system using
3749 ``nonlinsolve``, because :func:`~.linsolve` is better for general linear systems.
3751 >>> nonlinsolve([x + 2*y -z - 3, x - y - 4*z + 9, y + z - 4], [x, y, z])
3752 {(3*z - 5, 4 - z, z)}
3754 5. System having polynomial equations and only real solution is
3755 solved using :func:`~.solve_poly_system`:
3757 >>> e1 = sqrt(x**2 + y**2) - 10
3758 >>> e2 = sqrt(y**2 + (-x + 10)**2) - 3
3759 >>> nonlinsolve((e1, e2), (x, y))
3760 {(191/20, -3*sqrt(391)/20), (191/20, 3*sqrt(391)/20)}
3761 >>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [x, y])
3762 {(1, 2), (1 - sqrt(5), 2 + sqrt(5)), (1 + sqrt(5), 2 - sqrt(5))}
3763 >>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [y, x])
3764 {(2, 1), (2 - sqrt(5), 1 + sqrt(5)), (2 + sqrt(5), 1 - sqrt(5))}
3766 6. It is better to use symbols instead of trigonometric functions or
3767 :class:`~.Function`. For example, replace $\sin(x)$ with a symbol, replace
3768 $f(x)$ with a symbol and so on. Get a solution from ``nonlinsolve`` and then
3769 use :func:`~.solveset` to get the value of $x$.
3771 How nonlinsolve is better than old solver ``_solve_system`` :
3772 =============================================================
3774 1. A positive dimensional system solver: nonlinsolve can return
3775 solution for positive dimensional system. It finds the
3776 Groebner Basis of the positive dimensional system(calling it as
3777 basis) then we can start solving equation(having least number of
3778 variable first in the basis) using solveset and substituting that
3779 solved solutions into other equation(of basis) to get solution in
3780 terms of minimum variables. Here the important thing is how we
3781 are substituting the known values and in which equations.
3783 2. Real and complex solutions: nonlinsolve returns both real
3784 and complex solution. If all the equations in the system are polynomial
3785 then using :func:`~.solve_poly_system` both real and complex solution is returned.
3786 If all the equations in the system are not polynomial equation then goes to
3787 ``substitution`` method with this polynomial and non polynomial equation(s),
3788 to solve for unsolved variables. Here to solve for particular variable
3789 solveset_real and solveset_complex is used. For both real and complex
3790 solution ``_solve_using_known_values`` is used inside ``substitution``
3791 (``substitution`` will be called when any non-polynomial equation is present).
3792 If a solution is valid its general solution is added to the final result.
3794 3. :class:`~.Complement` and :class:`~.Intersection` will be added:
3795 nonlinsolve maintains dict for complements and intersections. If solveset
3796 find complements or/and intersections with any interval or set during the
3797 execution of ``substitution`` function, then complement or/and
3798 intersection for that variable is added before returning final solution.
3800 """
3801 if not system:
3802 return S.EmptySet
3804 if not symbols:
3805 msg = ('Symbols must be given, for which solution of the '
3806 'system is to be found.')
3807 raise ValueError(filldedent(msg))
3809 if hasattr(symbols[0], '__iter__'):
3810 symbols = symbols[0]
3812 if not is_sequence(symbols) or not symbols:
3813 msg = ('Symbols must be given, for which solution of the '
3814 'system is to be found.')
3815 raise IndexError(filldedent(msg))
3817 symbols = list(map(_sympify, symbols))
3818 system, symbols, swap = recast_to_symbols(system, symbols)
3819 if swap:
3820 soln = nonlinsolve(system, symbols)
3821 return FiniteSet(*[tuple(i.xreplace(swap) for i in s) for s in soln])
3823 if len(system) == 1 and len(symbols) == 1:
3824 return _solveset_work(system, symbols)
3826 # main code of def nonlinsolve() starts from here
3828 polys, polys_expr, nonpolys, denominators, unrad_changed = \
3829 _separate_poly_nonpoly(system, symbols)
3831 poly_eqs = []
3832 poly_sol = [{}]
3834 if polys:
3835 poly_sol, poly_eqs = _handle_poly(polys, symbols)
3836 if poly_sol and poly_sol[0]:
3837 poly_syms = set().union(*(eq.free_symbols for eq in polys))
3838 unrad_syms = set().union(*(eq.free_symbols for eq in unrad_changed))
3839 if unrad_syms == poly_syms and unrad_changed:
3840 # if all the symbols have been solved by _handle_poly
3841 # and unrad has been used then check solutions
3842 poly_sol = [sol for sol in poly_sol if checksol(unrad_changed, sol)]
3844 # Collect together the unsolved polynomials with the non-polynomial
3845 # equations.
3846 remaining = poly_eqs + nonpolys
3848 # to_tuple converts a solution dictionary to a tuple containing the
3849 # value for each symbol
3850 to_tuple = lambda sol: tuple(sol[s] for s in symbols)
3852 if not remaining:
3853 # If there is nothing left to solve then return the solution from
3854 # solve_poly_system directly.
3855 return FiniteSet(*map(to_tuple, poly_sol))
3856 else:
3857 # Here we handle:
3858 #
3859 # 1. The Groebner basis if solve_poly_system failed.
3860 # 2. The Groebner basis in the positive-dimensional case.
3861 # 3. Any non-polynomial equations
3862 #
3863 # If solve_poly_system did succeed then we pass those solutions in as
3864 # preliminary results.
3865 subs_res = substitution(remaining, symbols, result=poly_sol, exclude=denominators)
3867 if not isinstance(subs_res, FiniteSet):
3868 return subs_res
3870 # check solutions produced by substitution. Currently, checking is done for
3871 # only those solutions which have non-Set variable values.
3872 if unrad_changed:
3873 result = [dict(zip(symbols, sol)) for sol in subs_res.args]
3874 correct_sols = [sol for sol in result if any(isinstance(v, Set) for v in sol)
3875 or checksol(unrad_changed, sol) != False]
3876 return FiniteSet(*map(to_tuple, correct_sols))
3877 else:
3878 return subs_res