Repository: sympy/sympy

Issue:
solve_univariate_inequality gives wrong solution to the rational inequality abs(1/(x-1)) > 1
In current SymPy (commit c2ffd9), I see
```python
In [1]: x = symbols('x', real=True)
In [2]: from sympy.solvers.solveset import solveset
In [3]: solveset(abs(1/(x-1)) > 1)
Out[3]: (1, 2)
```
while the result should be ```(0, 1) ∪ (1, 2)```.

The following similar case works:
```python
In [4]: solveset(abs(1/x) > 1)
Out[4]: (-1, 0) ∪ (0, 1)
```
@aktech - can you please take a look?

Discussion:
asmeurer: Does this fix the original issue in #9567 as well? 
asmeurer: I don't understand what solve_linear does, but if it fixes the issue then +1.
smichr: I showed on the issue page that the reported issue returns correctly. 

As for solve_linear, it is trying to solve for a linear variable and also report (early in the solve process) if there is no solution possible (like when a  variable appears in an expression but actually cancels after simplification).

I can't believe we didn't see this error earlier.
asmeurer: I guess what I don't understand is why it sometimes returns a solution and sometimes returns the numerator and denominator. 
smichr: > why it sometimes returns a solution and sometimes returns the numerator and denominator

it should probably have been a private function of solve; I think the original idea was to not throw away the extracted numerator which is needed for the more extended work of solve. Since the function was returning 2 values when it found a solution it made sense to return the denominator with the numerator, too, so 2 values would be returned when there wasn't a linear solution.

Patch:
diff --git a/sympy/core/mod.py b/sympy/core/mod.py
index f31e75d4f8..4d58ebde75 100644
--- a/sympy/core/mod.py
+++ b/sympy/core/mod.py
@@ -36,7 +36,7 @@ def doit(p, q):
             to be less than or equal q.
             """
 
-            if p.is_infinite or q.is_infinite or p is nan or q is nan:
+            if p.is_infinite or q.is_infinite:
                 return nan
             if (p == q or p == -q or
                     p.is_Pow and p.exp.is_Integer and p.base == q or
diff --git a/sympy/geometry/polygon.py b/sympy/geometry/polygon.py
index fbc93fbaf0..24a4e90960 100644
--- a/sympy/geometry/polygon.py
+++ b/sympy/geometry/polygon.py
@@ -516,7 +516,7 @@ def encloses_point(self, p):
         References
         ==========
 
-        [1] http://paulbourke.net/geometry/polygonmesh/#insidepoly
+        [1] http://www.ariel.com.au/a/python-point-int-poly.html
 
         """
         p = Point(p)
@@ -559,8 +559,8 @@ def encloses_point(self, p):
                     if 0 <= max(p1x, p2x):
                         if p1y != p2y:
                             xinters = (-p1y)*(p2x - p1x)/(p2y - p1y) + p1x
-                            if p1x == p2x or 0 <= xinters:
-                                hit_odd = not hit_odd
+                        if p1x == p2x or 0 <= xinters:
+                            hit_odd = not hit_odd
             p1x, p1y = p2x, p2y
         return hit_odd
 
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index f9283083ed..f2c3fd2446 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1271,17 +1271,18 @@ def _solve(f, *symbols, **flags):
         got_s = set([])
         result = []
         for s in symbols:
-            n, d = solve_linear(f, symbols=[s])
-            if n.is_Symbol:
+            xi, v = solve_linear(f, symbols=[s])
+            if xi == s:
                 # no need to check but we should simplify if desired
                 if flags.get('simplify', True):
-                    d = simplify(d)
-                if got_s and any([ss in d.free_symbols for ss in got_s]):
+                    v = simplify(v)
+                vfree = v.free_symbols
+                if got_s and any([ss in vfree for ss in got_s]):
                     # sol depends on previously solved symbols: discard it
                     continue
-                got_s.add(n)
-                result.append({n: d})
-            elif n and d:  # otherwise there was no solution for s
+                got_s.add(xi)
+                result.append({xi: v})
+            elif xi:  # there might be a non-linear solution if xi is not 0
                 failed.append(s)
         if not failed:
             return result
@@ -1329,7 +1330,7 @@ def _solve(f, *symbols, **flags):
     elif f.is_Piecewise:
         result = set()
         for n, (expr, cond) in enumerate(f.args):
-            candidates = _solve(expr, *symbols, **flags)
+            candidates = _solve(expr, symbol, **flags)
             for candidate in candidates:
                 if candidate in result:
                     continue
