chore: 添加虚拟环境到仓库

- 添加 backend_service/venv 虚拟环境
- 包含所有Python依赖包
- 注意：虚拟环境约393MB，包含12655个文件

backend_service/venv/lib/python3.13/site-packages/sympy/physics/hep/gamma_matrices.py

@@ -0,0 +1,716 @@
+"""
+    Module to handle gamma matrices expressed as tensor objects.
+
+    Examples
+    ========
+
+    >>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex
+    >>> from sympy.tensor.tensor import tensor_indices
+    >>> i = tensor_indices('i', LorentzIndex)
+    >>> G(i)
+    GammaMatrix(i)
+
+    Note that there is already an instance of GammaMatrixHead in four dimensions:
+    GammaMatrix, which is simply declare as
+
+    >>> from sympy.physics.hep.gamma_matrices import GammaMatrix
+    >>> from sympy.tensor.tensor import tensor_indices
+    >>> i = tensor_indices('i', LorentzIndex)
+    >>> GammaMatrix(i)
+    GammaMatrix(i)
+
+    To access the metric tensor
+
+    >>> LorentzIndex.metric
+    metric(LorentzIndex,LorentzIndex)
+
+"""
+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
+
+
+# DiracSpinorIndex = TensorIndexType('DiracSpinorIndex', dim=4, dummy_name="S")
+
+
+LorentzIndex = TensorIndexType('LorentzIndex', dim=4, dummy_name="L")
+
+
+GammaMatrix = TensorHead("GammaMatrix", [LorentzIndex],
+                         TensorSymmetry.no_symmetry(1), comm=None)
+
+
+def extract_type_tens(expression, component):
+    """
+    Extract from a ``TensExpr`` all tensors with `component`.
+
+    Returns two tensor expressions:
+
+    * the first contains all ``Tensor`` of having `component`.
+    * the second contains all remaining.
+
+
+    """
+    if isinstance(expression, Tensor):
+        sp = [expression]
+    elif isinstance(expression, TensMul):
+        sp = expression.args
+    else:
+        raise ValueError('wrong type')
+
+    # 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 i.component == component:
+            new_expr *= i
+        else:
+            residual_expr *= i
+    return new_expr, residual_expr
+
+
+def simplify_gamma_expression(expression):
+    extracted_expr, residual_expr = extract_type_tens(expression, GammaMatrix)
+    res_expr = _simplify_single_line(extracted_expr)
+    return res_expr * residual_expr
+
+
+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
+
+
+def gamma_trace(t):
+    """
+    trace of a single line of gamma matrices
+
+    Examples
+    ========
+
+    >>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
+        gamma_trace, LorentzIndex
+    >>> 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)
+    >>> gamma_trace(G(i0)*G(i1))
+    4*metric(i0, i1)
+    >>> gamma_trace(ps*ps) - 4*p(i0)*p(-i0)
+    0
+    >>> gamma_trace(ps*qs + ps*ps) - 4*p(i0)*p(-i0) - 4*p(i0)*q(-i0)
+    0
+
+    """
+    if isinstance(t, TensAdd):
+        res = TensAdd(*[gamma_trace(x) for x in t.args])
+        return res
+    t = _simplify_single_line(t)
+    res = _trace_single_line(t)
+    return res
+
+
+def _simplify_single_line(expression):
+    """
+    Simplify single-line product of gamma matrices.
+
+    Examples
+    ========
+
+    >>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
+        LorentzIndex, _simplify_single_line
+    >>> from sympy.tensor.tensor import tensor_indices, TensorHead
+    >>> p = TensorHead('p', [LorentzIndex])
+    >>> i0,i1 = tensor_indices('i0:2', LorentzIndex)
+    >>> _simplify_single_line(G(i0)*G(i1)*p(-i1)*G(-i0)) + 2*G(i0)*p(-i0)
+    0
+
+    """
+    t1, t2 = extract_type_tens(expression, GammaMatrix)
+    if t1 != 1:
+        t1 = kahane_simplify(t1)
+    res = t1*t2
+    return res
+
+
+def _trace_single_line(t):
+    """
+    Evaluate the trace of a single gamma matrix line inside a ``TensExpr``.
+
+    Notes
+    =====
+
+    If there are ``DiracSpinorIndex.auto_left`` and ``DiracSpinorIndex.auto_right``
+    indices trace over them; otherwise traces are not implied (explain)
+
+
+    Examples
+    ========
+
+    >>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, \
+        LorentzIndex, _trace_single_line
+    >>> from sympy.tensor.tensor import tensor_indices, TensorHead
+    >>> p = TensorHead('p', [LorentzIndex])
+    >>> i0,i1,i2,i3,i4,i5 = tensor_indices('i0:6', LorentzIndex)
+    >>> _trace_single_line(G(i0)*G(i1))
+    4*metric(i0, i1)
+    >>> _trace_single_line(G(i0)*p(-i0)*G(i1)*p(-i1)) - 4*p(i0)*p(-i0)
+    0
+
+    """
+    def _trace_single_line1(t):
+        t = t.sorted_components()
+        components = t.components
+        ncomps = len(components)
+        g = LorentzIndex.metric
+        # gamma matirices are in a[i:j]
+        hit = 0
+        for i in range(ncomps):
+            if components[i] == GammaMatrix:
+                hit = 1
+                break
+
+        for j in range(i + hit, ncomps):
+            if components[j] != GammaMatrix:
+                break
+        else:
+            j = ncomps
+        numG = j - i
+        if numG == 0:
+            tcoeff = t.coeff
+            return t.nocoeff if tcoeff else t
+        if numG % 2 == 1:
+            return TensMul.from_data(S.Zero, [], [], [])
+        elif numG > 4:
+            # find the open matrix indices and connect them:
+            a = t.split()
+            ind1 = a[i].get_indices()[0]
+            ind2 = a[i + 1].get_indices()[0]
+            aa = a[:i] + a[i + 2:]
+            t1 = tensor_mul(*aa)*g(ind1, ind2)
+            t1 = t1.contract_metric(g)
+            args = [t1]
+            sign = 1
+            for k in range(i + 2, j):
+                sign = -sign
+                ind2 = a[k].get_indices()[0]
+                aa = a[:i] + a[i + 1:k] + a[k + 1:]
+                t2 = sign*tensor_mul(*aa)*g(ind1, ind2)
+                t2 = t2.contract_metric(g)
+                t2 = simplify_gpgp(t2, False)
+                args.append(t2)
+            t3 = TensAdd(*args)
+            t3 = _trace_single_line(t3)
+            return t3
+        else:
+            a = t.split()
+            t1 = _gamma_trace1(*a[i:j])
+            a2 = a[:i] + a[j:]
+            t2 = tensor_mul(*a2)
+            t3 = t1*t2
+            if not t3:
+                return t3
+            t3 = t3.contract_metric(g)
+            return t3
+
+    t = t.expand()
+    if isinstance(t, TensAdd):
+        a = [_trace_single_line1(x)*x.coeff for x in t.args]
+        return TensAdd(*a)
+    elif isinstance(t, (Tensor, TensMul)):
+        r = t.coeff*_trace_single_line1(t)
+        return r
+    else:
+        return trace(t)
+
+
+def _gamma_trace1(*a):
+    gctr = 4  # FIXME specific for d=4
+    g = LorentzIndex.metric
+    if not a:
+        return gctr
+    n = len(a)
+    if n%2 == 1:
+        #return TensMul.from_data(S.Zero, [], [], [])
+        return S.Zero
+    if n == 2:
+        ind0 = a[0].get_indices()[0]
+        ind1 = a[1].get_indices()[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]
+
+        return gctr*(g(ind0, ind1)*g(ind2, ind3) - \
+           g(ind0, ind2)*g(ind1, ind3) + g(ind0, ind3)*g(ind1, ind2))
+
+
+def kahane_simplify(expression):
+    r"""
+    This function cancels contracted elements in a product of four
+    dimensional gamma matrices, resulting in an expression equal to the given
+    one, without the contracted gamma matrices.
+
+    Parameters
+    ==========
+
+    `expression`    the tensor expression containing the gamma matrices to simplify.
+
+    Notes
+    =====
+
+    If spinor indices are given, the matrices must be given in
+    the order given in the product.
+
+    Algorithm
+    =========
+
+    The idea behind the algorithm is to use some well-known identities,
+    i.e., for contractions enclosing an even number of `\gamma` matrices
+
+    `\gamma^\mu \gamma_{a_1} \cdots \gamma_{a_{2N}} \gamma_\mu = 2 (\gamma_{a_{2N}} \gamma_{a_1} \cdots \gamma_{a_{2N-1}} + \gamma_{a_{2N-1}} \cdots \gamma_{a_1} \gamma_{a_{2N}} )`
+
+    for an odd number of `\gamma` matrices
+
+    `\gamma^\mu \gamma_{a_1} \cdots \gamma_{a_{2N+1}} \gamma_\mu = -2 \gamma_{a_{2N+1}} \gamma_{a_{2N}} \cdots \gamma_{a_{1}}`
+
+    Instead of repeatedly applying these identities to cancel out all contracted indices,
+    it is possible to recognize the links that would result from such an operation,
+    the problem is thus reduced to a simple rearrangement of free gamma matrices.
+
+    Examples
+    ========
+
+    When using, always remember that the original expression coefficient
+    has to be handled separately
+
+    >>> from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex
+    >>> from sympy.physics.hep.gamma_matrices import kahane_simplify
+    >>> from sympy.tensor.tensor import tensor_indices
+    >>> i0, i1, i2 = tensor_indices('i0:3', LorentzIndex)
+    >>> ta = G(i0)*G(-i0)
+    >>> kahane_simplify(ta)
+    Matrix([
+    [4, 0, 0, 0],
+    [0, 4, 0, 0],
+    [0, 0, 4, 0],
+    [0, 0, 0, 4]])
+    >>> tb = G(i0)*G(i1)*G(-i0)
+    >>> kahane_simplify(tb)
+    -2*GammaMatrix(i1)
+    >>> t = G(i0)*G(-i0)
+    >>> kahane_simplify(t)
+    Matrix([
+    [4, 0, 0, 0],
+    [0, 4, 0, 0],
+    [0, 0, 4, 0],
+    [0, 0, 0, 4]])
+    >>> t = G(i0)*G(-i0)
+    >>> kahane_simplify(t)
+    Matrix([
+    [4, 0, 0, 0],
+    [0, 4, 0, 0],
+    [0, 0, 4, 0],
+    [0, 0, 0, 4]])
+
+    If there are no contractions, the same expression is returned
+
+    >>> tc = G(i0)*G(i1)
+    >>> kahane_simplify(tc)
+    GammaMatrix(i0)*GammaMatrix(i1)
+
+    References
+    ==========
+
+    [1] Algorithm for Reducing Contracted Products of gamma Matrices,
+    Joseph Kahane, Journal of Mathematical Physics, Vol. 9, No. 10, October 1968.
+    """
+
+    if isinstance(expression, Mul):
+        return expression
+    if isinstance(expression, TensAdd):
+        return TensAdd(*[kahane_simplify(arg) for arg in expression.args])
+
+    if isinstance(expression, Tensor):
+        return expression
+
+    assert isinstance(expression, TensMul)
+
+    gammas = expression.args
+
+    for gamma in gammas:
+        assert gamma.component == GammaMatrix
+
+    free = expression.free
+    # spinor_free = [_ for _ in expression.free_in_args if _[1] != 0]
+
+    # if len(spinor_free) == 2:
+    #     spinor_free.sort(key=lambda x: x[2])
+    #     assert spinor_free[0][1] == 1 and spinor_free[-1][1] == 2
+    #     assert spinor_free[0][2] == 0
+    # elif spinor_free:
+    #     raise ValueError('spinor indices do not match')
+
+    dum = []
+    for dum_pair in expression.dum:
+        if expression.index_types[dum_pair[0]] == LorentzIndex:
+            dum.append((dum_pair[0], dum_pair[1]))
+
+    dum = sorted(dum)
+
+    if len(dum) == 0:  # or GammaMatrixHead:
+        # no contractions in `expression`, just return it.
+        return expression
+
+    # find the `first_dum_pos`, i.e. the position of the first contracted
+    # gamma matrix, Kahane's algorithm as described in his paper requires the
+    # gamma matrix expression to start with a contracted gamma matrix, this is
+    # a workaround which ignores possible initial free indices, and re-adds
+    # them later.
+
+    first_dum_pos = min(map(min, dum))
+
+    # for p1, p2, a1, a2 in expression.dum_in_args:
+    #     if p1 != 0 or p2 != 0:
+    #         # only Lorentz indices, skip Dirac indices:
+    #         continue
+    #     first_dum_pos = min(p1, p2)
+    #     break
+
+    total_number = len(free) + len(dum)*2
+    number_of_contractions = len(dum)
+
+    free_pos = [None]*total_number
+    for i in free:
+        free_pos[i[1]] = i[0]
+
+    # `index_is_free` is a list of booleans, to identify index position
+    # and whether that index is free or dummy.
+    index_is_free = [False]*total_number
+
+    for i, indx in enumerate(free):
+        index_is_free[indx[1]] = True
+
+    # `links` is a dictionary containing the graph described in Kahane's paper,
+    # to every key correspond one or two values, representing the linked indices.
+    # All values in `links` are integers, negative numbers are used in the case
+    # where it is necessary to insert gamma matrices between free indices, in
+    # order to make Kahane's algorithm work (see paper).
+    links = {i: [] for i in range(first_dum_pos, total_number)}
+
+    # `cum_sign` is a step variable to mark the sign of every index, see paper.
+    cum_sign = -1
+    # `cum_sign_list` keeps storage for all `cum_sign` (every index).
+    cum_sign_list = [None]*total_number
+    block_free_count = 0
+
+    # multiply `resulting_coeff` by the coefficient parameter, the rest
+    # of the algorithm ignores a scalar coefficient.
+    resulting_coeff = S.One
+
+    # initialize a list of lists of indices. The outer list will contain all
+    # additive tensor expressions, while the inner list will contain the
+    # free indices (rearranged according to the algorithm).
+    resulting_indices = [[]]
+
+    # start to count the `connected_components`, which together with the number
+    # of contractions, determines a -1 or +1 factor to be multiplied.
+    connected_components = 1
+
+    # First loop: here we fill `cum_sign_list`, and draw the links
+    # among consecutive indices (they are stored in `links`). Links among
+    # non-consecutive indices will be drawn later.
+    for i, is_free in enumerate(index_is_free):
+        # if `expression` starts with free indices, they are ignored here;
+        # they are later added as they are to the beginning of all
+        # `resulting_indices` list of lists of indices.
+        if i < first_dum_pos:
+            continue
+
+        if is_free:
+            block_free_count += 1
+            # if previous index was free as well, draw an arch in `links`.
+            if block_free_count > 1:
+                links[i - 1].append(i)
+                links[i].append(i - 1)
+        else:
+            # Change the sign of the index (`cum_sign`) if the number of free
+            # indices preceding it is even.
+            cum_sign *= 1 if (block_free_count % 2) else -1
+            if block_free_count == 0 and i != first_dum_pos:
+                # check if there are two consecutive dummy indices:
+                # in this case create virtual indices with negative position,
+                # these "virtual" indices represent the insertion of two
+                # gamma^0 matrices to separate consecutive dummy indices, as
+                # Kahane's algorithm requires dummy indices to be separated by
+                # free indices. The product of two gamma^0 matrices is unity,
+                # so the new expression being examined is the same as the
+                # original one.
+                if cum_sign == -1:
+                    links[-1-i] = [-1-i+1]
+                    links[-1-i+1] = [-1-i]
+            if (i - cum_sign) in links:
+                if i != first_dum_pos:
+                    links[i].append(i - cum_sign)
+                if block_free_count != 0:
+                    if i - cum_sign < len(index_is_free):
+                        if index_is_free[i - cum_sign]:
+                            links[i - cum_sign].append(i)
+            block_free_count = 0
+
+        cum_sign_list[i] = cum_sign
+
+    # The previous loop has only created links between consecutive free indices,
+    # it is necessary to properly create links among dummy (contracted) indices,
+    # according to the rules described in Kahane's paper. There is only one exception
+    # to Kahane's rules: the negative indices, which handle the case of some
+    # consecutive free indices (Kahane's paper just describes dummy indices
+    # separated by free indices, hinting that free indices can be added without
+    # altering the expression result).
+    for i in dum:
+        # get the positions of the two contracted indices:
+        pos1 = i[0]
+        pos2 = i[1]
+
+        # create Kahane's upper links, i.e. the upper arcs between dummy
+        # (i.e. contracted) indices:
+        links[pos1].append(pos2)
+        links[pos2].append(pos1)
+
+        # create Kahane's lower links, this corresponds to the arcs below
+        # the line described in the paper:
+
+        # first we move `pos1` and `pos2` according to the sign of the indices:
+        linkpos1 = pos1 + cum_sign_list[pos1]
+        linkpos2 = pos2 + cum_sign_list[pos2]
+
+        # otherwise, perform some checks before creating the lower arcs:
+
+        # make sure we are not exceeding the total number of indices:
+        if linkpos1 >= total_number:
+            continue
+        if linkpos2 >= total_number:
+            continue
+
+        # make sure we are not below the first dummy index in `expression`:
+        if linkpos1 < first_dum_pos:
+            continue
+        if linkpos2 < first_dum_pos:
+            continue
+
+        # check if the previous loop created "virtual" indices between dummy
+        # indices, in such a case relink `linkpos1` and `linkpos2`:
+        if (-1-linkpos1) in links:
+            linkpos1 = -1-linkpos1
+        if (-1-linkpos2) in links:
+            linkpos2 = -1-linkpos2
+
+        # move only if not next to free index:
+        if linkpos1 >= 0 and not index_is_free[linkpos1]:
+            linkpos1 = pos1
+
+        if linkpos2 >=0 and not index_is_free[linkpos2]:
+            linkpos2 = pos2
+
+        # create the lower arcs:
+        if linkpos2 not in links[linkpos1]:
+            links[linkpos1].append(linkpos2)
+        if linkpos1 not in links[linkpos2]:
+            links[linkpos2].append(linkpos1)
+
+    # This loop starts from the `first_dum_pos` index (first dummy index)
+    # walks through the graph deleting the visited indices from `links`,
+    # it adds a gamma matrix for every free index in encounters, while it
+    # completely ignores dummy indices and virtual indices.
+    pointer = first_dum_pos
+    previous_pointer = 0
+    while True:
+        if pointer in links:
+            next_ones = links.pop(pointer)
+        else:
+            break
+
+        if previous_pointer in next_ones:
+            next_ones.remove(previous_pointer)
+
+        previous_pointer = pointer
+
+        if next_ones:
+            pointer = next_ones[0]
+        else:
+            break
+
+        if pointer == previous_pointer:
+            break
+        if pointer >=0 and free_pos[pointer] is not None:
+            for ri in resulting_indices:
+                ri.append(free_pos[pointer])
+
+    # The following loop removes the remaining connected components in `links`.
+    # If there are free indices inside a connected component, it gives a
+    # contribution to the resulting expression given by the factor
+    # `gamma_a gamma_b ... gamma_z + gamma_z ... gamma_b gamma_a`, in Kahanes's
+    # paper represented as  {gamma_a, gamma_b, ... , gamma_z},
+    # virtual indices are ignored. The variable `connected_components` is
+    # increased by one for every connected component this loop encounters.
+
+    # If the connected component has virtual and dummy indices only
+    # (no free indices), it contributes to `resulting_indices` by a factor of two.
+    # The multiplication by two is a result of the
+    # factor {gamma^0, gamma^0} = 2 I, as it appears in Kahane's paper.
+    # Note: curly brackets are meant as in the paper, as a generalized
+    # multi-element anticommutator!
+
+    while links:
+        connected_components += 1
+        pointer = min(links.keys())
+        previous_pointer = pointer
+        # the inner loop erases the visited indices from `links`, and it adds
+        # all free indices to `prepend_indices` list, virtual indices are
+        # ignored.
+        prepend_indices = []
+        while True:
+            if pointer in links:
+                next_ones = links.pop(pointer)
+            else:
+                break
+
+            if previous_pointer in next_ones:
+                if len(next_ones) > 1:
+                    next_ones.remove(previous_pointer)
+
+            previous_pointer = pointer
+
+            if next_ones:
+                pointer = next_ones[0]
+
+            if pointer >= first_dum_pos and free_pos[pointer] is not None:
+                prepend_indices.insert(0, free_pos[pointer])
+        # if `prepend_indices` is void, it means there are no free indices
+        # in the loop (and it can be shown that there must be a virtual index),
+        # loops of virtual indices only contribute by a factor of two:
+        if len(prepend_indices) == 0:
+            resulting_coeff *= 2
+        # otherwise, add the free indices in `prepend_indices` to
+        # the `resulting_indices`:
+        else:
+            expr1 = prepend_indices
+            expr2 = list(reversed(prepend_indices))
+            resulting_indices = [expri + ri for ri in resulting_indices for expri in (expr1, expr2)]
+
+    # sign correction, as described in Kahane's paper:
+    resulting_coeff *= -1 if (number_of_contractions - connected_components + 1) % 2 else 1
+    # power of two factor, as described in Kahane's paper:
+    resulting_coeff *= 2**(number_of_contractions)
+
+    # If `first_dum_pos` is not zero, it means that there are trailing free gamma
+    # matrices in front of `expression`, so multiply by them:
+    resulting_indices = [ free_pos[0:first_dum_pos] + ri for ri in resulting_indices ]
+
+    resulting_expr = S.Zero
+    for i in resulting_indices:
+        temp_expr = S.One
+        for j in i:
+            temp_expr *= GammaMatrix(j)
+        resulting_expr += temp_expr
+
+    t = resulting_coeff * resulting_expr
+    t1 = None
+    if isinstance(t, TensAdd):
+        t1 = t.args[0]
+    elif isinstance(t, TensMul):
+        t1 = t
+    if t1:
+        pass
+    else:
+        t = eye(4)*t
+    return t

