Repository: sympy/sympy

Issue:
solve(x < oo) and solve(x > -oo) returns False
```python
>>> x = symbols('x')
>>> solve(x < oo)
>>> False
>>> solve(x > -oo)
>>> False
```
It should return `And(x > -oo, x < oo)`.

Discussion:
smichr: +1
smichr: Will commit in 24h if neither @hargup nor @skirpichev  have further comments.
hargup: +1, but I'm not very sure if the concerns by @skirpichev in https://github.com/sympy/sympy/pull/8259 are addressed here. I think we should wait for him to give a +1.
skirpichev: We have funny things in this pr:
```python
In [4]: solve(x >= oo)
Out[4]: x = ∞
In [3]: solve(x <= oo)
Out[3]: x = ∞
```
pelegm: > We have funny things in this pr

I'm ok with these "funny" things, as long as the assumptions on `x` allow it. For example, if `x` is real (but not `extended_real`), then the solution should be `[]` (or `False`, I'm not sure what is the convention).

However, as these assumptions are not properly handled in the real assumption system, I wouldn't bother to handle this now.
pelegm: Oh, I apologize: I was talking only about the **first** "funny" thing. The second funny thing is indeed funny, in the sense that we should consider it as a bug...
aktech: I think the first funny thing is related to this #8260 , which is a part of solvers, since inequality solver is dependent on solvers, that should be good for now.

For second funny thing `x <= oo`: which is read as x  is less OR equal to infinity for that solution should be:  `And(-oo < x, x < oo)` , I have FIXED that.
skirpichev: > For second funny thing x <= oo: which is read as x is less OR equal to infinity for that solution should be: And(-oo < x, x < oo) , I have FIXED that.

I don't think so:
```python
In [1]: solve(x <= oo)
Out[1]: -∞ < x ∧ x < ∞
```

Maybe it's a documentation bug, i.e. we disallow oo as solution (i.e. we search for solutions over real/complex domains, not in extended reals).

Patch:
diff --git a/AUTHORS b/AUTHORS
index bbe66d9108..f90115794e 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -375,4 +375,3 @@ Cody Herbst <cyherbst@gmail.com>
 AMiT Kumar <dtu.amit@gmail.com>
 Nishith Shah <nishithshah.2211@gmail.com>
 Guillaume Gay <contact@damcb.com>
-Ray Cathcart <github@cathcart.us>
diff --git a/doc/src/aboutus.rst b/doc/src/aboutus.rst
index affa4e79de..55d17a1dcf 100644
--- a/doc/src/aboutus.rst
+++ b/doc/src/aboutus.rst
@@ -380,7 +380,6 @@ want to be mentioned here, so see our repository history for a full list).
 #. AMiT Kumar: fix sign error in unrad
 #. Nishith Shah: fix solving of Piecewise functions
 #. Guillaume Gay: bugfix for LagrangesMethod
-#. Ray Cathcart: improve error handling in _random
 
 Up-to-date list in the order of the first contribution is given in the `AUTHORS
 <https://github.com/sympy/sympy/blob/master/AUTHORS>`_ file.
diff --git a/doc/src/install.rst b/doc/src/install.rst
index 7a67b2d296..e281939c48 100644
--- a/doc/src/install.rst
+++ b/doc/src/install.rst
@@ -3,12 +3,12 @@
 Installation
 ------------
 
-The SymPy CAS can be installed on virtually any computer with Python
-2.6 or above. SymPy does require `mpmath`_ Python library to be
-installed first.  The current recommended method of installation is
-directly from the source files.  Alternatively, executables are
-available for Windows, and some Linux distributions have SymPy
-packages available.
+The SymPy CAS can be installed on virtually any computer with Python 2.6 or
+above. SymPy does not require any special Python modules: let us know if you
+have any problems with SymPy on a standard Python install. The current
+recommended method of installation is directly from the source files.
+Alternatively, executables are available for Windows, and some Linux
+distributions have SymPy packages available.
 
 SymPy officially supports Python 2.6, 2.7, 3.2, 3.3, 3.4, and PyPy.
 
@@ -120,4 +120,3 @@ an `issue ticket`_.
 .. _Gitter: https://gitter.im/sympy/sympy
 .. _issue ticket: https://github.com/sympy/sympy/issues
 .. _mailing list: http://groups.google.com/group/sympy
-.. _mpmath: http://mpmath.org/
diff --git a/doc/src/modules/functions/combinatorial.rst b/doc/src/modules/functions/combinatorial.rst
index 21148aec87..d5c4a5d756 100644
--- a/doc/src/modules/functions/combinatorial.rst
+++ b/doc/src/modules/functions/combinatorial.rst
@@ -41,12 +41,6 @@ factorial
 .. autoclass:: sympy.functions.combinatorial.factorials.factorial
    :members:
 
-subfactorial
-------------
-
-.. autoclass:: sympy.functions.combinatorial.factorials.subfactorial
-   :members:
-
 factorial2 / double factorial
 -----------------------------
 
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 2b6875fabd..c6f4f9fd84 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -369,7 +369,7 @@ def _random(self, n=None, re_min=-1, im_min=-1, re_max=1, im_max=1):
                            for zi in free])))
             try:
                 nmag = abs(self.evalf(2, subs=reps))
-            except (ValueError, TypeError):
+            except TypeError:
                 # if an out of range value resulted in evalf problems
                 # then return None -- XXX is there a way to know how to
                 # select a good random number for a given expression?
diff --git a/sympy/core/function.py b/sympy/core/function.py
index 900cb31161..be6798e036 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -1556,7 +1556,9 @@ def _hashable_content(self):
 
     def _eval_subs(self, old, new):
         if old in self.variables:
-            return self
+            pts = list(self.point.args)
+            pts[self.variables.index(old)] = new
+            return self.func(self.expr, self.variables, pts)
 
     def _eval_derivative(self, s):
         if s not in self.free_symbols:
diff --git a/sympy/functions/combinatorial/factorials.py b/sympy/functions/combinatorial/factorials.py
index 185ae0d56c..20ffd3718c 100644
--- a/sympy/functions/combinatorial/factorials.py
+++ b/sympy/functions/combinatorial/factorials.py
@@ -1,6 +1,6 @@
 from __future__ import print_function, division
 
-from sympy.core import S, C, sympify, Dummy
+from sympy.core import S, C, sympify
 from sympy.core.function import Function, ArgumentIndexError
 from sympy.core.logic import fuzzy_and
 from sympy.ntheory import sieve
@@ -164,11 +164,6 @@ def eval(cls, n):
     def _eval_rewrite_as_gamma(self, n):
         return C.gamma(n + 1)
 
-    def _eval_rewrite_as_Product(self, n):
-        if n.is_nonnegative and n.is_integer:
-            i = Dummy('i', integer=True)
-            return C.Product(i, (i, 1, n))
-
     def _eval_is_integer(self):
         if self.args[0].is_integer and self.args[0].is_nonnegative:
             return True
@@ -205,15 +200,9 @@ class subfactorial(CombinatorialFunction):
     It can also be written as int(round(n!/exp(1))) but the recursive
     definition with caching is implemented for this function.
 
-    This function is generalized to noninteger arguments [2]_ as
-
-    .. math:: !x = \Gamma(x + 1, -1)/e
-
     References
     ==========
-
     .. [1] http://en.wikipedia.org/wiki/Subfactorial
-    .. [2] http://mathworld.wolfram.com/Subfactorial.html
 
     Examples
     ========
@@ -227,35 +216,32 @@ class subfactorial(CombinatorialFunction):
 
     See Also
     ========
-
-    sympy.functions.combinatorial.factorials.factorial,
-    sympy.utilities.iterables.generate_derangements,
-    sympy.functions.special.gamma_functions.uppergamma
+    factorial, sympy.utilities.iterables.generate_derangements
     """
 
     @classmethod
     @cacheit
     def _eval(self, n):
         if not n:
-            return S.One
+            return 1
         elif n == 1:
-            return S.Zero
+            return 0
         return (n - 1)*(self._eval(n - 1) + self._eval(n - 2))
 
     @classmethod
     def eval(cls, arg):
-        if arg.is_Number:
-            if arg.is_Integer and arg.is_nonnegative:
-                return cls._eval(arg)
-            elif arg is S.Infinity:
-                return arg
+        try:
+            arg = as_int(arg)
+            if arg < 0:
+                raise ValueError
+            return C.Integer(cls._eval(arg))
+        except ValueError:
+            if sympify(arg).is_Number:
+                raise ValueError("argument must be a nonnegative integer")
 
     def _eval_is_integer(self):
-        if self.args[0].is_integer and self.args[0].is_nonnegative:
-            return True
-
-    def _eval_rewrite_as_uppergamma(self, arg):
-        return C.uppergamma(arg + 1, -1)/S.Exp1
+        return fuzzy_and((self.args[0].is_integer,
+                          self.args[0].is_nonnegative))
 
 
 class factorial2(CombinatorialFunction):
diff --git a/sympy/physics/hep/gamma_matrices.py b/sympy/physics/hep/gamma_matrices.py
index 013c771f0e..2fe6756854 100644
--- a/sympy/physics/hep/gamma_matrices.py
+++ b/sympy/physics/hep/gamma_matrices.py
@@ -1,7 +1,7 @@
 from sympy import S
 from sympy.tensor.tensor import TensorIndexType, TensorIndex,\
     TensMul, TensorHead, tensorsymmetry, TensorType,\
-    TensAdd, tensor_mul, get_lines, Tensor
+    TensAdd, tensor_mul, get_lines
 from sympy.core.containers import Tuple
 
 
@@ -86,18 +86,13 @@ def extract_type_tens(expression):
 
 
         """
-        if isinstance(expression, Tensor):
-            sp = [expression]
-        elif isinstance(expression, TensMul):
-            sp = expression.args
-        else:
-            raise ValueError('wrong type')
+        sp = expression.split()
 
         # Collect all gamma matrices of the same dimension
         new_expr = S.One
         residual_expr = S.One
         for i in sp:
-            if isinstance(i, Tensor) and isinstance(i.args[0], GammaMatrixHead):
+            if isinstance(i.args[1][0], GammaMatrixHead):
                 new_expr *= i
             else:
                 residual_expr *= i
@@ -167,8 +162,8 @@ def _simplify_gpgp(ex):
                     if not ta:
                         ta = ex.split()
                         mu = TensorIndex('mu', GammaMatrix.LorentzIndex)
-                    ind1 = ta[ai[0]].get_indices()[1]
-                    ind2 = ta[ai[0] + 1].get_indices()[2]
+                    ind1 = ta[ai[0]].args[-1][1]
+                    ind2 = ta[ai[0] + 1].args[-1][2]
                     hit = True
                     if i == 0:
                         coeff = ex.coeff
@@ -377,7 +372,7 @@ def _trace_single_line1(t):
         if isinstance(t, TensAdd):
             a = [x.coeff*_trace_single_line1(x) for x in t.args]
             return TensAdd(*a)
-        elif isinstance(t, (Tensor, TensMul)):
+        elif isinstance(t, TensMul):
             r = t.coeff*_trace_single_line1(t)
             return r
         else:
