code: def simplify_gpgp(ex, sort=True):
      """
      simplify products ``G(i)*p(-i)*G(j)*p(-j) -> p(i)*p(-i)``
  
      Examples
      ========
  
      >>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
          LorentzIndex, simplify_gpgp
      >>> from sympy.tensor.tensor import tensor_indices, tensor_heads
      >>> p, q = tensor_heads('p, q', [LorentzIndex])
      >>> i0,i1,i2,i3,i4,i5 = tensor_indices('i0:6', LorentzIndex)
      >>> ps = p(i0)*G(-i0)
      >>> qs = q(i0)*G(-i0)
      >>> simplify_gpgp(ps*qs*qs)
      GammaMatrix(-L_0)*p(L_0)*q(L_1)*q(-L_1)
      """
  
      def _simplify_gpgp(ex):
          components = ex.components
          a = []
          comp_map = []
          for i, comp in enumerate(components):
              comp_map.extend([i] * comp.rank)
          dum = [
              (i[0], i[1], comp_map[i[0]], comp_map[i[1]])
              for i in ex.dum
          ]
          for i in range(len(components)):
              if components[i] != GammaMatrix:
                  continue
              for dx in dum:
                  if dx[2] == i:
                      p_pos1 = dx[3]
                  elif dx[3] == i:
                      p_pos1 = dx[2]
                  else:
                      continue
                  comp1 = components[p_pos1]
                  if comp1.comm == 0 and comp1.rank == 1:
                      a.append((i, p_pos1))
          if not a:
              return ex
          elim = set()
          tv = []
          hit = True
          coeff = S.One
          ta = None
          while hit:
              hit = False
              for i, ai in enumerate(a[:-1]):
                  if ai[0] in elim:
                      continue
                  if ai[0] != a[i + 1][0] - 1:
                      continue
                  if components[ai[1]] != components[a[i + 1][1]]:
                      continue
                  elim.add(ai[0])
                  elim.add(ai[1])
                  elim.add(a[i + 1][0])
                  elim.add(a[i + 1][1])
                  if not ta:
                      ta = ex.split()
                      mu = TensorIndex("mu", LorentzIndex)
                  hit = True
                  if i == 0:
                      coeff = ex.coeff
                  tx = components[ai[1]](mu) * components[ai[1]](-mu)
                  if len(a) == 2:
                      tx *= 4  # eye(4)
                  tv.append(tx)
                  break
  
          if tv:
              a = [x for j, x in enumerate(ta) if j not in elim]
              a.extend(tv)
              t = tensor_mul(*a) * coeff
              # t = t.replace(lambda x: x.is_Matrix, lambda x: 1)
              return t
          else:
              return ex
  
      if sort:
          ex = ex.sorted_components()
      # this would be better off with pattern matching
      while 1:
          t = _simplify_gpgp(ex)
          if t != ex:
              ex = t
          else:
              return t
  
file_path: sympy/sympy/physics/hep/gamma_matrices.py
imports: from sympy.core.mul import Mul
  from sympy.core.singleton import S
  from sympy.matrices.dense import eye
  from sympy.matrices.expressions.trace import trace
  from sympy.tensor.tensor import TensorIndexType, TensorIndex, TensMul, TensAdd, tensor_mul, Tensor, TensorHead, TensorSymmetry
label_desc: Use image node_id 3 for calling a global function with example usage: simplify_gpgp(ex, sort) and returns: ex, t, ex, t