backend_service/venv/lib/python3.13/site-packages/sympy/physics/hep/tests/test_gamma_matrices.py

@@ -0,0 +1,427 @@
+from sympy.matrices.dense import eye, Matrix
+from sympy.tensor.tensor import tensor_indices, TensorHead, tensor_heads, \
+    TensExpr, canon_bp
+from sympy.physics.hep.gamma_matrices import GammaMatrix as G, LorentzIndex, \
+    kahane_simplify, gamma_trace, _simplify_single_line, simplify_gamma_expression
+from sympy import Symbol
+
+
+def _is_tensor_eq(arg1, arg2):
+    arg1 = canon_bp(arg1)
+    arg2 = canon_bp(arg2)
+    if isinstance(arg1, TensExpr):
+        return arg1.equals(arg2)
+    elif isinstance(arg2, TensExpr):
+        return arg2.equals(arg1)
+    return arg1 == arg2
+
+def execute_gamma_simplify_tests_for_function(tfunc, D):
+    """
+    Perform tests to check if sfunc is able to simplify gamma matrix expressions.
+
+    Parameters
+    ==========
+
+    `sfunc`     a function to simplify a `TIDS`, shall return the simplified `TIDS`.
+    `D`         the number of dimension (in most cases `D=4`).
+
+    """
+
+    mu, nu, rho, sigma = tensor_indices("mu, nu, rho, sigma", LorentzIndex)
+    a1, a2, a3, a4, a5, a6 = tensor_indices("a1:7", LorentzIndex)
+    mu11, mu12, mu21, mu31, mu32, mu41, mu51, mu52 = tensor_indices("mu11, mu12, mu21, mu31, mu32, mu41, mu51, mu52", LorentzIndex)
+    mu61, mu71, mu72 = tensor_indices("mu61, mu71, mu72", LorentzIndex)
+    m0, m1, m2, m3, m4, m5, m6 = tensor_indices("m0:7", LorentzIndex)
+
+    def g(xx, yy):
+        return (G(xx)*G(yy) + G(yy)*G(xx))/2
+
+    # Some examples taken from Kahane's paper, 4 dim only:
+    if D == 4:
+        t = (G(a1)*G(mu11)*G(a2)*G(mu21)*G(-a1)*G(mu31)*G(-a2))
+        assert _is_tensor_eq(tfunc(t), -4*G(mu11)*G(mu31)*G(mu21) - 4*G(mu31)*G(mu11)*G(mu21))
+
+        t = (G(a1)*G(mu11)*G(mu12)*\
+                              G(a2)*G(mu21)*\
+                              G(a3)*G(mu31)*G(mu32)*\
+                              G(a4)*G(mu41)*\
+                              G(-a2)*G(mu51)*G(mu52)*\
+                              G(-a1)*G(mu61)*\
+                              G(-a3)*G(mu71)*G(mu72)*\
+                              G(-a4))
+        assert _is_tensor_eq(tfunc(t), \
+            16*G(mu31)*G(mu32)*G(mu72)*G(mu71)*G(mu11)*G(mu52)*G(mu51)*G(mu12)*G(mu61)*G(mu21)*G(mu41) + 16*G(mu31)*G(mu32)*G(mu72)*G(mu71)*G(mu12)*G(mu51)*G(mu52)*G(mu11)*G(mu61)*G(mu21)*G(mu41) + 16*G(mu71)*G(mu72)*G(mu32)*G(mu31)*G(mu11)*G(mu52)*G(mu51)*G(mu12)*G(mu61)*G(mu21)*G(mu41) + 16*G(mu71)*G(mu72)*G(mu32)*G(mu31)*G(mu12)*G(mu51)*G(mu52)*G(mu11)*G(mu61)*G(mu21)*G(mu41))
+
+    # Fully Lorentz-contracted expressions, these return scalars:
+
+    def add_delta(ne):
+        return ne * eye(4)  # DiracSpinorIndex.delta(DiracSpinorIndex.auto_left, -DiracSpinorIndex.auto_right)
+
+    t = (G(mu)*G(-mu))
+    ts = add_delta(D)
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(mu)*G(nu)*G(-mu)*G(-nu))
+    ts = add_delta(2*D - D**2)  # -8
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(mu)*G(nu)*G(-nu)*G(-mu))
+    ts = add_delta(D**2)  # 16
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(mu)*G(nu)*G(-rho)*G(-nu)*G(-mu)*G(rho))
+    ts = add_delta(4*D - 4*D**2 + D**3)  # 16
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(mu)*G(nu)*G(rho)*G(-rho)*G(-nu)*G(-mu))
+    ts = add_delta(D**3)  # 64
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(a1)*G(a2)*G(a3)*G(a4)*G(-a3)*G(-a1)*G(-a2)*G(-a4))
+    ts = add_delta(-8*D + 16*D**2 - 8*D**3 + D**4)  # -32
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(-mu)*G(-nu)*G(-rho)*G(-sigma)*G(nu)*G(mu)*G(sigma)*G(rho))
+    ts = add_delta(-16*D + 24*D**2 - 8*D**3 + D**4)  # 64
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(-mu)*G(nu)*G(-rho)*G(sigma)*G(rho)*G(-nu)*G(mu)*G(-sigma))
+    ts = add_delta(8*D - 12*D**2 + 6*D**3 - D**4)  # -32
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(-a3)*G(-a2)*G(-a1)*G(-a5)*G(-a4))
+    ts = add_delta(64*D - 112*D**2 + 60*D**3 - 12*D**4 + D**5)  # 256
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(-a3)*G(-a1)*G(-a2)*G(-a4)*G(-a5))
+    ts = add_delta(64*D - 120*D**2 + 72*D**3 - 16*D**4 + D**5)  # -128
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(a6)*G(-a3)*G(-a2)*G(-a1)*G(-a6)*G(-a5)*G(-a4))
+    ts = add_delta(416*D - 816*D**2 + 528*D**3 - 144*D**4 + 18*D**5 - D**6)  # -128
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    t = (G(a1)*G(a2)*G(a3)*G(a4)*G(a5)*G(a6)*G(-a2)*G(-a3)*G(-a1)*G(-a6)*G(-a4)*G(-a5))
+    ts = add_delta(416*D - 848*D**2 + 584*D**3 - 172*D**4 + 22*D**5 - D**6)  # -128
+    assert _is_tensor_eq(tfunc(t), ts)
+
+    # Expressions with free indices:
+
+    t = (G(mu)*G(nu)*G(rho)*G(sigma)*G(-mu))
+    assert _is_tensor_eq(tfunc(t), (-2*G(sigma)*G(rho)*G(nu) + (4-D)*G(nu)*G(rho)*G(sigma)))
+
+    t = (G(mu)*G(nu)*G(-mu))
+    assert _is_tensor_eq(tfunc(t), (2-D)*G(nu))
+
+    t = (G(mu)*G(nu)*G(rho)*G(-mu))
+    assert _is_tensor_eq(tfunc(t), 2*G(nu)*G(rho) + 2*G(rho)*G(nu) - (4-D)*G(nu)*G(rho))
+
+    t = 2*G(m2)*G(m0)*G(m1)*G(-m0)*G(-m1)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, (D*(-2*D + 4))*G(m2))
+
+    t = G(m2)*G(m0)*G(m1)*G(-m0)*G(-m2)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, ((-D + 2)**2)*G(m1))
+
+    t = G(m0)*G(m1)*G(m2)*G(m3)*G(-m1)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, (D - 4)*G(m0)*G(m2)*G(m3) + 4*G(m0)*g(m2, m3))
+
+    t = G(m0)*G(m1)*G(m2)*G(m3)*G(-m1)*G(-m0)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, ((D - 4)**2)*G(m2)*G(m3) + (8*D - 16)*g(m2, m3))
+
+    t = G(m2)*G(m0)*G(m1)*G(-m2)*G(-m0)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, ((-D + 2)*(D - 4) + 4)*G(m1))
+
+    t = G(m3)*G(m1)*G(m0)*G(m2)*G(-m3)*G(-m0)*G(-m2)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, (-4*D + (-D + 2)**2*(D - 4) + 8)*G(m1))
+
+    t = 2*G(m0)*G(m1)*G(m2)*G(m3)*G(-m0)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, ((-2*D + 8)*G(m1)*G(m2)*G(m3) - 4*G(m3)*G(m2)*G(m1)))
+
+    t = G(m5)*G(m0)*G(m1)*G(m4)*G(m2)*G(-m4)*G(m3)*G(-m0)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, (((-D + 2)*(-D + 4))*G(m5)*G(m1)*G(m2)*G(m3) + (2*D - 4)*G(m5)*G(m3)*G(m2)*G(m1)))
+
+    t = -G(m0)*G(m1)*G(m2)*G(m3)*G(-m0)*G(m4)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, ((D - 4)*G(m1)*G(m2)*G(m3)*G(m4) + 2*G(m3)*G(m2)*G(m1)*G(m4)))
+
+    t = G(-m5)*G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)*G(m5)
+    st = tfunc(t)
+
+    result1 = ((-D + 4)**2 + 4)*G(m1)*G(m2)*G(m3)*G(m4) +\
+        (4*D - 16)*G(m3)*G(m2)*G(m1)*G(m4) + (4*D - 16)*G(m4)*G(m1)*G(m2)*G(m3)\
+        + 4*G(m2)*G(m1)*G(m4)*G(m3) + 4*G(m3)*G(m4)*G(m1)*G(m2) +\
+        4*G(m4)*G(m3)*G(m2)*G(m1)
+
+    # Kahane's algorithm yields this result, which is equivalent to `result1`
+    # in four dimensions, but is not automatically recognized as equal:
+    result2 = 8*G(m1)*G(m2)*G(m3)*G(m4) + 8*G(m4)*G(m3)*G(m2)*G(m1)
+
+    if D == 4:
+        assert _is_tensor_eq(st, (result1)) or _is_tensor_eq(st, (result2))
+    else:
+        assert _is_tensor_eq(st, (result1))
+
+    # and a few very simple cases, with no contracted indices:
+
+    t = G(m0)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, t)
+
+    t = -7*G(m0)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, t)
+
+    t = 224*G(m0)*G(m1)*G(-m2)*G(m3)
+    st = tfunc(t)
+    assert _is_tensor_eq(st, t)
+
+
+def test_kahane_algorithm():
+    # Wrap this function to convert to and from TIDS:
+
+    def tfunc(e):
+        return _simplify_single_line(e)
+
+    execute_gamma_simplify_tests_for_function(tfunc, D=4)
+
+
+def test_kahane_simplify1():
+    i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15 = tensor_indices('i0:16', LorentzIndex)
+    mu, nu, rho, sigma = tensor_indices("mu, nu, rho, sigma", LorentzIndex)
+    D = 4
+    t = G(i0)*G(i1)
+    r = kahane_simplify(t)
+    assert r.equals(t)
+
+    t = G(i0)*G(i1)*G(-i0)
+    r = kahane_simplify(t)
+    assert r.equals(-2*G(i1))
+    t = G(i0)*G(i1)*G(-i0)
+    r = kahane_simplify(t)
+    assert r.equals(-2*G(i1))
+
+    t = G(i0)*G(i1)
+    r = kahane_simplify(t)
+    assert r.equals(t)
+    t = G(i0)*G(i1)
+    r = kahane_simplify(t)
+    assert r.equals(t)
+    t = G(i0)*G(-i0)
+    r = kahane_simplify(t)
+    assert r.equals(4*eye(4))
+    t = G(i0)*G(-i0)
+    r = kahane_simplify(t)
+    assert r.equals(4*eye(4))
+    t = G(i0)*G(-i0)
+    r = kahane_simplify(t)
+    assert r.equals(4*eye(4))
+    t = G(i0)*G(i1)*G(-i0)
+    r = kahane_simplify(t)
+    assert r.equals(-2*G(i1))
+    t = G(i0)*G(i1)*G(-i0)*G(-i1)
+    r = kahane_simplify(t)
+    assert r.equals((2*D - D**2)*eye(4))
+    t = G(i0)*G(i1)*G(-i0)*G(-i1)
+    r = kahane_simplify(t)
+    assert r.equals((2*D - D**2)*eye(4))
+    t = G(i0)*G(-i0)*G(i1)*G(-i1)
+    r = kahane_simplify(t)
+    assert r.equals(16*eye(4))
+    t = (G(mu)*G(nu)*G(-nu)*G(-mu))
+    r = kahane_simplify(t)
+    assert r.equals(D**2*eye(4))
+    t = (G(mu)*G(nu)*G(-nu)*G(-mu))
+    r = kahane_simplify(t)
+    assert r.equals(D**2*eye(4))
+    t = (G(mu)*G(nu)*G(-nu)*G(-mu))
+    r = kahane_simplify(t)
+    assert r.equals(D**2*eye(4))
+    t = (G(mu)*G(nu)*G(-rho)*G(-nu)*G(-mu)*G(rho))
+    r = kahane_simplify(t)
+    assert r.equals((4*D - 4*D**2 + D**3)*eye(4))
+    t = (G(-mu)*G(-nu)*G(-rho)*G(-sigma)*G(nu)*G(mu)*G(sigma)*G(rho))
+    r = kahane_simplify(t)
+    assert r.equals((-16*D + 24*D**2 - 8*D**3 + D**4)*eye(4))
+    t = (G(-mu)*G(nu)*G(-rho)*G(sigma)*G(rho)*G(-nu)*G(mu)*G(-sigma))
+    r = kahane_simplify(t)
+    assert r.equals((8*D - 12*D**2 + 6*D**3 - D**4)*eye(4))
+
+    # Expressions with free indices:
+    t = (G(mu)*G(nu)*G(rho)*G(sigma)*G(-mu))
+    r = kahane_simplify(t)
+    assert r.equals(-2*G(sigma)*G(rho)*G(nu))
+    t = (G(mu)*G(-mu)*G(rho)*G(sigma))
+    r = kahane_simplify(t)
+    assert r.equals(4*G(rho)*G(sigma))
+    t = (G(rho)*G(sigma)*G(mu)*G(-mu))
+    r = kahane_simplify(t)
+    assert r.equals(4*G(rho)*G(sigma))
+
+def test_gamma_matrix_class():
+    i, j, k = tensor_indices('i,j,k', LorentzIndex)
+
+    # define another type of TensorHead to see if exprs are correctly handled:
+    A = TensorHead('A', [LorentzIndex])
+
+    t = A(k)*G(i)*G(-i)
+    ts = simplify_gamma_expression(t)
+    assert _is_tensor_eq(ts, Matrix([
+        [4, 0, 0, 0],
+        [0, 4, 0, 0],
+        [0, 0, 4, 0],
+        [0, 0, 0, 4]])*A(k))
+
+    t = G(i)*A(k)*G(j)
+    ts = simplify_gamma_expression(t)
+    assert _is_tensor_eq(ts, A(k)*G(i)*G(j))
+
+    execute_gamma_simplify_tests_for_function(simplify_gamma_expression, D=4)
+
+
+def test_gamma_matrix_trace():
+    g = LorentzIndex.metric
+
+    m0, m1, m2, m3, m4, m5, m6 = tensor_indices('m0:7', LorentzIndex)
+    n0, n1, n2, n3, n4, n5 = tensor_indices('n0:6', LorentzIndex)
+
+    # working in D=4 dimensions
+    D = 4
+
+    # traces of odd number of gamma matrices are zero:
+    t = G(m0)
+    t1 = gamma_trace(t)
+    assert t1.equals(0)
+
+    t = G(m0)*G(m1)*G(m2)
+    t1 = gamma_trace(t)
+    assert t1.equals(0)
+
+    t = G(m0)*G(m1)*G(-m0)
+    t1 = gamma_trace(t)
+    assert t1.equals(0)
+
+    t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)
+    t1 = gamma_trace(t)
+    assert t1.equals(0)
+
+    # traces without internal contractions:
+    t = G(m0)*G(m1)
+    t1 = gamma_trace(t)
+    assert _is_tensor_eq(t1, 4*g(m0, m1))
+
+    t = G(m0)*G(m1)*G(m2)*G(m3)
+    t1 = gamma_trace(t)
+    t2 = -4*g(m0, m2)*g(m1, m3) + 4*g(m0, m1)*g(m2, m3) + 4*g(m0, m3)*g(m1, m2)
+    assert _is_tensor_eq(t1, t2)
+
+    t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(m5)
+    t1 = gamma_trace(t)
+    t2 = t1*g(-m0, -m5)
+    t2 = t2.contract_metric(g)
+    assert _is_tensor_eq(t2, D*gamma_trace(G(m1)*G(m2)*G(m3)*G(m4)))
+
+    # traces of expressions with internal contractions:
+    t = G(m0)*G(-m0)
+    t1 = gamma_trace(t)
+    assert t1.equals(4*D)
+
+    t = G(m0)*G(m1)*G(-m0)*G(-m1)
+    t1 = gamma_trace(t)
+    assert t1.equals(8*D - 4*D**2)
+
+    t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)
+    t1 = gamma_trace(t)
+    t2 = (-4*D)*g(m1, m3)*g(m2, m4) + (4*D)*g(m1, m2)*g(m3, m4) + \
+                 (4*D)*g(m1, m4)*g(m2, m3)
+    assert _is_tensor_eq(t1, t2)
+
+    t = G(-m5)*G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(-m0)*G(m5)
+    t1 = gamma_trace(t)
+    t2 = (32*D + 4*(-D + 4)**2 - 64)*(g(m1, m2)*g(m3, m4) - \
+            g(m1, m3)*g(m2, m4) + g(m1, m4)*g(m2, m3))
+    assert _is_tensor_eq(t1, t2)
+
+    t = G(m0)*G(m1)*G(-m0)*G(m3)
+    t1 = gamma_trace(t)
+    assert t1.equals((-4*D + 8)*g(m1, m3))
+
+#    p, q = S1('p,q')
+#    ps = p(m0)*G(-m0)
+#    qs = q(m0)*G(-m0)
+#    t = ps*qs*ps*qs
+#    t1 = gamma_trace(t)
+#    assert t1 == 8*p(m0)*q(-m0)*p(m1)*q(-m1) - 4*p(m0)*p(-m0)*q(m1)*q(-m1)
+
+    t = G(m0)*G(m1)*G(m2)*G(m3)*G(m4)*G(m5)*G(-m0)*G(-m1)*G(-m2)*G(-m3)*G(-m4)*G(-m5)
+    t1 = gamma_trace(t)
+    assert t1.equals(-4*D**6 + 120*D**5 - 1040*D**4 + 3360*D**3 - 4480*D**2 + 2048*D)
+
+    t = G(m0)*G(m1)*G(n1)*G(m2)*G(n2)*G(m3)*G(m4)*G(-n2)*G(-n1)*G(-m0)*G(-m1)*G(-m2)*G(-m3)*G(-m4)
+    t1 = gamma_trace(t)
+    tresu = -7168*D + 16768*D**2 - 14400*D**3 + 5920*D**4 - 1232*D**5 + 120*D**6 - 4*D**7
+    assert t1.equals(tresu)
+
+    # checked with Mathematica
+    # In[1]:= <<Tracer.m
+    # In[2]:= Spur[l];
+    # In[3]:= GammaTrace[l, {m0},{m1},{n1},{m2},{n2},{m3},{m4},{n3},{n4},{m0},{m1},{m2},{m3},{m4}]
+    t = G(m0)*G(m1)*G(n1)*G(m2)*G(n2)*G(m3)*G(m4)*G(n3)*G(n4)*G(-m0)*G(-m1)*G(-m2)*G(-m3)*G(-m4)
+    t1 = gamma_trace(t)
+#    t1 = t1.expand_coeff()
+    c1 = -4*D**5 + 120*D**4 - 1200*D**3 + 5280*D**2 - 10560*D + 7808
+    c2 = -4*D**5 + 88*D**4 - 560*D**3 + 1440*D**2 - 1600*D + 640
+    assert _is_tensor_eq(t1, c1*g(n1, n4)*g(n2, n3) + c2*g(n1, n2)*g(n3, n4) + \
+            (-c1)*g(n1, n3)*g(n2, n4))
+
+    p, q = tensor_heads('p,q', [LorentzIndex])
+    ps = p(m0)*G(-m0)
+    qs = q(m0)*G(-m0)
+    p2 = p(m0)*p(-m0)
+    q2 = q(m0)*q(-m0)
+    pq = p(m0)*q(-m0)
+    t = ps*qs*ps*qs
+    r = gamma_trace(t)
+    assert _is_tensor_eq(r, 8*pq*pq - 4*p2*q2)
+    t = ps*qs*ps*qs*ps*qs
+    r = gamma_trace(t)
+    assert _is_tensor_eq(r, -12*p2*pq*q2 + 16*pq*pq*pq)
+    t = ps*qs*ps*qs*ps*qs*ps*qs
+    r = gamma_trace(t)
+    assert _is_tensor_eq(r, -32*pq*pq*p2*q2 + 32*pq*pq*pq*pq + 4*p2*p2*q2*q2)
+
+    t = 4*p(m1)*p(m0)*p(-m0)*q(-m1)*q(m2)*q(-m2)
+    assert _is_tensor_eq(gamma_trace(t), t)
+    t = ps*ps*ps*ps*ps*ps*ps*ps
+    r = gamma_trace(t)
+    assert r.equals(4*p2*p2*p2*p2)
+
+
+def test_bug_13636():
+    """Test issue 13636 regarding handling traces of sums of products
+    of GammaMatrix mixed with other factors."""
+    pi, ki, pf = tensor_heads("pi, ki, pf", [LorentzIndex])
+    i0, i1, i2, i3, i4 = tensor_indices("i0:5", LorentzIndex)
+    x = Symbol("x")
+    pis = pi(i2) * G(-i2)
+    kis = ki(i3) * G(-i3)
+    pfs = pf(i4) * G(-i4)
+
+    a = pfs * G(i0) * kis * G(i1) * pis * G(-i1) * kis * G(-i0)
+    b = pfs * G(i0) * kis * G(i1) * pis * x * G(-i0) * pi(-i1)
+    ta = gamma_trace(a)
+    tb = gamma_trace(b)
+    t_a_plus_b = gamma_trace(a + b)
+    assert ta == 4 * (
+        -4 * ki(i0) * ki(-i0) * pf(i1) * pi(-i1)
+        + 8 * ki(i0) * ki(i1) * pf(-i0) * pi(-i1)
+    )
+    assert tb == -8 * x * ki(i0) * pf(-i0) * pi(i1) * pi(-i1)
+    assert t_a_plus_b == ta + tb