@@ -393,14 +388,14 @@ def _gamma_trace1(self, *a):
             #return TensMul.from_data(S.Zero, [], [], [])
             return S.Zero
         if n == 2:
-            ind0 = a[0].get_indices()[0]
-            ind1 = a[1].get_indices()[0]
+            ind0 = a[0].args[-1][0]
+            ind1 = a[1].args[-1][0]
             return gctr*g(ind0, ind1)
         if n == 4:
-            ind0 = a[0].get_indices()[0]
-            ind1 = a[1].get_indices()[0]
-            ind2 = a[2].get_indices()[0]
-            ind3 = a[3].get_indices()[0]
+            ind0 = a[0].args[-1][0]
+            ind1 = a[1].args[-1][0]
+            ind2 = a[2].args[-1][0]
+            ind3 = a[3].args[-1][0]
 
             return gctr*(g(ind0, ind1)*g(ind2, ind3) - \
                g(ind0, ind2)*g(ind1, ind3) + g(ind0, ind3)*g(ind1, ind2))
@@ -467,7 +462,7 @@ def _kahane_simplify(coeff, tids):
 
         >>> tc = 3*G(i0)*G(i1)
         >>> G._kahane_simplify(tc.coeff, tc._tids)
-        3*gamma(i0, auto_left, S_0)*gamma(i1, -S_0, -auto_right)
+        3*gamma(i0, auto_left, -S_0)*gamma(i1, S_0, -auto_right)
 
         References
         ==========
diff --git a/sympy/polys/rootoftools.py b/sympy/polys/rootoftools.py
index f889495708..d3d9a03457 100644
--- a/sympy/polys/rootoftools.py
+++ b/sympy/polys/rootoftools.py
@@ -554,41 +554,39 @@ def _eval_evalf(self, prec):
 
             while True:
                 if self.is_real:
-                    a = mpf(str(interval.a))
-                    b = mpf(str(interval.b))
-                    if a == b:
-                        root = a
-                        break
                     x0 = mpf(str(interval.center))
                 else:
-                    ax = mpf(str(interval.ax))
-                    bx = mpf(str(interval.bx))
-                    ay = mpf(str(interval.ay))
-                    by = mpf(str(interval.by))
-                    if ax == bx and ay == by:
-                        root = ax + S.ImaginaryUnit*by
-                        break
                     x0 = mpc(*map(str, interval.center))
-
                 try:
-                    root = findroot(func, x0)
+                    root = findroot(func, x0, verify=False)
                     # If the (real or complex) root is not in the 'interval',
                     # then keep refining the interval. This happens if findroot
                     # accidentally finds a different root outside of this
                     # interval because our initial estimate 'x0' was not close
-                    # enough. It is also possible that the secant method will
-                    # get trapped by a max/min in the interval; the root
-                    # verification by findroot will raise a ValueError in this
-                    # case and the interval will then be tightened -- and
-                    # eventually the root will be found.
+                    # enough.
                     if self.is_real:
-                        if (a < root < b):
+                        a = mpf(str(interval.a))
+                        b = mpf(str(interval.b))
+                        if a == b:
+                            root = a
+                            break
+                        if not (a < root < b):
+                            raise ValueError("Root not in the interval.")
+                    else:
+                        ax = mpf(str(interval.ax))
+                        bx = mpf(str(interval.bx))
+                        ay = mpf(str(interval.ay))
+                        by = mpf(str(interval.by))
+                        if ax == bx and ay == by:
+                            root = ax + S.ImaginaryUnit*by
                             break
-                    elif (ax < root.real < bx and ay < root.imag < by):
-                        break
+                        if not (ax < root.real < bx and ay < root.imag < by):
+                            raise ValueError("Root not in the interval.")
                 except ValueError:
-                    pass
-                interval = interval.refine()
+                    interval = interval.refine()
+                    continue
+                else:
+                    break
 
         return Float._new(root.real._mpf_, prec) + I*Float._new(root.imag._mpf_, prec)
 
diff --git a/sympy/printing/str.py b/sympy/printing/str.py
index b7658005de..4c0b8bfa65 100644
--- a/sympy/printing/str.py
+++ b/sympy/printing/str.py
@@ -364,9 +364,6 @@ def _print_TensorIndex(self, expr):
     def _print_TensorHead(self, expr):
         return expr._print()
 
-    def _print_Tensor(self, expr):
-        return expr._print()
-
     def _print_TensMul(self, expr):
         return expr._print()
 
diff --git a/sympy/solvers/inequalities.py b/sympy/solvers/inequalities.py
index 33afd9742a..102e38e279 100644
--- a/sympy/solvers/inequalities.py
+++ b/sympy/solvers/inequalities.py
@@ -432,6 +432,12 @@ def valid(x):
         raise NotImplementedError('sorting of these roots is not supported')
     for x in reals:
         end = x
+
+        if ((end is S.NegativeInfinity and expr.rel_op in ['>', '>=']) or
+           (end is S.Infinity and expr.rel_op in ['<', '<='])):
+            sol_sets.append(Interval(start, S.Infinity, True, True))
+            break
+
         if valid((start + end)/2 if start != S.NegativeInfinity else end - 1):
             sol_sets.append(Interval(start, end, True, True))
 
diff --git a/sympy/tensor/tensor.py b/sympy/tensor/tensor.py
index b28aa45651..70bbc2b256 100644
--- a/sympy/tensor/tensor.py
+++ b/sympy/tensor/tensor.py
@@ -49,7 +49,7 @@
 
 class TIDS(CantSympify):
     """
-    Tensor-index data structure. This contains internal data structures about
+    Tensor internal data structure. This contains internal data structures about
     components of a tensor expression, its free and dummy indices.
 
     To create a ``TIDS`` object via the standard constructor, the required
@@ -107,22 +107,6 @@ def __init__(self, components, free, dum):
         self.free = free
         self.dum = dum
         self._ext_rank = len(self.free) + 2*len(self.dum)
-        self.dum.sort(key=lambda x: (x[2], x[0]))
-
-    def get_tensors(self):
-        """
-        Get a list of ``Tensor`` objects having the same ``TIDS`` if multiplied
-        by one another.
-        """
-        indices = self.get_indices()
-        components = self.components
-        tensors = [None for i in components]  # pre-allocate list
-        ind_pos = 0
-        for i, component in enumerate(components):
-            prev_pos = ind_pos
-            ind_pos += component.rank
-            tensors[i] = Tensor(component, indices[prev_pos:ind_pos])
-        return tensors
 
     def get_components_with_free_indices(self):
         """
@@ -220,9 +204,28 @@ def from_components_and_indices(components, indices):
 
         return tids
 
-    @deprecated(useinstead="get_indices")
     def to_indices(self):
-        return self.get_indices()
+        """
+        Get a list of indices, creating new tensor indices to complete dummy indices.
+        """
+        component_indices = []
+        for i in self.components:
+            component_indices.append([None]*i.rank)
+
+        for i in self.free:
+            component_indices[i[2]][i[1]] = i[0]
+
+        for i, dummy_pos in enumerate(self.dum):
+            tensor_index_type = self.components[dummy_pos[2]].args[1].args[0][0]
+            dummy_index = TensorIndex('dummy_index_{0}'.format(i), tensor_index_type)
+            component_indices[dummy_pos[2]][dummy_pos[0]] = dummy_index
+            component_indices[dummy_pos[3]][dummy_pos[1]] = -dummy_index
+
+        indices = []
+        for i in component_indices:
+            indices.extend(i)
+
+        return indices
 
     @staticmethod
     def free_dum_from_indices(*indices):
@@ -283,59 +286,6 @@ def free_dum_from_indices(*indices):
         free.sort()
         return free, dum
 
-    @staticmethod
-    def _check_matrix_indices(f_free, g_free, nc1):
-        # This "private" method checks matrix indices.
-        # Matrix indices are special as there are only two, and observe
-        # anomalous substitution rules to determine contractions.
-
-        dum = []
-        # make sure that free indices appear in the same order as in their component:
-        f_free.sort(key=lambda x: (x[2], x[1]))
-        g_free.sort(key=lambda x: (x[2], x[1]))
-        matrix_indices_storage = {}
-        transform_right_to_left = {}
-        f_pop_pos = []
-        g_pop_pos = []
-        for free_pos, (ind, i, c) in enumerate(f_free):
-            index_type = ind._tensortype
-            if ind not in (index_type.auto_left, -index_type.auto_right):
-                continue
-            matrix_indices_storage[ind] = (free_pos, i, c)
-
-        for free_pos, (ind, i, c) in enumerate(g_free):
-            index_type = ind._tensortype
-            if ind not in (index_type.auto_left, -index_type.auto_right):
-                continue
-
-            if ind == index_type.auto_left:
-                if -index_type.auto_right in matrix_indices_storage:
-                    other_pos, other_i, other_c = matrix_indices_storage.pop(-index_type.auto_right)
-                    dum.append((other_i, i, other_c, c + nc1))
-                    # mark to remove other_pos and free_pos from free:
-                    g_pop_pos.append(free_pos)
-                    f_pop_pos.append(other_pos)
-                    continue
-                if ind in matrix_indices_storage:
-                    other_pos, other_i, other_c = matrix_indices_storage.pop(ind)
-                    dum.append((other_i, i, other_c, c + nc1))
-                    # mark to remove other_pos and free_pos from free:
-                    g_pop_pos.append(free_pos)
-                    f_pop_pos.append(other_pos)
-                    transform_right_to_left[-index_type.auto_right] = c
-                    continue
-
-            if ind in transform_right_to_left:
-                other_c = transform_right_to_left.pop(ind)
-                if c == other_c:
-                    g_free[free_pos] = (index_type.auto_left, i, c)
-
-        for i in reversed(sorted(f_pop_pos)):
-            f_free.pop(i)
-        for i in reversed(sorted(g_pop_pos)):
-            g_free.pop(i)
-        return dum
-
     @staticmethod
     def mul(f, g):
         """
@@ -379,24 +329,18 @@ def mul(f, g):
         """
         index_up = lambda u: u if u.is_up else -u
 
-        # lambda returns True is index is not a matrix index:
-        notmat = lambda i: i not in (i._tensortype.auto_left, -i._tensortype.auto_right)
-        f_free = f.free[:]
-        g_free = g.free[:]
-        nc1 = len(f.components)
-        dum = TIDS._check_matrix_indices(f_free, g_free, nc1)
-
         # find out which free indices of f and g are contracted
-        free_dict1 = dict([(i if i.is_up else -i, (pos, cpos, i)) for i, pos, cpos in f_free])
-        free_dict2 = dict([(i if i.is_up else -i, (pos, cpos, i)) for i, pos, cpos in g_free])
+        free_dict1 = dict([(i if i.is_up else -i, (pos, cpos, i)) for i, pos, cpos in f.free])
+        free_dict2 = dict([(i if i.is_up else -i, (pos, cpos, i)) for i, pos, cpos in g.free])
+
         free_names = set(free_dict1.keys()) & set(free_dict2.keys())
         # find the new `free` and `dum`
-
+        nc1 = len(f.components)
         dum2 = [(i1, i2, c1 + nc1, c2 + nc1) for i1, i2, c1, c2 in g.dum]
