Repository: sympy/sympy

Issue:
Tensor module: automatic relabelling of indices causes conflicts when using subs or xreplace on TensMul
```
>>> from sympy.tensor.tensor import TensorIndexType, tensor_indices, TensorHead
>>> R3 = TensorIndexType('R3', dim=3)
>>> p, q, r = tensor_indices("p q r", R3)
>>> K = TensorHead("K", [R3])
>>> V = TensorHead("V", [R3])
>>> ( K(p)*V(r)*K(-p) ).xreplace({V(r): K(q)*K(-q)})
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.10/site-packages/sympy/core/basic.py", line 1195, in xreplace
    value, _ = self._xreplace(rule)
  File "/usr/lib/python3.10/site-packages/sympy/core/basic.py", line 1217, in _xreplace
    return self.func(*args), True
  File "/usr/lib/python3.10/site-packages/sympy/tensor/tensor.py", line 3260, in __new__
    args, indices, free, dum = TensMul._tensMul_contract_indices(args, replace_indices=False)
  File "/usr/lib/python3.10/site-packages/sympy/tensor/tensor.py", line 3344, in _tensMul_contract_indices
    indices, free, free_names, dummy_data = TensMul._indices_to_free_dum(args_indices)
  File "/usr/lib/python3.10/site-packages/sympy/tensor/tensor.py", line 3320, in _indices_to_free_dum
    raise ValueError("Repeated index: %s" % index)
ValueError: Repeated index: R_0
```

The same error is also triggered by `subs`.

This happens because the tensor module automatically relabels dummy indices when a TensMul is created using the `*` operator (its `__mul__` method calls `doit` on the returned expression):
```
>>> K(q)*K(-q)
K(R_0)*K(-R_0)
>>> K(p)*V(r)*K(-p)
K(R_0)*V(r)*K(-R_0)
```
When a new TensMul (with the requested substitutions performed) is created in `subs` or `xreplace`, what happens is then that the internally created index `R_0` repeats, which is not allowed in a tensor multiplication when using the Einstein summation convention.

#### Software versions
Sympy: built from commit 5ba43252efd64bcda5aa0f6bc27debdb6b0fefc2
Python: 3.10.8

#### Attempted fix
I attempted to fix this in https://github.com/sympy/sympy/pull/24333

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index 12ebd404ed..419d5a2fec 100644
--- a/.mailmap
+++ b/.mailmap
@@ -781,6 +781,7 @@ Kirill Smelkov <kirr@landau.phys.spbu.ru> convert-repo <devnull@localhost>
 Kirill Smelkov <kirr@landau.phys.spbu.ru> kirill.smelkov <devnull@localhost>
 Kirtan Mali <kirtanmali555@gmail.com> KIRTAN  MALI <43683545+kmm555@users.noreply.github.com>
 Kirtan Mali <kirtanmali555@gmail.com> kmm555 <kirtanmali555@gmail.com>
+Kishore Gopalakrishnan <kishore96@gmail.com>
 Kiyohito Yamazaki <kyamaz@openql.org>
 Klaus Rettinghaus <klaus.rettinghaus@enote.com>
 Konrad Meyer <konrad.meyer@gmail.com>
diff --git a/sympy/tensor/tensor.py b/sympy/tensor/tensor.py
index d0decae116..a1c822ec26 100644
--- a/sympy/tensor/tensor.py
+++ b/sympy/tensor/tensor.py
@@ -3254,6 +3254,27 @@ def __new__(cls, *args, **kw_args):
         is_canon_bp = kw_args.get('is_canon_bp', False)
         args = list(map(_sympify, args))
 
+        """
+        If the internal dummy indices in one arg conflict with the free indices
+        of the remaining args, we need to rename those internal dummy indices.
+        """
+        free = [get_free_indices(arg) for arg in args]
+        free = set(itertools.chain(*free)) #flatten free
+        newargs = []
+        for arg in args:
+            dum_this = set(get_dummy_indices(arg))
+            dum_other = [get_dummy_indices(a) for a in newargs]
+            dum_other = set(itertools.chain(*dum_other)) #flatten dum_other
+            free_this = set(get_free_indices(arg))
+            if len(dum_this.intersection(free)) > 0:
+                exclude = free_this.union(free, dum_other)
+                newarg = TensMul._dedupe_indices(arg, exclude)
+            else:
+                newarg = arg
+            newargs.append(newarg)
+
+        args = newargs
+
         # Flatten:
         args = [i for arg in args for i in (arg.args if isinstance(arg, (TensMul, Mul)) else [arg])]
 
@@ -3401,6 +3422,17 @@ def doit(self, **hints):
         deep = hints.get('deep', True)
         if deep:
             args = [arg.doit(**hints) for arg in self.args]
+
+            """
+            There may now be conflicts between dummy indices of different args
+            (each arg's doit method does not have any information about which
+            dummy indices are already used in the other args), so we
+            deduplicate them.
+            """
+            rule = dict(zip(self.args, args))
+            rule = self._dedupe_indices_in_rule(rule)
+            args = [rule[a] for a in self.args]
+
         else:
             args = self.args
 
