From 43e0a8f5ca5e1c48cef5daaf014acdbca4e44568 Mon Sep 17 00:00:00 2001
From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Thu, 26 Apr 2001 19:36:07 +0000
Subject: [PATCH 1/1] - dirac_trace() handles arbitrary expressions (even
 unexpanded), including   gamma5 - implemented color_trace() - eps0123()
 creates a 4-dimensional epsilon tensor without checking the   indices (so
 they may be D-dimensional); this may become a different class   soon, but for
 now all the implemented properties of the standard epsilon   tensor also
 apply to eps0123

---
 ginac/clifford.cpp | 55 ++++++++++++++++++++++-------
 ginac/clifford.h   | 13 +++++--
 ginac/color.cpp    | 87 ++++++++++++++++++++++++++++++++++++++++++++++
 ginac/color.h      |  6 ++++
 ginac/tensor.cpp   | 12 +++++++
 ginac/tensor.h     | 12 +++++++
 6 files changed, 170 insertions(+), 15 deletions(-)

diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp
index d346c44a..8a8bfe09 100644
--- a/ginac/clifford.cpp
+++ b/ginac/clifford.cpp
@@ -24,6 +24,7 @@
 #include "ex.h"
 #include "idx.h"
 #include "ncmul.h"
+#include "symbol.h"
 #include "print.h"
 #include "archive.h"
 #include "debugmsg.h"
@@ -301,7 +302,20 @@ ex dirac_gamma5(unsigned char rl)
 	return clifford(diracgamma5(), rl);
 }
 
-ex dirac_trace(const ex & e, unsigned char rl = 0)
+ex dirac_slash(const ex & e, const ex & dim, unsigned char rl)
+{
+	varidx mu((new symbol)->setflag(status_flags::dynallocated), dim);
+	return indexed(e, mu.toggle_variance()) * dirac_gamma(mu, rl);
+}
+
+/** Check whether a given tinfo key (as returned by return_type_tinfo()
+ *  is that of a clifford object with the specified representation label. */
+static bool is_clifford_tinfo(unsigned ti, unsigned char rl)
+{
+	return ti == (TINFO_clifford + rl);
+}
+
+ex dirac_trace(const ex & e, unsigned char rl)
 {
 	if (is_ex_of_type(e, clifford)) {
 
@@ -325,12 +339,8 @@ ex dirac_trace(const ex & e, unsigned char rl = 0)
 		ex prod = _ex1();
 		for (unsigned i=0; i<e.nops(); i++) {
 			const ex &o = e.op(i);
-			if (is_ex_of_type(o, clifford)
-			 && ex_to_clifford(o).get_representation_label() == rl)
-				prod *= dirac_trace(o, rl);
-			else if (is_ex_of_type(o, ncmul)
-			 && is_ex_of_type(o.op(0), clifford)
-			 && ex_to_clifford(o.op(0)).get_representation_label() == rl)
+			unsigned ti = o.return_type_tinfo();
+			if (is_clifford_tinfo(o.return_type_tinfo(), rl))
 				prod *= dirac_trace(o, rl);
 			else
 				prod *= o;
@@ -339,10 +349,14 @@ ex dirac_trace(const ex & e, unsigned char rl = 0)
 
 	} else if (is_ex_exactly_of_type(e, ncmul)) {
 
-		if (!is_ex_of_type(e.op(0), clifford)
-		 || ex_to_clifford(e.op(0)).get_representation_label() != rl)
+		if (!is_clifford_tinfo(e.return_type_tinfo(), rl))
 			return _ex0();
 
+		// Expand product, if necessary
+		ex e_expanded = e.expand();
+		if (!is_ex_of_type(e_expanded, ncmul))
+			return dirac_trace(e_expanded, rl);
+
 		// gamma5 gets moved to the front so this check is enough
 		bool has_gamma5 = is_ex_of_type(e.op(0).op(0), diracgamma5);
 		unsigned num = e.nops();
@@ -351,9 +365,28 @@ ex dirac_trace(const ex & e, unsigned char rl = 0)
 
 			// Trace of gamma5 * odd number of gammas and trace of
 			// gamma5 * gamma_mu * gamma_nu are zero
-			if ((num & 1) == 0 || num == 2)
+			if ((num & 1) == 0 || num == 3)
 				return _ex0();
 
+			// Tr gamma5 S_2k =
+			//   epsilon0123_mu1_mu2_mu3_mu4 * Tr gamma_mu1 gamma_mu2 gamma_mu3 gamma_mu4 S_2k
+			ex dim = ex_to_idx(e.op(1).op(1)).get_dim();
+			varidx mu1((new symbol)->setflag(status_flags::dynallocated), dim),
+			       mu2((new symbol)->setflag(status_flags::dynallocated), dim),
+			       mu3((new symbol)->setflag(status_flags::dynallocated), dim),
+			       mu4((new symbol)->setflag(status_flags::dynallocated), dim);
+			exvector v;
+			v.reserve(num + 3);
+			v.push_back(dirac_gamma(mu1, rl));
+			v.push_back(dirac_gamma(mu2, rl));
+			v.push_back(dirac_gamma(mu3, rl));
+			v.push_back(dirac_gamma(mu4, rl));
+			for (int i=1; i<num; i++)
+				v.push_back(e.op(i));
+
+			return (eps0123(mu1.toggle_variance(), mu2.toggle_variance(), mu3.toggle_variance(), mu4.toggle_variance()) *
+			        dirac_trace(ncmul(v), rl)).simplify_indexed() / 24;
+
 		} else { // no gamma5
 
 			// Trace of odd number of gammas is zero
@@ -386,8 +419,6 @@ ex dirac_trace(const ex & e, unsigned char rl = 0)
 			}
 			return result;
 		}
-
-		throw (std::logic_error("dirac_trace: don't know how to compute trace"));
 	}
 
 	return _ex0();