-        free1 = [(ind, i, c) for ind, i, c in f_free if index_up(ind) not in free_names]
-        free2 = [(ind, i, c + nc1) for ind, i, c in g_free if index_up(ind) not in free_names]
+        free1 = [(ind, i, c) for ind, i, c in f.free if index_up(ind) not in free_names]
+        free2 = [(ind, i, c + nc1) for ind, i, c in g.free if index_up(ind) not in free_names]
         free = free1 + free2
-        dum.extend(f.dum + dum2)
+        dum = f.dum + dum2
         for name in free_names:
             ipos1, cpos1, ind1 = free_dict1[name]
             ipos2, cpos2, ind2 = free_dict2[name]
@@ -452,14 +396,6 @@ def sorted_components(self):
 
         return TIDS(components, free, dum), sign
 
-    def _get_sorted_free_indices_for_canon(self):
-        sorted_free = self.free[:]
-        sorted_free.sort(key=lambda x: x[0])
-        return sorted_free
-
-    def _get_sorted_dum_indices_for_canon(self):
-        return sorted(self.dum, key=lambda x: (x[2], x[0]))
-
     def canon_args(self):
         """
         Returns ``(g, dummies, msym, v)``, the entries of ``canonicalize``
@@ -482,7 +418,7 @@ def canon_args(self):
         # then the dummy indices, ordered by types and contravariant before
         # covariant
         # g[position in tensor] = position in ordered indices
-        for i, (indx, ipos, cpos) in enumerate(self._get_sorted_free_indices_for_canon()):
+        for i, (indx, ipos, cpos) in enumerate(self.free):
             pos = vpos[cpos] + ipos
             g[pos] = i
         pos = len(self.free)
@@ -491,7 +427,7 @@ def canon_args(self):
         prev = None
         a = []
         msym = []
-        for ipos1, ipos2, cpos1, cpos2 in self._get_sorted_dum_indices_for_canon():
+        for ipos1, ipos2, cpos1, cpos2 in self.dum:
             pos1 = vpos[cpos1] + ipos1
             pos2 = vpos[cpos2] + ipos2
             g[pos1] = j
@@ -542,7 +478,8 @@ def perm2tensor(self, g, canon_bp=False):
         for t in components:
             vpos.append(pos)
             pos += t._rank
-        sorted_free = [i[0] for i in self._get_sorted_free_indices_for_canon()]
+        sorted_free = [x[0] for x in self.free]
+        sorted_free.sort()
         nfree = len(sorted_free)
         rank = self._ext_rank
         dum = [[None]*4 for i in range((rank - nfree)//2)]
@@ -570,175 +507,6 @@ def perm2tensor(self, g, canon_bp=False):
 
         return TIDS(components, free, dum)
 
-    def get_indices(self):
-        """
-        Get a list of indices, creating new tensor indices to complete dummy indices.
-        """
-        components = self.components
-        free = self.free
-        dum = self.dum
-        indices = [None]*self._ext_rank
-        start = 0
-        pos = 0
-        vpos = []
-        for t in components:
-            vpos.append(pos)
-            pos += t.rank
-        cdt = defaultdict(int)
-        # if the free indices have names with dummy_fmt, start with an
-        # index higher than those for the dummy indices
-        # to avoid name collisions
-        for indx, ipos, cpos in free:
-            if indx._name.split('_')[0] == indx._tensortype._dummy_fmt[:-3]:
-                cdt[indx._tensortype] = max(cdt[indx._tensortype], int(indx._name.split('_')[1]) + 1)
-            start = vpos[cpos]
-            indices[start + ipos] = indx
-        for ipos1, ipos2, cpos1, cpos2 in dum:
-            start1 = vpos[cpos1]
-            start2 = vpos[cpos2]
-            typ1 = components[cpos1].index_types[ipos1]
-            assert typ1 == components[cpos2].index_types[ipos2]
-            fmt = typ1._dummy_fmt
-            nd = cdt[typ1]
-            indices[start1 + ipos1] = TensorIndex(fmt % nd, typ1)
-            indices[start2 + ipos2] = TensorIndex(fmt % nd, typ1, False)
-            cdt[typ1] += 1
-        return indices
-
-    def contract_metric(self, g):
-        """
-        Returns new TIDS and sign.
-
-        Sign is either 1 or -1, to correct the sign after metric contraction
-        (for spinor indices).
-        """
-        components = self.components
-        antisym = g.index_types[0].metric_antisym
-        #if not any(x == g for x in components):
-        #    return self
-        # list of positions of the metric ``g``
-        gpos = [i for i, x in enumerate(components) if x == g]
-        if not gpos:
-            return self, 1
-        sign = 1
-        dum = self.dum[:]
-        free = self.free[:]
-        elim = set()
-        for gposx in gpos:
-            if gposx in elim:
-                continue
-            free1 = [x for x in free if x[-1] == gposx]
-            dum1 = [x for x in dum if x[-2] == gposx or x[-1] == gposx]
-            if not dum1:
-                continue
-            elim.add(gposx)
-            if len(dum1) == 2:
-                if not antisym:
-                    dum10, dum11 = dum1
-                    if dum10[3] == gposx:
-                        # the index with pos p0 and component c0 is contravariant
-                        c0 = dum10[2]
-                        p0 = dum10[0]
-                    else:
-                        # the index with pos p0 and component c0 is covariant
-                        c0 = dum10[3]
-                        p0 = dum10[1]
-                    if dum11[3] == gposx:
-                        # the index with pos p1 and component c1 is contravariant
-                        c1 = dum11[2]
-                        p1 = dum11[0]
-                    else:
-                        # the index with pos p1 and component c1 is covariant
-                        c1 = dum11[3]
-                        p1 = dum11[1]
-                    dum.append((p0, p1, c0, c1))
-                else:
-                    dum10, dum11 = dum1
-                    # change the sign to bring the indices of the metric to contravariant
-                    # form; change the sign if dum10 has the metric index in position 0
-                    if dum10[3] == gposx:
-                        # the index with pos p0 and component c0 is contravariant
-                        c0 = dum10[2]
-                        p0 = dum10[0]
-                        if dum10[1] == 1:
-                            sign = -sign
-                    else:
-                        # the index with pos p0 and component c0 is covariant
-                        c0 = dum10[3]
-                        p0 = dum10[1]
-                        if dum10[0] == 0:
-                            sign = -sign
-                    if dum11[3] == gposx:
-                        # the index with pos p1 and component c1 is contravariant
-                        c1 = dum11[2]
-                        p1 = dum11[0]
-                        sign = -sign
-                    else:
-                        # the index with pos p1 and component c1 is covariant
-                        c1 = dum11[3]
-                        p1 = dum11[1]
-                    dum.append((p0, p1, c0, c1))
-
-            elif len(dum1) == 1:
-                if not antisym:
-                    dp0, dp1, dc0, dc1 = dum1[0]
-                    if dc0 == dc1:
-                        # g(i, -i)
-                        typ = g.index_types[0]
-                        if typ._dim is None:
-                            raise ValueError('dimension not assigned')
-                        sign = sign*typ._dim
-
-                    else:
-                        # g(i0, i1)*p(-i1)
-                        if dc0 == gposx:
-                            p1 = dp1
-                            c1 = dc1
-                        else:
-                            p1 = dp0
-                            c1 = dc0
-                        ind, p, c = free1[0]
-                        free.append((ind, p1, c1))
-                else:
-                    dp0, dp1, dc0, dc1 = dum1[0]
-                    if dc0 == dc1:
-                        # g(i, -i)
-                        typ = g.index_types[0]
-                        if typ._dim is None:
-                            raise ValueError('dimension not assigned')
-                        sign = sign*typ._dim
-
-                        if dp0 < dp1:
-                            # g(i, -i) = -D with antisymmetric metric
-                            sign = -sign
-                    else:
-                        # g(i0, i1)*p(-i1)
-                        if dc0 == gposx:
-                            p1 = dp1
-                            c1 = dc1
-                            if dp0 == 0:
-                                sign = -sign
-                        else:
-                            p1 = dp0
-                            c1 = dc0
-                        ind, p, c = free1[0]
-                        free.append((ind, p1, c1))
-            dum = [x for x in dum if x not in dum1]
-            free = [x for x in free if x not in free1]
-
-        shift = 0
-        shifts = [0]*len(components)
-        for i in range(len(components)):
-            if i in elim:
-                shift += 1
-                continue
-            shifts[i] = shift
-        free = [(ind, p, c - shifts[c]) for (ind, p, c) in free if c not in elim]
-        dum = [(p0, p1, c0 - shifts[c0], c1 - shifts[c1]) for  i, (p0, p1, c0, c1) in enumerate(dum) if c0 not in elim and c1 not in elim]
-        components = [c for i, c in enumerate(components) if i not in elim]
-        tids = TIDS(components, free, dum)
-        return tids, sign
-
 
 class VTIDS(TIDS):
     """
@@ -826,16 +594,6 @@ def _get(self, key):
         if isinstance(key, TensorHead):
             return None
 
-        if isinstance(key, Tensor):
-            # special case to handle metrics. Metric tensors cannot be
-            # constructed through contraction by the metric, their
-            # components show if they are a matrix or its inverse.
-            signature = tuple([i.is_up for i in key.get_indices()])
-            srch = (key.component,) + signature
-            if srch in self._substitutions_dict_tensmul:
-                return self._substitutions_dict_tensmul[srch]
-            return self.data_tensmul_from_tensorhead(key, key.component)
-
         if isinstance(key, TensMul):
             tensmul_list = key.split()
             if len(tensmul_list) == 1 and len(tensmul_list[0].components) == 1:
@@ -1676,7 +1434,6 @@ def __neg__(self):
                 (not self._is_up))
         return t1
 
-
 def tensor_indices(s, typ):
     """
     Returns list of tensor indices given their names and their types
@@ -2227,7 +1984,7 @@ def _check_auto_matrix_indices_in_call(self, *indices):
             if not self._matrix_behavior:
                 raise ValueError('wrong number of indices')
 
-            # Take the last one or two missing
+            # _matrix_behavior is True, so take the last one or two missing
             # indices as auto-matrix indices:
             ldiff = len(self.index_types) - len(indices)
             if ldiff > 2:
@@ -2265,7 +2022,7 @@ def _check_auto_matrix_indices_in_call(self, *indices):
 
         return indices, matrix_behavior_kinds
 
-    def __call__(self, *indices, **kw_args):
+    def __call__(self, *indices):
         """
         Returns a tensor with indices.
 
@@ -2314,7 +2071,7 @@ def __call__(self, *indices, **kw_args):
         Auto-matrix indices are automatically contracted upon multiplication,
 
         >>> r*s
-        A(auto_left, L_0)*B(-L_0, -auto_right, auto_left, -auto_right)
+        A(auto_left, -L_0)*B(L_0, -auto_right, auto_left, -auto_right)
 
         The multiplication algorithm has found an ``auto_right`` index in ``A``
         and an ``auto_left`` index in ``B`` referring to the same