@@ -3988,6 +4020,73 @@ def __iter__(self):
                 raise ValueError("No iteration on abstract tensors")
             return self.data.__iter__()
 
+    @staticmethod
+    def _dedupe_indices(new, exclude):
+        """
+        exclude: set
+        new: TensExpr
+
+        If ``new`` has any dummy indices that are in ``exclude``, return a version
+        of new with those indices replaced. If no replacements are needed,
+        return None
+
+        """
+        exclude = set(exclude)
+        dums_new = set(get_dummy_indices(new))
+        free_new = set(get_free_indices(new))
+
+        conflicts = dums_new.intersection(exclude)
+        if len(conflicts) == 0:
+            return None
+
+        """
+        ``exclude_for_gen`` is to be passed to ``_IndexStructure._get_generator_for_dummy_indices()``.
+        Since the latter does not use the index position for anything, we just
+        set it as ``None`` here.
+        """
+        exclude.update(dums_new)
+        exclude.update(free_new)
+        exclude_for_gen = [(i, None) for i in exclude]
+        gen = _IndexStructure._get_generator_for_dummy_indices(exclude_for_gen)
+        repl = {}
+        for d in conflicts:
+            if -d in repl.keys():
+                continue
+            newname = gen(d.tensor_index_type)
+            new_d = d.func(newname, *d.args[1:])
+            repl[d] = new_d
+            repl[-d] = -new_d
+
+        if len(repl) == 0:
+            return None
+
+        new_renamed = new._replace_indices(repl)
+        return new_renamed
+
+    def _dedupe_indices_in_rule(self, rule):
+        """
+        rule: dict
+
+        This applies TensMul._dedupe_indices on all values of rule.
+
+        """
+        index_rules = {k:v for k,v in rule.items() if isinstance(k, TensorIndex)}
+        other_rules = {k:v for k,v in rule.items() if k not in index_rules.keys()}
+        exclude = set(self.get_indices())
+
+        newrule = {}
+        newrule.update(index_rules)
+        exclude.update(index_rules.keys())
+        exclude.update(index_rules.values())
+        for old, new in other_rules.items():
+            new_renamed = TensMul._dedupe_indices(new, exclude)
+            if old == new or new_renamed is None:
+                newrule[old] = new
+            else:
+                newrule[old] = new_renamed
+                exclude.update(get_indices(new_renamed))
+        return newrule
+
     def _eval_rewrite_as_Indexed(self, *args):
         from sympy.concrete.summations import Sum
         index_symbols = [i.args[0] for i in self.get_indices()]
@@ -4305,6 +4404,12 @@ def get_indices(t):
         return ()
     return t.get_indices()
 
+def get_dummy_indices(t):
+    if not isinstance(t, TensExpr):
+        return ()
+    inds = t.get_indices()
+    free = t.get_free_indices()
+    return [i for i in inds if i not in free]
 
 def get_index_structure(t):
     if isinstance(t, TensExpr):


Test patch:
diff --git a/sympy/tensor/tests/test_tensor.py b/sympy/tensor/tests/test_tensor.py
index 94d6bd7252..7ab076b3b0 100644
--- a/sympy/tensor/tests/test_tensor.py
+++ b/sympy/tensor/tests/test_tensor.py
@@ -1961,6 +1961,23 @@ def test_rewrite_tensor_to_Indexed():
     assert b2.rewrite(Indexed) == Sum(Indexed(Symbol("B"), L_1)*Indexed(Symbol("A"), L_0, L_0, i2, L_1), (L_0, 0, 3), (L_1, 0, 3))
 
 
+def test_TensMul_subs():
+    """
+    Test subs and xreplace in TensMul. See bug #24337
+    """
+    R3 = TensorIndexType('R3', dim=3)
+    p, q, r = tensor_indices("p q r", R3)
+    K = TensorHead("K", [R3])
+    V = TensorHead("V", [R3])
+    A = TensorHead("A", [R3, R3])
+    C0 = TensorIndex(R3.dummy_name + "_0", R3, True)
+
+    assert ( K(p)*V(r)*K(-p) ).subs({V(r): K(q)*K(-q)}) == K(p)*K(q)*K(-q)*K(-p)
+    assert ( K(p)*V(r)*K(-p) ).xreplace({V(r): K(q)*K(-q)}) == K(p)*K(q)*K(-q)*K(-p)
+    assert ( K(p)*V(r) ).xreplace({p: C0, V(r): K(q)*K(-q)}) == K(C0)*K(q)*K(-q)
+    assert ( K(p)*A(q,-q)*K(-p) ).doit() == K(p)*A(q,-q)*K(-p)
+
+
 def test_tensorsymmetry():
     with warns_deprecated_sympy():
         tensorsymmetry([1]*2)