Repository: sympy/sympy

Issue:
solveset needs to clip solution with constraining set
When a Piecewise is being solved by solveset, a check is not made that the solution is within the conditions under which the expression is constrained. This leads to spurious solutions as shown below. The solution should be clipped by the constraining set:

```diff
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 1ce7e20..0595502 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -797,7 +797,7 @@ def _solveset(f, symbol, domain, _check=False):
                 in_set = in_set.as_set()
             if in_set.is_Interval:
                 dom -= in_set
-            solns = solver(expr, symbol, in_set)
+            solns = solver(expr, symbol, in_set) & in_set
             result += solns
     else:
         lhs, rhs_s = inverter(f, 0, symbol)
```

current behavior

```python
>>> solveset((Abs(x+Min(x,2))-2).rewrite(Piecewise),x,S.Reals)
{-1, 0, 1}
```
after fix
```python
>>> solveset((Abs(x+Min(x,2))-2).rewrite(Piecewise),x,S.Reals)
{-1, 1}
```

On review of the situation, there should be no need for clipping the result since the domain passed to `solver` was already identified as `in_set`. Perhaps the following is better:

```diff
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 1ce7e20..ca3f6ad 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -848,13 +848,20 @@ def _solveset(f, symbol, domain, _check=False):
             # -- leave it alone
             return result

-        # whittle away all but the symbol-containing core
-        # to use this for testing
-        fx = orig_f.as_independent(symbol, as_Add=True)[1]
-        fx = fx.as_independent(symbol, as_Add=False)[1]
+        ### XXX improve this. result & domain should just work
+        # but consider the current state where an ImageSet doesn't
+        # know how to intersect with Complexes:
+        # >>> solveset(exp(x) - 1, x)
+        # Intersection(ImageSet(Lambda(_n, 2*_n*I*pi), S.Integers), Complexes(S.Reals x S.Reals, False))
+        if isinstance(result, (Interval, FiniteSet)):
+            result = result & domain

         if isinstance(result, FiniteSet):
             # check the result for invalid solutions
+            # whittle away all but the symbol-containing core
+            # to use this for testing
+            fx = orig_f.as_independent(symbol, as_Add=True)[1]
+            fx = fx.as_independent(symbol, as_Add=False)[1]
             result = FiniteSet(*[s for s in result
                       if isinstance(s, RootOf)
                       or domain_check(fx, symbol, s)])
```

Discussion:
Yathartha22: > invert_real does and simply passes it along to _invert so if one uses a domain other than S.Reals it doesn't matter if invert_complex (alias for _invert) or invert_real is called with that domain.

+1 to the point you are trying to make, I think `invert_real` and `invert_complex` should just have `S.Reals` and `S.Complexes` as the domain respectively.

