Repository: sympy/sympy

Issue:
where possible, use python fractions for Rational
Since we have dropped support for 2.5 we can use the [fractions module](https://docs.python.org/2/library/fractions.html?highlight=rational) for Rational internals. The `gcd` method could also replace `igcd`:

```python
>>> timeit('gcd(12,3)','from fractions import gcd')
0.5296611342553904
>>> timeit('gcd(12,3)','from sympy import igcd as gcd')
6.766568780548827
```

Discussion:
smichr: ping @isuruf -- please take a look at the mods here in contrast to #11089
isuruf: > I wonder if it might be better to use Fraction as the underlying representation and then update methods like Rational.p and Rational.args to return the correct args of the Fraction

Implementation of Fractions doesn't do the `gcd` hint trick introduced in #11019, therefore will be slower to use. Also, arithmetic operators would use `gcd` instead of `igcd` and therefore caching will not be used in that case.
smichr: ping @isuruf - ok, have a look now. I benchmark with

```
for I in range(10**4):
 for J in range(10**4):
   g = igcd(I,j)
```

I think it is now more similar to the original `igcd` now. 

Patch:
diff --git a/doc/src/modules/stats.rst b/doc/src/modules/stats.rst
index 7b301a4ba5..d53f8b44d1 100644
--- a/doc/src/modules/stats.rst
+++ b/doc/src/modules/stats.rst
@@ -64,16 +64,11 @@ Interface
 ^^^^^^^^^
 
 .. autofunction:: P
-.. autoclass:: Probability
 .. autofunction:: E
-.. autoclass:: Expectation
 .. autofunction:: density
 .. autofunction:: given
 .. autofunction:: where
 .. autofunction:: variance
-.. autoclass:: Variance
-.. autofunction:: covariance
-.. autoclass:: Covariance
 .. autofunction:: std
 .. autofunction:: sample
 .. autofunction:: sample_iter
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index aeff883ef7..dde9bb8713 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -140,7 +140,7 @@ def _literal_float(f):
 
 
 def igcd(*args):
-    """Computes positive integer greatest common divisor.
+    """Computes nonnegative integer greatest common divisor.
 
     The algorithm is based on the well known Euclid's algorithm. To
     improve speed, igcd() has its own caching mechanism implemented.
@@ -155,25 +155,37 @@ def igcd(*args):
     5
 
     """
-    a = args[0]
-    for b in args[1:]:
-        try:
-            a = _gcdcache[(a, b)]
-        except KeyError:
-            a, b = as_int(a), as_int(b)
-
-            if a and b:
+    if len(args) < 2:
+        raise TypeError(
+            'igcd() takes at least 2 arguments (%s given)' % len(args))
+    if 1 in args:
+        a = 1
+        k = 0
+    else:
+        a = abs(as_int(args[0]))
+        k = 1
+    if a != 1:
+        while k < len(args):
+            b = args[k]
+            k += 1
+            try:
+                a = _gcdcache[(a, b)]
+            except KeyError:
+                b = as_int(b)
+                if not b:
+                    continue
+                if b == 1:
+                    a = 1
+                    break
                 if b < 0:
                     b = -b
-
+                t = a, b
                 while b:
                     a, b = b, a % b
-            else:
-                a = abs(a or b)
-
-            _gcdcache[(a, b)] = a
-        if a == 1 or b == 1:
-            return 1
+                _gcdcache[t] = _gcdcache[t[1], t[0]] = a
+    while k < len(args):
+        ok = as_int(args[k])
+        k += 1
     return a
 
 
@@ -192,6 +204,9 @@ def ilcm(*args):
     30
 
     """
+    if len(args) < 2:
+        raise TypeError(
+            'ilcm() takes at least 2 arguments (%s given)' % len(args))
     if 0 in args:
         return 0
     a = args[0]
@@ -1233,37 +1248,36 @@ def __new__(cls, p, q=None, gcd=None):
                 return p
 
             if isinstance(p, string_types):
+                if p.count('/') > 1:
+                    raise TypeError('invalid input: %s' % p)
+                pq = p.rsplit('/', 1)
+                if len(pq) == 2:
+                    p, q = pq
+                    fp = fractions.Fraction(p)
+                    fq = fractions.Fraction(q)
+                    f = fp/fq
+                    return Rational(f.numerator, f.denominator, 1)
                 p = p.replace(' ', '')
                 try:
-                    # we might have a Float
-                    neg_pow, digits, expt = decimal.Decimal(p).as_tuple()
-                    p = [1, -1][neg_pow]*int("".join(str(x) for x in digits))
-                    if expt > 0:
-                        return Rational(p*Pow(10, expt), 1)
-                    return Rational(p, Pow(10, -expt))
-                except decimal.InvalidOperation:
-                    f = regex.match('^([-+]?[0-9]+)/([0-9]+)$', p)
-                    if f:
-                        n, d = f.groups()
-                        return Rational(int(n), int(d))
-                    elif p.count('/') == 1:
-                        p, q = p.split('/')
-                        return Rational(Rational(p), Rational(q))
-                    else:
-                        pass  # error will raise below
-            else:
+                    p = fractions.Fraction(p)
+                except ValueError:
+                    pass  # error will raise below
+            elif isinstance(p, float):
+                p = fractions.Fraction(p)
+
+            if not isinstance(p, string_types):
                 try:
                     if isinstance(p, fractions.Fraction):
                         return Rational(p.numerator, p.denominator, 1)
                 except NameError:
                     pass  # error will raise below
 
-                if isinstance(p, (float, Float)):
+                if isinstance(p, Float):
                     return Rational(*float(p).as_integer_ratio())
 
             if not isinstance(p, SYMPY_INTS + (Rational,)):
                 raise TypeError('invalid input: %s' % p)
-            q = S.One
+            q = q or S.One
             gcd = 1
         else:
             p = Rational(p)
@@ -1311,49 +1325,8 @@ def limit_denominator(self, max_denominator=1000000):
         311/99
 
         """
-        # Algorithm notes: For any real number x, define a *best upper
-        # approximation* to x to be a rational number p/q such that:
-        #
-        #   (1) p/q >= x, and
-        #   (2) if p/q > r/s >= x then s > q, for any rational r/s.
-        #
-        # Define *best lower approximation* similarly.  Then it can be
-        # proved that a rational number is a best upper or lower
-        # approximation to x if, and only if, it is a convergent or
-        # semiconvergent of the (unique shortest) continued fraction
-        # associated to x.
-        #
-        # To find a best rational approximation with denominator <= M,
-        # we find the best upper and lower approximations with
-        # denominator <= M and take whichever of these is closer to x.
-        # In the event of a tie, the bound with smaller denominator is
-        # chosen.  If both denominators are equal (which can happen
-        # only when max_denominator == 1 and self is midway between
-        # two integers) the lower bound---i.e., the floor of self, is
-        # taken.
-
-        if max_denominator < 1:
-            raise ValueError("max_denominator should be at least 1")
-        if self.q <= max_denominator:
-            return self
-
-        p0, q0, p1, q1 = 0, 1, 1, 0
-        n, d = self.p, self.q
-        while True:
-            a = n//d
-            q2 = q0 + a*q1
-            if q2 > max_denominator:
-                break
-            p0, q0, p1, q1 = p1, q1, p0 + a*p1, q2
-            n, d = d, n - a*d
-
-        k = (max_denominator - q0)//q1
-        bound1 = Rational(p0 + k*p1, q0 + k*q1)
-        bound2 = Rational(p1, q1)
-        if abs(bound2 - self) <= abs(bound1 - self):
-            return bound2
-        else:
-            return bound1
+        f = fractions.Fraction(self.p, self.q)
+        return Rational(f.limit_denominator(fractions.Fraction(int(max_denominator))))
 
     def __getnewargs__(self):
         return (self.p, self.q)
diff --git a/sympy/core/power.py b/sympy/core/power.py
index ac07ca460f..91835e8679 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -12,19 +12,9 @@
 from .logic import fuzzy_bool, fuzzy_not
 from .compatibility import as_int, range
 from .evaluate import global_evaluate
-from sympy.utilities.iterables import sift
 
 from mpmath.libmp import sqrtrem as mpmath_sqrtrem
-
-from math import sqrt as _sqrt
-
-
-
-def isqrt(n):
-    """Return the largest integer less than or equal to sqrt(n)."""
-    if n < 17984395633462800708566937239552:
-        return int(_sqrt(n))
-    return integer_nthroot(int(n), 2)[0]
+from sympy.utilities.iterables import sift
 
 
 def integer_nthroot(y, n):
@@ -96,7 +86,7 @@ def integer_nthroot(y, n):
     while t > y:
         x -= 1
         t = x**n
-    return int(x), t == y  # int converts long to int if possible
+    return x, t == y
 
 
 class Pow(Expr):
diff --git a/sympy/functions/elementary/miscellaneous.py b/sympy/functions/elementary/miscellaneous.py
index 82cdb8b7bc..43dc9affc1 100644
--- a/sympy/functions/elementary/miscellaneous.py
+++ b/sympy/functions/elementary/miscellaneous.py
@@ -454,7 +454,7 @@ def _eval_derivative(self, s):
         return Add(*l)
 
     def evalf(self, prec=None, **options):
-        return self.func(*[a.evalf(prec, **options) for a in self.args])
+        return self.func(*[a.evalf(prec, options) for a in self.args])
     n = evalf
 
     _eval_is_algebraic = lambda s: _torf(i.is_algebraic for i in s.args)
diff --git a/sympy/integrals/deltafunctions.py b/sympy/integrals/deltafunctions.py
index fc53923b31..098179ae74 100644
--- a/sympy/integrals/deltafunctions.py
+++ b/sympy/integrals/deltafunctions.py
@@ -3,7 +3,6 @@
 from sympy.core import Mul
 from sympy.functions import DiracDelta, Heaviside
 from sympy.core.compatibility import default_sort_key
-from sympy.core.singleton import S
 
 
 def change_mul(node, x):
@@ -41,6 +40,9 @@ def change_mul(node, x):
        sympy.functions.special.delta_functions.DiracDelta
        deltaintegrate
     """
+    if not (node.is_Mul or node.is_Pow):
+        return node
+
     new_args = []
     dirac = None
 
@@ -54,7 +56,8 @@ def change_mul(node, x):
         if arg.is_Pow and arg.base.func is DiracDelta:
             new_args.append(arg.func(arg.base, arg.exp - 1))
             arg = arg.base
-        if dirac is None and (arg.func is DiracDelta and arg.is_simple(x)):
+        if dirac is None and (arg.func is DiracDelta and arg.is_simple(x)
+                and (len(arg.args) <= 1 or arg.args[1] == 0)):
             dirac = arg
         else:
             new_args.append(arg)
@@ -66,7 +69,7 @@ def change_mul(node, x):
             elif arg.is_Pow and arg.base.func is DiracDelta:
                 new_args.append(arg.func(arg.base.simplify(x), arg.exp))
             else:
-                new_args.append(arg)
+                new_args.append(change_mul(arg, x))
         if new_args != sorted_args:
             nnode = Mul(*new_args).expand()
         else:  # if the node didn't change there is nothing to do
@@ -157,41 +160,17 @@ def deltaintegrate(f, x):
                 return fh
         else:
             # no expansion performed, try to extract a simple DiracDelta term
-            deltaterm, rest_mult = change_mul(f, x)
+            dg, rest_mult = change_mul(f, x)
 
-            if not deltaterm:
+            if not dg:
                 if rest_mult:
                     fh = integrate(rest_mult, x)
                     return fh
             else:
-                deltaterm = deltaterm.simplify(x)
-                if deltaterm.is_Mul:  # Take out any extracted factors
-                    deltaterm, rest_mult_2 = change_mul(deltaterm, x)
+                dg = dg.simplify(x)
+                if dg.is_Mul:  # Take out any extracted factors
+                    dg, rest_mult_2 = change_mul(dg, x)
                     rest_mult = rest_mult*rest_mult_2
-                point = solve(deltaterm.args[0], x)[0]
-
-                # Return the largest hyperreal term left after
-                # repeated integration by parts.  For example,
-                #
-                #   integrate(y*DiracDelta(x, 1),x) == y*DiracDelta(x,0),  not 0
-                #
-                # This is so Integral(y*DiracDelta(x).diff(x),x).doit()
-                # will return y*DiracDelta(x) instead of 0 or DiracDelta(x),
-                # both of which are correct everywhere the value is defined
-                # but give wrong answers for nested integration.
-                n = (0 if len(deltaterm.args)==1 else deltaterm.args[1])
-                m = 0
-                while n >= 0:
-                    r = (-1)**n*rest_mult.diff(x, n).subs(x, point)
-                    if r is S.Zero:
-                        n -= 1
-                        m += 1
-                    else:
-                        if m == 0:
-                            return r*Heaviside(x - point)
-                        else:
-                            return r*DiracDelta(x,m-1)
-                # In some very weak sense, x=0 is still a singularity,
-                # but we hope will not be of any practial consequence.
-                return S.Zero
+                point = solve(dg.args[0], x)[0]
+                return (rest_mult.subs(x, point)*Heaviside(x - point))
     return None
diff --git a/sympy/matrices/densesolve.py b/sympy/matrices/densesolve.py
index bfbfc31660..9bac3f5587 100644
--- a/sympy/matrices/densesolve.py
+++ b/sympy/matrices/densesolve.py
@@ -7,14 +7,10 @@
 
 from sympy.matrices.densetools import col, eye, augment
 from sympy.matrices.densetools import rowadd, rowmul, conjugate_transpose
-from sympy.core.symbol import symbols
+from sympy import sqrt, var
 from sympy.core.compatibility import range
-from sympy.core.power import isqrt
-
 import copy
 
-
-
 def row_echelon(matlist, K):
     """
     Returns the row echelon form of a matrix with diagonal elements
@@ -139,7 +135,7 @@ def cholesky(matlist, K):
             for k in range(j):
                 a += L[i][k]*L[j][k]
             if i == j:
-                L[i][j] = isqrt(new_matlist[i][j] - a)
+                L[i][j] = int(sqrt(new_matlist[i][j] - a))
             else:
                 L[i][j] = (new_matlist[i][j] - a)/L[j][j]
     return L, conjugate_transpose(L, K)
@@ -308,8 +304,10 @@ def LU_solve(matlist, variable, constant, K):
     """
     new_matlist = copy.deepcopy(matlist)
     nrow = len(new_matlist)
+    y = []
     L, U = LU(new_matlist, K)
-    y = [[i] for i in symbols('y:%i' % nrow)]
+    for i in range(nrow):
+        y.append([var('y' + str(i))])
     forward_substitution(L, y, constant, K)
     backward_substitution(U, variable, y, K)
     return variable
@@ -351,8 +349,10 @@ def cholesky_solve(matlist, variable, constant, K):
     """
     new_matlist = copy.deepcopy(matlist)
     nrow = len(new_matlist)
+    y = []
     L, Lstar = cholesky(new_matlist, K)
-    y = [[i] for i in symbols('y:%i' % nrow)]
+    for i in range(nrow):
+        y.append([var('y' + str(i))])
     forward_substitution(L, y, constant, K)
     backward_substitution(Lstar, variable, y, K)
     return variable
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 7c8f2b03a3..706c578d74 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -338,7 +338,6 @@ class MatrixElement(Expr):
     i = property(lambda self: self.args[1])
     j = property(lambda self: self.args[2])
     _diff_wrt = True
-    is_commutative = True
 
     def doit(self, **kwargs):
         deep = kwargs.get('deep', True)
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 41314fc251..a44a70f567 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -505,11 +505,7 @@ def __mul__(self, other):
                 raise ShapeError("Matrices size mismatch.")
             if A.cols == 0:
                 return classof(A, B)._new(A.rows, B.cols, lambda i, j: 0)
-            try:
-                blst = B.T.tolist()
-            except AttributeError:
-                # If B is a MatrixSymbol, B.T.tolist does not exist
-                return NotImplemented
+            blst = B.T.tolist()
             alst = A.tolist()
             return classof(A, B)._new(A.rows, B.cols, lambda i, j:
                 reduce(lambda k, l: k + l,
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index cff676835b..8a50b13f87 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -911,8 +911,6 @@ def _union(self, other):
 
         See Set._union for docstring
         """
-        if other.is_UniversalSet:
-            return S.UniversalSet
         if other.is_Interval and self._is_comparable(other):
             from sympy.functions.elementary.miscellaneous import Min, Max
             # Non-overlapping intervals
diff --git a/sympy/simplify/trigsimp.py b/sympy/simplify/trigsimp.py
index 69f03cc98e..81f179f685 100644
--- a/sympy/simplify/trigsimp.py
+++ b/sympy/simplify/trigsimp.py
@@ -15,7 +15,7 @@
 from sympy.strategies.core import identity
 from sympy.strategies.tree import greedy
 
-from sympy.polys import Poly
+
 from sympy.polys.polyerrors import PolificationFailed
 from sympy.polys.polytools import groebner
 from sympy.polys.domains import ZZ
diff --git a/sympy/solvers/diophantine.py b/sympy/solvers/diophantine.py
index 2b30ad0f88..a05e2e50f4 100644
--- a/sympy/solvers/diophantine.py
+++ b/sympy/solvers/diophantine.py
@@ -7,7 +7,7 @@
 from sympy.core.mul import Mul
 from sympy.core.numbers import Rational
 from sympy.core.numbers import igcdex, ilcm, igcd
-from sympy.core.power import integer_nthroot, isqrt
+from sympy.core.power import integer_nthroot
 from sympy.core.relational import Eq
 from sympy.core.singleton import S
 from sympy.core.symbol import Symbol, symbols
@@ -29,7 +29,6 @@
 from sympy.utilities.misc import filldedent
 
 
-
 # these are imported with 'from sympy.solvers.diophantine import *
 __all__ = ['diophantine', 'classify_diop']
 
@@ -856,27 +855,26 @@ def _diop_quadratic(var, coeff, t):
             c = C // g
             e = sign(B/A)
 
-            sqa = isqrt(a)
-            sqc = isqrt(c)
-            _c = e*sqc*D - sqa*E
-            if not _c:
+
+            _c = e*sqrt(c)*D - sqrt(a)*E
+            if _c == 0:
                 z = symbols("z", real=True)
-                eq = sqa*g*z**2 + D*z + sqa*F
+                eq = sqrt(a)*g*z**2 + D*z + sqrt(a)*F
                 roots = solveset_real(eq, z).intersect(S.Integers)
                 for root in roots:
-                    ans = diop_solve(sqa*x + e*sqc*y - root)
+                    ans = diop_solve(sqrt(a)*x + e*sqrt(c)*y - root)
                     sol.add((ans[0], ans[1]))
 
-            elif _is_int(c):
-                solve_x = lambda u: -e*sqc*g*_c*t**2 - (E + 2*e*sqc*g*u)*t\
-                    - (e*sqc*g*u**2 + E*u + e*sqc*F) // _c
+            elif _is_int(_c):
+                solve_x = lambda u: -e*sqrt(c)*g*_c*t**2 - (E + 2*e*sqrt(c)*g*u)*t\
+                    - (e*sqrt(c)*g*u**2 + E*u + e*sqrt(c)*F) // _c
 
-                solve_y = lambda u: sqa*g*_c*t**2 + (D + 2*sqa*g*u)*t \
-                    + (sqa*g*u**2 + D*u + sqa*F) // _c
+                solve_y = lambda u: sqrt(a)*g*_c*t**2 + (D + 2*sqrt(a)*g*u)*t \
+                    + (sqrt(a)*g*u**2 + D*u + sqrt(a)*F) // _c
 
                 for z0 in range(0, abs(_c)):
                     if divisible(
-                            sqa*g*z0**2 + D*z0 + sqa*F,
+                            sqrt(a)*g*z0**2 + D*z0 + sqrt(a)*F,
                             _c):
                         sol.add((solve_x(z0), solve_y(z0)))
 
@@ -2548,7 +2546,7 @@ def diop_general_sum_of_squares(eq, limit=1):
     >>> from sympy.solvers.diophantine import diop_general_sum_of_squares
     >>> from sympy.abc import a, b, c, d, e, f
     >>> diop_general_sum_of_squares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345)
