Repository: sympy/sympy

Issue:
`(S.Naturals).contains(oo)` hangs
```
>>> sympy.__version__
'0.7.7.dev'
>>> python3 --version
Python 3.4.0
```


Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v147). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* sets
  - Rationals added to fancysets ([#16896](https://github.com/sympy/sympy/pull/16896) by [@smichr](https://github.com/smichr))

  - ImageSet._contains solver improvements ([#16896](https://github.com/sympy/sympy/pull/16896) by [@smichr](https://github.com/smichr))

  - ImageSet._contains considers all regions in a ComplexRegion ([#16896](https://github.com/sympy/sympy/pull/16896) by [@smichr](https://github.com/smichr))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.5.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234". See
    https://github.com/blog/1506-closing-issues-via-pull-requests . Please also
    write a comment on that issue linking back to this pull request once it is
    open. -->


    #### Brief description of what is fixed or changed


    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. See
    https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
    on how to write release notes. The bot will check your release notes
    automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    - sets
        - Rationals added to fancysets
        - ImageSet._contains solver improvements
        - ImageSet._contains considers all regions in a ComplexRegion
    <!-- END RELEASE NOTES -->


</details><p>

oscarbenjamin: Rationals here doesn't print nicely:
```
In [1]: S.Rationals                                                                                                               
Out[1]: Rationals

In [2]: S.Reals                                                                                                                   
Out[2]: ℝ

In [3]: S.Integers                                                                                                                
Out[3]: ℤ
```
What is `Rationals.__iter__` useful for? It currently only gives rational numbers between +-1. There is a way to enumerate all rationals in the diagram here though:
https://en.wikipedia.org/wiki/Rational_number#Properties


smichr: This result is driving me crazy...I will have to step through and see why 
`S.Naturals.intersect(S.Reals) -> Range(1, oo, 1)` but `S.Naturals0` gives itself.

smichr: ```console
def e2d(s):
    if type(s) is str:
        s = s.replace('==', '=')
        return [tuple(map(lambda i:i.strip(), li.split('='))) for li in s.strip().splitlines()]
    elif isinstance(s, Eq):
        return s.args
    elif type(s) is list:
        return list(map(e2d, s))

e2d('v_f=x\nsigma_f=f\nsigma_c=c\nsigma_p=p')

s2d = lambda s: [tuple(map(lambda i:i.strip(), li.split('='))) for li in s.strip().splitlines()]
```
codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/16896?src=pr&el=h1) Report
> Merging [#16896](https://codecov.io/gh/sympy/sympy/pull/16896?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/5158ecf2cb5b926f0f6421aab7b0279280305554?src=pr&el=desc) will **increase** coverage by `0.032%`.
> The diff coverage is `76.282%`.

```diff
@@              Coverage Diff              @@
##            master    #16896       +/-   ##
=============================================
+ Coverage   73.906%   73.939%   +0.032%     
=============================================
  Files          619       619               
  Lines       160001    160067       +66     
  Branches     37554     37575       +21     
=============================================
+ Hits        118251    118352      +101     
+ Misses       36269     36234       -35     
  Partials      5481      5481
```

smichr: @oscarbenjamin , do you want to give this a final check?
smichr: I am checking coverage.
codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/16896?src=pr&el=h1) Report
> Merging [#16896](https://codecov.io/gh/sympy/sympy/pull/16896?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/bfe26eae1cee10a5d221be028e43fe6bc1a49341?src=pr&el=desc) will **increase** coverage by `0.01%`.
> The diff coverage is `92.473%`.

```diff
@@             Coverage Diff              @@
##            master    #16896      +/-   ##
============================================
+ Coverage   73.921%   73.932%   +0.01%     
============================================
  Files          620       620              
  Lines       160369    160329      -40     
  Branches     37632     37622      -10     
============================================
- Hits        118547    118535      -12     
+ Misses       36335     36296      -39     
- Partials      5487      5498      +11
```


Patch:
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index 4c8a2559d0..90fce09242 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -280,8 +280,8 @@ def not_empty_in(finset_intersection, *syms):
     if len(syms) == 0:
         raise ValueError("One or more symbols must be given in syms.")
 
-    if finset_intersection.is_EmptySet:
-        return EmptySet()
+    if finset_intersection is S.EmptySet:
+        return S.EmptySet
 
     if isinstance(finset_intersection, Union):
         elm_in_sets = finset_intersection.args[0]
diff --git a/sympy/core/function.py b/sympy/core/function.py
index 2ddee5df54..7399654f14 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -1864,6 +1864,9 @@ def __new__(cls, variables, expr):
         for i in v:
             if not getattr(i, 'is_symbol', False):
                 raise TypeError('variable is not a symbol: %s' % i)
+        if len(v) != len(set(v)):
+            x = [i for i in v if v.count(i) > 1][0]
+            raise SyntaxError("duplicate argument '%s' in Lambda args" % x)
         if len(v) == 1 and v[0] == expr:
             return S.IdentityFunction
 
diff --git a/sympy/ntheory/continued_fraction.py b/sympy/ntheory/continued_fraction.py
index 1c38d3bf24..848f30e9c0 100644
--- a/sympy/ntheory/continued_fraction.py
+++ b/sympy/ntheory/continued_fraction.py
@@ -171,8 +171,8 @@ def continued_fraction_periodic(p, q, d=0, s=1):
     if (d - p**2)%q:
         d *= q**2
         sd *= q
-        p *= abs(q)
-        q *= abs(q)
+        p *= q
+        q *= q
 
     terms = []
     pq = {}
diff --git a/sympy/printing/pretty/pretty.py b/sympy/printing/pretty/pretty.py
index e070adf3ac..d67ff0b199 100644
--- a/sympy/printing/pretty/pretty.py
+++ b/sympy/printing/pretty/pretty.py
@@ -163,6 +163,7 @@ def _print_Atom(self, e):
     _print_Naturals = _print_Atom
     _print_Naturals0 = _print_Atom
     _print_Integers = _print_Atom
+    _print_Rationals = _print_Atom
     _print_Complexes = _print_Atom
 
     def _print_Reals(self, e):
diff --git a/sympy/printing/pretty/pretty_symbology.py b/sympy/printing/pretty/pretty_symbology.py
index d27eb51653..348a59bcf7 100644
--- a/sympy/printing/pretty/pretty_symbology.py
+++ b/sympy/printing/pretty/pretty_symbology.py
@@ -501,6 +501,7 @@ def xsym(sym):
                                 (U('DOUBLE-STRUCK CAPITAL N') +
                                  U('SUBSCRIPT ZERO'))),
     'Integers':                U('DOUBLE-STRUCK CAPITAL Z'),
+    'Rationals':               U('DOUBLE-STRUCK CAPITAL Q'),
     'Reals':                   U('DOUBLE-STRUCK CAPITAL R'),
     'Complexes':               U('DOUBLE-STRUCK CAPITAL C'),
     'Union':                   U('UNION'),
diff --git a/sympy/printing/str.py b/sympy/printing/str.py
index ee560ca0a8..0ce283b18b 100644
--- a/sympy/printing/str.py
+++ b/sympy/printing/str.py
@@ -582,6 +582,9 @@ def _print_Naturals(self, expr):
     def _print_Naturals0(self, expr):
         return 'Naturals0'
 
+    def _print_Rationals(self, expr):
+        return 'Rationals'
+
     def _print_Reals(self, expr):
         return 'Reals'
 
diff --git a/sympy/sets/__init__.py b/sympy/sets/__init__.py
index 453fc449e9..a66ce07138 100644
--- a/sympy/sets/__init__.py
+++ b/sympy/sets/__init__.py
@@ -10,4 +10,5 @@
 Naturals0 = S.Naturals0
 UniversalSet = S.UniversalSet
 Integers = S.Integers
+Rationals = S.Rationals
 del S
diff --git a/sympy/sets/conditionset.py b/sympy/sets/conditionset.py
index 3a4f3a8282..6437f471a8 100644
--- a/sympy/sets/conditionset.py
+++ b/sympy/sets/conditionset.py
@@ -175,7 +175,15 @@ def free_symbols(self):
         s, c, b = self.args
         return (c.free_symbols - s.free_symbols) | b.free_symbols
 
-    def contains(self, other):
+    def _contains(self, other):
+        d = Dummy()
+        try:
+            return self.as_relational(d).subs(d, other)
+        except TypeError:
+            # couldn't do the substitution without error
+            return False
+
+    def as_relational(self, other):
         return And(Lambda(self.sym, self.condition)(
             other), self.base_set.contains(other))
 
diff --git a/sympy/sets/fancysets.py b/sympy/sets/fancysets.py
index 3da73cf460..0d4804359b 100644
--- a/sympy/sets/fancysets.py
+++ b/sympy/sets/fancysets.py
@@ -4,14 +4,67 @@
 from sympy.core.compatibility import as_int, with_metaclass, range, PY3
 from sympy.core.expr import Expr
 from sympy.core.function import Lambda
+from sympy.core.numbers import oo
+from sympy.core.relational import Eq
 from sympy.core.singleton import Singleton, S
 from sympy.core.symbol import Dummy, symbols
 from sympy.core.sympify import _sympify, sympify, converter
 from sympy.logic.boolalg import And
-from sympy.sets.sets import Set, Interval, Union, FiniteSet, ProductSet
+from sympy.sets.sets import (Set, Interval, Union, FiniteSet,
+    ProductSet, Intersection)
+from sympy.sets.contains import Contains
+from sympy.sets.conditionset import ConditionSet
+from sympy.utilities.iterables import flatten
 from sympy.utilities.misc import filldedent
 
 
+class Rationals(with_metaclass(Singleton, Set)):
+    """
+    Represents the rational numbers. This set is also available as
+    the Singleton, S.Rationals.
+
+    Examples
+    ========
+
+    >>> from sympy import S
+    >>> S.Half in S.Rationals
+    True
+    >>> iterable = iter(S.Rationals)
+    >>> [next(iterable) for i in range(12)]
+    [0, 1, -1, 1/2, 2, -1/2, -2, 1/3, 3, -1/3, -3, 2/3]
+    """
+
+    is_iterable = True
+    _inf = S.NegativeInfinity
+    _sup = S.Infinity
+
+    def _contains(self, other):
+        if not isinstance(other, Expr):
+            return False
+        if other.is_Number:
+            return other.is_Rational
+        return other.is_rational
+
+    def __iter__(self):
+        from sympy.core.numbers import igcd, Rational
+        yield S.Zero
+        yield S.One
+        yield S.NegativeOne
+        d = 2
+        while True:
+            for n in range(d):
+                if igcd(n, d) == 1:
+                    yield Rational(n, d)
+                    yield Rational(d, n)
+                    yield Rational(-n, d)
+                    yield Rational(-d, n)
+            d += 1
+
+    @property
+    def _boundary(self):
+        return self
+
+
 class Naturals(with_metaclass(Singleton, Set)):
     """
     Represents the natural numbers (or counting numbers) which are all
@@ -47,11 +100,11 @@ class Naturals(with_metaclass(Singleton, Set)):
 
     def _contains(self, other):
         if not isinstance(other, Expr):
-            return S.false
+            return False
         elif other.is_positive and other.is_integer:
-            return S.true
+            return True
         elif other.is_integer is False or other.is_positive is False:
-            return S.false
+            return False
 
     def __iter__(self):
         i = self._inf
@@ -63,6 +116,10 @@ def __iter__(self):
     def _boundary(self):
         return self
 
+    def as_relational(self, x):
+        from sympy.functions.elementary.integers import floor
+        return And(Eq(floor(x), x), x >= self.inf, x < oo)
+
 
 class Naturals0(Naturals):
     """Represents the whole numbers which are all the non-negative integers,
@@ -121,10 +178,7 @@ class Integers(with_metaclass(Singleton, Set)):
     def _contains(self, other):
         if not isinstance(other, Expr):
             return S.false
-        elif other.is_integer:
-            return S.true
-        elif other.is_integer is False:
-            return S.false
+        return other.is_integer
 
     def __iter__(self):
         yield S.Zero
@@ -146,6 +200,10 @@ def _sup(self):
     def _boundary(self):
         return self
 
+    def as_relational(self, x):
+        from sympy.functions.elementary.integers import floor
+        return And(Eq(floor(x), x), -oo < x, x < oo)
+
 
 class Reals(with_metaclass(Singleton, Interval)):
     """
@@ -246,9 +304,13 @@ class ImageSet(Set):
     def __new__(cls, flambda, *sets):
         if not isinstance(flambda, Lambda):
             raise ValueError('first argument must be a Lambda')
-        if flambda is S.IdentityFunction and len(sets) == 1:
+
+        if flambda is S.IdentityFunction:
+            if len(sets) != 1:
+                raise ValueError('identify function requires a single set')
             return sets[0]
-        if not flambda.expr.free_symbols or not flambda.expr.args:
+
+        if not set(flambda.variables) & flambda.expr.free_symbols:
             return FiniteSet(flambda.expr)
 
         return Basic.__new__(cls, flambda, *sets)
@@ -275,17 +337,11 @@ def _contains(self, other):
         from sympy.solvers.solvers import solve
         from sympy.utilities.iterables import is_sequence, iterable, cartes
         L = self.lamda
-        if is_sequence(other):
-            if not is_sequence(L.expr):
-                return S.false
-            if len(L.expr) != len(other):
-                raise ValueError(filldedent('''
-    Dimensions of other and output of Lambda are different.'''))
-        elif iterable(other):
-                raise ValueError(filldedent('''
-    `other` should be an ordered object like a Tuple.'''))
+        if is_sequence(other) != is_sequence(L.expr):
+            return False
+        elif is_sequence(other) and len(L.expr) != len(other):
+            return False
 
-        solns = None
         if self._is_multivariate():
             if not is_sequence(L.expr):
                 # exprs -> (numer, denom) and check again
@@ -333,6 +389,12 @@ def _contains(self, other):
                         # for x, y and getting (x, 0), (0, y), (0, 0)
                         solns = [i for i in solns if not any(
                             s in i for s in variables)]
+                        if not solns:
+                            return False
+                    else:
+                        # not sure if [] means no solution or
+                        # couldn't find one
+                        return
         else:
             x = L.variables[0]
             if isinstance(L.expr, Expr):
@@ -344,31 +406,29 @@ def _contains(self, other):
                     msgset = solnsSet
             else:
                 # scalar -> vector
+                # note: it is not necessary for components of other
+                # to be in the corresponding base set unless the
+                # computed component is always in the corresponding
+                # domain. e.g. 1/2 is in imageset(x, x/2, Integers)
+                # while it cannot be in imageset(x, x + 2, Integers).
+                # So when the base set is comprised of integers or reals
+                # perhaps a pre-check could be done to see if the computed
+                # values are still in the set.
+                dom = self.base_set
                 for e, o in zip(L.expr, other):
-                    solns = solveset(e - o, x)
-                    if solns is S.EmptySet:
-                        return S.false
-                    for soln in solns:
-                        try:
-                            if soln in self.base_set:
-                                break  # check next pair
-                        except TypeError:
-                            if self.base_set.contains(soln.evalf()):
-                                break
-                    else:
-                        return S.false  # never broke so there was no True
-                return S.true
-
-        if solns is None:
-            raise NotImplementedError(filldedent('''
-            Determining whether %s contains %s has not
-            been implemented.''' % (msgset, other)))
+                    msgset = dom
+                    other = e - o
+                    dom = dom.intersection(solveset(e - o, x, domain=dom))
+                    if not dom:
+                        # there is no solution in common
+                        return False
+                return not isinstance(dom, Intersection)
         for soln in solns:
             try:
                 if soln in self.base_set:
-                    return S.true
+                    return True
             except TypeError:
-                return self.base_set.contains(soln.evalf())
+                return
         return S.false
 
     @property
@@ -913,7 +973,6 @@ def __new__(cls, sets, polar=False):
                 new_sets.append(sets)
             # Normalize input theta
             for k, v in enumerate(new_sets):
-                from sympy.sets import ProductSet
                 new_sets[k] = ProductSet(v.args[0],
                                          normalize_theta_set(v.args[1]))
             sets = Union(*new_sets)
diff --git a/sympy/sets/handlers/functions.py b/sympy/sets/handlers/functions.py
index d811e87b6d..c4b6423c62 100644
--- a/sympy/sets/handlers/functions.py
+++ b/sympy/sets/handlers/functions.py
@@ -1,4 +1,4 @@
-from sympy import Set, symbols, exp, log, S, Wild
+from sympy import Set, symbols, exp, log, S, Wild, Dummy, oo
 from sympy.core import Expr, Add
 from sympy.core.function import Lambda, _coeff_isneg, FunctionClass
 from sympy.core.mod import Mod
@@ -6,7 +6,7 @@
 from sympy.multipledispatch import dispatch
 from sympy.sets import (imageset, Interval, FiniteSet, Union, ImageSet,
                         EmptySet, Intersection, Range)
-from sympy.sets.fancysets import Integers, Naturals
+from sympy.sets.fancysets import Integers, Naturals, Reals
 
 
 _x, _y = symbols("x y")
@@ -37,6 +37,9 @@ def _set_function(f, x):
     if len(expr.free_symbols) > 1 or len(f.variables) != 1:
         return
     var = f.variables[0]
+    if not var.is_real:
+        if expr.subs(var, Dummy(real=True)).is_real is False:
+            return
 
     if expr.is_Piecewise:
         result = S.EmptySet
@@ -56,7 +59,7 @@ def _set_function(f, x):
 
             # remove the part which has been `imaged`
             domain_set = Complement(domain_set, intrvl)
-            if domain_set.is_EmptySet:
+            if domain_set is S.EmptySet:
                 break
         return result
 
@@ -169,6 +172,8 @@ def _set_function(f, self):
         return
 
     n = f.variables[0]
+    if expr == abs(n):
+        return S.Naturals0
 
     # f(x) + c and f(-x) + c cover the same integers
     # so choose the form that has the fewest negatives
@@ -212,11 +217,28 @@ def _set_function(f, self):
 
     x = f.variables[0]
     if not expr.free_symbols - {x}:
+        if expr == abs(x):
+            if self is S.Naturals:
+                return self
+            return S.Naturals0
         step = expr.coeff(x)
         c = expr.subs(x, 0)
         if c.is_Integer and step.is_Integer and expr == step*x + c:
             if self is S.Naturals:
                 c += step
             if step > 0:
-                return Range(c, S.Infinity, step)
-            return Range(c, S.NegativeInfinity, step)
+                if step == 1:
+                    if c == 0:
+                        return S.Naturals0
+                    elif c == 1:
+                        return S.Naturals
+                return Range(c, oo, step)
+            return Range(c, -oo, step)
+
+
+@dispatch(FunctionUnion, Reals)
+def _set_function(f, self):
+    expr = f.expr
+    if not isinstance(expr, Expr):
+        return
+    return _set_function(f, Interval(-oo, oo))
diff --git a/sympy/sets/handlers/intersection.py b/sympy/sets/handlers/intersection.py
index 34b9a94cf6..ea4108fae3 100644
--- a/sympy/sets/handlers/intersection.py
+++ b/sympy/sets/handlers/intersection.py
@@ -3,7 +3,7 @@
 from sympy.multipledispatch import dispatch
 from sympy.sets.conditionset import ConditionSet
 from sympy.sets.fancysets import (Integers, Naturals, Reals, Range,
-    ImageSet, Naturals0)
+    ImageSet, Naturals0, Rationals)
 from sympy.sets.sets import UniversalSet, imageset, ProductSet
 
 
@@ -25,27 +25,12 @@ def intersection_sets(a, b):
 
 @dispatch(Naturals, Naturals)
 def intersection_sets(a, b):
-    return a if a is S.Naturals0 else b
-
-@dispatch(Naturals, Interval)
-def intersection_sets(a, b):
-    return Intersection(S.Integers, b, Interval(a._inf, S.Infinity))
+    return a if a is S.Naturals else b
 
 @dispatch(Interval, Naturals)
 def intersection_sets(a, b):
     return intersection_sets(b, a)
 
-@dispatch(Integers, Interval)
-def intersection_sets(a, b):
-    try:
-        from sympy.functions.elementary.integers import floor, ceiling
-        if b._inf is S.NegativeInfinity and b._sup is S.Infinity:
-            return a
-        s = Range(ceiling(b.left), floor(b.right) + 1)
-        return intersection_sets(s, b)  # take out endpoints if open interval
-    except ValueError:
-        return None
-
 @dispatch(ComplexRegion, Set)
 def intersection_sets(self, other):
     if other.is_ComplexRegion:
@@ -157,7 +142,7 @@ def intersection_sets(a, b):
     # we want to know when the two equations might
     # have integer solutions so we use the diophantine
     # solver
-    va, vb = diop_linear(eq(r1, Dummy()) - eq(r2, Dummy()))
+    va, vb = diop_linear(eq(r1, Dummy('a')) - eq(r2, Dummy('b')))
 
     # check for no solution
     no_solution = va is None and vb is None
@@ -446,3 +431,33 @@ def intersection_sets(a, b):
 @dispatch(Set, Set)
 def intersection_sets(a, b):
     return None
+
+@dispatch(Integers, Rationals)
+def intersection_sets(a, b):
+    return a
+
+@dispatch(Naturals, Rationals)
+def intersection_sets(a, b):
+    return a
+
+@dispatch(Rationals, Reals)
+def intersection_sets(a, b):
+    return a
+
+def _intlike_interval(a, b):
+    try:
+        from sympy.functions.elementary.integers import floor, ceiling
+        if b._inf is S.NegativeInfinity and b._sup is S.Infinity:
+            return a
+        s = Range(max(a.inf, ceiling(b.left)), floor(b.right) + 1)
+        return intersection_sets(s, b)  # take out endpoints if open interval
+    except ValueError:
+        return None
+
+@dispatch(Integers, Interval)
+def intersection_sets(a, b):
+    return _intlike_interval(a, b)
+
+@dispatch(Naturals, Interval)
+def intersection_sets(a, b):
+    return _intlike_interval(a, b)
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index b85427e793..03cf57a30b 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -1,17 +1,18 @@
 from __future__ import print_function, division
 
 from itertools import product
+from collections import defaultdict
 import inspect
 
 from sympy.core.basic import Basic
 from sympy.core.compatibility import (iterable, with_metaclass,
-    ordered, range, PY3)
+    ordered, range, PY3, is_sequence)
 from sympy.core.cache import cacheit
 from sympy.core.evalf import EvalfMixin
 from sympy.core.evaluate import global_evaluate
 from sympy.core.expr import Expr
 from sympy.core.function import FunctionClass
-from sympy.core.logic import fuzzy_bool
+from sympy.core.logic import fuzzy_bool, fuzzy_or
 from sympy.core.mul import Mul
 from sympy.core.numbers import Float
 from sympy.core.operations import LatticeOp
@@ -27,6 +28,13 @@
 
 from mpmath import mpi, mpf
 
+
+tfn = defaultdict(lambda: None, {
+    True: S.true,
+    S.true: S.true,
+    False: S.false,
+    S.false: S.false})
+
 class Set(Basic):
     """
     The base class for any kind of set.
@@ -270,28 +278,57 @@ def _sup(self):
 
     def contains(self, other):
         """
-        Returns True if 'other' is contained in 'self' as an element.
-
-        As a shortcut it is possible to use the 'in' operator:
+        Returns a SymPy value indicating whether ``other`` is contained
+        in ``self``: ``true`` if it is, ``false`` if it isn't, else
+        an unevaluated ``Contains`` expression (or, as in the case of
+        ConditionSet and a union of FiniteSet/Intervals, an expression
+        indicating the conditions for containment).
 
         Examples
         ========
 
-        >>> from sympy import Interval
+        >>> from sympy import Interval, S
+        >>> from sympy.abc import x
+
         >>> Interval(0, 1).contains(0.5)
         True
-        >>> 0.5 in Interval(0, 1)
-        True
 
+        As a shortcut it is possible to use the 'in' operator, but that
+        will raise an error unless an affirmative true or false is not
+        obtained.
+
+        >>> Interval(0, 1).contains(x)
+        (0 <= x) & (x <= 1)
+        >>> x in Interval(0, 1)
+        Traceback (most recent call last):
+        ...
+        TypeError: did not evaluate to a bool: None
+
+        The result of 'in' is a bool, not a SymPy value
+
+        >>> 1 in Interval(0, 2)
+        True
+        >>> _ is S.true
+        False
         """
         other = sympify(other, strict=True)
-        ret = sympify(self._contains(other))
-        if ret is None:
-            ret = Contains(other, self, evaluate=False)
-        return ret
+        c = self._contains(other)
+        if c is None:
+            return Contains(other, self, evaluate=False)
+        b = tfn[c]
+        if b is None:
+            return c
+        return b
 
     def _contains(self, other):
-        raise NotImplementedError("(%s)._contains(%s)" % (self, other))
+        raise NotImplementedError(filldedent('''
+            (%s)._contains(%s) is not defined. This method, when
+            defined, will receive a sympified object. The method
+            should return True, False, None or something that
+            expresses what must be true for the containment of that
+            object in self to be evaluated. If None is returned
+            then a generic Contains object will be returned
+            by the ``contains`` method.''' % (self, other)))
 
     def is_subset(self, other):
         """
@@ -308,10 +345,12 @@ def is_subset(self, other):
 
         """
         if isinstance(other, Set):
-            # XXX issue 16873
-            # self might be an unevaluated form of self
-            # so the equality test will fail
-            return self.intersect(other) == self
+            s_o = self.intersect(other)
+            if s_o == self:
+                return True
+            elif not isinstance(other, Intersection):
+                return False
+            return s_o
         else:
             raise ValueError("Unknown argument '%s'" % other)
 
@@ -556,10 +595,12 @@ def __sub__(self, other):
         return Complement(self, other)
 
     def __contains__(self, other):
-        symb = sympify(self.contains(other))
-        if not (symb is S.true or symb is S.false):
-            raise TypeError('contains did not evaluate to a bool: %r' % symb)
-        return bool(symb)
+        other = sympify(other)
+        c = self._contains(other)
+        b = tfn[c]
+        if b is None:
+            raise TypeError('did not evaluate to a bool: %r' % c)
+        return b
 
 
 class ProductSet(Set):
@@ -648,13 +689,21 @@ def _contains(self, element):
 
         Passes operation on to constituent sets
         """
-        try:
+        if is_sequence(element):
             if len(element) != len(self.args):
-                return false
-        except TypeError:  # maybe element isn't an iterable
-            return false
-        return And(*
-            [set.contains(item) for set, item in zip(self.sets, element)])
+                return False
+        elif len(self.args) > 1:
+            return False
+        d = [Dummy() for i in element]
+        reps = dict(zip(d, element))
+        return tfn[self.as_relational(*d).xreplace(reps)]
+
+    def as_relational(self, *symbols):
+        if len(symbols) != len(self.args) or not all(
+                i.is_Symbol for i in symbols):
+            raise ValueError(
+                'number of symbols must match the number of sets')
+        return And(*[s.contains(i) for s, i in zip(self.args, symbols)])
 
     @property
     def sets(self):
@@ -915,17 +964,21 @@ def _contains(self, other):
             if not other.is_extended_real is None:
                 return other.is_extended_real
 
-        if self.left_open:
-            expr = other > self.start
-        else:
-            expr = other >= self.start
+        d = Dummy()
+        return self.as_relational(d).subs(d, other)
 
+    def as_relational(self, x):
+        """Rewrite an interval in terms of inequalities and logic operators."""
+        x = sympify(x)
         if self.right_open:
-            expr = And(expr, other < self.end)
+            right = x < self.end
         else:
-            expr = And(expr, other <= self.end)
-
-        return _sympify(expr)
+            right = x <= self.end
+        if self.left_open:
+            left = self.start < x
+        else:
+            left = self.start <= x
+        return And(left, right)
 
     @property
     def _measure(self):
@@ -958,19 +1011,6 @@ def is_right_unbounded(self):
         """Return ``True`` if the right endpoint is positive infinity. """
         return self.right is S.Infinity or self.right == Float("+inf")
 
-    def as_relational(self, x):
-        """Rewrite an interval in terms of inequalities and logic operators."""
-        x = sympify(x)
-        if self.right_open:
-            right = x < self.end
-        else:
-            right = x <= self.end
-        if self.left_open:
-            left = self.start < x
-        else:
-            left = self.start <= x
-        return And(left, right)
-
     def _eval_Eq(self, other):
         if not isinstance(other, Interval):
             if isinstance(other, FiniteSet):
@@ -1062,9 +1102,6 @@ def _sup(self):
         from sympy.functions.elementary.miscellaneous import Max
         return Max(*[set.sup for set in self.args])
 
-    def _contains(self, other):
-        return Or(*[set.contains(other) for set in self.args])
-
     @property
     def _measure(self):
         # Measure of a union is the sum of the measures of the sets minus
@@ -1120,14 +1157,28 @@ def boundary_of_set(i):
             return b
         return Union(*map(boundary_of_set, range(len(self.args))))
 
+    def _contains(self, other):
+        try:
+            d = Dummy()
+            r = self.as_relational(d).subs(d, other)
+            b = tfn[r]
+            if b is None and not any(isinstance(i.contains(other), Contains)
+                    for i in self.args):
+                return r
+            return b
+        except (TypeError, NotImplementedError):
+            return Or(*[s.contains(other) for s in self.args])
+
     def as_relational(self, symbol):
         """Rewrite a Union in terms of equalities and logic operators. """
-        if len(self.args) == 2:
-            a, b = self.args
-            if (a.sup == b.inf and a.inf is S.NegativeInfinity
-                    and b.sup is S.Infinity):
-                return And(Ne(symbol, a.sup), symbol < b.sup, symbol > a.inf)
-        return Or(*[set.as_relational(symbol) for set in self.args])
+        if all(isinstance(i, (FiniteSet, Interval)) for i in self.args):
+            if len(self.args) == 2:
+                a, b = self.args
+                if (a.sup == b.inf and a.inf is S.NegativeInfinity
+                        and b.sup is S.Infinity):
+                    return And(Ne(symbol, a.sup), symbol < b.sup, symbol > a.inf)
+            return Or(*[set.as_relational(symbol) for set in self.args])
+        raise NotImplementedError('relational of Union with non-Intervals')
 
     @property
     def is_iterable(self):
@@ -1613,16 +1664,9 @@ def _contains(self, other):
         False
 
         """
-        r = false
-        for e in self.args:
-            # override global evaluation so we can use Eq to do
-            # do the evaluation
-            t = Eq(e, other, evaluate=True)
-            if t is true:
-                return t
-            elif t is not false:
-                r = None
-        return r
+        # evaluate=True is needed to override evaluate=False context;
+        # we need Eq to do the evaluation
+        return fuzzy_or([tfn[Eq(e, other, evaluate=True)] for e in self.args])
 
     @property
     def _boundary(self):
@@ -1840,9 +1884,11 @@ def imageset(*args):
 
         if isinstance(set, ImageSet):
             if len(set.lamda.variables) == 1 and len(f.variables) == 1:
-                return imageset(Lambda(set.lamda.variables[0],
-                                       f.expr.subs(f.variables[0], set.lamda.expr)),
-                                set.base_set)
+                x = set.lamda.variables[0]
+                y = f.variables[0]
+                return imageset(
+                    Lambda(x, f.expr.subs(y, set.lamda.expr)),
+                    set.base_set)
 
         if r is not None:
             return r
@@ -1942,7 +1988,7 @@ def simplify_intersection(args):
             raise TypeError("Input args to Union must be Sets")
 
     # If any EmptySets return EmptySet
-    if any(s.is_EmptySet for s in args):
+    if S.EmptySet in args:
         return S.EmptySet
 
     # Handle Finite sets
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index c3f81669db..9eb4ba861b 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -2364,8 +2364,8 @@ def substitution(system, symbols, result=[{}], known_symbols=[],
     complements = {}
     intersections = {}
 
-    # when total_solveset_call is equals to total_conditionset
-    # means solvest fail to solve all the eq.
+    # when total_solveset_call equals total_conditionset
+    # it means that solveset failed to solve all eqs.
     total_conditionset = -1
     total_solveset_call = -1
 


Test patch:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index ebb6a98cf2..b4911a41de 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -898,14 +898,21 @@ def test_sympy__sets__setexpr__SetExpr():
     assert _test_args(SetExpr(Interval(0, 1)))
 
 
+def test_sympy__sets__fancysets__Rationals():
+    from sympy.sets.fancysets import Rationals
+    assert _test_args(Rationals())
+
+
 def test_sympy__sets__fancysets__Naturals():
     from sympy.sets.fancysets import Naturals
     assert _test_args(Naturals())
 
+
 def test_sympy__sets__fancysets__Naturals0():
     from sympy.sets.fancysets import Naturals0
     assert _test_args(Naturals0())
 
+
 def test_sympy__sets__fancysets__Integers():
     from sympy.sets.fancysets import Integers
     assert _test_args(Integers())
diff --git a/sympy/core/tests/test_function.py b/sympy/core/tests/test_function.py
index ff2d4a6ed2..a7bf674bb8 100644
--- a/sympy/core/tests/test_function.py
+++ b/sympy/core/tests/test_function.py
@@ -207,6 +207,9 @@ def test_Lambda():
     raises(TypeError, lambda: Lambda(1, x))
     assert Lambda(x, 1)(1) is S.One
 
+    raises(SyntaxError, lambda: Lambda((x, x), x + 2))
+
+
 
 def test_IdentityFunction():
     assert Lambda(x, x) is Lambda(y, y) is S.IdentityFunction
diff --git a/sympy/sets/tests/test_conditionset.py b/sympy/sets/tests/test_conditionset.py
index b976c2c919..16e34c9de1 100644
--- a/sympy/sets/tests/test_conditionset.py
+++ b/sympy/sets/tests/test_conditionset.py
@@ -167,4 +167,5 @@ def test_contains():
     assert ConditionSet(x, y > 5, Interval(1, 7)
         ).contains(8) is S.false
     assert ConditionSet(x, y > 5, Interval(1, 7)
-        ).contains(w) == And(w >= 1, w <= 7, y > 5)
+        ).contains(w) == And(S(1) <= w, w <= 7, y > 5)
+    assert 0 not in ConditionSet(x, 1/x >= 0, S.Reals)
diff --git a/sympy/sets/tests/test_contains.py b/sympy/sets/tests/test_contains.py
index 4efd8aee97..5e8d87d175 100644
--- a/sympy/sets/tests/test_contains.py
+++ b/sympy/sets/tests/test_contains.py
@@ -1,4 +1,5 @@
 from sympy import Symbol, Contains, S, Interval, FiniteSet, oo, Eq
+from sympy.core.expr import unchanged
 from sympy.utilities.pytest import raises
 
 def test_contains_basic():
@@ -12,8 +13,8 @@ def test_contains_basic():
 
 def test_issue_6194():
     x = Symbol('x')
-    assert Contains(x, Interval(0, 1)) != (x >= 0) & (x <= 1)
-    assert Interval(0, 1).contains(x) == (x >= 0) & (x <= 1)
+    assert unchanged(Contains, x, Interval(0, 1))
+    assert Interval(0, 1).contains(x) == (S(0) <= x) & (x <= 1)
     assert Contains(x, FiniteSet(0)) != S.false
     assert Contains(x, Interval(1, 1)) != S.false
     assert Contains(x, S.Integers) != S.false
diff --git a/sympy/sets/tests/test_fancysets.py b/sympy/sets/tests/test_fancysets.py
index 2764c1b7c8..87748fb5e8 100644
--- a/sympy/sets/tests/test_fancysets.py
+++ b/sympy/sets/tests/test_fancysets.py
@@ -1,5 +1,4 @@
 from sympy.core.compatibility import range, PY3
-from sympy.core.mod import Mod
 from sympy.sets.fancysets import (ImageSet, Range, normalize_theta_set,
                                   ComplexRegion)
 from sympy.sets.sets import (FiniteSet, Interval, imageset, Union,
@@ -32,6 +31,10 @@ def test_naturals():
 
     assert N.inf == 1
     assert N.sup == oo
+    assert not N.contains(oo)
+    for s in (S.Naturals0, S.Naturals):
+        assert s.intersection(S.Reals) is s
+        assert s.is_subset(S.Reals)
 
 
 def test_naturals0():
@@ -39,6 +42,7 @@ def test_naturals0():
     assert 0 in N
     assert -1 not in N
     assert next(iter(N)) == 0
+    assert not N.contains(oo)
 
 
 def test_integers():
@@ -46,6 +50,9 @@ def test_integers():
     assert 5 in Z
     assert -5 in Z
     assert 5.5 not in Z
+    assert not Z.contains(oo)
+    assert not Z.contains(-oo)
+
     zi = iter(Z)
     a, b, c, d = next(zi), next(zi), next(zi), next(zi)
     assert (a, b, c, d) == (0, 1, -1, 2)
@@ -64,7 +71,8 @@ def test_integers():
 
 def test_ImageSet():
     assert ImageSet(Lambda(x, 1), S.Integers) == FiniteSet(1)
-    assert ImageSet(Lambda(x, y), S.Integers) == FiniteSet(y)
+    assert ImageSet(Lambda(x, y), S.Integers
+        ) == {y}
     squares = ImageSet(Lambda(x, x**2), S.Naturals)
     assert 4 in squares
     assert 5 not in squares
@@ -385,6 +393,7 @@ def test_Reals():
     assert sqrt(-1) not in S.Reals
     assert S.Reals == Interval(-oo, oo)
     assert S.Reals != Interval(0, oo)
+    assert S.Reals.is_subset(Interval(-oo, oo))
 
 
 def test_Complex():
@@ -802,3 +811,31 @@ def test_issue_16871b():
 def test_no_mod_on_imaginary():
     assert imageset(Lambda(x, 2*x + 3*I), S.Integers
         ) == ImageSet(Lambda(x, 2*x + I), S.Integers)
+
+
+def test_Rationals():
+    assert S.Integers.is_subset(S.Rationals)
+    assert S.Naturals.is_subset(S.Rationals)
+    assert S.Naturals0.is_subset(S.Rationals)
+    assert S.Rationals.is_subset(S.Reals)
+    assert S.Rationals.inf == -oo
+    assert S.Rationals.sup == oo
+    it = iter(S.Rationals)
+    assert [next(it) for i in range(12)] == [
+        0, 1, -1, S(1)/2, 2, -S(1)/2, -2,
+        S(1)/3, 3, -S(1)/3, -3, S(2)/3]
+    assert Basic() not in S.Rationals
+    assert S.Half in S.Rationals
+    assert 1.0 not in S.Rationals
+    assert 2 in S.Rationals
+    r = symbols('r', rational=True)
+    assert r in S.Rationals
+    raises(TypeError, lambda: x in S.Rationals)
+
+
+def test_imageset_intersection():
+    n = Dummy()
+    s = ImageSet(Lambda(n, -I*(I*(2*pi*n - pi/4) +
+        log(Abs(sqrt(-I))))), S.Integers)
+    assert s.intersect(S.Reals) == ImageSet(
+        Lambda(n, 2*pi*n + 7*pi/4), S.Integers)
diff --git a/sympy/sets/tests/test_sets.py b/sympy/sets/tests/test_sets.py
index c839ceb4ed..67761ce004 100644
--- a/sympy/sets/tests/test_sets.py
+++ b/sympy/sets/tests/test_sets.py
@@ -3,7 +3,7 @@
     FiniteSet, Intersection, imageset, I, true, false, ProductSet, E,
     sqrt, Complement, EmptySet, sin, cos, Lambda, ImageSet, pi,
     Eq, Pow, Contains, Sum, rootof, SymmetricDifference, Piecewise,
-    Matrix, signsimp, Range, Add, symbols)
+    Matrix, signsimp, Range, Add, symbols, zoo)
 from mpmath import mpi
 
 from sympy.core.compatibility import range
@@ -14,9 +14,22 @@
 
 def test_imageset():
     ints = S.Integers
+    assert imageset(x, x - 1, S.Naturals) is S.Naturals0
+    assert imageset(x, x + 1, S.Naturals0) is S.Naturals
+    assert imageset(x, abs(x), S.Naturals0) is S.Naturals0
+    assert imageset(x, abs(x), S.Naturals) is S.Naturals
+    assert imageset(x, abs(x), S.Integers) is S.Naturals0
+    # issue 16878a
+    r = symbols('r', real=True)
+    assert (1, r) not in imageset(x, (x, x), S.Reals)
+    assert (r, r) in imageset(x, (x, x), S.Reals)
+    assert 1 + I in imageset(x, x + I, S.Reals)
+    assert {1} not in imageset(x, (x,), S.Reals)
+    assert (1, 1) not in imageset(x, (x,) , S.Reals)
     raises(TypeError, lambda: imageset(x, ints))
     raises(ValueError, lambda: imageset(x, y, z, ints))
     raises(ValueError, lambda: imageset(Lambda(x, cos(x)), y))
+    raises(ValueError, lambda: imageset(Lambda(x, x), ints, ints))
     assert imageset(cos, ints) == ImageSet(Lambda(x, cos(x)), ints)
     def f(x):
         return cos(x)
@@ -24,7 +37,8 @@ def f(x):
     f = lambda x: cos(x)
     assert imageset(f, ints) == ImageSet(Lambda(x, cos(x)), ints)
     assert imageset(x, 1, ints) == FiniteSet(1)
-    assert imageset(x, y, ints) == FiniteSet(y)
+    assert imageset(x, y, ints) == {y}
+    assert imageset((x, y), (1, z), ints*S.Reals) == {(1, z)}
     clash = Symbol('x', integer=true)
     assert (str(imageset(lambda x: x + clash, Interval(-2, 1)).lamda.expr)
         in ('_x + x', 'x + _x'))
@@ -250,8 +264,8 @@ def test_intersect1():
     assert all(i.intersection(S.Integers) is i for i in
         (S.Naturals, S.Naturals0))
     s =  S.Naturals0
-    assert S.Naturals.intersection(s) is s
-    assert s.intersection(S.Naturals) is s
+    assert S.Naturals.intersection(s) is S.Naturals
+    assert s.intersection(S.Naturals) is S.Naturals
     x = Symbol('x')
     assert Interval(0, 2).intersect(Interval(1, 2)) == Interval(1, 2)
     assert Interval(0, 2).intersect(Interval(1, 2, True)) == \
@@ -508,9 +522,9 @@ def test_contains():
 
     # non-bool results
     assert Union(Interval(1, 2), Interval(3, 4)).contains(x) == \
-        Or(And(x <= 2, x >= 1), And(x <= 4, x >= 3))
+        Or(And(S(1) <= x, x <= 2), And(S(3) <= x, x <= 4))
     assert Intersection(Interval(1, x), Interval(2, 3)).contains(y) == \
-        And(y <= 3, y <= x, y >= 1, y >= 2)
+        And(y <= 3, y <= x, S(1) <= y, S(2) <= y)
 
     assert (S.Complexes).contains(S.ComplexInfinity) == S.false
 
@@ -518,10 +532,10 @@ def test_contains():
 def test_interval_symbolic():
     x = Symbol('x')
     e = Interval(0, 1)
-    assert e.contains(x) == And(0 <= x, x <= 1)
+    assert e.contains(x) == And(S(0) <= x, x <= 1)
     raises(TypeError, lambda: x in e)
     e = Interval(0, 1, True, True)
-    assert e.contains(x) == And(0 < x, x < 1)
+    assert e.contains(x) == And(S(0) < x, x < 1)
 
 
 def test_union_contains():
@@ -529,9 +543,10 @@ def test_union_contains():
     i1 = Interval(0, 1)
     i2 = Interval(2, 3)
     i3 = Union(i1, i2)
+    assert i3.as_relational(x) == Or(And(S(0) <= x, x <= 1), And(S(2) <= x, x <= 3))
     raises(TypeError, lambda: x in i3)
     e = i3.contains(x)
-    assert e == Or(And(0 <= x, x <= 1), And(2 <= x, x <= 3))
+    assert e == i3.as_relational(x)
     assert e.subs(x, -0.5) is false
     assert e.subs(x, 0.5) is true
     assert e.subs(x, 1.5) is false
@@ -1000,7 +1015,7 @@ def test_issue_10113():
 
 def test_issue_10248():
     assert list(Intersection(S.Reals, FiniteSet(x))) == [
-        And(x < oo, x > -oo)]
+        (-oo < x) & (x < oo)]
 
 
 def test_issue_9447():
@@ -1114,3 +1129,16 @@ def test_union_intersection_constructor():
 
     assert Union({1}, {2}) == FiniteSet(1, 2)
     assert Intersection({1, 2}, {2, 3}) == FiniteSet(2)
+
+
+def test_Union_contains():
+    assert zoo not in Union(
+        Interval.open(-oo, 0), Interval.open(0, oo))
+
+
+@XFAIL
+def test_issue_16878b():
+    # in intersection_sets for (ImageSet, Set) there is no code
+    # that handles the base_set of S.Reals like there is
+    # for Integers
+    assert imageset(x, (x, x), S.Reals).is_subset(S.Reals**2) is True
diff --git a/sympy/vector/tests/test_coordsysrect.py b/sympy/vector/tests/test_coordsysrect.py
index 5f749c10e0..3f601a0f4c 100644
--- a/sympy/vector/tests/test_coordsysrect.py
+++ b/sympy/vector/tests/test_coordsysrect.py
@@ -436,7 +436,7 @@ def test_check_orthogonality():
     a = CoordSys3D('a', transformation=((u, v, z), (cosh(u) * cos(v), sinh(u) * sin(v), z)))
     assert a._check_orthogonality(a._transformation) is True
 
-    raises(ValueError, lambda: CoordSys3D('a', transformation=((x, x, z), (x, y, z))))
+    raises(ValueError, lambda: CoordSys3D('a', transformation=((x, y, z), (x, x, z))))
     raises(ValueError, lambda: CoordSys3D('a', transformation=(
         (x, y, z), (x*sin(y/2)*cos(z), x*sin(y)*sin(z), x*cos(y)))))