Patch:
diff --git a/examples/notebooks/sho1d_example.ipynb b/examples/notebooks/sho1d_example.ipynb
index 426b413cad..52d62a62b0 100644
--- a/examples/notebooks/sho1d_example.ipynb
+++ b/examples/notebooks/sho1d_example.ipynb
@@ -340,7 +340,7 @@
      "cell_type": "markdown",
      "metadata": {},
      "source": [
-      "Take the dagger of the raising and lowering operators. They should return each other."
+      "Take the dagger of the raising and lowering operators. They should return eachother."
      ]
     },
     {
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index 0ce7def251..fe2a0280ce 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -345,7 +345,7 @@ def calc_part(re_im, nexpr):
         nint = int(to_int(nexpr, rnd))
         if is_int:
             # make sure that we had enough precision to distinguish
-            # between nint and the re or im part (re_im) of expr that
+            # bewteen nint and the re or im part (re_im) of expr that
             # was passed to calc_part
             ire, iim, ire_acc, iim_acc = evalf(
                 re_im - nint, 10, options)  # don't need much precision
diff --git a/sympy/core/power.py b/sympy/core/power.py
index bd4bb7ae15..02ec6480cb 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -404,21 +404,10 @@ def _eval_Mod(self, q):
 
             '''
             For unevaluated Integer power, use built-in pow modular
-            exponentiation, if powers are not too large wrt base.
+            exponentiation.
             '''
             if self.base.is_Integer and self.exp.is_Integer and q.is_Integer:
-                b, e, m = int(self.base), int(self.exp), int(q)
-                # For very large powers, use totient reduction if e >= lg(m).
-                # Bound on m, is for safe factorization memory wise ie m^(1/4).
-                # For pollard-rho to be faster than built-in pow lg(e) > m^(1/4)
-                # check is added.
-                mb = m.bit_length()
-                if mb <= 80  and e >= mb and e.bit_length()**4 >= m:
-                    from sympy.ntheory import totient
-                    phi = totient(m)
-                    return pow(b, phi + e%phi, m)
-                else:
-                    return pow(b, e, m)
+                return pow(int(self.base), int(self.exp), int(q))
 
     def _eval_is_even(self):
         if self.exp.is_integer and self.exp.is_positive:
diff --git a/sympy/functions/elementary/complexes.py b/sympy/functions/elementary/complexes.py
index 791d5c6779..77bc6da274 100644
--- a/sympy/functions/elementary/complexes.py
+++ b/sympy/functions/elementary/complexes.py
@@ -480,12 +480,6 @@ def eval(cls, arg):
                 if base.is_negative:
                     return (-base)**re(exponent)*exp(-S.Pi*im(exponent))
                 return
-            elif not base.has(Symbol): # complex base
-                # express base**exponent as exp(exponent*log(base))
-                a, b = log(base).as_real_imag()
-                z = a + I*b
-                return exp(re(exponent*z))
-
         if isinstance(arg, exp):
             return exp(re(arg.args[0]))
         if isinstance(arg, AppliedUndef):
@@ -505,7 +499,7 @@ def eval(cls, arg):
                 return arg2
         # reject result if all new conjugates are just wrappers around
         # an expression that was already in the arg
-        conj = arg.conjugate()
+        conj = signsimp(arg.conjugate(), evaluate=False)
         new_conj = conj.atoms(conjugate) - arg.atoms(conjugate)
         if new_conj and all(arg.has(i.args[0]) for i in new_conj):
             return
diff --git a/sympy/functions/elementary/piecewise.py b/sympy/functions/elementary/piecewise.py
index 7cb3cec826..adb6db521a 100644
--- a/sympy/functions/elementary/piecewise.py
+++ b/sympy/functions/elementary/piecewise.py
@@ -3,7 +3,8 @@
 from sympy.core import Basic, S, Function, diff, Tuple, Dummy, Number
 from sympy.core.basic import as_Basic
 from sympy.core.sympify import SympifyError
-from sympy.core.relational import Equality, Relational, _canonical
+from sympy.core.relational import (Equality, Unequality, Relational,
+    _canonical)
 from sympy.functions.elementary.miscellaneous import Max, Min
 from sympy.logic.boolalg import (And, Boolean, distribute_and_over_or,
     true, false, Not, Or, ITE, simplify_logic)
@@ -857,28 +858,42 @@ def __eval_cond(cls, cond):
             except TypeError:
                 pass
 
-    def as_expr_set_pairs(self):
+    def as_expr_set_pairs(self, domain=S.Reals):
         """Return tuples for each argument of self that give
-        the expression and the interval in which it is valid.
+        the expression and the interval in which it is valid
+        which is contained within the given domain.
         If a condition cannot be converted to a set, an error
         will be raised. The variable of the conditions is
         assumed to be real; sets of real values are returned.
 
         Examples
         ========
-        >>> from sympy import Piecewise
+
+        >>> from sympy import Piecewise, Interval
         >>> from sympy.abc import x
-        >>> Piecewise(
+        >>> p = Piecewise(
         ...     (1, x < 2),
         ...     (2,(x > 0) & (x < 4)),
-        ...     (3, True)).as_expr_set_pairs()
+        ...     (3, True))
+        >>> p.as_expr_set_pairs()
         [(1, Interval.open(-oo, 2)),
          (2, Interval.Ropen(2, 4)),
          (3, Interval(4, oo))]
+        >>> p.as_expr_set_pairs(Interval(0, 3))
+        [(1, Interval.Ropen(0, 2)),
+         (2, Interval(2, 3)), (3, EmptySet())]
         """
         exp_sets = []
-        U = S.Reals
+        U = domain
+        complex = not domain.is_subset(S.Reals)
         for expr, cond in self.args:
+            if complex:
+                for i in cond.atoms(Relational):
+                    if not isinstance(i, (Equality, Unequality)):
+                        raise ValueError(filldedent('''
+                            Inequalities in the complex domain are
+                            not supported. Try the real domain by
+                            setting domain=S.Reals'''))
             cond_int = U.intersect(cond.as_set())
             U = U - cond_int
             exp_sets.append((expr, cond_int))
diff --git a/sympy/geometry/line.py b/sympy/geometry/line.py
index 1411bdd5ac..447cbdb8a2 100644
--- a/sympy/geometry/line.py
+++ b/sympy/geometry/line.py
@@ -234,7 +234,7 @@ def smallest_angle_between(l1, l2):
 
         See Also
         ========
-        angle_between, Ray2D.closing_angle
+        angle_bewteen, Ray2D.closing_angle
         """
         if not isinstance(l1, LinearEntity) and not isinstance(l2, LinearEntity):
             raise TypeError('Must pass only LinearEntity objects')
diff --git a/sympy/integrals/rubi/parsetools/parse.py b/sympy/integrals/rubi/parsetools/parse.py
index b2581062ba..bc20114972 100644
--- a/sympy/integrals/rubi/parsetools/parse.py
+++ b/sympy/integrals/rubi/parsetools/parse.py
@@ -293,7 +293,7 @@ def replaceWith(s, symbols, i):
     if type(s) == Function('With') or type(s) == Function('Module'):
         constraints = ', '
         result = '    def With{}({}):'.format(i, ', '.join(symbols))
-        if type(s.args[0]) == Function('List'): # get all local variables of With and Module
+        if type(s.args[0]) == Function('List'): # get all local varibles of With and Module
             L = list(s.args[0].args)
         else:
             L = [s.args[0]]
diff --git a/sympy/physics/mechanics/functions.py b/sympy/physics/mechanics/functions.py
index cee733eca8..5d7806b212 100644
--- a/sympy/physics/mechanics/functions.py
+++ b/sympy/physics/mechanics/functions.py
@@ -341,45 +341,6 @@ def potential_energy(*body):
             raise TypeError('*body must have only Particle or RigidBody')
     return pe_sys
 
-def gravity(acceleration, *bodies):
-    """
-    Returns a list of gravity forces given the acceleration
-    due to gravity and any number of particles or rigidbodies.
-
-    Examples
-    ========
-
-    >>> from sympy.physics.mechanics import ReferenceFrame, Point, Particle, outer, RigidBody
-    >>> from sympy.physics.mechanics.functions import gravity
-    >>> from sympy import symbols
-    >>> N = ReferenceFrame('N')
-    >>> m, M, g = symbols('m M g')
-    >>> F1, F2 = symbols('F1 F2')
-    >>> po = Point('po')
-    >>> pa = Particle('pa', po, m)
-    >>> A = ReferenceFrame('A')
-    >>> P = Point('P')
-    >>> I = outer(A.x, A.x)
-    >>> B = RigidBody('B', P, A, M, (I, P))
-    >>> forceList = [(po, F1), (P, F2)]
-    >>> forceList.extend(gravity(g*N.y, pa, B))
-    >>> forceList
-    [(po, F1), (P, F2), (po, g*m*N.y), (P, M*g*N.y)]
-    """
-
-    gravity_force = []
-    if not bodies:
-        raise TypeError("No bodies(instances of Particle or Rigidbody) were passed.")
-
-    for e in bodies:
-        try:
-            point = e.masscenter
-        except AttributeError:
-            point = e.point
-
-        gravity_force.append((point, e.mass*acceleration))
-
-    return gravity_force
 
 def Lagrangian(frame, *body):
     """Lagrangian of a multibody system.
diff --git a/sympy/series/formal.py b/sympy/series/formal.py
index c82b30ca5e..0678633296 100644
--- a/sympy/series/formal.py
+++ b/sympy/series/formal.py
@@ -70,7 +70,7 @@ def rational_algorithm(f, x, k, order=4, full=False):
 
     By setting ``full=True``, range of admissible functions to be solved using
     ``rational_algorithm`` can be increased. This option should be used
-    carefully as it can significantly slow down the computation as ``doit`` is
+    carefully as it can signifcantly slow down the computation as ``doit`` is
     performed on the :class:`RootSum` object returned by the ``apart`` function.
     Use ``full=False`` whenever possible.
 
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 1ce7e20f5c..fc89050605 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -11,19 +11,22 @@
 
 from sympy.core.sympify import sympify
 from sympy.core import S, Pow, Dummy, pi, Expr, Wild, Mul, Equality
+from sympy.core.containers import Tuple
+from sympy.core.facts import InconsistentAssumptions
 from sympy.core.numbers import I, Number, Rational, oo
-from sympy.core.function import (Lambda, expand_complex, AppliedUndef)
+from sympy.core.function import (Lambda, expand_complex, AppliedUndef, Function)
 from sympy.core.relational import Eq
 from sympy.core.symbol import Symbol
 from sympy.simplify.simplify import simplify, fraction, trigsimp
 from sympy.functions import (log, Abs, tan, cot, sin, cos, sec, csc, exp,
                              acos, asin, acsc, asec, arg,
-                             piecewise_fold)
+                             piecewise_fold, Piecewise)
 from sympy.functions.elementary.trigonometric import (TrigonometricFunction,
                                                       HyperbolicFunction)