@@ -2340,8 +2097,14 @@ def __call__(self, *indices, **kw_args):
         """
 
         indices, matrix_behavior_kinds = self._check_auto_matrix_indices_in_call(*indices)
-        tensor = Tensor._new_with_dummy_replacement(self, indices, **kw_args)
-        return tensor
+
+        components = [self]
+        tids = TIDS.from_components_and_indices(components, indices)
+
+        tmul = TensMul.from_TIDS(S.One, tids)
+        tmul._matrix_behavior_kinds = matrix_behavior_kinds
+
+        return tmul
 
     def __pow__(self, other):
         if self.data is None:
@@ -2613,26 +2376,20 @@ def __new__(cls, *args, **kw_args):
         if not args:
             return S.Zero
 
-        if len(args) == 1 and not isinstance(args[0], TensExpr):
-            return args[0]
-
         # replace auto-matrix indices so that they are the same in all addends
         args = TensAdd._tensAdd_check_automatrix(args)
 
         # now check that all addends have the same indices:
         TensAdd._tensAdd_check(args)
+        args = Tuple(*args)
 
         # if TensAdd has only 1 TensMul element in its `args`:
         if len(args) == 1 and isinstance(args[0], TensMul):
             obj = Basic.__new__(cls, *args, **kw_args)
             return obj
 
-        # TODO: do not or do canonicalize by default?
-        # Technically, one may wish to have additions of non-canonicalized
-        # tensors. This feature should be removed in the future.
-        # Unfortunately this would require to rewrite a lot of tests.
         # canonicalize all TensMul
-        args = [canon_bp(x) for x in args if x]
+        args = [x.canon_bp() for x in args if x]
         args = [x for x in args if x]
 
         # if there are no more args (i.e. have cancelled out),
@@ -2640,22 +2397,18 @@ def __new__(cls, *args, **kw_args):
         if not args:
             return S.Zero
 
-        if len(args) == 1:
-            return args[0]
-
         # collect canonicalized terms
-        def sort_key(t):
-            x = get_tids(t)
-            return (x.components, x.free, x.dum)
-        args.sort(key=sort_key)
-        args = TensAdd._tensAdd_collect_terms(args)
-        if not args:
+        args.sort(key=lambda x: (x.components, x.free, x.dum))
+        a = TensAdd._tensAdd_collect_terms(args)
+        if not a:
             return S.Zero
         # it there is only a component tensor return it
-        if len(args) == 1:
-            return args[0]
+        if len(a) == 1:
+            return a[0]
 
+        args = Tuple(*args)
         obj = Basic.__new__(cls, *args, **kw_args)
+        obj._args = tuple(a)
         return obj
 
     @staticmethod
@@ -2670,7 +2423,7 @@ def _tensAdd_flatten(args):
                         args1.extend(list(x.args))
                     else:
                         args1.append(x)
-            args1 = [x for x in args1 if isinstance(x, TensExpr) and x.coeff]
+            args1 = [x for x in args1 if isinstance(x, TensExpr) and x._coeff]
             args2 = [x for x in args if not isinstance(x, TensExpr)]
             t1 = TensMul.from_data(Add(*args2), [], [], [])
             args = [t1] + args1
@@ -2680,7 +2433,7 @@ def _tensAdd_flatten(args):
                 a.extend(list(x.args))
             else:
                 a.append(x)
-        args = [x for x in a if x.coeff]
+        args = [x for x in a if x._coeff]
         return args
 
     @staticmethod
@@ -2699,7 +2452,7 @@ def _tensAdd_check_automatrix(args):
         for i, arg in enumerate(args):
             arg_auto_left_types = set([])
             arg_auto_right_types = set([])
-            for index in get_indices(arg):
+            for index in arg.get_indices():
                 # @type index: TensorIndex
                 if index in (index._tensortype.auto_left, -index._tensortype.auto_left):
                     auto_left_types.add(index._tensortype)
@@ -2723,8 +2476,8 @@ def _tensAdd_check_automatrix(args):
     @staticmethod
     def _tensAdd_check(args):
         # check that all addends have the same free indices
-        indices0 = set([x[0] for x in get_tids(args[0]).free])
-        list_indices = [set([y[0] for y in get_tids(x).free]) for x in args[1:]]
+        indices0 = set([x[0] for x in args[0].free])
+        list_indices = [set([y[0] for y in x.free]) for x in args[1:]]
         if not all(x == indices0 for x in list_indices):
             raise ValueError('all tensors must have the same indices')
 
@@ -2733,16 +2486,14 @@ def _tensAdd_collect_terms(args):
         # collect TensMul terms differing at most by their coefficient
         a = []
         prev = args[0]
-        prev_coeff = get_coeff(prev)
+        prev_coeff = prev._coeff
         changed = False
 
         for x in args[1:]:
             # if x and prev have the same tensor, update the coeff of prev
-            x_tids = get_tids(x)
-            prev_tids = get_tids(prev)
-            if x_tids.components == prev_tids.components \
-                    and x_tids.free == prev_tids.free and x_tids.dum == prev_tids.dum:
-                prev_coeff = prev_coeff + get_coeff(x)
+            if x.components == prev.components \
+                    and x.free == prev.free and x.dum == prev.dum:
+                prev_coeff = prev_coeff + x._coeff
                 changed = True
                 op = 0
             else:
@@ -2753,18 +2504,18 @@ def _tensAdd_collect_terms(args):
                 else:
                     # get a tensor from prev with coeff=prev_coeff and store it
                     if prev_coeff:
-                        t = TensMul.from_data(prev_coeff, prev_tids.components,
-                            prev_tids.free, prev_tids.dum)
+                        t = TensMul.from_data(prev_coeff, prev.components,
+                            prev.free, prev.dum)
                         a.append(t)
                 # move x to prev
                 op = 1
                 pprev, prev = prev, x
-                pprev_coeff, prev_coeff = prev_coeff, get_coeff(x)
+                pprev_coeff, prev_coeff = prev_coeff, x._coeff
                 changed = False
         # if the case op=0 prev was not stored; store it now
         # in the case op=1 x was not stored; store it now (as prev)
         if op == 0 and prev_coeff:
-            prev = TensMul.from_data(prev_coeff, prev_tids.components, prev_tids.free, prev_tids.dum)
+            prev = TensMul.from_data(prev_coeff, prev.components, prev.free, prev.dum)
             a.append(prev)
         elif op == 1:
             a.append(prev)
@@ -2809,7 +2560,7 @@ def __call__(self, *indices):
         if indices == free_args:
             return self
         index_tuples = list(zip(free_args, indices))
-        a = [x.func(*x.fun_eval(*index_tuples).args) for x in self.args]
+        a = [TensMul(*x.fun_eval(*index_tuples).args) for x in self.args]
         res = TensAdd(*a)
 
         return res
@@ -2899,7 +2650,7 @@ def contract_metric(self, g):
         see the ``TensorIndexType`` docstring for the contraction conventions
         """
 
-        args = [contract_metric(x, g) for x in self.args]
+        args = [x.contract_metric(g) for x in self.args]
         t = TensAdd(*args)
         return canon_bp(t)
 
@@ -3024,369 +2775,98 @@ def __iter__(self):
 
 
 @doctest_depends_on(modules=('numpy',))
-class Tensor(TensExpr):
+class TensMul(TensExpr):
     """
-    Base tensor class, i.e. this represents a tensor, the single unit to be
-    put into an expression.
+    Product of tensors
 
-    This object is usually created from a ``TensorHead``, by attaching indices
-    to it. Indices preceded by a minus sign are considered contravariant,
-    otherwise covariant.
+    Parameters
+    ==========
 
-    Examples
-    ========
+    coeff : SymPy coefficient of the tensor
+    args
 
-    >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead
-    >>> Lorentz = TensorIndexType("Lorentz", dummy_fmt="L")
-    >>> mu, nu = tensor_indices('mu nu', Lorentz)
-    >>> A = tensorhead("A", [Lorentz, Lorentz], [[1], [1]])
-    >>> A(mu, -nu)
-    A(mu, -nu)
-    >>> A(mu, -mu)
-    A(L_0, -L_0)
+    Attributes
+    ==========
+
+    ``components`` : list of ``TensorHead`` of the component tensors
+    ``types`` : list of nonrepeated ``TensorIndexType``
+    ``free`` : list of ``(ind, ipos, icomp)``, see Notes
+    ``dum`` : list of ``(ipos1, ipos2, icomp1, icomp2)``, see Notes
+    ``ext_rank`` : rank of the tensor counting the dummy indices
+    ``rank`` : rank of the tensor
+    ``coeff`` : SymPy coefficient of the tensor
+    ``free_args`` : list of the free indices in sorted order
+    ``is_canon_bp`` : ``True`` if the tensor in in canonical form
+
+    Notes
+    =====
+
+    ``args[0]``   list of ``TensorHead`` of the component tensors.
+
+    ``args[1]``   list of ``(ind, ipos, icomp)``
+    where ``ind`` is a free index, ``ipos`` is the slot position
+    of ``ind`` in the ``icomp``-th component tensor.
+
+    ``args[2]`` list of tuples representing dummy indices.
+    ``(ipos1, ipos2, icomp1, icomp2)`` indicates that the contravariant
+    dummy index is the ``ipos1``-th slot position in the ``icomp1``-th
+    component tensor; the corresponding covariant index is
+    in the ``ipos2`` slot position in the ``icomp2``-th component tensor.
 
     """
 
-    is_commutative = False
+    def __new__(cls, coeff, *args, **kw_args):
+        coeff = sympify(coeff)
+
+        if len(args) == 2:
+            components = args[0]
+            indices = args[1]
+            tids = TIDS.from_components_and_indices(components, indices)
+        elif len(args) == 1:
+            tids = args[0]
+            components = tids.components
+            indices = tids.to_indices()
+        else:
+            raise TypeError("wrong construction")
 
-    def __new__(cls, tensor_head, indices, **kw_args):
-        tids = TIDS.from_components_and_indices((tensor_head,), indices)
-        obj = Basic.__new__(cls, tensor_head, Tuple(*indices), **kw_args)
+        for i in indices:
+            if not isinstance(i, TensorIndex):
+                raise TypeError("i should be of type TensorIndex")
+
+        t_components = Tuple(*components)
+        t_indices = Tuple(*indices)
+
+        obj = Basic.__new__(cls, coeff, t_components, t_indices)
+        obj._types = []
+        for t in tids.components:
+            obj._types.extend(t._types)
         obj._tids = tids
-        obj._indices = indices
+        obj._ext_rank = len(obj._tids.free) + 2*len(obj._tids.dum)
+        obj._coeff = coeff
         obj._is_canon_bp = kw_args.get('is_canon_bp', False)
+        obj._matrix_behavior_kinds = dict()
         return obj
 
     @staticmethod
-    def _new_with_dummy_replacement(tensor_head, indices, **kw_args):
-        tids = TIDS.from_components_and_indices((tensor_head,), indices)
-        indices = tids.get_indices()
-        return Tensor(tensor_head, indices, **kw_args)
+    def from_data(coeff, components, free, dum, **kw_args):
+        tids = TIDS(components, free, dum)
+        return TensMul.from_TIDS(coeff, tids, **kw_args)
+
+    @staticmethod
+    def from_TIDS(coeff, tids, **kw_args):
+        return TensMul(coeff, tids, **kw_args)
 
     @property
