Repository: sympy/sympy

Issue:
digamma and trigamma should be Function subclasses
Right now `digamma` and `trigamma` are Python functions. They should be Function subclasses, so they can be unevaluated. 

Discussion:

Patch:
diff --git a/sympy/core/containers.py b/sympy/core/containers.py
index 933f3202cd..e805179cf7 100644
--- a/sympy/core/containers.py
+++ b/sympy/core/containers.py
@@ -10,7 +10,6 @@
 
 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
@@ -144,24 +143,6 @@ 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 894bd54222..8750d4708e 100644
--- a/sympy/core/relational.py
+++ b/sympy/core/relational.py
@@ -455,7 +455,6 @@ 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/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index 5a94e17adc..a05d3ee6f9 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -992,7 +992,7 @@ def _sage_(self):
         return sage.log_gamma(self.args[0]._sage_())
 
 
-def digamma(x):
+class digamma(Function):
     r"""
     The digamma function is the first derivative of the loggamma function i.e,
 
@@ -1020,10 +1020,52 @@ def digamma(x):
     .. [2] http://mathworld.wolfram.com/DigammaFunction.html
     .. [3] http://functions.wolfram.com/GammaBetaErf/PolyGamma2/
     """
-    return polygamma(0, x)
+    def _eval_evalf(self, prec):
+        z = self.args[0]
+        return polygamma(0, z).evalf(prec)
+
+    def fdiff(self, argindex=1):
+        z = self.args[0]
+        return polygamma(0, z).fdiff()
+
+    def _eval_is_real(self):
+        z = self.args[0]
+        return polygamma(0, z).is_real
 
+    def _eval_is_positive(self):
+        z = self.args[0]
+        return polygamma(0, z).is_positive
+
+    def _eval_is_negative(self):
+        z = self.args[0]
+        return polygamma(0, z).is_negative
+
+    def _eval_aseries(self, n, args0, x, logx):
+        as_polygamma = self.rewrite(polygamma)
+        args0 = [S.Zero,] + args0
+        return as_polygamma._eval_aseries(n, args0, x, logx)
 
-def trigamma(x):
+    @classmethod
+    def eval(cls, z):
+        return polygamma(0, z)
+
+    def _eval_expand_func(self, **hints):
+        z = self.args[0]
+        return polygamma(0, z).expand(func=True)
+
+    def _eval_rewrite_as_harmonic(self, z, **kwargs):
+        return harmonic(z - 1) - S.EulerGamma
+
+    def _eval_rewrite_as_polygamma(self, z, **kwargs):
+        return polygamma(0, z)
+
+    def _eval_as_leading_term(self, x):
+        z = self.args[0]
+        return polygamma(0, z).as_leading_term(x)
+
+
+
+class trigamma(Function):
     r"""
     The trigamma function is the second derivative of the loggamma function i.e,
 
@@ -1050,7 +1092,52 @@ def trigamma(x):
     .. [2] http://mathworld.wolfram.com/TrigammaFunction.html
     .. [3] http://functions.wolfram.com/GammaBetaErf/PolyGamma2/
     """