-from sympy.functions.elementary.miscellaneous import real_root
+from sympy.functions.elementary.miscellaneous import real_root, Application
 from sympy.sets import (FiniteSet, EmptySet, imageset, Interval, Intersection,
                         Union, ConditionSet, ImageSet, Complement)
+from sympy.sets.sets import Set
 from sympy.matrices import Matrix
 from sympy.polys import (roots, Poly, degree, together, PolynomialError,
                          RootOf)
@@ -32,12 +35,54 @@
 from sympy.solvers.polysys import solve_poly_system
 from sympy.solvers.inequalities import solve_univariate_inequality
 from sympy.utilities import filldedent
+from sympy.utilities.iterables import numbered_symbols
 from sympy.calculus.util import periodicity, continuous_domain
 from sympy.core.compatibility import ordered, default_sort_key, is_sequence
 
 from types import GeneratorType
 
 
+def _masked(f, *atoms):
+    """Return ``f``, with all objects given by ``atoms`` replaced with
+    Dummy symbols, ``d``, and the list of replacements, ``(d, e)``,
+    where ``e`` is an object of type given by ``atoms`` in which
+    any other instances of atoms have been recursively replaced with
+    Dummy symbols, too. The tuples are ordered so that if they are
+    applied in sequence, the orgin ``f`` will be restored.
+
+    Examples
+    ========
+
+    >>> from sympy import cos
+    >>> from sympy.abc import x
+    >>> from sympy.solvers.solveset import _masked
+
+    >>> f = cos(cos(x) + 1)
+    >>> f, reps = _masked(cos(1 + cos(x)), cos)
+    >>> f
+    _a1
+    >>> reps
+    [(_a1, cos(_a0 + 1)), (_a0, cos(x))]
+    >>> for d, e in reps:
+    ...     f = f.xreplace({d: e})
+    >>> f
+    cos(cos(x) + 1)
+    """
+    sym = numbered_symbols('a', cls=Dummy)
+    mask = []
+    for a in ordered(f.atoms(*atoms)):
+        for i in mask:
+            a = a.replace(*i)
+        mask.append((a, next(sym)))
+    mask = list(reversed(mask))
+    for i, (o, n) in enumerate(mask):
+        f = f.replace(o, n)
+        mask[i] = (n, o)
+    if mask and f == mask[0][1]:
+        f = mask[0][0]
+    return f, mask
+
+
 def _invert(f_x, y, x, domain=S.Complexes):
     r"""
     Reduce the complex valued equation ``f(x) = y`` to a set of equations
@@ -105,7 +150,7 @@ def _invert(f_x, y, x, domain=S.Complexes):
     else:
         x1, s = _invert_complex(f_x, FiniteSet(y), x)
 
-    if not isinstance(s, FiniteSet) or x1 == f_x or x1 != x:
+    if not isinstance(s, FiniteSet) or x1 != x:
         return x1, s
 
     return x1, s.intersection(domain)
@@ -757,15 +802,9 @@ def _solveset(f, symbol, domain, _check=False):
                               S.NegativeInfinity]):
             f = a/m + h  # XXX condition `m != 0` should be added to soln
 
-    f = piecewise_fold(f)
-
     # assign the solvers to use
     solver = lambda f, x, domain=domain: _solveset(f, x, domain)
-    if domain.is_subset(S.Reals):
-        inverter_func = invert_real
-    else:
-        inverter_func = invert_complex
-    inverter = lambda f, rhs, symbol: inverter_func(f, rhs, symbol, domain)
+    inverter = lambda f, rhs, symbol: _invert(f, rhs, symbol, domain)
 
     result = EmptySet()
 
@@ -789,16 +828,28 @@ def _solveset(f, symbol, domain, _check=False):
         a = f.args[0]
         result = solveset_real(a > 0, symbol)
     elif f.is_Piecewise:
-        dom = domain
         result = EmptySet()
-        expr_set_pairs = f.as_expr_set_pairs()
+        expr_set_pairs = f.as_expr_set_pairs(domain)
         for (expr, in_set) in expr_set_pairs:
             if in_set.is_Relational:
                 in_set = in_set.as_set()
-            if in_set.is_Interval:
-                dom -= in_set
             solns = solver(expr, symbol, in_set)
             result += solns
+    elif isinstance(f, Eq):
+        from sympy.core import Add
+        result = solver(Add(f.lhs, - f.rhs, evaluate=False), symbol, domain)
+    elif f.is_Relational:
+        if not domain.is_subset(S.Reals):
+            raise NotImplementedError(filldedent('''
+                Inequalities in the complex domain are
+                not supported. Try the real domain by
+                setting domain=S.Reals'''))
+        try:
+            result = solve_univariate_inequality(
+            f, symbol, domain=domain, relational=False)
+        except NotImplementedError:
+            result = ConditionSet(symbol, f, domain)
+        return result
     else:
         lhs, rhs_s = inverter(f, 0, symbol)
         if lhs == symbol:
@@ -961,6 +1012,7 @@ def solveset(f, symbol=None, domain=S.Complexes):
 
     """
     f = sympify(f)