@@ -1363,9 +1364,9 @@ def _solve(f, *symbols, **flags):
         check = False
     else:
         # first see if it really depends on symbol and whether there
-        # is a linear solution
+        # is only a linear solution
         f_num, sol = solve_linear(f, symbols=symbols)
-        if not symbol in f_num.free_symbols:
+        if f_num is S.Zero:
             return []
         elif f_num.is_Symbol:
             # no need to check but simplify if desired
@@ -1821,36 +1822,52 @@ def _ok_syms(e, sort=False):
 
 
 def solve_linear(lhs, rhs=0, symbols=[], exclude=[]):
-    r""" Return a tuple derived from f = lhs - rhs that is either:
+    r""" Return a tuple derived from f = lhs - rhs that is one of
+    the following:
+
+        (0, 1) meaning that ``f`` is independent of the symbols in
+        ``symbols`` that aren't in ``exclude``, e.g::
+
+            >>> from sympy.solvers.solvers import solve_linear
+            >>> from sympy.abc import x, y, z
+            >>> from sympy import cos, sin
+            >>> eq = y*cos(x)**2 + y*sin(x)**2 - y  # = y*(1 - 1) = 0
+            >>> solve_linear(eq)
+            (0, 1)
+            >>> eq = cos(x)**2 + sin(x)**2  # = 1
+            >>> solve_linear(eq)
+            (0, 1)
+            >>> solve_linear(x, exclude=[x])
+            (0, 1)
+
+        (0, 0) meaning that there is no solution to the equation
+        amongst the symbols given.
+
+            (If the first element of the tuple is not zero then
+            the function is guaranteed to be dependent on a symbol
+            in ``symbols``.)
+
+        (symbol, solution) where symbol appears linearly in the
+        numerator of ``f``, is in ``symbols`` (if given) and is
+        not in ``exclude`` (if given). No simplification is done
+        to ``f`` other than a ``mul=True`` expansion, so the
+        solution will correspond strictly to a unique solution.
+
+        ``(n, d)`` where ``n`` and ``d`` are the numerator and
+        denominator of ``f`` when the numerator was not linear
+        in any symbol of interest; ``n`` will never be a symbol
+        unless a solution for that symbol was found (in which case
+        the second element is the solution, not the denominator).
 
-        (numerator, denominator) of ``f``
-            If this comes back as (0, 1) it means
-            that ``f`` is independent of the symbols in ``symbols``, e.g::
-
-                y*cos(x)**2 + y*sin(x)**2 - y = y*(0) = 0
-                cos(x)**2 + sin(x)**2 = 1
-
-            If it comes back as (0, 0) there is no solution to the equation
-            amongst the symbols given.
-
-            If the numerator is not zero then the function is guaranteed
-            to be dependent on a symbol in ``symbols``.
-
-        or
-
-        (symbol, solution) where symbol appears linearly in the numerator of
-        ``f``, is in ``symbols`` (if given) and is not in ``exclude`` (if given).
-
-        No simplification is done to ``f`` other than and mul=True expansion,
-        so the solution will correspond strictly to a unique solution.
 
     Examples
     ========
 
-    >>> from sympy.solvers.solvers import solve_linear
-    >>> from sympy.abc import x, y, z
+    >>> from sympy.core.power import Pow
+    >>> from sympy.polys.polytools import cancel
 
-    These are linear in x and 1/x:
+    The variable ``x`` appears as a linear variable in each of the
+    following:
 
     >>> solve_linear(x + y**2)
     (x, -y**2)
@@ -1862,34 +1879,53 @@ def solve_linear(lhs, rhs=0, symbols=[], exclude=[]):
     >>> solve_linear(x**2/y**2 - 3)
     (x**2 - 3*y**2, y**2)
 
-    If the numerator is a symbol then (0, 0) is returned if the solution for
-    that symbol would have set any denominator to 0:
+    If the numerator of the expression is a symbol then (0, 0) is
+    returned if the solution for that symbol would have set any
+    denominator to 0:
 
-    >>> solve_linear(1/(1/x - 2))
+    >>> eq = 1/(1/x - 2)
+    >>> eq.as_numer_denom()
+    (x, -2*x + 1)
+    >>> solve_linear(eq)
     (0, 0)
-    >>> 1/(1/x) # to SymPy, this looks like x ...
+
+    But automatic rewriting may cause a symbol in the denominator to
+    appear in the numerator so a solution will be returned:
+
+    >>> (1/x)**-1
     x