-    return polygamma(1, x)
+    def _eval_evalf(self, prec):
+        z = self.args[0]
+        return polygamma(1, z).evalf(prec)
+
+    def fdiff(self, argindex=1):
+        z = self.args[0]
+        return polygamma(1, z).fdiff()
+
+    def _eval_is_real(self):
+        z = self.args[0]
+        return polygamma(1, z).is_real
+
+    def _eval_is_positive(self):
+        z = self.args[0]
+        return polygamma(1, z).is_positive
+
+    def _eval_is_negative(self):
+        z = self.args[0]
+        return polygamma(1, z).is_negative
+
+    def _eval_aseries(self, n, args0, x, logx):
+        as_polygamma = self.rewrite(polygamma)
+        args0 = [S.One,] + args0
+        return as_polygamma._eval_aseries(n, args0, x, logx)
+
+    @classmethod
+    def eval(cls, z):
+        return polygamma(1, z)
+
+    def _eval_expand_func(self, **hints):
+        z = self.args[0]
+        return polygamma(1, z).expand(func=True)
+
+    def _eval_rewrite_as_zeta(self, z, **kwargs):
+        return zeta(2, z)
+
+    def _eval_rewrite_as_polygamma(self, z, **kwargs):
+        return polygamma(1, z)
+
+    def _eval_rewrite_as_harmonic(self, z, **kwargs):
+        return -harmonic(z - 1, 2) + S.Pi**2 / 6
+
+    def _eval_as_leading_term(self, x):
+        z = self.args[0]
+        return polygamma(1, z).as_leading_term(x)
+
 
 ###############################################################################
 ##################### COMPLETE MULTIVARIATE GAMMA FUNCTION ####################
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index bd38273d2f..be129d7571 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 a44f1a6619..d79b5bc06d 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 5c6986a756..d80f7bff2b 100644
--- a/sympy/sets/fancysets.py
+++ b/sympy/sets/fancysets.py
@@ -305,20 +305,13 @@ class ImageSet(Set):
     """
     def __new__(cls, flambda, *sets):
         if not isinstance(flambda, Lambda):
-            raise ValueError('First argument must be a Lambda')
-
-        sets = [_sympify(s) for s in sets]
+            raise ValueError('first argument must be a Lambda')
 
         if flambda is S.IdentityFunction:
             if len(sets) != 1:
-                raise ValueError('Identity function requires a single set')
+                raise ValueError('identify 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:
@@ -329,14 +322,7 @@ def __new__(cls, flambda, *sets):
         return Basic.__new__(cls, flambda, *sets)
 
     lamda = property(lambda self: self.args[0])
-
-    @property
-    def base_set(self):
-        sets = self.args[1:]
-        if len(sets) == 1:
-            return sets[0]
-        else:
-            return ProductSet(*self.args[1:]).flatten()
+    base_set = property(lambda self: ProductSet(self.args[1:]))
 
     def __iter__(self):
         already_seen = set()
diff --git a/sympy/sets/handlers/intersection.py b/sympy/sets/handlers/intersection.py
index 6c8e0d7399..3cee005286 100644
--- a/sympy/sets/handlers/intersection.py
+++ b/sympy/sets/handlers/intersection.py
@@ -350,8 +350,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 e84befe068..743f8755f3 100644
--- a/sympy/sets/handlers/union.py
+++ b/sympy/sets/handlers/union.py
@@ -41,15 +41,14 @@ def union_sets(a, b):
 def union_sets(a, b):
     if b.is_subset(a):
         return a
-    if len(b.sets) != len(a.sets):
+    if len(b.args) != len(a.args):
         return None
-    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
+    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]
     return None
 
 @dispatch(ProductSet, Set)
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index db9b7f2d06..13096ac933 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -8,7 +8,6 @@
 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,7 +23,6 @@
 from sympy.logic.boolalg import And, Or, Not, Xor, 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
 
@@ -192,19 +190,21 @@ def complement(self, universe):
 
     def _complement(self, other):
         # this behaves as other - self
-        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)
+        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))
 
         elif isinstance(other, Interval):
             if isinstance(self, Interval) or isinstance(self, FiniteSet):
@@ -611,9 +611,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)
@@ -654,17 +654,12 @@ 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
     ==========
@@ -674,52 +669,33 @@ class ProductSet(Set):
     is_ProductSet = True
 
     def __new__(cls, *sets, **assumptions):
-        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])
-
-        sets = [sympify(s) for s in sets]
-
-        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(())
+        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 EmptySet() in sets:
+        if EmptySet() in sets or len(sets) == 0:
             return EmptySet()
 
-        return Basic.__new__(cls, *sets, **assumptions)
-
-    @property
-    def sets(self):
-        return self.args
+        if len(sets) == 1:
+            return sets[0]
 
-    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))
+        return Basic.__new__(cls, *sets, **assumptions)
 
     def _eval_Eq(self, other):
         if not other.is_ProductSet:
             return
 
-        if len(self.sets) != len(other.sets):
+        if len(self.args) != len(other.args):
             return false
 
-        eqs = (Eq(x, y) for x, y in zip(self.sets, other.sets))
-        return tfn[fuzzy_and(map(fuzzy_bool, eqs))]
+        return And(*(Eq(x, y) for x, y in zip(self.args, other.args)))
 
     def _contains(self, element):
         """