+    symbol = sympify(symbol)
 
     if f is S.true:
         return domain
@@ -969,10 +1021,26 @@ def solveset(f, symbol=None, domain=S.Complexes):
         return S.EmptySet
 
     if not isinstance(f, (Expr, Number)):
-        raise ValueError("%s is not a valid SymPy expression" % (f))
+        raise ValueError("%s is not a valid SymPy expression" % f)
+
+    if not isinstance(symbol, Expr) and  symbol is not None:
+        raise ValueError("%s is not a valid SymPy symbol" % symbol)
+
+    if not isinstance(domain, Set):
+        raise ValueError("%s is not a valid domain" %(domain))
 
     free_symbols = f.free_symbols
 
+    if symbol is None and not free_symbols:
+        b = Eq(f, 0)
+        if b is S.true:
+            return domain
+        elif b is S.false:
+            return S.EmptySet
+        else:
+            raise NotImplementedError(filldedent('''
+                relationship between value and 0 is unknown: %s''' % b))
+
     if symbol is None:
         if len(free_symbols) == 1:
             symbol = free_symbols.pop()
@@ -985,31 +1053,29 @@ def solveset(f, symbol=None, domain=S.Complexes):
         # the xreplace will be needed if a ConditionSet is returned
         return solveset(f[0], s[0], domain).xreplace(swap)
 