-    set([(15, 22, 22, 24, 24)])
+    set([(2, 4, 4, 10, 47)])
 
     Reference
     =========
@@ -2810,23 +2808,26 @@ def sum_of_three_squares(n):
     if n % 8 == 3:
         s = s if _odd(s) else s - 1
 
-        for x in range(s, -1, -2):
-            N = (n - x**2) // 2
-            if isprime(N):
-                y, z = prime_as_sum_of_two_squares(N)
-                return _sorted_tuple(2**v*x, 2**v*(y + z), 2**v*abs(y - z))
-        return
+        for i in range(s, -1, -2):
+            if isprime((n - i**2) // 2):
+                x = i
+                break
+
+        y, z = prime_as_sum_of_two_squares((n - x**2) // 2)
+        return _sorted_tuple(2**v*x, 2**v*(y + z), 2**v*abs(y - z))
 
     if n % 8 == 2 or n % 8 == 6:
         s = s if _odd(s) else s - 1
     else:
         s = s - 1 if _odd(s) else s
 
-    for x in range(s, -1, -2):
-        N = n - x**2
-        if isprime(N):
-            y, z = prime_as_sum_of_two_squares(N)
-            return _sorted_tuple(2**v*x, 2**v*y, 2**v*z)
+    for i in range(s, -1, -2):
+        if isprime(n - i**2):
+            x = i
+            break
+
+    y, z = prime_as_sum_of_two_squares(n - x**2)
+    return _sorted_tuple(2**v*x, 2**v*y, 2**v*z)
 
 
 def sum_of_four_squares(n):
@@ -2901,16 +2902,16 @@ def power_representation(n, p, k, zeros=False):
 
     >>> f = power_representation(1729, 3, 2)
     >>> next(f)
-    (9, 10)
-    >>> next(f)
     (1, 12)
+    >>> next(f)
+    (9, 10)
 
     If the flag `zeros` is True, the solution may contain tuples with
     zeros; any such solutions will be generated after the solutions
     without zeros:
 
     >>> list(power_representation(125, 2, 3, zeros=True))
-    [(5, 6, 8), (3, 4, 10), (0, 5, 10), (0, 2, 11)]
+    [(3, 4, 10), (5, 6, 8), (0, 2, 11), (0, 5, 10)]
 
     For even `p` the `permute_sign` function can be used to get all
     signed values:
@@ -2981,13 +2982,13 @@ def power_representation(n, p, k, zeros=False):
     if n >= k:
         a = integer_nthroot(n - (k - 1), p)[0]
         for t in pow_rep_recursive(a, k, n, [], p):
-            yield tuple(reversed(t))
+            yield t
 
     if zeros:
         a = integer_nthroot(n, p)[0]
         for i in range(1, k):
             for t in pow_rep_recursive(a, i, n, [], p):
-                yield tuple(reversed(t + (0,) * (k - i)))
+                yield (0,) * (k - i) + t
 
 
 sum_of_powers = power_representation
@@ -2996,15 +2997,15 @@ def power_representation(n, p, k, zeros=False):
 def pow_rep_recursive(n_i, k, n_remaining, terms, p):
 
     if k == 0 and n_remaining == 0:
-        yield tuple(terms)
+        yield _sorted_tuple(*terms)
     else:
-        if n_i >= 1 and k > 0:
+        if n_i >= 1 and k > 0 and n_remaining >= 0:
+            if n_i**p <= n_remaining:
+                for t in pow_rep_recursive(n_i, k - 1, n_remaining - n_i**p, terms + [n_i], p):
+                    yield _sorted_tuple(*t)
+
             for t in pow_rep_recursive(n_i - 1, k, n_remaining, terms, p):
-                yield t
-            residual = n_remaining - pow(n_i, p)
-            if residual >= 0:
-                for t in pow_rep_recursive(n_i, k - 1, residual, terms + [n_i], p):
-                    yield t
+                yield _sorted_tuple(*t)
 
 
 def sum_of_squares(n, k, zeros=False):
diff --git a/sympy/stats/__init__.py b/sympy/stats/__init__.py
index f426633dd0..d34f8f93a1 100644
--- a/sympy/stats/__init__.py
+++ b/sympy/stats/__init__.py
@@ -68,7 +68,3 @@
 from . import drv_types
 from .drv_types import (Geometric, Poisson)
 __all__.extend(drv_types.__all__)
-
-from . import symbolic_probability
-from .symbolic_probability import Probability, Expectation, Variance, Covariance
-__all__.extend(symbolic_probability.__all__)
diff --git a/sympy/stats/symbolic_probability.py b/sympy/stats/symbolic_probability.py
index 4f8de12702..8ac27c4696 100644
--- a/sympy/stats/symbolic_probability.py
+++ b/sympy/stats/symbolic_probability.py
@@ -9,24 +9,10 @@
 from sympy.stats import variance, covariance
 from sympy.stats.rv import RandomSymbol, probability, expectation
 
-__all__ = ['Probability', 'Expectation', 'Variance', 'Covariance']
-
 
 class Probability(Expr):
     """
     Symbolic expression for the probability.
-
-    Examples
-    ========
-
-    >>> from sympy.stats import Probability, Normal
-    >>> X = Normal("X", 0, 1)
-    >>> prob = Probability(X > 1)
-    >>> prob
-    Probability(X > 1)
-    >>> prob.doit()
-    sqrt(2)*(-sqrt(2)*sqrt(pi)*erf(sqrt(2)/2) + sqrt(2)*sqrt(pi))/(4*sqrt(pi))
-
     """
     def __new__(cls, prob, condition=None, **kwargs):
         prob = _sympify(prob)
@@ -45,37 +31,6 @@ def doit(self, **kwargs):
 class Expectation(Expr):
     """
     Symbolic expression for the expectation.
-
-    Examples
-    ========
-
-    >>> from sympy.stats import Expectation, Normal
-    >>> from sympy import symbols
-    >>> mu = symbols("mu")
-    >>> sigma = symbols("sigma", positive=True)
-    >>> X = Normal("X", mu, sigma)
-    >>> Expectation(X)
-    Expectation(X)
-    >>> Expectation(X).doit().simplify()
-    mu
-    >>> Expectation(X).doit(evaluate=False)
-    Integral(sqrt(2)*X*exp(-(X - mu)**2/(2*sigma**2))/(2*sqrt(pi)*sigma), (X, -oo, oo))
-
-    This class is aware of some properties of the expectation:
-
-    >>> from sympy.abc import a
-    >>> Expectation(a*X)
-    a*Expectation(X)
-    >>> Y = Normal("Y", 0, 1)
-    >>> Expectation(X + Y)
-    Expectation(X) + Expectation(Y)
-
-    To prevent this kind of automatic evaluation, use ``evaluate=False``:
-
-    >>> Expectation(a*X + Y)
-    a*Expectation(X) + Expectation(Y)
-    >>> Expectation(a*X + Y, evaluate=False)
-    Expectation(a*X + Y)
     """
     def __new__(cls, expr, condition=None, **kwargs):
         expr = _sympify(expr)
@@ -120,39 +75,6 @@ def doit(self, **kwargs):
 class Variance(Expr):
     """
     Symbolic expression for the variance.
-
-    Examples
-    ========
-
-    >>> from sympy import symbols
-    >>> from sympy.stats import Normal, Expectation, Variance
-    >>> mu = symbols("mu", positive=True)
-    >>> sigma = symbols("sigma", positive=True)
-    >>> X = Normal("X", mu, sigma)
-    >>> Variance(X)
-    Variance(X)
-    >>> Variance(X).doit()
-    sigma**2
-
-    Rewrite the variance in terms of the expectation
-
-    >>> Variance(X).rewrite(Expectation)
-    -Expectation(X)**2 + Expectation(X**2)
-
-    Some transformations based on the properties of the variance may happen:
-
-    >>> from sympy.abc import a
-    >>> Y = Normal("Y", 0, 1)
-    >>> Variance(a*X)
-    a**2*Variance(X)
-
-    To prevent automatic evaluation, use ``evaluate=False``:
-
-    >>> Variance(X + Y)
-    2*Covariance(X, Y) + Variance(X) + Variance(Y)
-    >>> Variance(X + Y, evaluate=False)
-    Variance(X + Y)
-
     """
     def __new__(cls, arg, condition=None, **kwargs):
         arg = _sympify(arg)
@@ -208,11 +130,6 @@ def __new__(cls, arg, condition=None, **kwargs):
             # this expression contains a RandomSymbol somehow:
             return Variance(arg, condition, evaluate=False)
 
-    def _eval_rewrite_as_Expectation(self, arg, condition=None):
-            e1 = Expectation(arg**2, condition)
-            e2 = Expectation(arg, condition)**2
-            return e1 - e2
-
     def doit(self, **kwargs):
         return variance(self.args[0], self._condition, **kwargs)
 
@@ -220,57 +137,23 @@ def doit(self, **kwargs):
 class Covariance(Expr):
     """
     Symbolic expression for the covariance.
-
-    Examples
-    ========
-
-    >>> from sympy.stats import Covariance
-    >>> from sympy.stats import Normal
-    >>> X = Normal("X", 3, 2)
-    >>> Y = Normal("Y", 0, 1)
-    >>> Z = Normal("Z", 0, 1)
-    >>> W = Normal("W", 0, 1)
-    >>> cexpr = Covariance(X, Y)
-    >>> cexpr
-    Covariance(X, Y)
-
-    Evaluate the covariance, `X` and `Y` are independent,
-    therefore zero is the result:
-
-    >>> cexpr.doit()
-    0
-
-    Rewrite the covariance expression in terms of expectations:
-
-    >>> from sympy.stats import Expectation
-    >>> cexpr.rewrite(Expectation)
-    Expectation(X*Y) - Expectation(X)*Expectation(Y)
-
-    This class is aware of some properties of the covariance:
-
-    >>> Covariance(X, X)
-    Variance(X)
-    >>> from sympy.abc import a, b, c, d
-    >>> Covariance(a*X, b*Y)
-    a*b*Covariance(X, Y)
-
-    In order to prevent such transformations, use ``evaluate=False``:
-
-    >>> Covariance(a*X + b*Y, c*Z + d*W)
-    a*c*Covariance(X, Z) + a*d*Covariance(W, X) + b*c*Covariance(Y, Z) + b*d*Covariance(W, Y)
-    >>> Covariance(a*X + b*Y, c*Z + d*W, evaluate=False)
-    Covariance(a*X + b*Y, c*Z + d*W)
     """
     def __new__(cls, arg1, arg2, condition=None, **kwargs):
         arg1 = _sympify(arg1)
         arg2 = _sympify(arg2)
 
+        if not kwargs.pop('evaluate', global_evaluate[0]):
+            if condition is None:
+                obj = Expr.__new__(cls, arg1, arg2)
+            else:
+                condition = _sympify(condition)
+                obj = Expr.__new__(cls, arg1, arg2, condition)
+            obj._condition = condition
+            return obj
+
         if condition is not None:
             condition = _sympify(condition)
 
-        if not kwargs.pop('evaluate', global_evaluate[0]):
-            return cls._call_super_constructor(arg1, arg2, condition)
-
         if arg1 == arg2:
             return Variance(arg1, condition)
 
@@ -282,24 +165,15 @@ def __new__(cls, arg1, arg2, condition=None, **kwargs):
         arg1, arg2 = sorted([arg1, arg2], key=default_sort_key)
 
         if isinstance(arg1, RandomSymbol) and isinstance(arg2, RandomSymbol):
-            return cls._call_super_constructor(arg1, arg2, condition)
+            return Expr.__new__(cls, arg1, arg2)
 
         coeff_rv_list1 = cls._expand_single_argument(arg1.expand())
         coeff_rv_list2 = cls._expand_single_argument(arg2.expand())
 
-        addends = [a*b*Covariance(*sorted([r1, r2], key=default_sort_key), condition=condition, evaluate=False)
+        addends = [a*b*Covariance(*sorted([r1, r2], key=default_sort_key), evaluate=False)
                    for (a, r1) in coeff_rv_list1 for (b, r2) in coeff_rv_list2]
         return Add(*addends)
 
-    @classmethod
-    def _call_super_constructor(cls, arg1, arg2, condition):
-        if condition is not None:
-            obj = Expr.__new__(cls, arg1, arg2, condition)
-        else:
-            obj = Expr.__new__(cls, arg1, arg2)
-        obj._condition = condition
-        return obj
-
     @classmethod
     def _expand_single_argument(cls, expr):
         # return (coefficient, random_symbol) pairs:
@@ -330,10 +204,5 @@ def _get_mul_nonrv_rv_tuple(cls, m):
                 nonrv.append(a)
         return (Mul(*nonrv), Mul(*rv))
 
-    def _eval_rewrite_as_Expectation(self, arg1, arg2, condition=None):
-        e1 = Expectation(arg1*arg2, condition)
-        e2 = Expectation(arg1, condition)*Expectation(arg2, condition)
-        return e1 - e2
-
     def doit(self, **kwargs):
         return covariance(self.args[0], self.args[1], self._condition, **kwargs)
diff --git a/sympy/tensor/tensor.py b/sympy/tensor/tensor.py
index ceb19944c3..1f14351c3b 100644
--- a/sympy/tensor/tensor.py
+++ b/sympy/tensor/tensor.py
@@ -828,16 +828,8 @@ def _get(self, key):
             if any([i is None for i in data_list]):
                 raise ValueError("Mixing tensors with associated components "\
                                  "data with tensors without components data")
-
-            free_args_list = [[x[0] for x in arg.free] for arg in key.args]
-            numpy = import_module("numpy")
-            for data, free_args in zip(data_list, free_args_list):
-                if len(free_args) < 2:
-                    sumvar += data
-                else:
-                    free_args_pos = {y: x for x, y in enumerate(free_args)}
-                    axes = [free_args_pos[arg] for arg in key.free_args]
-                    sumvar += numpy.transpose(data, axes)
+            for i in data_list:
+                sumvar += i
             return sumvar
 
         return None


Test patch:
diff --git a/sympy/core/tests/test_numbers.py b/sympy/core/tests/test_numbers.py
index 0df107ec5d..1ebb1c31a9 100644
--- a/sympy/core/tests/test_numbers.py
+++ b/sympy/core/tests/test_numbers.py
@@ -4,14 +4,18 @@
                    Number, zoo, log, Mul, Pow, Tuple, latex, Gt, Lt, Ge, Le,
                    AlgebraicNumber, simplify, sin)
 from sympy.core.compatibility import long, u
-from sympy.core.power import integer_nthroot, isqrt
+from sympy.core.power import integer_nthroot
 from sympy.core.logic import fuzzy_not
 from sympy.core.numbers import (igcd, ilcm, igcdex, seterr, _intcache,
     mpf_norm, comp, mod_inverse)
-from mpmath import mpf
+from sympy.utilities.iterables import permutations
 from sympy.utilities.pytest import XFAIL, raises
+
+from mpmath import mpf
 import mpmath
 
+
+
 t = Symbol('t', real=False)
 
 def same_and_same_prec(a, b):
@@ -203,8 +207,14 @@ def test_igcd():
     assert igcd(-7, 3) == 1
     assert igcd(-7, -3) == 1
     assert igcd(*[10, 20, 30]) == 10
-    raises(ValueError, lambda: igcd(45.1, 30))
-    raises(ValueError, lambda: igcd(45, 30.1))
+    raises(TypeError, lambda: igcd())
+    raises(TypeError, lambda: igcd(2))
+    raises(ValueError, lambda: igcd(0, None))
+    raises(ValueError, lambda: igcd(1, 2.2))
+    for args in permutations((45.1, 1, 30)):
+        raises(ValueError, lambda: igcd(*args))
+    for args in permutations((1, 2, None)):
+        raises(ValueError, lambda: igcd(*args))
 
 
 def test_ilcm():
@@ -855,24 +865,12 @@ def test_powers():
     # Test that this is fast
     assert integer_nthroot(2, 10**10) == (1, False)
 
-    # output should be int if possible
-    assert type(integer_nthroot(2**61, 2)[0]) is int
-
 
 def test_integer_nthroot_overflow():
     assert integer_nthroot(10**(50*50), 50) == (10**50, True)
     assert integer_nthroot(10**100000, 10000) == (10**10, True)
 
 
-def test_isqrt():
-    from math import sqrt as _sqrt
-    limit = 17984395633462800708566937239551
-    assert int(_sqrt(limit)) == integer_nthroot(limit, 2)[0]
-    assert int(_sqrt(limit + 1)) != integer_nthroot(limit + 1, 2)[0]
-    assert isqrt(limit + 1) == integer_nthroot(limit + 1, 2)[0]
-    assert isqrt(limit + 1 + S.Half) == integer_nthroot(limit + 1, 2)[0]
-
-
 def test_powers_Integer():
     """Test Integer._eval_power"""
     # check infinity
diff --git a/sympy/functions/elementary/tests/test_miscellaneous.py b/sympy/functions/elementary/tests/test_miscellaneous.py
index 1fb2ccd184..6e4b301f2c 100644
--- a/sympy/functions/elementary/tests/test_miscellaneous.py
+++ b/sympy/functions/elementary/tests/test_miscellaneous.py
@@ -351,21 +351,3 @@ def test_rewrite_MaxMin_as_Heaviside():
         x*Heaviside(-2*x)*Heaviside(-x - 2) - \
         x*Heaviside(2*x)*Heaviside(x - 2) \
         - 2*Heaviside(-x + 2)*Heaviside(x + 2)
-
-
-def test_issue_11099():
-    from sympy.abc import x, y
-    # some fixed value tests
-    fixed_test_data = {x: -2, y: 3}
-    assert Min(x, y).evalf(subs=fixed_test_data) == \
-        Min(x, y).subs(fixed_test_data).evalf()
-    assert Max(x, y).evalf(subs=fixed_test_data) == \
-        Max(x, y).subs(fixed_test_data).evalf()
-    # randomly generate some test data
-    from random import randint
-    for i in range(20):
-        random_test_data = {x: randint(-100, 100), y: randint(-100, 100)}
-        assert Min(x, y).evalf(subs=random_test_data) == \
-            Min(x, y).subs(random_test_data).evalf()
-        assert Max(x, y).evalf(subs=random_test_data) == \
-            Max(x, y).subs(random_test_data).evalf()
diff --git a/sympy/functions/special/tests/test_delta_functions.py b/sympy/functions/special/tests/test_delta_functions.py
index e120380a3d..dd14bbdf7e 100644
--- a/sympy/functions/special/tests/test_delta_functions.py
+++ b/sympy/functions/special/tests/test_delta_functions.py
@@ -18,10 +18,6 @@ def test_DiracDelta():
     assert DiracDelta(nan) == nan
     assert DiracDelta(0).func is DiracDelta
     assert DiracDelta(x).func is DiracDelta
-    # FIXME: this is generally undefined @ x=0
-    #         But then limit(Delta(c)*Heaviside(x),x,-oo)
-    #         need's to be implemented.
-    #assert 0*DiracDelta(x) == 0
 
     assert adjoint(DiracDelta(x)) == DiracDelta(x)
     assert adjoint(DiracDelta(x - y)) == DiracDelta(x - y)
diff --git a/sympy/integrals/tests/test_deltafunctions.py b/sympy/integrals/tests/test_deltafunctions.py
index e1977280ca..8b3789aca0 100644
--- a/sympy/integrals/tests/test_deltafunctions.py
+++ b/sympy/integrals/tests/test_deltafunctions.py
@@ -7,7 +7,7 @@
 
 
 def test_change_mul():
-    assert change_mul(x, x) == (None, None)
+    assert change_mul(x, x) == x
     assert change_mul(x*y, x) == (None, None)
     assert change_mul(x*y*DiracDelta(x), x) == (DiracDelta(x), x*y)
     assert change_mul(x*y*DiracDelta(x)*DiracDelta(y), x) == \
@@ -35,11 +35,8 @@ def test_deltaintegrate():
     assert deltaintegrate(DiracDelta(x)**2, x) == DiracDelta(0)*Heaviside(x)
     assert deltaintegrate(y*DiracDelta(x)**2, x) == \
         y*DiracDelta(0)*Heaviside(x)
-    assert deltaintegrate(DiracDelta(x, 1), x) == DiracDelta(x, 0)
-    assert deltaintegrate(y*DiracDelta(x, 1), x) == y*DiracDelta(x, 0)
-    assert deltaintegrate(DiracDelta(x, 1)**2, x) == -DiracDelta(0, 2)*Heaviside(x)
-    assert deltaintegrate(y*DiracDelta(x, 1)**2, x) == -y*DiracDelta(0, 2)*Heaviside(x)
-
+    assert deltaintegrate(DiracDelta(x, 1)**2, x) is None
+    assert deltaintegrate(y*DiracDelta(x, 1)**2, x) is None
 
     assert deltaintegrate(DiracDelta(x) * f(x), x) == f(0) * Heaviside(x)
     assert deltaintegrate(DiracDelta(-x) * f(x), x) == f(0) * Heaviside(x)
diff --git a/sympy/matrices/expressions/tests/test_matrix_exprs.py b/sympy/matrices/expressions/tests/test_matrix_exprs.py
index 851df5a622..0dc25a0dc7 100644
--- a/sympy/matrices/expressions/tests/test_matrix_exprs.py
+++ b/sympy/matrices/expressions/tests/test_matrix_exprs.py
@@ -204,24 +204,6 @@ def test_single_indexing():
     B = MatrixSymbol('B', n, m)
     raises(IndexError, lambda: B[1])
 
-def test_MatrixElement_commutative():
-    assert A[0, 1]*A[1, 0] == A[1, 0]*A[0, 1]
-
-def test_MatrixSymbol_determinant():
-    A = MatrixSymbol('A', 4, 4)
-    assert A.as_explicit().det() == A[0, 0]*A[1, 1]*A[2, 2]*A[3, 3] - \
-        A[0, 0]*A[1, 1]*A[2, 3]*A[3, 2] - A[0, 0]*A[1, 2]*A[2, 1]*A[3, 3] + \
-        A[0, 0]*A[1, 2]*A[2, 3]*A[3, 1] + A[0, 0]*A[1, 3]*A[2, 1]*A[3, 2] - \
-        A[0, 0]*A[1, 3]*A[2, 2]*A[3, 1] - A[0, 1]*A[1, 0]*A[2, 2]*A[3, 3] + \
-        A[0, 1]*A[1, 0]*A[2, 3]*A[3, 2] + A[0, 1]*A[1, 2]*A[2, 0]*A[3, 3] - \
-        A[0, 1]*A[1, 2]*A[2, 3]*A[3, 0] - A[0, 1]*A[1, 3]*A[2, 0]*A[3, 2] + \
-        A[0, 1]*A[1, 3]*A[2, 2]*A[3, 0] + A[0, 2]*A[1, 0]*A[2, 1]*A[3, 3] - \
-        A[0, 2]*A[1, 0]*A[2, 3]*A[3, 1] - A[0, 2]*A[1, 1]*A[2, 0]*A[3, 3] + \
-        A[0, 2]*A[1, 1]*A[2, 3]*A[3, 0] + A[0, 2]*A[1, 3]*A[2, 0]*A[3, 1] - \
-        A[0, 2]*A[1, 3]*A[2, 1]*A[3, 0] - A[0, 3]*A[1, 0]*A[2, 1]*A[3, 2] + \
-        A[0, 3]*A[1, 0]*A[2, 2]*A[3, 1] + A[0, 3]*A[1, 1]*A[2, 0]*A[3, 2] - \
-        A[0, 3]*A[1, 1]*A[2, 2]*A[3, 0] - A[0, 3]*A[1, 2]*A[2, 0]*A[3, 1] + \
-        A[0, 3]*A[1, 2]*A[2, 1]*A[3, 0]
 
 def test_MatrixElement_diff():
     assert (A[3, 0]*A[0, 0]).diff(A[0, 0]) == A[3, 0]
diff --git a/sympy/matrices/tests/test_interactions.py b/sympy/matrices/tests/test_interactions.py
index 7bb0d7991c..bdccd438f6 100644
--- a/sympy/matrices/tests/test_interactions.py
+++ b/sympy/matrices/tests/test_interactions.py
@@ -14,7 +14,6 @@
 from sympy.utilities.pytest import raises
 
 SM = MatrixSymbol('X', 3, 3)
-SV = MatrixSymbol('v', 3, 1)
 MM = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
 IM = ImmutableMatrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
 meye = eye(3)
@@ -50,16 +49,6 @@ def test_matrix_symbol_MM():
     assert Y[1, 1] == 1 + X[1, 1]
 
 
-def test_matrix_symbol_vector_matrix_multiplication():
-    A = MM * SV
-    B = IM * SV
-    assert A == B
-    C = (SV.T * MM.T).T
-    assert B == C
-    D = (SV.T * IM.T).T
-    assert C == D
-
-
 def test_indexing_interactions():
     assert (a * IM)[1, 1] == 5*a
     assert (SM + IM)[1, 1] == SM[1, 1] + IM[1, 1]
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index b9145f01f6..e2c796fd9a 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -58,7 +58,6 @@
 from sympy.core.compatibility import iterable
 from sympy.core.mul import _keep_coeff
 from sympy.utilities.pytest import raises, XFAIL
-from sympy.simplify import simplify
 
 from sympy.abc import a, b, c, d, p, q, t, w, x, y, z
 from sympy import MatrixSymbol
@@ -2882,7 +2881,8 @@ def test_cancel():
     M = MatrixSymbol('M', 5, 5)
     assert cancel(M[0,0] + 7) == M[0,0] + 7
     expr = sin(M[1, 4] + M[2, 1] * 5 * M[4, 0]) - 5 * M[1, 2] / z
-    assert cancel(expr) == (z*sin(M[1, 4] + M[2, 1] * 5 * M[4, 0]) - 5 * M[1, 2]) / z
+    assert cancel(expr) == expr
+
 
 def test_reduced():
     f = 2*x**4 + y**2 - x**2 + y**3
diff --git a/sympy/printing/tests/test_ccode.py b/sympy/printing/tests/test_ccode.py
index 26062c507e..b97368a2f8 100644
--- a/sympy/printing/tests/test_ccode.py
+++ b/sympy/printing/tests/test_ccode.py
@@ -446,8 +446,8 @@ def test_Matrix_printing():
         "M[3] = q[1] + q[2];\n"
         "M[4] = q[3];\n"
         "M[5] = 5;\n"
-        "M[6] = 2*q[4]/q[1];\n"
-        "M[7] = sqrt(q[0]) + 4;\n"
+        "M[6] = 2*q[4]*1.0/q[1];\n"
+        "M[7] = 4 + sqrt(q[0]);\n"
         "M[8] = 0;")
 
 
diff --git a/sympy/printing/tests/test_fcode.py b/sympy/printing/tests/test_fcode.py
index bdb96a69d2..c077a6b8c0 100644
--- a/sympy/printing/tests/test_fcode.py
+++ b/sympy/printing/tests/test_fcode.py
@@ -658,10 +658,10 @@ def test_Matrix_printing():
     assert fcode(m, M) == (
         "      M(1, 1) = sin(q(2, 1))\n"
         "      M(2, 1) = q(2, 1) + q(3, 1)\n"
-        "      M(3, 1) = 2*q(5, 1)/q(2, 1)\n"
+        "      M(3, 1) = 2*q(5, 1)*1.0/q(2, 1)\n"
         "      M(1, 2) = 0\n"
         "      M(2, 2) = q(4, 1)\n"
-        "      M(3, 2) = sqrt(q(1, 1)) + 4\n"
+        "      M(3, 2) = 4 + sqrt(q(1, 1))\n"
         "      M(1, 3) = cos(q(3, 1))\n"
         "      M(2, 3) = 5\n"
         "      M(3, 3) = 0")
diff --git a/sympy/printing/tests/test_jscode.py b/sympy/printing/tests/test_jscode.py
index bff26227e9..588039f440 100644
--- a/sympy/printing/tests/test_jscode.py
+++ b/sympy/printing/tests/test_jscode.py
@@ -361,6 +361,6 @@ def test_Matrix_printing():
         "M[3] = q[1] + q[2];\n"
         "M[4] = q[3];\n"
         "M[5] = 5;\n"
-        "M[6] = 2*q[4]/q[1];\n"
-        "M[7] = Math.sqrt(q[0]) + 4;\n"
+        "M[6] = 2*q[4]*1/q[1];\n"
+        "M[7] = 4 + Math.sqrt(q[0]);\n"
         "M[8] = 0;")