-    >>> solve_linear(1/(1/x)) # so a solution is given
+    >>> solve_linear((1/x)**-1)
     (x, 0)
 
-    If x is allowed to cancel, then this appears linear, but this sort of
-    cancellation is not done so the solution will always satisfy the original
-    expression without causing a division by zero error.
+    Use an unevaluated expression to avoid this:
 
-    >>> solve_linear(x**2*(1/x - z**2/x))
+    >>> solve_linear(Pow(1/x, -1, evaluate=False))
+    (0, 0)
+
+    If ``x`` is allowed to cancel in the following expression, then it
+    appears to be linear in ``x``, but this sort of cancellation is not
+    done by ``solve_linear`` so the solution will always satisfy the
+    original expression without causing a division by zero error.
+
+    >>> eq = x**2*(1/x - z**2/x)
+    >>> solve_linear(cancel(eq))
+    (x, 0)
+    >>> solve_linear(eq)
     (x**2*(-z**2 + 1), x)
 
-    You can give a list of what you prefer for x candidates:
+    A list of symbols for which a solution is desired may be given:
 
     >>> solve_linear(x + y + z, symbols=[y])
     (y, -x - z)
 
-    You can also indicate what variables you don't want to consider:
+    A list of symbols to ignore may also be given:
 
-    >>> solve_linear(x + y + z, exclude=[x, z])
+    >>> solve_linear(x + y + z, exclude=[x])
     (y, -x - z)
 
-    If only x was excluded then a solution for y or z might be obtained.
+    (A solution for ``y`` is obtained because it is the first variable
+    from the canonically sorted list of symbols that had a linear
+    solution.)
 
     """
     if isinstance(lhs, Equality):
@@ -1922,48 +1958,50 @@ def solve_linear(lhs, rhs=0, symbols=[], exclude=[]):
                              ''' % (bad, eg)))
         symbols = free.intersection(symbols)
     symbols = symbols.difference(exclude)
+    if not symbols:
+        return S.Zero, S.One
     dfree = d.free_symbols
 
-    # derivatives are easy to do but tricky to analyze to see if they are going
-    # to disallow a linear solution, so for simplicity we just evaluate the
-    # ones that have the symbols of interest
+    # derivatives are easy to do but tricky to analyze to see if they
+    # are going to disallow a linear solution, so for simplicity we
+    # just evaluate the ones that have the symbols of interest
     derivs = defaultdict(list)
     for der in n.atoms(Derivative):
         csym = der.free_symbols & symbols
         for c in csym:
             derivs[c].append(der)
 
-    if symbols:
-        all_zero = True
-        for xi in symbols:
-            # if there are derivatives in this var, calculate them now
-            if type(derivs[xi]) is list:
-                derivs[xi] = dict([(der, der.doit()) for der in derivs[xi]])
-            nn = n.subs(derivs[xi])
-            dn = nn.diff(xi)
-            if dn:
-                all_zero = False
-                if dn is S.NaN:
-                    break
-                if not xi in dn.free_symbols:
-                    vi = -(nn.subs(xi, 0))/dn
-                    if dens is None:
-                        dens = _simple_dens(eq, symbols)
-                    if not any(checksol(di, {xi: vi}, minimal=True) is True
-                              for di in dens):
-                        # simplify any trivial integral
-                        irep = [(i, i.doit()) for i in vi.atoms(Integral) if
-                                i.function.is_number]
-                        # do a slight bit of simplification
-                        vi = expand_mul(vi.subs(irep))
-                        if not d.has(xi) or not (d/xi).has(xi):
-                            return xi, vi
-
-        if all_zero:
-            return S.Zero, S.One
-    if n.is_Symbol:  # there was no valid solution
-        n = d = S.Zero
-    return n, d  # should we cancel now?
+    all_zero = True
+    for xi in sorted(symbols, key=default_sort_key):  # canonical order
+        # if there are derivatives in this var, calculate them now
+        if type(derivs[xi]) is list:
+            derivs[xi] = dict([(der, der.doit()) for der in derivs[xi]])
+        newn = n.subs(derivs[xi])
+        dnewn_dxi = newn.diff(xi)
+        # dnewn_dxi can be nonzero if it survives differentation by any
+        # of its free symbols
+        free = dnewn_dxi.free_symbols
+        if dnewn_dxi and (not free or any(dnewn_dxi.diff(s) for s in free)):
+            all_zero = False
+            if dnewn_dxi is S.NaN:
+                break
+            if xi not in dnewn_dxi.free_symbols:
+                vi = -(newn.subs(xi, 0))/dnewn_dxi
+                if dens is None:
+                    dens = _simple_dens(eq, symbols)
+                if not any(checksol(di, {xi: vi}, minimal=True) is True
+                          for di in dens):
+                    # simplify any trivial integral
+                    irep = [(i, i.doit()) for i in vi.atoms(Integral) if
+                            i.function.is_number]
+                    # do a slight bit of simplification
+                    vi = expand_mul(vi.subs(irep))
+                    return xi, vi
+    if all_zero:
+        return S.Zero, S.One
+    if n.is_Symbol: # no solution for this symbol was found
+        return S.Zero, S.Zero
+    return n, d
 
 
 def minsolve_linear_system(system, *symbols, **flags):
@@ -1971,9 +2009,8 @@ def minsolve_linear_system(system, *symbols, **flags):
     Find a particular solution to a linear system.
 
     In particular, try to find a solution with the minimal possible number
-    of non-zero variables. This is a very computationally hard prolem.
-    If ``quick=True``, a heuristic is used. Otherwise a naive algorithm with
-    exponential complexity is used.
+    of non-zero variables using a naive algorithm with exponential complexity.
+    If ``quick=True``, a heuristic is used.
     """
     quick = flags.get('quick', False)
     # Check if there are any non-zero solutions at all