-    elif not free_symbols:
-        b = Eq(f, 0)
-        if b is S.true:
-            return domain
-        elif b is S.false:
-            return S.EmptySet
-        else:
-            raise NotImplementedError(filldedent('''
-                relationship between value and 0 is unknown: %s''' % b))
-
-    if isinstance(f, Eq):
-        from sympy.core import Add
-        f = Add(f.lhs, - f.rhs, evaluate=False)
-    elif f.is_Relational:
-        if not domain.is_subset(S.Reals):
-            raise NotImplementedError(filldedent('''
-                Inequalities in the complex domain are
-                not supported. Try the real domain by
-                setting domain=S.Reals'''))
-        try:
-            result = solve_univariate_inequality(
-            f, symbol, domain=domain, relational=False)
-        except NotImplementedError:
-            result = ConditionSet(symbol, f, domain)
-        return result
+    if domain.is_subset(S.Reals):
+        if not symbol.is_real:
+            assumptions = symbol.assumptions0
+            assumptions['real'] = True
+            try:
+                r = Dummy('r', **assumptions)
+                return solveset(f.xreplace({symbol: r}), r, domain
+                    ).xreplace({r: symbol})
+            except InconsistentAssumptions:
+                pass
+    # Abs has its own handling method which avoids the
+    # rewriting property that the first piece of abs(x)
+    # is for x >= 0 and the 2nd piece for x < 0 -- solutions
+    # can look better if the 2nd condition is x <= 0. Since
+    # the solution is a set, duplication of results is not
+    # an issue, e.g. {y, -y} when y is 0 will be {0}
+    f, mask = _masked(f, Abs)
+    f = f.rewrite(Piecewise) # everything that's not an Abs
+    for d, e in mask:
+        # everything *in* an Abs
+        e = e.func(e.args[0].rewrite(Piecewise))
+        f = f.xreplace({d: e})
+    f = piecewise_fold(f)
 
     return _solveset(f, symbol, domain, _check=True)
 


Test patch:
diff --git a/sympy/core/tests/test_arit.py b/sympy/core/tests/test_arit.py
index 6e905451e5..9cc84a725e 100644
--- a/sympy/core/tests/test_arit.py
+++ b/sympy/core/tests/test_arit.py
@@ -1636,9 +1636,6 @@ def test_Mod():
     # modular exponentiation
     assert Mod(Pow(4, 13, evaluate=False), 497) == Mod(Pow(4, 13), 497)
     assert Mod(Pow(2, 10000000000, evaluate=False), 3) == 1
-    assert Mod(Pow(32131231232, 9**10**6, evaluate=False),10**12) == pow(32131231232,9**10**6,10**12)
-    assert Mod(Pow(33284959323, 123**999, evaluate=False),11**13) == pow(33284959323,123**999,11**13)
-    assert Mod(Pow(78789849597, 333**555, evaluate=False),12**9) == pow(78789849597,333**555,12**9)
 
     # Wilson's theorem
     factorial(18042, evaluate=False) % 18043 == 18042