diff --git a/sympy/sets/tests/test_sets.py b/sympy/sets/tests/test_sets.py
index 4d4a5175f6..f1871c3022 100644
--- a/sympy/sets/tests/test_sets.py
+++ b/sympy/sets/tests/test_sets.py
@@ -1009,17 +1009,6 @@ def test_issue_10326():
     assert Interval(-oo, oo).contains(-oo) is S.false
 
 
-def test_issue_2799():
-    U = S.UniversalSet
-    a = Symbol('a', real=True)
-    inf_interval = Interval(a, oo)
-    R = S.Reals
-
-    assert U + inf_interval == inf_interval + U
-    assert U + R == R + U
-    assert R + inf_interval == inf_interval + R
-
-
 def test_issue_9706():
     assert Interval(-oo, 0).closure == Interval(-oo, 0, True, False)
     assert Interval(0, oo).closure == Interval(0, oo, False, True)
diff --git a/sympy/simplify/tests/test_trigsimp.py b/sympy/simplify/tests/test_trigsimp.py
index 641ed67bed..d753c844a8 100644
--- a/sympy/simplify/tests/test_trigsimp.py
+++ b/sympy/simplify/tests/test_trigsimp.py
@@ -316,7 +316,7 @@ def test_trigsimp_groebner():
     results = [resnum/resdenom, (-resnum)/(-resdenom)]
     assert trigsimp_groebner(ex) in results
     assert trigsimp_groebner(s/c, hints=[tan]) == tan(x)
