Repository: sympy/sympy

Issue:
ProductSet of 0 or 1 sets
The following inconsistency demonstrates a potential flaw in the definition of ProductSet:
```julia
In [10]: [(1,)*n in ProductSet(*(Reals,)*n) for n in range(10)]                                                                   
Out[10]: [False, False, True, True, True, True, True, True, True, True]
```
There are two problems:

1. The empty ProductSet should not be the empty set. It should be the set containing only the empty tuple `FiniteSet(())`:
https://en.wikipedia.org/wiki/Empty_product#Nullary_Cartesian_product
2. A ProductSet of one set evaluates to the underlying set meaning that there is no way to represent the set of 1-tuples of elements from the underlying set. I think this evaluation should be disallowed so that ProductSet(Reals) is the set of 1-tuples of real numbers. This problem doesn't have a direct mathematical analogue since this ambiguity around 1-tuples is common in mathematical literature. Working in Python though we need to have well-defined behaviour for this - looking through ImageSet and ProductSet it is clear that they cannot be fixed without having more precise programmatic definitions. One implication of this change would be that `Reals != Reals**1`.


Discussion:
sympy-bot: 

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v148). 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).

:x: There was an issue with the release notes. **Please do not close this pull request;** instead edit the description after reading the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

* No release notes were found. Please edit the PR description and write the release notes under `<!-- BEGIN RELEASE NOTES -->`.

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 -->

    <!-- END RELEASE NOTES -->


</details><p>