diff --git a/ginac/clifford.h b/ginac/clifford.h
index f1783396..af3cb788 100644
--- a/ginac/clifford.h
+++ b/ginac/clifford.h
@@ -88,7 +88,8 @@ public:
 };
 
 
-/** This class represents the Dirac gamma5 object. */
+/** This class represents the Dirac gamma5 object which anticommutes with
+ *  all other gammas. */
 class diracgamma5 : public tensor
 {
 	GINAC_DECLARE_REGISTERED_CLASS(diracgamma5, tensor)
@@ -125,11 +126,17 @@ ex dirac_gamma(const ex & mu, unsigned char rl = 0);
  *  @return newly constructed object */
 ex dirac_gamma5(unsigned char rl = 0);
 
+/** Create a term of the form e_mu * gamma~mu with a unique index mu.
+ *
+ *  @param dim Dimension of index
+ *  @param rl Representation label */
+ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
+
 /** Calculate the trace of an expression containing gamma objects with
  *  a specified representation label. The computed trace is a linear
  *  functional that is equal to the usual trace only in D = 4 dimensions.
- *  In particular, the functional is non-cyclic in D != 4 dimensions when
- *  gamma5 is involved.
+ *  In particular, the functional is not always cyclic in D != 4 dimensions
+ *  when gamma5 is involved.
  *
  *  @param rl Representation label */
 ex dirac_trace(const ex & e, unsigned char rl = 0);
diff --git a/ginac/color.cpp b/ginac/color.cpp
index f8150c77..e36c1cbe 100644
--- a/ginac/color.cpp
+++ b/ginac/color.cpp
@@ -29,6 +29,7 @@
 #include "ncmul.h"
 #include "numeric.h"
 #include "power.h" // for sqrt()
+#include "symbol.h"
 #include "print.h"
 #include "archive.h"
 #include "debugmsg.h"
@@ -409,4 +410,90 @@ ex color_h(const ex & a, const ex & b, const ex & c)
 	return color_d(a, b, c) + I * color_f(a, b, c);
 }
 