-    def is_canon_bp(self):
-        return self._is_canon_bp
+    def free_args(self):
+        return sorted([x[0] for x in self.free])
 
     @property
-    def indices(self):
-        return self._indices
+    def components(self):
+        return self._tids.components[:]
 
     @property
     def free(self):
-        return self._tids.free
-
-    @property
-    def dum(self):
-        return self._tids.dum
-
-    @property
-    def rank(self):
-        return len(self.free)
-
-    @property
-    def free_args(self):
-        return sorted([x[0] for x in self.free])
-
-    def perm2tensor(self, g, canon_bp=False):
-        """
-        Returns the tensor corresponding to the permutation ``g``
-
-        For further details, see the method in ``TIDS`` with the same name.
-        """
-        return perm2tensor(self, g, canon_bp)
-
-    def canon_bp(self):
-        if self._is_canon_bp:
-            return self
-        g, dummies, msym, v = self._tids.canon_args()
-        can = canonicalize(g, dummies, msym, *v)
-        if can == 0:
-            return S.Zero
-        tensor = self.perm2tensor(can, True)
-        return tensor
-
-    @property
-    def types(self):
-        return get_tids(self).components[0]._types
-
-    @property
-    def coeff(self):
-        return S.One
-
-    @property
-    def component(self):
-        return self.args[0]
-
-    @property
-    def components(self):
-        return [self.args[0]]
-
-    def split(self):
-        return [self]
-
-    def expand(self):
-        return self
-
-    def sorted_components(self):
-        return self
-
-    def get_indices(self):
-        """
-        Get a list of indices, corresponding to those of the tensor.
-        """
-        return self._tids.get_indices()
-
-    def as_base_exp(self):
-        return self, S.One
-
-    def substitute_indices(self, *index_tuples):
-        return substitute_indices(self, *index_tuples)
-
-    def __call__(self, *indices):
-        """Returns tensor with ordered free indices replaced by ``indices``
-
-        Examples
-        ========
-
-        >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead
-        >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
-        >>> i0,i1,i2,i3,i4 = tensor_indices('i0:5', Lorentz)
-        >>> A = tensorhead('A', [Lorentz]*5, [[1]*5])
-        >>> t = A(i2, i1, -i2, -i3, i4)
-        >>> t
-        A(L_0, i1, -L_0, -i3, i4)
-        >>> t(i1, i2, i3)
-        A(L_0, i1, -L_0, i2, i3)
-        """
-
-        free_args = self.free_args
-        indices = list(indices)
-        if [x._tensortype for x in indices] != [x._tensortype for x in free_args]:
-            raise ValueError('incompatible types')
-        if indices == free_args:
-            return self
-        t = self.fun_eval(*list(zip(free_args, indices)))
-
-        # object is rebuilt in order to make sure that all contracted indices
-        # get recognized as dummies, but only if there are contracted indices.
-        if len(set(i if i.is_up else -i for i in indices)) != len(indices):
-            return t.func(*t.args)
-        return t
-
-    def fun_eval(self, *index_tuples):
-        free = self.free
-        free1 = []
-        for j, ipos, cpos in free:
-            # search j in index_tuples
-            for i, v in index_tuples:
-                if i == j:
-                    free1.append((v, ipos, cpos))
-                    break
-            else:
-                free1.append((j, ipos, cpos))
-        return TensMul.from_data(self.coeff, self.components, free1, self.dum)
-
-    # TODO: put this into TensExpr?
-    def __iter__(self):
-        return self.data.flatten().__iter__()
-
-    # TODO: put this into TensExpr?
-    def __getitem__(self, item):
-        return self.data[item]
-
-    @property
-    def data(self):
-        return _tensor_data_substitution_dict[self]
-
-    @data.setter
-    def data(self, data):
-        # TODO: check data compatibility with properties of tensor.
-        _tensor_data_substitution_dict[self] = data
-
-    @data.deleter
-    def data(self):
-        if self in _tensor_data_substitution_dict:
-            del _tensor_data_substitution_dict[self]
-        if self.metric in _tensor_data_substitution_dict:
-            del _tensor_data_substitution_dict[self.metric]
-
-    def __mul__(self, other):
-        if isinstance(other, TensAdd):
-            return TensAdd(*[self*arg for arg in other.args])
-        tmul = TensMul(self, other)
-        return tmul
-
-    def __rmul__(self, other):
-        return TensMul(other, self)
-
-    def __div__(self, other):
-        if isinstance(other, TensExpr):
-            raise ValueError('cannot divide by a tensor')
-        return TensMul(self, S.One/other, is_canon_bp=self.is_canon_bp)
-
-    def __rdiv__(self, other):
-        raise ValueError('cannot divide by a tensor')
-
-    def __add__(self, other):
-        return TensAdd(self, other)
-
-    def __radd__(self, other):
-        return TensAdd(other, self)
-
-    def __sub__(self, other):
-        return TensAdd(self, -other)
-
-    def __rsub__(self, other):
-        return TensAdd(other, self)
-
-    __truediv__ = __div__
-    __rtruediv__ = __rdiv__
-
-    def __neg__(self):
-        return TensMul(S.NegativeOne, self)
-
-    def _print(self):
-        indices = [str(ind) for ind in self.indices]
-        component = self.component
-        if component.rank > 0:
-            return ('%s(%s)' % (component.name, ', '.join(indices)))
-        else:
-            return ('%s' % component.name)
-
-    def equals(self, other):
-        if other == 0:
-            return self.coeff == 0
-        other = sympify(other)
-        if not isinstance(other, TensExpr):
-            assert not self.components
-            return S.One == other
-
-        def _get_compar_comp(self):
-            t = self.canon_bp()
-            r = (t.coeff, tuple(t.components), \
-                    tuple(sorted(t.free)), tuple(sorted(t.dum)))
-            return r
-
-        return _get_compar_comp(self) == _get_compar_comp(other)
-
-    def contract_metric(self, metric):
-        tids, sign = get_tids(self).contract_metric(metric)
-        return TensMul.from_TIDS(sign, tids)
-
-    def contract_delta(self, metric):
-        return self.contract_metric(metric)
-
-
-@doctest_depends_on(modules=('numpy',))
-class TensMul(TensExpr):
-    """
-    Product of tensors
-
-    Parameters
-    ==========
-
-    coeff : SymPy coefficient of the tensor
-    args
-
-    Attributes
-    ==========
-
-    ``components`` : list of ``TensorHead`` of the component tensors
-    ``types`` : list of nonrepeated ``TensorIndexType``
-    ``free`` : list of ``(ind, ipos, icomp)``, see Notes
-    ``dum`` : list of ``(ipos1, ipos2, icomp1, icomp2)``, see Notes
-    ``ext_rank`` : rank of the tensor counting the dummy indices
-    ``rank`` : rank of the tensor
-    ``coeff`` : SymPy coefficient of the tensor
-    ``free_args`` : list of the free indices in sorted order
-    ``is_canon_bp`` : ``True`` if the tensor in in canonical form
-
-    Notes
-    =====
-
-    ``args[0]``   list of ``TensorHead`` of the component tensors.
-
-    ``args[1]``   list of ``(ind, ipos, icomp)``
-    where ``ind`` is a free index, ``ipos`` is the slot position
-    of ``ind`` in the ``icomp``-th component tensor.
-
-    ``args[2]`` list of tuples representing dummy indices.
-    ``(ipos1, ipos2, icomp1, icomp2)`` indicates that the contravariant
-    dummy index is the ``ipos1``-th slot position in the ``icomp1``-th
-    component tensor; the corresponding covariant index is
-    in the ``ipos2`` slot position in the ``icomp2``-th component tensor.
-
-    """
-
-    def __new__(cls, *args, **kw_args):
-        # make sure everything is sympified:
-        args = [sympify(arg) for arg in args]
-
-        # flatten:
-        args = TensMul._flatten(args)
-
-        is_canon_bp = kw_args.get('is_canon_bp', False)
-        if not any([isinstance(arg, TensExpr) for arg in args]):
-            tids = TIDS([], [], [])
-        else:
-            tids_list = [arg._tids for arg in args if isinstance(arg, (Tensor, TensMul))]
-            if len(tids_list) == 1:
-                for arg in args:
-                    if not isinstance(arg, Tensor):
-                        continue
-                    is_canon_bp = kw_args.get('is_canon_bp', arg._is_canon_bp)
-            tids = functools.reduce(lambda a, b: a*b, tids_list)
-
-        if any([isinstance(arg, TensAdd) for arg in args]):
-            add_args = TensAdd._tensAdd_flatten(args)
-            return TensAdd(*add_args)
-        coeff = functools.reduce(lambda a, b: a*b, [S.One] + [arg for arg in args if not isinstance(arg, TensExpr)])
-        args = tids.get_tensors()
-        if coeff != 1:
-            args = [coeff] + args
-        if len(args) == 1:
-            return args[0]
-
-        obj = Basic.__new__(cls, *args)
-        obj._types = []
-        for t in tids.components:
-            obj._types.extend(t._types)
-        obj._tids = tids
-        obj._ext_rank = len(obj._tids.free) + 2*len(obj._tids.dum)
-        obj._coeff = coeff
-        obj._is_canon_bp = is_canon_bp
-        return obj
-
-    @staticmethod
-    def _flatten(args):
-        a = []
-        for arg in args:
-            if isinstance(arg, TensMul):
-                a.extend(arg.args)
-            else:
-                a.append(arg)
-        return a
-
-    @staticmethod
-    def from_data(coeff, components, free, dum, **kw_args):
-        tids = TIDS(components, free, dum)
-        return TensMul.from_TIDS(coeff, tids, **kw_args)
-
-    @staticmethod
-    def from_TIDS(coeff, tids, **kw_args):
-        return TensMul(coeff, *tids.get_tensors(), **kw_args)
-
-    @property
-    def free_args(self):
-        return sorted([x[0] for x in self.free])
-
-    @property
-    def components(self):
-        return self._tids.components[:]
-
-    @property
-    def free(self):
-        return self._tids.free[:]
+        return self._tids.free[:]
 
     @property
     def coeff(self):
@@ -3414,7 +2894,7 @@ def equals(self, other):
 
         def _get_compar_comp(self):
             t = self.canon_bp()