diff --git a/sympy/functions/elementary/tests/test_complexes.py b/sympy/functions/elementary/tests/test_complexes.py
index 069c2cc45a..866846d1d8 100644
--- a/sympy/functions/elementary/tests/test_complexes.py
+++ b/sympy/functions/elementary/tests/test_complexes.py
@@ -1,5 +1,5 @@
 from sympy import (
-    Abs, adjoint, arg, atan, atan2, conjugate, cos, DiracDelta, E, exp, expand,
+    Abs, adjoint, arg, atan2, conjugate, cos, DiracDelta, E, exp, expand,
     Expr, Function, Heaviside, I, im, log, nan, oo, pi, Rational, re, S,
     sign, sin, sqrt, Symbol, symbols, transpose, zoo, exp_polar, Piecewise,
     Interval, comp, Integral, Matrix, ImmutableMatrix, SparseMatrix,
@@ -399,12 +399,6 @@ def test_Abs():
     assert 1/Abs(x)**3 == 1/(x**2*Abs(x))
     assert Abs(x)**-3 == Abs(x)/(x**4)
     assert Abs(x**3) == x**2*Abs(x)
-    assert Abs(I**I) == exp(-pi/2)
-    assert Abs((4 + 5*I)**(6 + 7*I)) == 68921*exp(-7*atan(S(5)/4))
-    y = Symbol('y', real=True)
-    assert Abs(I**y) == 1
-    y = Symbol('y')
-    assert Abs(I**y) == exp(-pi*im(y)/2)
 
     x = Symbol('x', imaginary=True)
     assert Abs(x).diff(x) == -sign(x)
@@ -883,3 +877,9 @@ def test_issue_14216():
     A = MatrixSymbol("A", 2, 2)
     assert unpolarify(A[0, 0]) == A[0, 0]
     assert unpolarify(A[0, 0]*A[1, 0]) == A[0, 0]*A[1, 0]
+
+
+def test_issue_14238():
+    # doesn't cause recursion error
+    r = Symbol('r', real=True)
+    assert Abs(r + Piecewise((0, r > 0), (1 - r, True)))
diff --git a/sympy/physics/mechanics/tests/test_functions.py b/sympy/physics/mechanics/tests/test_functions.py
index dd94b908d2..db8d3ca149 100644
--- a/sympy/physics/mechanics/tests/test_functions.py
+++ b/sympy/physics/mechanics/tests/test_functions.py
@@ -6,8 +6,6 @@
                                      kinetic_energy, linear_momentum,
                                      outer, potential_energy, msubs,
                                      find_dynamicsymbols)
-
-from sympy.physics.mechanics.functions import gravity
 from sympy.physics.vector.vector import Vector
 from sympy.utilities.pytest import raises
 
@@ -180,21 +178,3 @@ def test_find_dynamicsymbols():
     assert find_dynamicsymbols(v, reference_frame=A) == sol
     # Test if a ValueError is raised on supplying only a vector as input
     raises(ValueError, lambda: find_dynamicsymbols(v))
-
-def test_gravity():
-    N = ReferenceFrame('N')
-    m, M, g = symbols('m M g')
-    F1, F2 = dynamicsymbols('F1 F2')
-    po = Point('po')
-    pa = Particle('pa', po, m)
-    A = ReferenceFrame('A')
-    P = Point('P')
-    I = outer(A.x, A.x)
-    B = RigidBody('B', P, A, M, (I, P))
-    forceList = [(po, F1), (P, F2)]
-    forceList.extend(gravity(g*N.y, pa, B))
-    l = [(po, F1), (P, F2), (po, g*m*N.y), (P, g*M*N.y)]
-
-    for i in range(len(l)):
-        for j in range(len(l[i])):
-            assert forceList[i][j] == l[i][j]
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 03eb18c3c5..00c82d6c93 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -9,7 +9,7 @@
 from sympy.functions.elementary.exponential import (LambertW, exp, log)
 from sympy.functions.elementary.hyperbolic import (HyperbolicFunction,
     atanh, sinh, tanh)