+/** Check whether a given tinfo key (as returned by return_type_tinfo()
+ *  is that of a color object with the specified representation label. */
+static bool is_color_tinfo(unsigned ti, unsigned char rl)
+{
+	return ti == (TINFO_color + rl);
+}
+
+ex color_trace(const ex & e, unsigned char rl)
+{
+	if (is_ex_of_type(e, color)) {
+
+		if (ex_to_color(e).get_representation_label() == rl
+		 && is_ex_of_type(e.op(0), su3one))
+			return _ex3();
+		else
+			return _ex0();
+
+	} else if (is_ex_exactly_of_type(e, add)) {
+
+		// Trace of sum = sum of traces
+		ex sum = _ex0();
+		for (unsigned i=0; i<e.nops(); i++)
+			sum += color_trace(e.op(i), rl);
+		return sum;
+
+	} else if (is_ex_exactly_of_type(e, mul)) {
+
+		// Trace of product: pull out non-color factors
+		ex prod = _ex1();
+		for (unsigned i=0; i<e.nops(); i++) {
+			const ex &o = e.op(i);
+			if (is_color_tinfo(o.return_type_tinfo(), rl))
+				prod *= color_trace(o, rl);
+			else
+				prod *= o;
+		}
+		return prod;
+
+	} else if (is_ex_exactly_of_type(e, ncmul)) {
+
+		if (!is_color_tinfo(e.return_type_tinfo(), rl))
+			return _ex0();
+
+		// Expand product, if necessary
+		ex e_expanded = e.expand();
+		if (!is_ex_of_type(e_expanded, ncmul))
+			return color_trace(e_expanded, rl);
+
+		unsigned num = e.nops();
+
+		if (num == 2) {
+
+			// Tr T_a T_b = 1/2 delta_a_b
+			return delta_tensor(e.op(0).op(1), e.op(1).op(1)) / 2;
+
+		} else if (num == 3) {
+
+			// Tr T_a T_b T_c = 1/4 h_a_b_c
+			return color_h(e.op(0).op(1), e.op(1).op(1), e.op(2).op(1)) / 4;
+
+		} else {
+
+			// Traces of 4 or more generators are computed recursively:
+			// Tr T_a1 .. T_an =
+			//     1/6 delta_a(n-1)_an Tr T_a1 .. T_a(n-2)
+			//   + 1/2 h_a(n-1)_an_k Tr T_a1 .. T_a(n-2) T_k
+			const ex &last_index = e.op(num - 1).op(1);
+			const ex &next_to_last_index = e.op(num - 2).op(1);
+			idx summation_index((new symbol)->setflag(status_flags::dynallocated), 8);
+
+			exvector v1;
+			v1.reserve(num - 2);
+			for (int i=0; i<num-2; i++)
+				v1.push_back(e.op(i));
+
+			exvector v2 = v1;
+			v2.push_back(color_T(summation_index, rl));
+
+			return delta_tensor(next_to_last_index, last_index) * color_trace(ncmul(v1), rl) / 6
+			       + color_h(next_to_last_index, last_index, summation_index) * color_trace(ncmul(v2), rl) / 2;
+		}
+	}
+
+	return _ex0();
+}
+
 } // namespace GiNaC
diff --git a/ginac/color.h b/ginac/color.h
index 9a89387c..70049325 100644
--- a/ginac/color.h
+++ b/ginac/color.h
@@ -153,6 +153,12 @@ ex color_d(const ex & a, const ex & b, const ex & c);
 /** This returns the linear combination d.a.b.c+I*f.a.b.c. */
 ex color_h(const ex & a, const ex & b, const ex & c);
 
+/** Calculate the trace of an expression containing color objects with a
+ *  specified representation label.
+ *
+ *  @param rl Representation label */
+ex color_trace(const ex & e, unsigned char rl = 0);
+
 
 } // namespace GiNaC
 
diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
index c14b50d7..d16bbfce 100644
--- a/ginac/tensor.cpp
+++ b/ginac/tensor.cpp
@@ -440,4 +440,16 @@ ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool
 	return indexed(tensepsilon(true, pos_sig), indexed::antisymmetric, i1, i2, i3, i4);
 }
 
+ex eps0123(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig)
+{
+	if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx) || !is_ex_of_type(i3, varidx) || !is_ex_of_type(i4, varidx))
+		throw(std::invalid_argument("indices of epsilon tensor must be of type varidx"));
+
+	ex dim = ex_to_idx(i1).get_dim();
+	if (dim.is_equal(4))
+		return lorentz_eps(i1, i2, i3, i4, pos_sig);
+	else
+		return indexed(tensepsilon(true, pos_sig), indexed::antisymmetric, i1, i2, i3, i4);
+}
+
 } // namespace GiNaC
diff --git a/ginac/tensor.h b/ginac/tensor.h
index 7396aca1..b077ee4d 100644
--- a/ginac/tensor.h
+++ b/ginac/tensor.h
@@ -183,6 +183,18 @@ ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
  *  @return newly constructed epsilon tensor */
 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig = false);
 
+/** Create an epsilon tensor in a 4-dimensional projection of a D-dimensional
+ *  Minkowski space. It vanishes whenever one of the indices is not in the
+ *  set {0, 1, 2, 3}.
+ *
+ *  @param i1 First index
+ *  @param i2 Second index
+ *  @param i3 Third index
+ *  @param i4 Fourth index
+ *  @param pos_sig Whether the signature of the metric is positive
+ *  @return newly constructed epsilon tensor */
+ex eps0123(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig = false);
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_TENSOR_H__
-- 
2.44.0