-            r = (get_coeff(t), tuple(t.components), \
+            r = (t._coeff, tuple(t.components), \
                     tuple(sorted(t.free)), tuple(sorted(t.dum)))
             return r
 
@@ -3441,7 +2921,34 @@ def get_indices(self):
         >>> t.get_indices()
         [m1, m0, m2]
         """
-        return self._tids.get_indices()
+        indices = [None]*self._ext_rank
+        start = 0
+        pos = 0
+        vpos = []
+        components = self.components
+        for t in components:
+            vpos.append(pos)
+            pos += t._rank
+        cdt = defaultdict(int)
+        # if the free indices have names with dummy_fmt, start with an
+        # index higher than those for the dummy indices
+        # to avoid name collisions
+        for indx, ipos, cpos in self.free:
+            if indx._name.split('_')[0] == indx._tensortype._dummy_fmt[:-3]:
+                cdt[indx._tensortype] = max(cdt[indx._tensortype], int(indx._name.split('_')[1]) + 1)
+            start = vpos[cpos]
+            indices[start + ipos] = indx
+        for ipos1, ipos2, cpos1, cpos2 in self.dum:
+            start1 = vpos[cpos1]
+            start2 = vpos[cpos2]
+            typ1 = components[cpos1].index_types[ipos1]
+            assert typ1 == components[cpos2].index_types[ipos2]
+            fmt = typ1._dummy_fmt
+            nd = cdt[typ1]
+            indices[start1 + ipos1] = TensorIndex(fmt % nd, typ1)
+            indices[start2 + ipos2] = TensorIndex(fmt % nd, typ1, False)
+            cdt[typ1] += 1
+        return indices
 
     def split(self):
         """
@@ -3464,17 +2971,18 @@ def split(self):
         >>> t.split()
         [A(a, L_0), B(-L_0, c)]
         """
-        if self.args == ():
-            return [self]
-        splitp = []
-        res = 1
-        for arg in self.args:
-            if isinstance(arg, Tensor):
-                splitp.append(res*arg)
-                res = 1
-            else:
-                res *= arg
-        return splitp
+        indices = self.get_indices()
+        pos = 0
+        components = self.components
+        if not components:
+            return [TensMul.from_data(self._coeff, [], [], [])]
+        res = []
+        for t in components:
+            t1 = t(*indices[pos:pos + t._rank])
+            pos += t._rank
+            res.append(t1)
+        res[0] = TensMul.from_data(self._coeff, res[0].components, res[0]._tids.free, res[0]._tids.dum, is_canon_bp=res[0]._is_canon_bp)
+        return res
 
     def __add__(self, other):
         return TensAdd(self, other)
@@ -3510,21 +3018,55 @@ def __mul__(self, other):
         """
         other = sympify(other)
         if not isinstance(other, TensExpr):
-            coeff = self.coeff*other
+            coeff = self._coeff*other
             tmul = TensMul.from_TIDS(coeff, self._tids, is_canon_bp=self._is_canon_bp)
+            tmul._matrix_behavior_kinds = self._matrix_behavior_kinds
             return tmul
         if isinstance(other, TensAdd):
             return TensAdd(*[self*x for x in other.args])
 
-        new_tids = self._tids*other._tids
-        coeff = self.coeff*other.coeff
+        matrix_behavior_kinds = dict()
+
+        self_matrix_behavior_kinds = self._matrix_behavior_kinds
+        other_matrix_behavior_kinds = other._matrix_behavior_kinds
+        self2 = self
+
+        for key, v1 in self_matrix_behavior_kinds.items():
+            if key in other_matrix_behavior_kinds:
+                v2 = other_matrix_behavior_kinds[key]
+                if len(v1) == 1:
+                    other = other.substitute_indices((v2[0], -v1[0]))
+                    if len(v2) == 2:
+                        matrix_behavior_kinds[key] = (v2[1],)
+                elif len(v1) == 2:
+                    auto_index = v1[1]._tensortype.auto_index
+                    self2 = self2.substitute_indices((v1[1], -auto_index))
+                    other = other.substitute_indices((v2[0], auto_index))
+                    if len(v2) == 1:
+                        matrix_behavior_kinds[key] = (v1[0],)
+                    elif len(v2) == 2:
+                        matrix_behavior_kinds[key] = (v1[0], v2[1])
+            else:
+                matrix_behavior_kinds[key] = v1
+
+        for key, v2 in other_matrix_behavior_kinds.items():
+            if key in self_matrix_behavior_kinds:
+                continue
+            matrix_behavior_kinds[key] = v2
+
+        new_tids = self2._tids*other._tids
+        coeff = self2._coeff*other._coeff
         tmul = TensMul.from_TIDS(coeff, new_tids)
+        if isinstance(tmul, TensExpr):
+            tmul._matrix_behavior_kinds = matrix_behavior_kinds
+
         return tmul
 
     def __rmul__(self, other):
         other = sympify(other)
         coeff = other*self._coeff
         tmul = TensMul.from_TIDS(coeff, self._tids)
+        tmul._matrix_behavior_kinds = self._matrix_behavior_kinds
         return tmul
 
     def __div__(self, other):
@@ -3533,6 +3075,7 @@ def __div__(self, other):
             raise ValueError('cannot divide by a tensor')
         coeff = self._coeff/other
         tmul = TensMul.from_TIDS(coeff, self._tids, is_canon_bp=self._is_canon_bp)
+        tmul._matrix_behavior_kinds = self._matrix_behavior_kinds
         return tmul
 
     def __rdiv__(self, other):
@@ -3550,7 +3093,7 @@ def sorted_components(self):
         calling the corresponding method in a ``TIDS`` object.
         """
         new_tids, sign = self._tids.sorted_components()
-        coeff = -self.coeff if sign == -1 else self.coeff
+        coeff = -self._coeff if sign == -1 else self._coeff
         t = TensMul.from_TIDS(coeff, new_tids)
         return t
 
@@ -3560,7 +3103,13 @@ def perm2tensor(self, g, canon_bp=False):
 
         For further details, see the method in ``TIDS`` with the same name.
         """
-        return perm2tensor(self, g, canon_bp)
+        new_tids = self._tids.perm2tensor(g, canon_bp)
+        coeff = self._coeff
+        if g[-1] != len(g) - 1:
+            coeff = -coeff
+        res = TensMul.from_TIDS(coeff, new_tids, is_canon_bp=canon_bp)
+        res._matrix_behavior_kinds = self._matrix_behavior_kinds
+        return res
 
     def canon_bp(self):
         """
@@ -3591,6 +3140,7 @@ def canon_bp(self):
         if can == 0:
             return S.Zero
         tmul = t.perm2tensor(can, True)
+        tmul._matrix_behavior_kinds = self._matrix_behavior_kinds
         return tmul
 
     def contract_delta(self, delta):
@@ -3625,12 +3175,173 @@ def contract_metric(self, g):
         >>> t.contract_metric(g).canon_bp()
         p(L_0)*q(-L_0)
         """
-        tids, sign = get_tids(self).contract_metric(g)
-        res = TensMul.from_TIDS(sign*self.coeff, tids)
+        components = self.components
+        antisym = g.index_types[0].metric_antisym
+        #if not any(x == g for x in components):
+        #    return self
+        # list of positions of the metric ``g``
+        gpos = [i for i, x in enumerate(components) if x == g]
+        if not gpos:
+            return self
+        coeff = self._coeff
+        tids = self._tids
+        dum = tids.dum[:]
+        free = tids.free[:]
+        elim = set()
+        for gposx in gpos:
+            if gposx in elim:
+                continue
+            free1 = [x for x in free if x[-1] == gposx]
+            dum1 = [x for x in dum if x[-2] == gposx or x[-1] == gposx]
+            if not dum1:
+                continue
+            elim.add(gposx)
+            if len(dum1) == 2:
+                if not antisym:
+                    dum10, dum11 = dum1
+                    if dum10[3] == gposx:
+                        # the index with pos p0 and component c0 is contravariant
+                        c0 = dum10[2]
+                        p0 = dum10[0]
+                    else:
+                        # the index with pos p0 and component c0 is covariant
+                        c0 = dum10[3]
+                        p0 = dum10[1]
+                    if dum11[3] == gposx:
+                        # the index with pos p1 and component c1 is contravariant
+                        c1 = dum11[2]
+                        p1 = dum11[0]
+                    else:
+                        # the index with pos p1 and component c1 is covariant
+                        c1 = dum11[3]
+                        p1 = dum11[1]
+                    dum.append((p0, p1, c0, c1))
+                else:
+                    dum10, dum11 = dum1
+                    # change the sign to bring the indices of the metric to contravariant
+                    # form; change the sign if dum10 has the metric index in position 0
+                    if dum10[3] == gposx:
+                        # the index with pos p0 and component c0 is contravariant
+                        c0 = dum10[2]
+                        p0 = dum10[0]
+                        if dum10[1] == 1:
+                            coeff = -coeff
+                    else:
+                        # the index with pos p0 and component c0 is covariant
+                        c0 = dum10[3]
+                        p0 = dum10[1]
+                        if dum10[0] == 0:
+                            coeff = -coeff
+                    if dum11[3] == gposx:
+                        # the index with pos p1 and component c1 is contravariant
+                        c1 = dum11[2]
+                        p1 = dum11[0]
+                        coeff = -coeff
+                    else:
+                        # the index with pos p1 and component c1 is covariant
+                        c1 = dum11[3]
+                        p1 = dum11[1]
+                    dum.append((p0, p1, c0, c1))
+
+            elif len(dum1) == 1:
+                if not antisym:
+                    dp0, dp1, dc0, dc1 = dum1[0]
+                    if dc0 == dc1:
+                        # g(i, -i)
+                        typ = g.index_types[0]
+                        if typ._dim is None:
+                            raise ValueError('dimension not assigned')
+                        coeff = coeff*typ._dim
+
+                    else:
+                        # g(i0, i1)*p(-i1)
+                        if dc0 == gposx:
+                            p1 = dp1
+                            c1 = dc1
+                        else:
+                            p1 = dp0
+                            c1 = dc0
+                        ind, p, c = free1[0]
+                        free.append((ind, p1, c1))
+                else:
+                    dp0, dp1, dc0, dc1 = dum1[0]
+                    if dc0 == dc1:
+                        # g(i, -i)
+                        typ = g.index_types[0]
+                        if typ._dim is None:
+                            raise ValueError('dimension not assigned')
+                        coeff = coeff*typ._dim
+
+                        if dp0 < dp1:
+                            # g(i, -i) = -D with antisymmetric metric
+                            coeff = -coeff
+                    else:
+                        # g(i0, i1)*p(-i1)
+                        if dc0 == gposx:
+                            p1 = dp1
+                            c1 = dc1
+                            if dp0 == 0:
+                                coeff = -coeff
+                        else:
+                            p1 = dp0
+                            c1 = dc0
+                        ind, p, c = free1[0]
+                        free.append((ind, p1, c1))
+            dum = [x for x in dum if x not in dum1]
+            free = [x for x in free if x not in free1]
+
+        shift = 0
+        shifts = [0]*len(components)
+        for i in range(len(components)):
+            if i in elim:
+                shift += 1
+                continue
+            shifts[i] = shift
+        free = [(ind, p, c - shifts[c]) for (ind, p, c) in free if c not in elim]
+        dum = [(p0, p1, c0 - shifts[c0], c1 - shifts[c1]) for  i, (p0, p1, c0, c1) in enumerate(dum) if c0 not in elim and c1 not in elim]
+        components = [c for i, c in enumerate(components) if i not in elim]
+        tids = TIDS(components, free, dum)
+        res = TensMul.from_TIDS(coeff, tids)
         return res
 
     def substitute_indices(self, *index_tuples):
-        return substitute_indices(self, *index_tuples)
+        """
+        Return a tensor with free indices substituted according to ``index_tuples``
+
+        ``index_types`` list of tuples ``(old_index, new_index)``
+
+        Note: this method will neither raise or lower the indices, it will just replace their symbol.
+
+        Examples
+        ========
+
+        >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead
+        >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
+        >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz)
+        >>> A, B = tensorhead('A,B', [Lorentz]*2, [[1]*2])
+        >>> t = A(i, k)*B(-k, -j); t
+        A(i, L_0)*B(-L_0, -j)
+        >>> t.substitute_indices((i,j), (j, k))
+        A(j, L_0)*B(-L_0, -k)
+        """
+        free = self.free
+        free1 = []
+        for j, ipos, cpos in free:
+            for i, v in index_tuples:
+                if i._name == j._name and i._tensortype == j._tensortype:
+                    if i._is_up == j._is_up:
+                        free1.append((v, ipos, cpos))
+                    else:
+                        free1.append((-v, ipos, cpos))
+                    break
+            else:
+                free1.append((j, ipos, cpos))
+
+        t = TensMul.from_data(self._coeff, self.components, free1, self.dum)
+
+        # object is rebuilt in order to make sure that all contracted indices get recognized as dummies.
+        # t2 = TensMul(*t.args)
+        return t
 
     def fun_eval(self, *index_tuples):
         """
@@ -3660,10 +3371,10 @@ def fun_eval(self, *index_tuples):
                     break
             else:
                 free1.append((j, ipos, cpos))
-        return TensMul.from_data(self.coeff, self.components, free1, self.dum)
+        return TensMul.from_data(self._coeff, self.components, free1, self.dum)
 
     def __call__(self, *indices):
-        """Returns tensor product with ordered free indices replaced by ``indices``
+        """Returns tensor with ordered free indices replaced by ``indices``
 
         Examples
         ========
@@ -3694,19 +3405,26 @@ def __call__(self, *indices):
         return t
 
     def _print(self):
-        args = self.args
-        get_str = lambda arg: str(arg) if arg.is_Atom or isinstance(arg, TensExpr) else ("(%s)" % str(arg))
-
-        if not args:
-            # no arguments is equivalent to "1", i.e. TensMul().
-            # If tensors are constructed correctly, this should never occur.
-            return "1"
-        if self.coeff == S.NegativeOne:
-            # expressions like "-A(a)"
-            return "-"+"*".join([get_str(arg) for arg in args[1:]])
-
-        # prints expressions like "A(a)", "3*A(a)", "(1+x)*A(a)"
-        return "*".join([get_str(arg) for arg in self.args])
+        if len(self.components) == 0:
+            return str(self._coeff)
+        indices = [str(ind) for ind in self.get_indices()]
+        pos = 0
+        a = []
+        for t in self.components:
+            if t._rank > 0:
+                a.append('%s(%s)' % (t.name, ', '.join(indices[pos:pos + t._rank])))
+            else:
+                a.append('%s' % t.name)
+            pos += t._rank
+        res = '*'. join(a)
+        if self._coeff == S.One:
+            return res
+        elif self._coeff == -S.One:
+            return '-%s' % res
+        if self._coeff.is_Atom:
+            return '%s*%s' % (self._coeff, res)
+        else:
+            return '(%s)*%s' %(self._coeff, res)
 
     @property
     def data(self):
@@ -3717,11 +3435,13 @@ def data(self):
 
     @data.setter
     def data(self, data):
-        raise ValueError("Not possible to set component data to a tensor expression")
+        # TODO: check data compatibility with properties of tensor.
+        _tensor_data_substitution_dict[self] = data
 
     @data.deleter
     def data(self):
-        raise ValueError("Not possible to delete component data to a tensor expression")
+        if self in _tensor_data_substitution_dict:
+            del _tensor_data_substitution_dict[self]
 
     def __iter__(self):
         if self.data is None:
@@ -3782,7 +3502,7 @@ def riemann_cyclic(t2):
     >>> riemann_cyclic(t)
     0
     """
-    if isinstance(t2, (TensMul, Tensor)):
+    if isinstance(t2, TensMul):
         args = [t2]
     else:
         args = t2.args
@@ -3807,7 +3527,6 @@ def _join_lines(a):
         while i < len(a):
             x = a[i]
             xend = x[-1]
-            xstart = x[0]
             hit = True
             while hit:
                 hit = False
@@ -3819,27 +3538,6 @@ def _join_lines(a):
                         x.extend(a[j][1:])
                         xend = x[-1]
                         a.pop(j)
-                        continue
-                    if a[j][0] == xstart:
-                        hit = True
-                        a[i] = reversed(a[j][1:]) + x
-                        x = a[i]
-                        xstart = a[i][0]
-                        a.pop(j)
-                        continue
-                    if a[j][-1] == xend:
-                        hit = True
-                        x.extend(reversed(a[j][:-1]))
-                        xend = x[-1]
-                        a.pop(j)
-                        continue
-                    if a[j][-1] == xstart:
-                        hit = True
-                        a[i] = a[j][:-1] + x
-                        x = a[i]
-                        xstart = x[0]
-                        a.pop(j)
-                        continue
             i += 1
         return a
 
@@ -3910,80 +3608,3 @@ def _join_lines(a):
     rest = [x for x in range(len(components)) if x not in rest]
 
     return lines, traces, rest
-
-def get_indices(t):
-    if not isinstance(t, TensExpr):
-        return ()
-    return t.get_indices()
-
-def get_tids(t):
-    if isinstance(t, TensExpr):
-        return t._tids
-    return TIDS([], [], [])
-
-def get_coeff(t):
-    if isinstance(t, Tensor):
-        return S.One
-    if isinstance(t, TensMul):
-        return t.coeff
-    if isinstance(t, TensExpr):
-        raise ValueError("no coefficient associated to this tensor expression")
-    return t
-
-def contract_metric(t, g):
-    if isinstance(t, TensExpr):
-        return t.contract_metric(g)
-    return t
-
-def perm2tensor(t, g, canon_bp=False):
-    """
-    Returns the tensor corresponding to the permutation ``g``
-
-    For further details, see the method in ``TIDS`` with the same name.
-    """
-    if not isinstance(t, TensExpr):
-        return t
-    new_tids = get_tids(t).perm2tensor(g, canon_bp)
-    coeff = get_coeff(t)
-    if g[-1] != len(g) - 1:
-        coeff = -coeff
-    res = TensMul.from_TIDS(coeff, new_tids, is_canon_bp=canon_bp)
-    return res
-
-def substitute_indices(t, *index_tuples):
-    """
-    Return a tensor with free indices substituted according to ``index_tuples``
-
-    ``index_types`` list of tuples ``(old_index, new_index)``
-
-    Note: this method will neither raise or lower the indices, it will just replace their symbol.
-
-    Examples
-    ========
-
-    >>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, tensorhead
-    >>> Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
-    >>> i, j, k, l = tensor_indices('i,j,k,l', Lorentz)
-    >>> A, B = tensorhead('A,B', [Lorentz]*2, [[1]*2])
-    >>> t = A(i, k)*B(-k, -j); t
-    A(i, L_0)*B(-L_0, -j)
-    >>> t.substitute_indices((i,j), (j, k))
-    A(j, L_0)*B(-L_0, -k)
-    """
-    if not isinstance(t, TensExpr):
-        return t
-    free = t.free
-    free1 = []
-    for j, ipos, cpos in free:
-        for i, v in index_tuples:
-            if i._name == j._name and i._tensortype == j._tensortype:
-                if i._is_up == j._is_up:
-                    free1.append((v, ipos, cpos))
-                else:
-                    free1.append((-v, ipos, cpos))
-                break
-        else:
-            free1.append((j, ipos, cpos))
-
-    t = TensMul.from_data(t.coeff, t.components, free1, t.dum)
-    return t


Test patch:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index c1bdff1441..93d35e1939 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -3039,17 +3039,6 @@ def test_sympy__tensor__tensor__TensAdd():
     assert _test_args(TensAdd(t1, t2))
 
 
-def test_sympy__tensor__tensor__Tensor():
-    from sympy.core import S
-    from sympy.tensor.tensor import TensorIndexType, TensorSymmetry, TensorType, get_symmetric_group_sgs, tensor_indices, TensMul, TIDS
-    Lorentz = TensorIndexType('Lorentz', dummy_fmt='L')
-    a, b = tensor_indices('a,b', Lorentz)
-    sym = TensorSymmetry(get_symmetric_group_sgs(1))
-    S1 = TensorType([Lorentz], sym)
-    p = S1('p')
-    assert _test_args(p(a))
-
-
 def test_sympy__tensor__tensor__TensMul():
     from sympy.core import S
     from sympy.tensor.tensor import TensorIndexType, TensorSymmetry, TensorType, get_symmetric_group_sgs, tensor_indices, TensMul, TIDS
@@ -3058,8 +3047,8 @@ def test_sympy__tensor__tensor__TensMul():
     sym = TensorSymmetry(get_symmetric_group_sgs(1))
     S1 = TensorType([Lorentz], sym)
     p = S1('p')
-    q = S1('q')
-    assert _test_args(3*p(a)*q(b))
+    free, dum = TIDS.free_dum_from_indices(a)
+    assert _test_args(TensMul.from_data(S.One, [p], free, dum))
 
 
 def test_as_coeff_add():
diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py
index 118bd92ddc..ee1be6bdfc 100644
--- a/sympy/core/tests/test_expr.py
+++ b/sympy/core/tests/test_expr.py
@@ -6,7 +6,7 @@
     Piecewise, Mul, Pow, nsimplify, ratsimp, trigsimp, radsimp, powsimp,
     simplify, together, collect, factorial, apart, combsimp, factor, refine,
     cancel, Tuple, default_sort_key, DiracDelta, gamma, Dummy, Sum, E,
-    exp_polar, Lambda, expand, diff, O, Heaviside, Si, Max)
+    exp_polar, Lambda, expand, diff, O, Heaviside, Si)
 from sympy.core.function import AppliedUndef
 from sympy.physics.secondquant import FockState
 from sympy.physics.units import meter
