+/** Check whether a given tinfo key (as returned by return_type_tinfo()
+ * is that of a color object (with an arbitrary representation label). */
+static bool is_color_tinfo(tinfo_t ti)
+{
+ p_int start_loc=(p_int)&color::return_type_tinfo_static;
+ return (p_int)ti>=start_loc && (p_int)ti<start_loc+256;
+}
+
+/** Extract representation label from tinfo key (as returned by
+ * return_type_tinfo()). */
+static unsigned char get_representation_label(tinfo_t ti)
+{
+ return (unsigned char)((p_int)ti-(p_int)&color::return_type_tinfo_static);
+}
+
+ex color_trace(const ex & e, const std::set<unsigned char> & rls)
+{
+ if (is_a<color>(e)) {
+
+ unsigned char rl = ex_to<color>(e).get_representation_label();
+
+ // Are we taking the trace over this object's representation label?
+ if (rls.find(rl) == rls.end())
+ return e;
+
+ // Yes, all generators are traceless, except for color_ONE
+ if (is_a<su3one>(e.op(0)))
+ return _ex3;
+ else
+ return _ex0;
+
+ } else if (is_exactly_a<mul>(e)) {
+
+ // Trace of product: pull out non-color factors
+ ex prod = _ex1;
+ for (size_t i=0; i<e.nops(); i++) {
+ const ex &o = e.op(i);
+ if (is_color_tinfo(o.return_type_tinfo()))
+ prod *= color_trace(o, rls);
+ else
+ prod *= o;
+ }
+ return prod;
+
+ } else if (is_exactly_a<ncmul>(e)) {
+
+ unsigned char rl = get_representation_label(e.return_type_tinfo());
+
+ // Are we taking the trace over this string's representation label?
+ if (rls.find(rl) == rls.end())
+ return e;
+
+ // Yes, expand product if necessary
+ ex e_expanded = e.expand();
+ if (!is_a<ncmul>(e_expanded))
+ return color_trace(e_expanded, rls);
+
+ size_t 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 (size_t 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;
+ }
+
+ } else if (e.nops() > 0) {
+
+ // Trace maps to all other container classes (this includes sums)
+ pointer_to_map_function_1arg<const std::set<unsigned char> &> fcn(color_trace, rls);
+ return e.map(fcn);
+
+ } else
+ return _ex0;
+}
+
+ex color_trace(const ex & e, const lst & rll)
+{
+ // Convert list to set
+ std::set<unsigned char> rls;
+ for (lst::const_iterator i = rll.begin(); i != rll.end(); ++i) {
+ if (i->info(info_flags::nonnegint))
+ rls.insert(ex_to<numeric>(*i).to_int());
+ }
+
+ return color_trace(e, rls);
+}
+
+ex color_trace(const ex & e, unsigned char rl)
+{
+ // Convert label to set
+ std::set<unsigned char> rls;
+ rls.insert(rl);
+
+ return color_trace(e, rls);
+}
+