oscarbenjamin: This PR tightens the definition of ProductSet which I think is useful. It is now ready for review.
codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/17557?src=pr&el=h1) Report
> Merging [#17557](https://codecov.io/gh/sympy/sympy/pull/17557?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/edd9f10f8f4862ae481a65730abaaaaf66317224?src=pr&el=desc) will **increase** coverage by `12.734%`.
> The diff coverage is `79.661%`.

```diff
@@              Coverage Diff               @@
##            master    #17557        +/-   ##
==============================================
+ Coverage   62.029%   74.764%   +12.734%     
==============================================
  Files          633       633                
  Lines       164783    164798        +15     
  Branches     38720     38725         +5     
==============================================
+ Hits        102214    123210     +20996     
+ Misses       56755     36175     -20580     
+ Partials      5814      5413       -401
```


Patch:
diff --git a/sympy/core/containers.py b/sympy/core/containers.py
index e805179cf7..933f3202cd 100644
--- a/sympy/core/containers.py
+++ b/sympy/core/containers.py
@@ -10,6 +10,7 @@
 
 from collections import OrderedDict, defaultdict
 
+from sympy.core import S
 from sympy.core.basic import Basic
 from sympy.core.compatibility import as_int, range, MutableSet
 from sympy.core.sympify import sympify, converter
@@ -143,6 +144,24 @@ def index(self, value, start=None, stop=None):
         else:
             return self.args.index(value, start, stop)
 
+    def _eval_Eq(self, other):
+        from sympy.core.function import AppliedUndef
+        from sympy.core.logic import fuzzy_and, fuzzy_bool
+        from sympy.core.relational import Eq
+
+        if other.is_Symbol or isinstance(other, AppliedUndef):
+            return None
+
+        if not isinstance(other, Tuple) or len(self) != len(other):
+            return S.false
+
+        r = fuzzy_and(fuzzy_bool(Eq(s, o)) for s, o in zip(self, other))
+        if r is True:
+            return S.true
+        elif r is False:
+            return S.false
+
+
 converter[tuple] = lambda tup: Tuple(*tup)
 
 
diff --git a/sympy/core/relational.py b/sympy/core/relational.py
index 6dac83a600..eb6e37d1c4 100644
--- a/sympy/core/relational.py
+++ b/sympy/core/relational.py
@@ -449,6 +449,7 @@ class Equality(Relational):
 
     def __new__(cls, lhs, rhs=None, **options):
         from sympy.core.add import Add
+        from sympy.core.containers import Tuple
         from sympy.core.logic import fuzzy_bool
         from sympy.core.expr import _n2
         from sympy.simplify.simplify import clear_coefficients
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index be129d7571..bd38273d2f 100644
--- a/sympy/printing/latex.py
+++ b/sympy/printing/latex.py
@@ -2031,7 +2031,7 @@ def _print_SymmetricDifference(self, u):
 
     def _print_ProductSet(self, p):
         prec = precedence_traditional(p)
-        if len(p.sets) > 1 and not has_variety(p.sets):
+        if len(p.sets) >= 1 and not has_variety(p.sets):
             return self.parenthesize(p.sets[0], prec) + "^{%d}" % len(p.sets)
         return r" \times ".join(
             self.parenthesize(set, prec) for set in p.sets)
diff --git a/sympy/printing/pretty/pretty.py b/sympy/printing/pretty/pretty.py
index d79b5bc06d..a44f1a6619 100644
--- a/sympy/printing/pretty/pretty.py
+++ b/sympy/printing/pretty/pretty.py
@@ -1892,7 +1892,7 @@ def _print_Fraction(self, expr):
             return self.emptyPrinter(expr)
 
     def _print_ProductSet(self, p):
-        if len(p.sets) > 1 and not has_variety(p.sets):
+        if len(p.sets) >= 1 and not has_variety(p.sets):
             from sympy import Pow
             return self._print(Pow(p.sets[0], len(p.sets), evaluate=False))
         else:
diff --git a/sympy/sets/fancysets.py b/sympy/sets/fancysets.py
index 08c1834364..bd6b4eb9bf 100644
--- a/sympy/sets/fancysets.py
+++ b/sympy/sets/fancysets.py
@@ -305,13 +305,20 @@ class ImageSet(Set):
     """
     def __new__(cls, flambda, *sets):
         if not isinstance(flambda, Lambda):
-            raise ValueError('first argument must be a Lambda')
+            raise ValueError('First argument must be a Lambda')
+
+        sets = [_sympify(s) for s in sets]
 
         if flambda is S.IdentityFunction:
             if len(sets) != 1:
-                raise ValueError('identify function requires a single set')
+                raise ValueError('Identity function requires a single set')
             return sets[0]
 
+        if not all(isinstance(s, Set) for s in sets):
+            raise TypeError("Set arguments to ImageSet should of type Set")
+
+        sets = [s.flatten() if s.is_ProductSet else s for s in sets]
+
         if not set(flambda.variables) & flambda.expr.free_symbols:
             emptyprod = fuzzy_or(s.is_empty for s in sets)
             if emptyprod == True:
@@ -322,7 +329,14 @@ def __new__(cls, flambda, *sets):
         return Basic.__new__(cls, flambda, *sets)
 
     lamda = property(lambda self: self.args[0])
-    base_set = property(lambda self: ProductSet(self.args[1:]))
+
+    @property
+    def base_set(self):
+        sets = self.args[1:]
+        if len(sets) == 1:
+            return sets[0]
+        else:
+            return ProductSet(*self.args[1:]).flatten()
 
     def __iter__(self):
         already_seen = set()
diff --git a/sympy/sets/handlers/intersection.py b/sympy/sets/handlers/intersection.py
index 81690e42d3..e36cb7b95c 100644
--- a/sympy/sets/handlers/intersection.py
+++ b/sympy/sets/handlers/intersection.py
@@ -351,8 +351,8 @@ def intersection_sets(self, other):
 def intersection_sets(a, b):
     if len(b.args) != len(a.args):
         return S.EmptySet
-    return ProductSet(i.intersect(j)
-            for i, j in zip(a.sets, b.sets))
+    return ProductSet(*(i.intersect(j) for i, j in zip(a.sets, b.sets)))
+
 
 @dispatch(Interval, Interval)
 def intersection_sets(a, b):
diff --git a/sympy/sets/handlers/union.py b/sympy/sets/handlers/union.py
index 743f8755f3..e84befe068 100644
--- a/sympy/sets/handlers/union.py
+++ b/sympy/sets/handlers/union.py
@@ -41,14 +41,15 @@ def union_sets(a, b):
 def union_sets(a, b):
     if b.is_subset(a):
         return a
-    if len(b.args) != len(a.args):
+    if len(b.sets) != len(a.sets):
         return None
-    if a.args[0] == b.args[0]:
-        return a.args[0] * Union(ProductSet(a.args[1:]),
-                                    ProductSet(b.args[1:]))
-    if a.args[-1] == b.args[-1]:
-        return Union(ProductSet(a.args[:-1]),
-                     ProductSet(b.args[:-1])) * a.args[-1]
+    if len(a.sets) == 2:
+        a1, a2 = a.sets
+        b1, b2 = b.sets
+        if a1 == b1:
+            return a1 * Union(a2, b2)
+        if a2 == b2:
+            return Union(a1, b1) * a2
     return None
 
 @dispatch(ProductSet, Set)
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index f6e0d8b898..2f6eb3f042 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -8,6 +8,7 @@
 from sympy.core.compatibility import (iterable, with_metaclass,
     ordered, range, PY3, is_sequence, reduce)
 from sympy.core.cache import cacheit
+from sympy.core.containers import Tuple
 from sympy.core.decorators import deprecated
 from sympy.core.evalf import EvalfMixin
 from sympy.core.evaluate import global_evaluate
@@ -24,6 +25,7 @@
 from sympy.logic.boolalg import And, Or, Not, true, false
 from sympy.sets.contains import Contains
 from sympy.utilities import subsets
+from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.utilities.iterables import sift
 from sympy.utilities.misc import func_name, filldedent
 
@@ -191,21 +193,19 @@ def complement(self, universe):
 
     def _complement(self, other):
         # this behaves as other - self
-        if isinstance(other, ProductSet):
-            # For each set consider it or it's complement
-            # We need at least one of the sets to be complemented
-            # Consider all 2^n combinations.
-            # We can conveniently represent these options easily using a
-            # ProductSet
-
-            # XXX: this doesn't work if the dimensions of the sets isn't same.
-            # A - B is essentially same as A if B has a different
-            # dimensionality than A
-            switch_sets = ProductSet(FiniteSet(o, o - s) for s, o in
-                                     zip(self.sets, other.sets))
-            product_sets = (ProductSet(*set) for set in switch_sets)
-            # Union of all combinations but this one
-            return Union(*(p for p in product_sets if p != other))
+        if isinstance(self, ProductSet) and isinstance(other, ProductSet):
+            # If self and other are disjoint then other - self == self
+            if len(self.sets) != len(other.sets):
+                return other
+
+            # There can be other ways to represent this but this gives:
+            # (A x B) - (C x D) = ((A - C) x B) U (A x (B - D))
+            overlaps = []
+            pairs = list(zip(self.sets, other.sets))
+            for n in range(len(pairs)):
+                sets = (o if i != n else o-s for i, (s, o) in enumerate(pairs))
+                overlaps.append(ProductSet(*sets))
+            return Union(*overlaps)
 
         elif isinstance(other, Interval):
             if isinstance(self, Interval) or isinstance(self, FiniteSet):
@@ -612,9 +612,9 @@ def __xor__(self, other):
         return SymmetricDifference(self, other)
 
     def __pow__(self, exp):
-        if not sympify(exp).is_Integer and exp >= 0:
+        if not (sympify(exp).is_Integer and exp >= 0):
             raise ValueError("%s: Exponent must be a positive Integer" % exp)
-        return ProductSet([self]*exp)
+        return ProductSet(*[self]*exp)
 
     def __sub__(self, other):
         return Complement(self, other)
@@ -655,12 +655,17 @@ class ProductSet(Set):
     >>> set(coin**2)
     {(H, H), (H, T), (T, H), (T, T)}
 
+    The Cartesian product is not commutative or associative e.g.:
+
+    >>> I*S == S*I
+    False
+    >>> (I*I)*I == I*(I*I)
+    False
 
     Notes
     =====
 
     - Passes most operations down to the argument sets
-    - Flattens Products of ProductSets
 
     References
     ==========
@@ -670,33 +675,52 @@ class ProductSet(Set):
     is_ProductSet = True
 
     def __new__(cls, *sets, **assumptions):
-        def flatten(arg):
-            if isinstance(arg, Set):
-                if arg.is_ProductSet:
-                    return sum(map(flatten, arg.args), [])
-                else:
-                    return [arg]
-            elif iterable(arg):
-                return sum(map(flatten, arg), [])
-            raise TypeError("Input must be Sets or iterables of Sets")
-        sets = flatten(list(sets))
+        if len(sets) == 1 and iterable(sets[0]) and not isinstance(sets[0], (Set, set)):
+            SymPyDeprecationWarning(
+                feature="ProductSet(iterable)",
+                useinstead="ProductSet(*iterable)",
+                issue=17557,
+                deprecated_since_version="1.5"
+            ).warn()
+            sets = tuple(sets[0])
 
-        if EmptySet() in sets or len(sets) == 0:
-            return EmptySet()
+        sets = [sympify(s) for s in sets]
 
-        if len(sets) == 1:
-            return sets[0]
+        if not all(isinstance(s, Set) for s in sets):
+            raise TypeError("Arguments to ProductSet should be of type Set")
+
+        # Nullary product of sets is *not* the empty set
+        if len(sets) == 0:
+            return FiniteSet(())
+
+        if EmptySet() in sets:
+            return EmptySet()
 
         return Basic.__new__(cls, *sets, **assumptions)
 
+    @property
+    def sets(self):
+        return self.args
+
+    def flatten(self):
+        def _flatten(sets):
+            for s in sets:
+                if s.is_ProductSet:
+                    for s2 in _flatten(s.sets):
+                        yield s2
+                else:
+                    yield s
+        return ProductSet(*_flatten(self.sets))
+
     def _eval_Eq(self, other):
         if not other.is_ProductSet:
             return
 
-        if len(self.args) != len(other.args):
+        if len(self.sets) != len(other.sets):
             return false
 
-        return And(*(Eq(x, y) for x, y in zip(self.args, other.args)))
+        eqs = (Eq(x, y) for x, y in zip(self.sets, other.sets))
+        return tfn[fuzzy_and(map(fuzzy_bool, eqs))]
 
     def _contains(self, element):
         """
@@ -714,30 +738,26 @@ def _contains(self, element):
 
         Passes operation on to constituent sets
         """
-        if is_sequence(element):
-            if len(element) != len(self.args):
-                return False
-        elif len(self.args) > 1:
+        if element.is_Symbol:
+            return None
+
+        if not isinstance(element, Tuple) or len(element) != len(self.sets):
             return False
-        d = [Dummy() for i in element]
-        reps = dict(zip(d, element))
-        return tfn[self.as_relational(*d).xreplace(reps)]
+
+        return fuzzy_and(s._contains(e) for s, e in zip(self.sets, element))
 
     def as_relational(self, *symbols):
-        if len(symbols) != len(self.args) or not all(
+        symbols = [_sympify(s) for s in symbols]
+        if len(symbols) != len(self.sets) 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):
-        return self.args
+        return And(*[s.as_relational(i) for s, i in zip(self.sets, symbols)])
 
     @property
     def _boundary(self):
-        return Union(*(ProductSet(b + b.boundary if i != j else b.boundary
-                                for j, b in enumerate(self.sets))
+        return Union(*(ProductSet(*(b + b.boundary if i != j else b.boundary
+                                for j, b in enumerate(self.sets)))
                                 for i, a in enumerate(self.sets)))
 
     @property
@@ -778,15 +798,15 @@ def is_empty(self):
     @property
     def _measure(self):
         measure = 1
-        for set in self.sets:
-            measure *= set.measure
+        for s in self.sets:
+            measure *= s.measure
         return measure
 
     def __len__(self):
-        return Mul(*[len(s) for s in self.args])
+        return Mul(*[len(s) for s in self.sets])
 
     def __bool__(self):
-        return all([bool(s) for s in self.args])
+        return all([bool(s) for s in self.sets])
 
     __nonzero__ = __bool__
 
@@ -1693,7 +1713,7 @@ def all_in_both():
             yield fuzzy_and(self._contains(e) for e in o_set - s_set)
             yield fuzzy_and(other._contains(e) for e in s_set - o_set)
 
-        return fuzzy_and(all_in_both())
+        return tfn[fuzzy_and(all_in_both())]
 
     def __iter__(self):
         return iter(self.args)
@@ -1756,7 +1776,7 @@ def _contains(self, other):
         """
         # 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])
+        return fuzzy_or(fuzzy_bool(Eq(e, other, evaluate=True)) for e in self.args)
 
     @property
     def _boundary(self):
diff --git a/sympy/stats/joint_rv_types.py b/sympy/stats/joint_rv_types.py
index a491e54102..3faba969a5 100644
--- a/sympy/stats/joint_rv_types.py
+++ b/sympy/stats/joint_rv_types.py
@@ -361,7 +361,7 @@ def set(self):
         prod_set = Range(0, self.n + 1)
         for i in range(2, self.n + 1):
             prod_set *= Range(0, self.n//i + 1)
-        return prod_set
+        return prod_set.flatten()
 
     def pdf(self, *syms):
         n, theta = self.n, self.theta
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index cef6633c67..78129519f0 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -510,7 +510,7 @@ def domains(self):
 
     @property
     def set(self):
-        return ProductSet(domain.set for domain in self.domains)
+        return ProductSet(*(domain.set for domain in self.domains))
 
     def __contains__(self, other):
         # Split event into each subdomain


Test patch:
diff --git a/sympy/core/tests/test_containers.py b/sympy/core/tests/test_containers.py
index af871de597..af20f181c8 100644
--- a/sympy/core/tests/test_containers.py
+++ b/sympy/core/tests/test_containers.py
@@ -1,8 +1,14 @@
 from collections import defaultdict
+
 from sympy import Matrix, Tuple, symbols, sympify, Basic, Dict, S, FiniteSet, Integer
+from sympy.core.compatibility import is_sequence, iterable, range
 from sympy.core.containers import tuple_wrapper
+from sympy.core.expr import unchanged
+from sympy.core.function import Function, Lambda
+from sympy.core.relational import Eq
 from sympy.utilities.pytest import raises
-from sympy.core.compatibility import is_sequence, iterable, range
+
+from sympy.abc import x, y
 
 
 def test_Tuple():
@@ -63,6 +69,21 @@ def test_Tuple_equality():
     assert (Tuple(1, 2) != Tuple(1, 3)) is True
 
 
+def test_Tuple_Eq():
+    assert Eq(Tuple(), Tuple()) is S.true
+    assert Eq(Tuple(1), 1) is S.false
+    assert Eq(Tuple(1, 2), Tuple(1)) is S.false
+    assert Eq(Tuple(1), Tuple(1)) is S.true
+    assert Eq(Tuple(1, 2), Tuple(1, 3)) is S.false
+    assert Eq(Tuple(1, 2), Tuple(1, 2)) is S.true
+    assert unchanged(Eq, Tuple(1, x), Tuple(1, 2))
+    assert Eq(Tuple(1, x), Tuple(1, 2)).subs(x, 2) is S.true
+    assert unchanged(Eq, Tuple(1, 2), x)
+    f = Function('f')
+    assert unchanged(Eq, Tuple(1), f(x))
+    assert Eq(Tuple(1), f(x)).subs(x, 1).subs(f, Lambda(y, (y,))) is S.true
+
+
 def test_Tuple_comparision():
     assert (Tuple(1, 3) >= Tuple(-10, 30)) is S.true
     assert (Tuple(1, 3) <= Tuple(-10, 30)) is S.false
diff --git a/sympy/core/tests/test_relational.py b/sympy/core/tests/test_relational.py
index cb9f70c3c9..f6e303b68b 100644
--- a/sympy/core/tests/test_relational.py
+++ b/sympy/core/tests/test_relational.py
@@ -129,6 +129,9 @@ def test_Eq():
     # issue 13348
     assert Eq(True, 1) is S.false
 
+    assert Eq((), 1) is S.false
+
+
 def test_rel_Infinity():
     # NOTE: All of these are actually handled by sympy.core.Number, and do
     # not create Relational objects.
diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index 479b1beda0..dcb49b94ca 100644
--- a/sympy/printing/pretty/tests/test_pretty.py
+++ b/sympy/printing/pretty/tests/test_pretty.py
@@ -3858,7 +3858,13 @@ def test_pretty_Intersection_issue_10414():
     assert upretty(Intersection(a, b)) == ucode_str
     assert pretty(Intersection(a, b)) == ascii_str
 
-def test_ProductSet_paranthesis():
+def test_ProductSet_exponent():
+    ucode_str = '      1\n[0, 1] '
+    assert upretty(Interval(0, 1)**1) == ucode_str
+    ucode_str = '      2\n[0, 1] '
+    assert upretty(Interval(0, 1)**2) == ucode_str
+
+def test_ProductSet_parenthesis():
     ucode_str = u'([4, 7] × {1, 2}) ∪ ([2, 3] × [4, 7])'
 
     a, b, c = Interval(2, 3), Interval(4, 7), Interval(1, 9)
diff --git a/sympy/printing/tests/test_latex.py b/sympy/printing/tests/test_latex.py
index 402fd7b1fd..fbff7ae2e0 100644
--- a/sympy/printing/tests/test_latex.py
+++ b/sympy/printing/tests/test_latex.py
@@ -876,7 +876,7 @@ def test_latex_productset():
     fset = FiniteSet(1, 2, 3)
     assert latex(line**2) == r"%s^{2}" % latex(line)
     assert latex(line**10) == r"%s^{10}" % latex(line)
-    assert latex(line * bigline * fset) == r"%s \times %s \times %s" % (
+    assert latex((line * bigline * fset).flatten()) == r"%s \times %s \times %s" % (
         latex(line), latex(bigline), latex(fset))
 
 
@@ -978,7 +978,7 @@ def test_set_operators_parenthesis():
         '\\left\\{d\\right\\}\\right)'
 
     # XXX This can be incorrect since cartesian product is not associative
-    assert latex(ProductSet(A, P2)) == \
+    assert latex(ProductSet(A, P2).flatten()) == \
         '\\left\\{a\\right\\} \\times \\left\\{c\\right\\} \\times ' \
         '\\left\\{d\\right\\}'
     assert latex(ProductSet(U1, U2)) == \
diff --git a/sympy/sets/tests/test_fancysets.py b/sympy/sets/tests/test_fancysets.py
index f049ce7c7d..f3159d25c3 100644
--- a/sympy/sets/tests/test_fancysets.py
+++ b/sympy/sets/tests/test_fancysets.py
@@ -121,6 +121,11 @@ def test_ImageSet():
     assert imageset(lambda x, y: x + y, S.Integers, S.Naturals
         ).base_set == ProductSet(S.Integers, S.Naturals)
 
+    # Passing a set instead of a FiniteSet shouldn't raise
+    assert unchanged(ImageSet, Lambda(x, x**2), {1, 2, 3})
+
+    raises(TypeError, lambda: ImageSet(Lambda(x, x**2), 1))
+
 
 def test_image_is_ImageSet():
     assert isinstance(imageset(x, sqrt(sin(x)), Range(5)), ImageSet)
diff --git a/sympy/sets/tests/test_sets.py b/sympy/sets/tests/test_sets.py
index 0b2fc10e29..85c00587f6 100644
--- a/sympy/sets/tests/test_sets.py
+++ b/sympy/sets/tests/test_sets.py
@@ -258,6 +258,15 @@ def test_Complement():
     assert Complement(FiniteSet(x, y, 2), Interval(-10, 10)) == \
             Complement(FiniteSet(x, y), Interval(-10, 10))
 
+    A = FiniteSet(*symbols('a:c'))
+    B = FiniteSet(*symbols('d:f'))
+    assert unchanged(Complement, ProductSet(A, A), B)
+
+    A2 = ProductSet(A, A)
+    B3 = ProductSet(B, B, B)
+    assert A2 - B3 == A2
+    assert B3 - A2 == B3
+
 
 def test_complement():
     assert Interval(0, 1).complement(S.Reals) == \
@@ -422,8 +431,85 @@ def test_is_disjoint():
     assert Interval(0, 2).is_disjoint(Interval(3, 4)) == True
 
 
-def test_ProductSet_of_single_arg_is_arg():
-    assert ProductSet(Interval(0, 1)) == Interval(0, 1)
+def test_ProductSet():
+    # ProductSet is always a set of Tuples
+    assert ProductSet(S.Reals) == S.Reals ** 1
+    assert ProductSet(S.Reals, S.Reals) == S.Reals ** 2
+    assert ProductSet(S.Reals, S.Reals, S.Reals) == S.Reals ** 3
+
+    assert ProductSet(S.Reals) != S.Reals
+    assert ProductSet(S.Reals, S.Reals) == S.Reals * S.Reals
+    assert ProductSet(S.Reals, S.Reals, S.Reals) != S.Reals * S.Reals * S.Reals
+    assert ProductSet(S.Reals, S.Reals, S.Reals) == (S.Reals * S.Reals * S.Reals).flatten()
+
+    assert 1 not in ProductSet(S.Reals)
+    assert (1,) in ProductSet(S.Reals)
+
+    assert 1 not in ProductSet(S.Reals, S.Reals)
+    assert (1, 2) in ProductSet(S.Reals, S.Reals)
+    assert (1, I) not in ProductSet(S.Reals, S.Reals)
+
+    assert (1, 2, 3) in ProductSet(S.Reals, S.Reals, S.Reals)
+    assert (1, 2, 3) in S.Reals ** 3
+    assert (1, 2, 3) not in S.Reals * S.Reals * S.Reals
+    assert ((1, 2), 3) in S.Reals * S.Reals * S.Reals
+    assert (1, (2, 3)) not in S.Reals * S.Reals * S.Reals
+    assert (1, (2, 3)) in S.Reals * (S.Reals * S.Reals)
+
+    assert ProductSet() == FiniteSet(())
+    assert ProductSet(S.Reals, S.EmptySet) == S.EmptySet
+
+    # See GH-17458
+
+    for n in range(5):
+        Rn = ProductSet(*(S.Reals,) * n)
+        assert (1,) * n in Rn
+        assert 1 not in Rn
+
+    assert (S.Reals * S.Reals) * S.Reals != S.Reals * (S.Reals * S.Reals)
+
+    S1 = S.Reals
+    S2 = S.Integers
+    x1 = pi
+    x2 = 3
+    assert x1 in S1
+    assert x2 in S2
+    assert (x1, x2) in S1 * S2
+    S3 = S1 * S2
+    x3 = (x1, x2)
+    assert x3 in S3
+    assert (x3, x3) in S3 * S3
+    assert x3 + x3 not in S3 * S3
+
+    raises(ValueError, lambda: S.Reals**-1)
+    with warns_deprecated_sympy():
+        ProductSet(FiniteSet(s) for s in range(2))
+    raises(TypeError, lambda: ProductSet(None))
+
+    S1 = FiniteSet(1, 2)
+    S2 = FiniteSet(3, 4)
+    S3 = ProductSet(S1, S2)
+    assert (S3.as_relational(x, y)
+            == And(S1.as_relational(x), S2.as_relational(y))
+            == And(Or(Eq(x, 1), Eq(x, 2)), Or(Eq(y, 3), Eq(y, 4))))
+    raises(ValueError, lambda: S3.as_relational(x))
+    raises(ValueError, lambda: S3.as_relational(x, 1))
+    raises(ValueError, lambda: ProductSet(Interval(0, 1)).as_relational(x, y))
+
+    Z2 = ProductSet(S.Integers, S.Integers)
+    assert Z2.contains((1, 2)) is S.true
+    assert Z2.contains((1,)) is S.false
+    assert Z2.contains(x) == Contains(x, Z2, evaluate=False)
+    assert Z2.contains(x).subs(x, 1) is S.false
+    assert Z2.contains((x, 1)).subs(x, 2) is S.true
+    assert Z2.contains((x, y)) == Contains((x, y), Z2, evaluate=False)
+    assert unchanged(Contains, (x, y), Z2)
+    assert Contains((1, 2), Z2) is S.true
+
+
+def test_ProductSet_of_single_arg_is_not_arg():
+    assert unchanged(ProductSet, Interval(0, 1))
+    assert ProductSet(Interval(0, 1)) != Interval(0, 1)
 
 
 def test_ProductSet_is_empty():
@@ -773,8 +859,10 @@ def test_product_basic():
     assert (0, 0) in square
     assert 0 not in square
     assert (H, T) in coin ** 2
-    assert (.5, .5, .5) in square * unit_line
-    assert (H, 3, 3) in coin * d6* d6
+    assert (.5, .5, .5) in (square * unit_line).flatten()
+    assert ((.5, .5), .5) in square * unit_line
+    assert (H, 3, 3) in (coin * d6 * d6).flatten()
+    assert ((H, 3), 3) in coin * d6 * d6
     HH, TT = sympify(H), sympify(T)
     assert set(coin**2) == set(((HH, HH), (HH, TT), (TT, HH), (TT, TT)))
 
@@ -1014,10 +1102,20 @@ def test_Eq():
     assert Eq(s1*s2, s2*s1) == False
 
     assert unchanged(Eq, FiniteSet({x, y}), FiniteSet({x}))
-    assert Eq(FiniteSet({x, y}).subs(y, x), FiniteSet({x})) is True
-    assert Eq(FiniteSet({x, y}), FiniteSet({x})).subs(y, x) is True
-    assert Eq(FiniteSet({x, y}).subs(y, x+1), FiniteSet({x})) is False
-    assert Eq(FiniteSet({x, y}), FiniteSet({x})).subs(y, x+1) is False
+    assert Eq(FiniteSet({x, y}).subs(y, x), FiniteSet({x})) is S.true
+    assert Eq(FiniteSet({x, y}), FiniteSet({x})).subs(y, x) is S.true
+    assert Eq(FiniteSet({x, y}).subs(y, x+1), FiniteSet({x})) is S.false
+    assert Eq(FiniteSet({x, y}), FiniteSet({x})).subs(y, x+1) is S.false
+
+    assert Eq(ProductSet({1}, {2}), Interval(1, 2)) not in (S.true, S.false)
+    assert Eq(ProductSet({1}), ProductSet({1}, {2})) is S.false
+
+    assert Eq(FiniteSet(()), FiniteSet(1)) is S.false
+    assert Eq(ProductSet(), FiniteSet(1)) is S.false
+
+    i1 = Interval(0, 1)
+    i2 = Interval(x, y)
+    assert unchanged(Eq, ProductSet(i1, i1), ProductSet(i2, i2))
 
 
 def test_SymmetricDifference():