@@ -1531,9 +1531,6 @@ def test_random():
     assert posify(x)[0]._random() is not None
     assert lucas(n)._random(2, -2, 0, -1, 1) is None
 
-    # issue 8662
-    assert Piecewise((Max(x, y), z))._random() is None
-
 
 def test_round():
     from sympy.abc import x
diff --git a/sympy/core/tests/test_function.py b/sympy/core/tests/test_function.py
index 304cf5f421..0645b8f890 100644
--- a/sympy/core/tests/test_function.py
+++ b/sympy/core/tests/test_function.py
@@ -183,7 +183,7 @@ def test_Lambda_equality():
 
 def test_Subs():
     assert Subs(x, x, 0) == Subs(y, y, 0)
-    assert Subs(x, x, 0).subs(x, 1) == Subs(x, x, 0)
+    assert Subs(x, x, 0).subs(x, 1) == Subs(x, x, 1)
     assert Subs(y, x, 0).subs(y, 1) == Subs(1, x, 0)
     assert Subs(f(x), x, 0).doit() == f(0)
     assert Subs(f(x**2), x**2, 0).doit() == f(0)
@@ -206,7 +206,7 @@ def test_Subs():
     assert Subs(f(x)*y, (x, y), (0, 1)) == Subs(f(y)*x, (y, x), (0, 1))
     assert Subs(f(x)*y, (x, y), (1, 1)) == Subs(f(y)*x, (x, y), (1, 1))
 