-    assert trigsimp_groebner(c*s) == c*s
+
     assert trigsimp((-s + 1)/c + c/(-s + 1),
                     method='groebner') == 2/c
     assert trigsimp((-s + 1)/c + c/(-s + 1),
diff --git a/sympy/solvers/tests/test_diophantine.py b/sympy/solvers/tests/test_diophantine.py
index c502425159..95ed44ad4b 100644
--- a/sympy/solvers/tests/test_diophantine.py
+++ b/sympy/solvers/tests/test_diophantine.py
@@ -794,7 +794,7 @@ def test_diop_sum_of_even_powers():
     raises(NotImplementedError, lambda: diop_general_sum_of_even_powers(-eq, 2))
     neg = symbols('neg', negative=True)
     eq = x**4 + y**4 + neg**4 - 2673
-    assert diop_general_sum_of_even_powers(eq) == set([(-3, 6, 6)])
+    assert diop_general_sum_of_even_powers(eq) == set([(-2, 4, 7)])
     assert diophantine(x**4 + y**4 + 2) == set()
     assert diop_general_sum_of_even_powers(x**4 + y**4 - 2, limit=0) == set()
 
@@ -816,7 +816,7 @@ def test_sum_of_squares_powers():
     assert list(sum_of_squares(0, 3)) == []
     assert list(sum_of_squares(4, 1)) == [(2,)]
     assert list(sum_of_squares(5, 1)) == []
-    assert list(sum_of_squares(50, 2)) == [(5, 5), (1, 7)]
+    assert list(sum_of_squares(50, 2)) == [(1, 7), (5, 5)]
     assert list(sum_of_squares(11, 5, True)) == [
         (1, 1, 1, 2, 2), (0, 0, 1, 1, 3)]
     assert list(sum_of_squares(8, 8)) == [(1, 1, 1, 1, 1, 1, 1, 1)]
diff --git a/sympy/stats/tests/test_symbolic_probability.py b/sympy/stats/tests/test_symbolic_probability.py
index e578a3e063..baabbd5d0f 100644
--- a/sympy/stats/tests/test_symbolic_probability.py
+++ b/sympy/stats/tests/test_symbolic_probability.py
@@ -1,7 +1,7 @@
 from sympy import symbols, Mul, sin
 from sympy.stats import Normal, Poisson, variance
 from sympy.stats.rv import probability, expectation
-from sympy.stats import Covariance, Variance, Probability, Expectation
+from sympy.stats.symbolic_probability import Covariance, Variance, Probability, Expectation
 
 
 def test_literal_probability():
@@ -56,22 +56,3 @@ def test_literal_probability():
     assert Covariance(X, X**2) == Covariance(X, X**2, evaluate=False)
     assert Covariance(X, sin(X)) == Covariance(sin(X), X, evaluate=False)
     assert Covariance(X**2, sin(X)*Y) == Covariance(sin(X)*Y, X**2, evaluate=False)
-
-
-def test_probability_rewrite():
-    X = Normal('X', 2, 3)
-    Y = Normal('Y', 3, 4)
-    Z = Poisson('Z', 4)
-    W = Poisson('W', 3)
-    x, y, w, z = symbols('x, y, w, z')
-
-    assert Variance(w).rewrite(Expectation) == 0
-    assert Variance(X).rewrite(Expectation) == Expectation(X ** 2) - Expectation(X) ** 2
-    assert Variance(X, condition=Y).rewrite(Expectation) == Expectation(X ** 2, Y) - Expectation(X, Y) ** 2
-    assert Variance(X, Y) != Expectation(X**2) - Expectation(X)**2
-    assert Variance(X + z).rewrite(Expectation) == Expectation(X ** 2) - Expectation(X) ** 2
-    assert Variance(X * Y).rewrite(Expectation) == Expectation(X ** 2 * Y ** 2) - Expectation(X * Y) ** 2
-
-    assert Covariance(w, X).rewrite(Expectation) == 0
-    assert Covariance(X, Y).rewrite(Expectation) == Expectation(X*Y) - Expectation(X)*Expectation(Y)
-    assert Covariance(X, Y, condition=W).rewrite(Expectation) == Expectation(X * Y, W) - Expectation(X, W) * Expectation(Y, W)
diff --git a/sympy/tensor/tests/test_tensor.py b/sympy/tensor/tests/test_tensor.py
index bb8fcc5927..5a171d63df 100644
--- a/sympy/tensor/tests/test_tensor.py
+++ b/sympy/tensor/tests/test_tensor.py
@@ -1777,60 +1777,3 @@ def test_TensMul_data():
     assert ((F(alpha, mu) * F(beta, nu) * u(-alpha) * u(-beta)).data ==
             (E(mu) * E(nu)).data
             ).all()
-
-    S2 = TensorType([Lorentz] * 2, tensorsymmetry([1] * 2))
-    g = S2('g')
-    g.data = Lorentz.data
-
-    # tensor 'perp' is orthogonal to vector 'u'
-    perp = u(mu) * u(nu) + g(mu, nu)
-
-    mul_1 = u(-mu) * perp(mu, nu)
-    assert (mul_1.data == numpy.array([0, 0, 0, 0])).all()
-
-    mul_2 = u(-mu) * perp(mu, alpha) * perp(nu, beta)
-    assert (mul_2.data == numpy.zeros(shape=(4, 4, 4))).all()
-
-    Fperp = perp(mu, alpha) * perp(nu, beta) * F(-alpha, -beta)
-    assert (Fperp.data[0, :] == numpy.array([0, 0, 0, 0])).all()
-    assert (Fperp.data[:, 0] == numpy.array([0, 0, 0, 0])).all()
-
-    mul_3 = u(-mu) * Fperp(mu, nu)
-    assert (mul_3.data == numpy.array([0, 0, 0, 0])).all()
-
-
-def test_issue_11020_TensAdd_data():
-    numpy = import_module('numpy')
-    if numpy is None:
-        return
-
-    Lorentz = TensorIndexType('Lorentz', metric=False, dummy_fmt='i', dim=2)
-    Lorentz.data = [-1, 1]
-
-    a, b, c, d = tensor_indices('a, b, c, d', Lorentz)
-    i0, i1 = tensor_indices('i_0:2', Lorentz)
-
-    Vec = TensorType([Lorentz], tensorsymmetry([1]))
-    S2 = TensorType([Lorentz] * 2, tensorsymmetry([1] * 2))
-
-    # metric tensor
-    g = S2('g')
-    g.data = Lorentz.data
-
-    u = Vec('u')
-    u.data = [1, 0]
-
-    add_1 = g(b, c) * g(d, i0) * u(-i0) - g(b, c) * u(d)
-    assert (add_1.data == numpy.zeros(shape=(2, 2, 2))).all()
-    # Now let us replace index `d` with `a`:
-    add_2 = g(b, c) * g(a, i0) * u(-i0) - g(b, c) * u(a)
-    assert (add_2.data == numpy.zeros(shape=(2, 2, 2))).all()
-
-    # some more tests
-    # perp is tensor orthogonal to u^\mu
-    perp = u(a) * u(b) + g(a, b)
-    mul_1 = u(-a) * perp(a, b)
-    assert (mul_1.data == numpy.array([0, 0])).all()
-
-    mul_2 = u(-c) * perp(c, a) * perp(d, b)
-    assert (mul_2.data == numpy.zeros(shape=(2, 2, 2))).all()