@@ -737,26 +713,30 @@ def _contains(self, element):
 
         Passes operation on to constituent sets
         """
-        if element.is_Symbol:
-            return None
-
-        if not isinstance(element, Tuple) or len(element) != len(self.sets):
+        if is_sequence(element):
+            if len(element) != len(self.args):
+                return False
+        elif len(self.args) > 1:
             return False
-
-        return fuzzy_and(s._contains(e) for s, e in zip(self.sets, element))
+        d = [Dummy() for i in element]
+        reps = dict(zip(d, element))
+        return tfn[self.as_relational(*d).xreplace(reps)]
 
     def as_relational(self, *symbols):
-        symbols = [_sympify(s) for s in symbols]
-        if len(symbols) != len(self.sets) or not all(
+        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.as_relational(i) for s, i in zip(self.sets, symbols)])
+        return And(*[s.contains(i) for s, i in zip(self.args, symbols)])
+
+    @property
+    def sets(self):
+        return self.args
 
     @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
@@ -797,15 +777,15 @@ def is_empty(self):
     @property
     def _measure(self):
         measure = 1
-        for s in self.sets:
-            measure *= s.measure
+        for set in self.sets:
+            measure *= set.measure
         return measure
 
     def __len__(self):
-        return Mul(*[len(s) for s in self.sets])
+        return Mul(*[len(s) for s in self.args])
 
     def __bool__(self):
-        return all([bool(s) for s in self.sets])
+        return all([bool(s) for s in self.args])
 
     __nonzero__ = __bool__
 
@@ -1722,7 +1702,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 tfn[fuzzy_and(all_in_both())]
+        return fuzzy_and(all_in_both())
 
     def __iter__(self):
         return iter(self.args)
@@ -1785,7 +1765,7 @@ def _contains(self, other):
         """
         # evaluate=True is needed to override evaluate=False context;
         # we need Eq to do the evaluation