-from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.functions.elementary.miscellaneous import sqrt, Min, Max
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.functions.elementary.trigonometric import (
     TrigonometricFunction, acos, acot, acsc, asec, asin, atan, atan2,
@@ -56,6 +56,9 @@ def test_invert_real():
     def ireal(x, s=S.Reals):
         return Intersection(s, x)
 
+    # issue 14223
+    assert invert_real(x, 0, x, Interval(1, 2)) == (x, S.EmptySet)
+
     minus_n = Intersection(Interval(-oo, 0), FiniteSet(-n))
     plus_n = Intersection(Interval(0, oo), FiniteSet(n))
     assert solveset(abs(x) - n, x, S.Reals) == Union(minus_n, plus_n)
@@ -189,10 +192,12 @@ def test_domain_check():
     assert domain_check(x, x, oo) is False
     assert domain_check(0, x, oo) is False
 
+
 def test_issue_11536():
     assert solveset(0**x - 100, x, S.Reals) == S.EmptySet
     assert solveset(0**x - 1, x, S.Reals) == FiniteSet(0)
 
+
 def test_is_function_class_equation():
     from sympy.abc import x, a
     assert _is_function_class_equation(TrigonometricFunction,
@@ -263,6 +268,10 @@ def test_garbage_input():
     raises(ValueError, lambda: solveset_complex([x], x))
     assert solveset_complex(x, pi) == S.EmptySet
 
+    raises(ValueError, lambda: solveset((x, y), x))
+    raises(ValueError, lambda: solveset(x + 1, S.Reals))
+    raises(ValueError, lambda: solveset(x + 1, x, 2))
+
 
 def test_solve_mul():
     assert solveset_real((a*x + b)*(exp(x) - 3), x) == \
@@ -278,7 +287,9 @@ def test_solve_invert():
     assert solveset_real(3**(x + 2), x) == FiniteSet()
     assert solveset_real(3**(2 - x), x) == FiniteSet()
 
-    assert solveset_real(y - b*exp(a/x), x) == Intersection(S.Reals, FiniteSet(a/log(y/b)))
+    assert solveset_real(y - b*exp(a/x), x) == Intersection(
+        S.Reals, FiniteSet(a/log(y/b)))
+
     # issue 4504
     assert solveset_real(2**x - 10, x) == FiniteSet(log(10)/log(2))
 
@@ -652,7 +663,7 @@ def test_atan2():
     assert solveset_real(atan2(x, 2) - pi/3, x) == FiniteSet(2*sqrt(3))
 
 
-def test_piecewise():
+def test_piecewise_solveset():
     eq = Piecewise((x - 2, Gt(x, 2)), (2 - x, True)) - 3
     assert set(solveset_real(eq, x)) == set(FiniteSet(-1, 5))
 
@@ -665,8 +676,11 @@ def test_piecewise():
     f = Piecewise(((x - 2)**2, x >= 0), (0, True))
     assert solveset(f, x, domain=S.Reals) == Union(FiniteSet(2), Interval(-oo, 0, True, True))
 
-    assert solveset(Piecewise((x + 1, x > 0), (I, True)) - I, x) == \
-        Interval(-oo, 0)
+    assert solveset(
+        Piecewise((x + 1, x > 0), (I, True)) - I, x, S.Reals
+        ) == Interval(-oo, 0)
+
+    assert solveset(Piecewise((x - 1, Ne(x, I)), (x, True)), x) == FiniteSet(1)
 
 
 def test_solveset_complex_polynomial():
@@ -793,10 +807,11 @@ def test_solve_trig():
                   9*sqrt(57))**(S(1)/3))/(3*(67 + 9*sqrt(57))**(S(1)/6))) +
                   2*pi), S.Integers))
 
-    assert solveset_real(2*tan(x)*sin(x) + 1, x) == \
-        Union(ImageSet(Lambda(n, 2*n*pi + 2*atan(sqrt(4
-                + sqrt(17)))), S.Integers), ImageSet(Lambda(n, 2*n*pi
-                - 2*atan(sqrt(4 + sqrt(17))) + 2*pi), S.Integers))
+    assert solveset_real(2*tan(x)*sin(x) + 1, x) == Union(
+        ImageSet(Lambda(n, 2*n*pi + atan(sqrt(2)*sqrt(-1 + sqrt(17))/
+            (-sqrt(17) + 1)) + pi), S.Integers),
+        ImageSet(Lambda(n, 2*n*pi - atan(sqrt(2)*sqrt(-1 + sqrt(17))/
+            (-sqrt(17) + 1)) + pi), S.Integers))
 
     assert solveset_real(cos(2*x)*cos(4*x) - 1, x) == \
                             ImageSet(Lambda(n, n*pi), S.Integers)
@@ -958,17 +973,17 @@ def test_conditionset():
     assert solveset(Eq(sin(x)**2 + cos(x)**2, 1), x, domain=S.Reals) == \
         ConditionSet(x, True, S.Reals)
 