-    assert Subs(f(x), x, 0).subs(x, 1).doit() == f(0)
+    assert Subs(f(x), x, 0).subs(x, 1).doit() == f(1)
     assert Subs(f(x), x, y).subs(y, 0) == Subs(f(x), x, 0)
     assert Subs(y*f(x), x, y).subs(y, 2) == Subs(2*f(x), x, 2)
     assert (2 * Subs(f(x), x, 0)).subs(Subs(f(x), x, 0), y) == 2*y
@@ -235,9 +235,6 @@ def test_Subs():
         Subs(f(x)*cos(y) + z, (x, y), (0, pi/3)).evalf(2) == \
         z + Rational('1/2').n(2)*f(0)
 
-    assert f(x).diff(x).subs(x, 0).subs(x, y) == f(x).diff(x).subs(x, 0)
-    assert (x*f(x).diff(x).subs(x, 0)).subs(x, y) == y*f(x).diff(x).subs(x, 0)
-
 
 @XFAIL
 def test_Subs2():
diff --git a/sympy/functions/combinatorial/tests/test_comb_factorials.py b/sympy/functions/combinatorial/tests/test_comb_factorials.py
index 727c5fc999..9761dd2d6c 100644
--- a/sympy/functions/combinatorial/tests/test_comb_factorials.py
+++ b/sympy/functions/combinatorial/tests/test_comb_factorials.py
@@ -1,6 +1,6 @@
 from sympy import (Symbol, symbols, factorial, factorial2, binomial,
                    rf, ff, gamma, polygamma, EulerGamma, O, pi, nan,
-                   oo, zoo, simplify, expand_func, C, S, Product)
+                   oo, zoo, simplify, expand_func)
 from sympy.functions.combinatorial.factorials import subfactorial
 from sympy.utilities.pytest import XFAIL, raises
 
@@ -146,10 +146,8 @@ def test_factorial_series():
 
 def test_factorial_rewrite():
     n = Symbol('n', integer=True)
-    k = Symbol('k', integer=True, nonnegative=True)
 
     assert factorial(n).rewrite(gamma) == gamma(n + 1)
-    assert str(factorial(k).rewrite(Product)) == 'Product(_i, (_i, 1, k))'
 
 
 def test_factorial2():
@@ -264,11 +262,8 @@ def test_factorial_simplify_fail():
 def test_subfactorial():
     assert all(subfactorial(i) == ans for i, ans in enumerate(
         [1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496]))
-    assert subfactorial(oo) == oo
-
-    x = Symbol('x')
-    assert subfactorial(x).rewrite(C.uppergamma) == \
-        C.uppergamma(x + 1, -1)/S.Exp1
+    raises(ValueError, lambda: subfactorial(0.1))
+    raises(ValueError, lambda: subfactorial(-2))
 
     tt = Symbol('tt', integer=True, nonnegative=True)
     tf = Symbol('tf', integer=True, nonnegative=False)
@@ -280,11 +275,11 @@ def test_subfactorial():
     nf = Symbol('nf', nonnegative=False)
     nn = Symbol('nf')
     assert subfactorial(tt).is_integer
-    assert subfactorial(tf).is_integer is None
+    assert subfactorial(tf).is_integer is False
     assert subfactorial(tn).is_integer is None
-    assert subfactorial(ft).is_integer is None
-    assert subfactorial(ff).is_integer is None
-    assert subfactorial(fn).is_integer is None
+    assert subfactorial(ft).is_integer is False
+    assert subfactorial(ff).is_integer is False
+    assert subfactorial(fn).is_integer is False
     assert subfactorial(nt).is_integer is None
-    assert subfactorial(nf).is_integer is None
+    assert subfactorial(nf).is_integer is False
     assert subfactorial(nn).is_integer is None
diff --git a/sympy/logic/tests/test_boolalg.py b/sympy/logic/tests/test_boolalg.py
index 64be048911..6e9fb597c4 100644
--- a/sympy/logic/tests/test_boolalg.py
+++ b/sympy/logic/tests/test_boolalg.py
@@ -667,3 +667,11 @@ def test_all_or_nothing():
 def test_canonical_atoms():
     assert true.canonical == true
     assert false.canonical == false
+
+
+def test_issue_8777():
+    x = symbols('x')
+    assert And(x > 2, x < oo).as_set() == Interval(2, oo, left_open=True)
+    assert And(x >= 1, x < oo).as_set() == Interval(1, oo)
+    assert (x < oo).as_set() == Interval(-oo, oo)
+    assert (x > -oo).as_set() == Interval(-oo, oo)
diff --git a/sympy/physics/hep/tests/test_gamma_matrices.py b/sympy/physics/hep/tests/test_gamma_matrices.py
index d0c018761c..dcbcc81c44 100644
--- a/sympy/physics/hep/tests/test_gamma_matrices.py
+++ b/sympy/physics/hep/tests/test_gamma_matrices.py
@@ -435,18 +435,18 @@ def test_get_lines():
         G(i4,s4,-s5)*G(i5,s6,-s7)*G(i10,s11,-s12)*G(i6,s7,-s8)*G(i7,s8,-s9)*\
         G(i8,s9,-s6)*G(i9,s10,-s0)
     r = get_lines(t, DiracSpinorIndex)
-    assert r == ([[7, 1], [0, 3, 4, 5]], [[6, 8, 9, 10], [2, 11]], [])
+    assert r == ([[0, 3, 4, 5], [7, 1]], [[6, 8, 9, 10], [2, 11]], [])
     t = G(i4,s4,-s5)*G(i5,s6,-s7)*G(i10,s11,-s12)*G(i6,s7,-s8)*G(i7,s8,-s9)*\
         G(i8,s9,-s6)*G(i9,s10,-s0)*\
         G(i1,s1,-s2)*G(i11,s12,-s13)*G(i0,s0,-s10)*G(i2,s2,-s3)*G(i3,s3,-s4)
     r = get_lines(t, DiracSpinorIndex)
-    assert r == ([[7, 10, 11, 0], [2, 8]], [[1, 3, 4, 5], [6, 9]], [])
+    assert r == ([[2, 8], [7, 10, 11, 0]], [[1, 3, 4, 5], [6, 9]], [])
     t = G(i8,s9,-s6)*G(i9,s10,-s0)*G(i4,s4,-s5)*G(i13,s14,-s15)*\
         G(i10,s11,-s12)*G(i1,s1,-s2)*G(i11,s12,-s13)*\
         G(i0,s0,-s10)*G(i6,s7,-s8)*G(i7,s8,-s9)*\
         G(i2,s2,-s3)*G(i12,s13,-s14)*G(i3,s3,-s4)*G(i5,s6,-s7)
     r = get_lines(t, DiracSpinorIndex)
-    assert r == ([[5, 10, 12, 2], [4, 6, 11, 3]], [[1, 7], [0, 13, 8, 9]], [])
+    assert r == ([[4, 6, 11, 3], [5, 10, 12, 2]], [[1, 7], [0, 13, 8, 9]], [])
 
 def test_simplify_lines():
     i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12 = tensor_indices('i0:13', G.LorentzIndex)
diff --git a/sympy/polys/tests/test_rootoftools.py b/sympy/polys/tests/test_rootoftools.py
index 60c7bedd99..40e2fa7e61 100644
--- a/sympy/polys/tests/test_rootoftools.py
+++ b/sympy/polys/tests/test_rootoftools.py
@@ -228,9 +228,6 @@ def test_RootOf_evalf():
     ans = [w.n(2) for w in solve(x**3 - x - 4)]
     assert RootOf(exp(x)**3 - exp(x) - 4, 0).n(2) in ans
 
-    # make sure verification is used in case a max/min traps the "root"
-    assert str(RootOf(4*x**5 + 16*x**3 + 12*x**2 + 7, 0).n(3)) == '-0.976'
-
 
 def test_RootOf_evalf_caching_bug():
     r = RootOf(x**5 - 5*x + 12, 1)
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 07c795444c..5b030ffaf2 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -3,7 +3,7 @@
 from sympy import (
     limit, exp, oo, log, sqrt, Limit, sin, floor, cos, ceiling,
     atan, gamma, Symbol, S, pi, Integral, cot, Rational, I, zoo,
-    tan, cot, integrate, Sum, sign, Function, subfactorial)
+    tan, cot, integrate, Sum, sign, Function)
 
 from sympy.series.limits import heuristics
 from sympy.series.order import Order
@@ -433,7 +433,3 @@ def test_issue_4503():
     dx = Symbol('dx')
     assert limit((sqrt(1 + exp(x + dx)) - sqrt(1 + exp(x)))/dx, dx, 0) == \
         exp(x)/(2*sqrt(exp(x) + 1))
-
-
-def test_issue_8730():
-    assert limit(subfactorial(x), x, oo) == oo
diff --git a/sympy/solvers/tests/test_inequalities.py b/sympy/solvers/tests/test_inequalities.py
index d4bcc8ea3b..e2902de685 100644
--- a/sympy/solvers/tests/test_inequalities.py
+++ b/sympy/solvers/tests/test_inequalities.py
@@ -304,3 +304,11 @@ def test_issue_8545():
     assert reduce_abs_inequality(eq, '<', x) == ans
     eq = 1 - x - sqrt((1 - x)**2)
     assert reduce_inequalities(eq < 0) == ans
+
+
+def test_issue_8783():
+    x = Symbol('x')
+    assert solve(x > -oo) == And(-oo < x, x < oo)
+    assert solve(x < oo) == And(-oo < x, x < oo)
+    assert solve(x > oo) == False
+    assert solve(x < -oo) == False
diff --git a/sympy/tensor/tests/test_tensor.py b/sympy/tensor/tests/test_tensor.py
index d893c90797..bfc0c5196d 100644
--- a/sympy/tensor/tests/test_tensor.py
+++ b/sympy/tensor/tests/test_tensor.py
@@ -655,15 +655,6 @@ def test_mul():
     assert sorted(t.free) == [(a, 0, 0), (b, 1, 0)]
     assert t.components == [A]
 
-    ts = A(a, b)
-    assert str(ts) == 'A(a, b)'
-    assert ts.types == [Lorentz]
-    assert ts.rank == 2
-    assert ts.dum == []
-    assert ts.coeff == 1
-    assert sorted(ts.free) == [(a, 0, 0), (b, 1, 0)]
-    assert ts.components == [A]
-
     t = A(-b, a)*B(-a, c)*A(-c, d)
     t1 = tensor_mul(*t.split())
     assert t == t(-b, d)
@@ -1196,7 +1187,7 @@ def test_hidden_indices_for_matrix_multiplication():
     B1 = B(m1)
 
     assert _is_equal((B1*A0*B0), B(m1, s0)*A(m0, -s0, s1)*B(-m0, -s1))
-    assert _is_equal((B0*A0), B(-m0, s0)*A(m0, -s0, S.auto_left))
+    assert _is_equal((B0*A0), B(-m0, s0)*A(m0, -s0, -S.auto_right))
     assert _is_equal((A0*B0), A(m0, S.auto_left, s0)*B(-m0, -s0))
 
     C = tensorhead('C', [L, L], [[1]*2])