-        return fuzzy_or(fuzzy_bool(Eq(e, other, evaluate=True)) for e in self.args)
+        return fuzzy_or([tfn[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 3faba969a5..a491e54102 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.flatten()
+        return prod_set
 
     def pdf(self, *syms):
         n, theta = self.n, self.theta
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index d8a0c50da0..9520ad70f2 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -509,7 +509,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_args.py b/sympy/core/tests/test_args.py
index a8e3d511bc..67ee3df3e9 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -2288,6 +2288,13 @@ def test_sympy__functions__special__gamma_functions__polygamma():
     from sympy.functions.special.gamma_functions import polygamma
     assert _test_args(polygamma(x, 2))
 
+def test_sympy__functions__special__gamma_functions__digamma():
+    from sympy.functions.special.gamma_functions import digamma
+    assert _test_args(digamma(x))
+
+def test_sympy__functions__special__gamma_functions__trigamma():
+    from sympy.functions.special.gamma_functions import trigamma
+    assert _test_args(trigamma(x))
 
 def test_sympy__functions__special__gamma_functions__uppergamma():
     from sympy.functions.special.gamma_functions import uppergamma
diff --git a/sympy/core/tests/test_containers.py b/sympy/core/tests/test_containers.py
index af20f181c8..af871de597 100644
--- a/sympy/core/tests/test_containers.py
+++ b/sympy/core/tests/test_containers.py
@@ -1,14 +1,8 @@
 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.abc import x, y
+from sympy.core.compatibility import is_sequence, iterable, range
 
 
 def test_Tuple():
@@ -69,21 +63,6 @@ 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 4b81ca1b4b..62a70dae85 100644
--- a/sympy/core/tests/test_relational.py
+++ b/sympy/core/tests/test_relational.py
@@ -129,9 +129,6 @@ 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/functions/special/tests/test_beta_functions.py b/sympy/functions/special/tests/test_beta_functions.py
index 3c23ac4e10..fc19f1ed35 100644
--- a/sympy/functions/special/tests/test_beta_functions.py
+++ b/sympy/functions/special/tests/test_beta_functions.py
@@ -1,4 +1,5 @@
 from sympy import (Symbol, gamma, expand_func, beta, digamma, diff, conjugate)
+from sympy.functions.special.gamma_functions import polygamma
 from sympy.core.function import ArgumentIndexError
 from sympy.utilities.pytest import raises
 
@@ -12,8 +13,8 @@ def test_beta():
     assert expand_func(beta(x, y) - beta(y, x)) == 0  # Symmetric
     assert expand_func(beta(x, y)) == expand_func(beta(x, y + 1) + beta(x + 1, y)).simplify()
 
-    assert diff(beta(x, y), x) == beta(x, y)*(digamma(x) - digamma(x + y))
-    assert diff(beta(x, y), y) == beta(x, y)*(digamma(y) - digamma(x + y))
+    assert diff(beta(x, y), x) == beta(x, y)*(polygamma(0, x) - polygamma(0, x + y))
+    assert diff(beta(x, y), y) == beta(x, y)*(polygamma(0, y) - polygamma(0, x + y))
 
     assert conjugate(beta(x, y)) == beta(conjugate(x), conjugate(y))
 
diff --git a/sympy/functions/special/tests/test_gamma_functions.py b/sympy/functions/special/tests/test_gamma_functions.py
index 5acabb1442..5cc70ea927 100644
--- a/sympy/functions/special/tests/test_gamma_functions.py
+++ b/sympy/functions/special/tests/test_gamma_functions.py
@@ -1,6 +1,6 @@
 from sympy import (
     Symbol, Dummy, gamma, I, oo, nan, zoo, factorial, sqrt, Rational,
-    multigamma, log, polygamma, EulerGamma, pi, uppergamma, S, expand_func,
+    multigamma, log, polygamma, digamma, trigamma, EulerGamma, pi, uppergamma, S, expand_func,
     loggamma, sin, cos, O, lowergamma, exp, erf, erfc, exp_polar, harmonic,
     zeta, conjugate, Ei, im, re, tanh, Abs)
 
@@ -114,7 +114,7 @@ def test_lowergamma():
     assert td(lowergamma(randcplx(), y), y)
     assert td(lowergamma(x, randcplx()), x)
     assert lowergamma(x, y).diff(x) == \
-        gamma(x)*polygamma(0, x) - uppergamma(x, y)*log(y) \
+        gamma(x)*digamma(x) - uppergamma(x, y)*log(y) \
         - meijerg([], [1, 1], [0, 0, x], [], y)
 
     assert lowergamma(S.Half, x) == sqrt(pi)*erf(sqrt(x))
@@ -367,6 +367,119 @@ def test_polygamma_expand_func():
     e = polygamma(3, x + y + S(3)/4)
     assert e.expand(func=True, basic=False) == e
 
+def test_digamma():
+    from sympy import I
+
+    assert digamma(nan) == nan
+
+    assert digamma(oo) == oo
+    assert digamma(-oo) == oo
+    assert digamma(I*oo) == oo
+    assert digamma(-I*oo) == oo
+
+    assert digamma(-9) == zoo
+
+    assert digamma(-9) == zoo
+    assert digamma(-1) == zoo
+
+    assert digamma(0) == zoo
+
+    assert digamma(1) == -EulerGamma
+    assert digamma(7) == Rational(49, 20) - EulerGamma
+
+    def t(m, n):
+        x = S(m)/n
+        r = digamma(x)
+        if r.has(digamma):
+            return False
+        return abs(digamma(x.n()).n() - r.n()).n() < 1e-10
+    assert t(1, 2)
+    assert t(3, 2)
+    assert t(-1, 2)
+    assert t(1, 4)
+    assert t(-3, 4)
+    assert t(1, 3)
+    assert t(4, 3)
+    assert t(3, 4)
+    assert t(2, 3)
+    assert t(123, 5)
+
+    assert digamma(x).rewrite(zeta) == polygamma(0, x)
+
+    assert digamma(x).rewrite(harmonic) == harmonic(x - 1) - EulerGamma
+
+    assert digamma(I).is_real is None
+
+    assert digamma(x,evaluate=False).fdiff() == polygamma(1, x)
+
+    assert digamma(x,evaluate=False).is_real is None
+
+    assert digamma(x,evaluate=False).is_positive is None
+
+    assert digamma(x,evaluate=False).is_negative is None
+
+    assert digamma(x,evaluate=False).rewrite(polygamma) == polygamma(0, x)
+
+
+def test_digamma_expand_func():
+    assert digamma(x).expand(func=True) == polygamma(0, x)
+    assert digamma(2*x).expand(func=True) == \
+        polygamma(0, x)/2 + polygamma(0, Rational(1, 2) + x)/2 + log(2)
+    assert digamma(-1 + x).expand(func=True) == \
+        polygamma(0, x) - 1/(x - 1)
+    assert digamma(1 + x).expand(func=True) == \
+        1/x + polygamma(0, x )
+    assert digamma(2 + x).expand(func=True) == \
+        1/x + 1/(1 + x) + polygamma(0, x)
+    assert digamma(3 + x).expand(func=True) == \
+        polygamma(0, x) + 1/x + 1/(1 + x) + 1/(2 + x)
+    assert digamma(4 + x).expand(func=True) == \
+        polygamma(0, x) + 1/x + 1/(1 + x) + 1/(2 + x) + 1/(3 + x)
+    assert digamma(x + y).expand(func=True) == \
+        polygamma(0, x + y)
+
+def test_trigamma():
+    from sympy import I
+
+    assert trigamma(nan) == nan
+
+    assert trigamma(oo) == 0
+
+    assert trigamma(1) == pi**2/6
+    assert trigamma(2) == pi**2/6 - 1
+    assert trigamma(3) == pi**2/6 - Rational(5, 4)
+
+    assert trigamma(x, evaluate=False).rewrite(zeta) == zeta(2, x)
+    assert trigamma(x, evaluate=False).rewrite(harmonic) == \
+        trigamma(x).rewrite(polygamma).rewrite(harmonic)
+
+    assert trigamma(x,evaluate=False).fdiff() == polygamma(2, x)
+
+    assert trigamma(x,evaluate=False).is_real is None
+
+    assert trigamma(x,evaluate=False).is_positive is None
+
+    assert trigamma(x,evaluate=False).is_negative is None
+
+    assert trigamma(x,evaluate=False).rewrite(polygamma) == polygamma(1, x)
+
+def test_trigamma_expand_func():
+    assert trigamma(2*x).expand(func=True) == \
+        polygamma(1, x)/4 + polygamma(1, Rational(1, 2) + x)/4
+    assert trigamma(1 + x).expand(func=True) == \
+        polygamma(1, x) - 1/x**2
+    assert trigamma(2 + x).expand(func=True, multinomial=False) == \
+        polygamma(1, x) - 1/x**2 - 1/(1 + x)**2
+    assert trigamma(3 + x).expand(func=True, multinomial=False) == \
+        polygamma(1, x) - 1/x**2 - 1/(1 + x)**2 - 1/(2 + x)**2
+    assert trigamma(4 + x).expand(func=True, multinomial=False) == \
+        polygamma(1, x) - 1/x**2 - 1/(1 + x)**2 - \
+        1/(2 + x)**2 - 1/(3 + x)**2
+    assert trigamma(x + y).expand(func=True) == \
+        polygamma(1, x + y)
+    assert trigamma(3 + 4*x + y).expand(func=True, multinomial=False) == \
+        polygamma(1, y + 4*x) - 1/(y + 4*x)**2 - \
+        1/(1 + y + 4*x)**2 - 1/(2 + y + 4*x)**2
 
 def test_loggamma():
     raises(TypeError, lambda: loggamma(2, 3))
diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index dcb49b94ca..479b1beda0 100644
--- a/sympy/printing/pretty/tests/test_pretty.py
+++ b/sympy/printing/pretty/tests/test_pretty.py
@@ -3858,13 +3858,7 @@ def test_pretty_Intersection_issue_10414():
     assert upretty(Intersection(a, b)) == ucode_str
     assert pretty(Intersection(a, b)) == ascii_str
 
-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():
+def test_ProductSet_paranthesis():
     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 fbff7ae2e0..402fd7b1fd 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).flatten()) == r"%s \times %s \times %s" % (
+    assert latex(line * bigline * fset) == 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).flatten()) == \
+    assert latex(ProductSet(A, P2)) == \
         '\\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 f3159d25c3..f049ce7c7d 100644
--- a/sympy/sets/tests/test_fancysets.py
+++ b/sympy/sets/tests/test_fancysets.py
@@ -121,11 +121,6 @@ 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 559cd41df2..278e987944 100644
--- a/sympy/sets/tests/test_sets.py
+++ b/sympy/sets/tests/test_sets.py
@@ -261,15 +261,6 @@ 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) == \
@@ -434,85 +425,8 @@ def test_is_disjoint():
     assert Interval(0, 2).is_disjoint(Interval(3, 4)) == True
 
 
-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_of_single_arg_is_arg():
+    assert ProductSet(Interval(0, 1)) == Interval(0, 1)
 
 
 def test_ProductSet_is_empty():
@@ -885,10 +799,8 @@ 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).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
+    assert (.5, .5, .5) in square * unit_line
+    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)))
 
@@ -1128,20 +1040,10 @@ 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 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))
+    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
 
 
 def test_SymmetricDifference():