-    assert solveset(Eq(x**2 + x*sin(x), 1), x, domain=S.Reals) == \
-        ConditionSet(x, Eq(x*(x + sin(x)) - 1, 0), S.Reals)
+    assert solveset(Eq(x**2 + x*sin(x), 1), x, domain=S.Reals
+        ) == ConditionSet(x, Eq(x*(x + sin(x)) - 1, 0), S.Reals)
 
-    assert solveset(Eq(sin(Abs(x)), x), x, domain=S.Reals) == \
-        ConditionSet(x, Eq(-x + sin(Abs(x)), 0), Interval(-oo, oo))
+    assert solveset(Eq(-I*(exp(I*x) - exp(-I*x))/2, 1), x
+        ) == imageset(Lambda(n, 2*n*pi + pi/2), S.Integers)
 
-    assert solveset(Eq(-I*(exp(I*x) - exp(-I*x))/2, 1), x) == \
-        imageset(Lambda(n, 2*n*pi + pi/2), S.Integers)
+    assert solveset(x + sin(x) > 1, x, domain=S.Reals
+        ) == ConditionSet(x, x + sin(x) > 1, S.Reals)
 
-    assert solveset(x + sin(x) > 1, x, domain=S.Reals) == \
-        ConditionSet(x, x + sin(x) > 1, S.Reals)
+    assert solveset(Eq(sin(Abs(x)), x), x, domain=S.Reals
+        ) == ConditionSet(x, Eq(-x + sin(Abs(x)), 0), S.Reals)
 
 
 @XFAIL
@@ -988,9 +1003,8 @@ def test_solveset_domain():
 def test_improve_coverage():
     from sympy.solvers.solveset import _has_rational_power
     x = Symbol('x')
-    y = exp(x+1/x**2)
-    solution = solveset(y**2+y, x, S.Reals)
-    unsolved_object = ConditionSet(x, Eq((exp((x**3 + 1)/x**2) + 1)*exp((x**3 + 1)/x**2), 0), S.Reals)
+    solution = solveset(exp(x) + sin(x), x, S.Reals)
+    unsolved_object = ConditionSet(x, Eq(exp(x) + sin(x), 0), S.Reals)
     assert solution == unsolved_object
 
     assert _has_rational_power(sin(x)*exp(x) + 1, x) == (False, S.One)
@@ -1676,6 +1690,7 @@ def test__is_finite_with_finite_vars():
         Dummy(), Dummy(real=True), Dummy(complex=True)))
     assert f(1/Dummy(real=False)) is True  # b/c it's finite but not 0
 
+
 def test_issue_13550():
     assert solveset(x**2 - 2*x - 15, symbol = x, domain = Interval(-oo, 0)) == FiniteSet(-3)
 
@@ -1683,3 +1698,30 @@ def test_issue_13550():
 def test_issue_13849():
     t = symbols('t')
     assert nonlinsolve((t*(sqrt(5) + sqrt(2)) - sqrt(2), t), t) == EmptySet()
+
+
+def test_issue_14223():
+    x = Symbol('x')
+    assert solveset((Abs(x + Min(x, 2)) - 2).rewrite(Piecewise), x,
+        S.Reals) == FiniteSet(-1, 1)
+    assert solveset((Abs(x + Min(x, 2)) - 2).rewrite(Piecewise), x,
+        Interval(0, 2)) == FiniteSet(1)
+
+
+def test_issue_10158():
+    x = Symbol('x')
+    dom = S.Reals
+    assert solveset(x*Max(x, 15) - 10, x, dom) == FiniteSet(2/S(3))
+    assert solveset(x*Min(x, 15) - 10, x, dom) == FiniteSet(-sqrt(10), sqrt(10))
+    assert solveset(Max(Abs(x - 3) - 1, x + 2) - 3, x, dom) == FiniteSet(-1, 1)
+    assert solveset(Abs(x - 1) - Abs(y), x, dom) == FiniteSet(-Abs(y) + 1, Abs(y) + 1)
+    assert solveset(Abs(x + 4*Abs(x + 1)), x, dom) == FiniteSet(-4/S(3), -4/S(5))
+    assert solveset(2*Abs(x + Abs(x + Max(3, x))) - 2, x, S.Reals
+        ) == FiniteSet(-1, -2)
+    dom = S.Complexes
+    raises(ValueError, lambda: solveset(x*Max(x, 15) - 10, x, dom))
+    raises(ValueError, lambda: solveset(x*Min(x, 15) - 10, x, dom))
+    raises(ValueError, lambda:
+        solveset(Max(Abs(x - 3) - 1, x + 2) - 3, x, dom))
+    raises(ValueError, lambda: solveset(Abs(x - 1) - Abs(y), x, dom))
+    raises(ValueError, lambda: solveset(Abs(x + 4*Abs(x + 1)), x, dom))