Test patch:
diff --git a/sympy/core/tests/test_arit.py b/sympy/core/tests/test_arit.py
index 8c8625ce62..cbc6f77ebd 100644
--- a/sympy/core/tests/test_arit.py
+++ b/sympy/core/tests/test_arit.py
@@ -1502,9 +1502,6 @@ def test_Mod():
     assert x % 5 == Mod(x, 5)
     assert x % y == Mod(x, y)
     assert (x % y).subs({x: 5, y: 3}) == 2
-    assert Mod(nan, 1) == nan
-    assert Mod(1, nan) == nan
-    assert Mod(nan, nan) == nan
 
     # Float handling
     point3 = Float(3.3) % 1
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index cfd5ca0c74..9aee8f389e 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -501,6 +501,8 @@ def test_issue_3870():
 def test_solve_linear():
     w = Wild('w')
     assert solve_linear(x, x) == (0, 1)
+    assert solve_linear(x, exclude=[x]) == (0, 1)
+    assert solve_linear(x, symbols=[w]) == (0, 1)
     assert solve_linear(x, y - 2*x) in [(x, y/3), (y, 3*x)]
     assert solve_linear(x, y - 2*x, exclude=[x]) == (y, 3*x)
     assert solve_linear(3*x - y, 0) in [(x, y/3), (y, 3*x)]
@@ -513,10 +515,15 @@ def test_solve_linear():
     assert solve_linear(cos(x)**2 + sin(x)**2 + 2 + y, symbols=[x]) == (0, 1)
     assert solve_linear(Eq(x, 3)) == (x, 3)
     assert solve_linear(1/(1/x - 2)) == (0, 0)
-    assert solve_linear((x + 1)*exp(-x), symbols=[x]) == (x + 1, exp(x))
+    assert solve_linear((x + 1)*exp(-x), symbols=[x]) == (x, -1)
     assert solve_linear((x + 1)*exp(x), symbols=[x]) == ((x + 1)*exp(x), 1)
-    assert solve_linear(x*exp(-x**2), symbols=[x]) == (0, 0)
+    assert solve_linear(x*exp(-x**2), symbols=[x]) == (x, 0)
     assert solve_linear(0**x - 1) == (0**x - 1, 1)
+    assert solve_linear(1 + 1/(x - 1)) == (x, 0)
+    eq = y*cos(x)**2 + y*sin(x)**2 - y  # = y*(1 - 1) = 0
+    assert solve_linear(eq) == (0, 1)
+    eq = cos(x)**2 + sin(x)**2  # = 1
+    assert solve_linear(eq) == (0, 1)
     raises(ValueError, lambda: solve_linear(Eq(x, 3), 3))
 
 
@@ -1771,3 +1778,7 @@ def test_issue_2840_8155():
     assert solve(2*sin(x) - 2*sin(2*x)) == [
         0, -pi, pi, -2*I*log(-sqrt(3)/2 - I/2), -2*I*log(-sqrt(3)/2 + I/2),
         -2*I*log(sqrt(3)/2 - I/2), -2*I*log(sqrt(3)/2 + I/2)]
+
+
+def test_issue_9567():
+    assert solve(1 + 1/(x - 1)) == [0]