/** @file color.cpp
*
- * Implementation of GiNaC's color objects.
- * No real implementation yet, to be done. */
+ * Implementation of GiNaC's color (SU(3) Lie algebra) objects. */
/*
- * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
-#include <string>
-#include <list>
#include <algorithm>
-#include <iostream>
#include <stdexcept>
#include "color.h"
#include "ex.h"
-#include "coloridx.h"
+#include "idx.h"
#include "ncmul.h"
#include "numeric.h"
-#include "relational.h"
+#include "power.h" // for sqrt()
+#include "symbol.h"
+#include "print.h"
+#include "archive.h"
+#include "debugmsg.h"
+#include "utils.h"
+
+namespace GiNaC {
+
+GINAC_IMPLEMENT_REGISTERED_CLASS(color, indexed)
+GINAC_IMPLEMENT_REGISTERED_CLASS(su3one, tensor)
+GINAC_IMPLEMENT_REGISTERED_CLASS(su3t, tensor)
+GINAC_IMPLEMENT_REGISTERED_CLASS(su3f, tensor)
+GINAC_IMPLEMENT_REGISTERED_CLASS(su3d, tensor)
//////////
// default constructor, destructor, copy constructor assignment operator and helpers
//////////
-// public
-
-color::color() : type(invalid), representation_label(0)
+color::color() : representation_label(0)
{
- debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
- tinfo_key=TINFO_color;
+ debugmsg("color default constructor", LOGLEVEL_CONSTRUCT);
+ tinfo_key = TINFO_color;
}
-color::~color()
+void color::copy(const color & other)
{
- debugmsg("color destructor",LOGLEVEL_DESTRUCT);
- destroy(0);
+ inherited::copy(other);
+ representation_label = other.representation_label;
}
-color::color(color const & other)
-{
- debugmsg("color copy constructor",LOGLEVEL_CONSTRUCT);
- copy (other);
-}
-
-color const & color::operator=(color const & other)
-{
- debugmsg("color operator=",LOGLEVEL_ASSIGNMENT);
- if (this != &other) {
- destroy(1);
- copy(other);
- }
- return *this;
-}
-
-// protected
-
-void color::copy(color const & other)
-{
- indexed::copy(other);
- type=other.type;
- representation_label=other.representation_label;
-}
-
-void color::destroy(bool call_parent)
-{
- if (call_parent) {
- indexed::destroy(call_parent);
- }
-}
+DEFAULT_DESTROY(color)
+DEFAULT_CTORS(su3one)
+DEFAULT_CTORS(su3t)
+DEFAULT_CTORS(su3f)
+DEFAULT_CTORS(su3d)
//////////
// other constructors
//////////
-// protected
-
-color::color(color_types const t, unsigned const rl) : type(t), representation_label(rl)
+/** Construct object without any color index. This constructor is for
+ * internal use only. Use the color_ONE() function instead.
+ * @see color_ONE */
+color::color(const ex & b, unsigned char rl) : inherited(b), representation_label(rl)
{
- debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
+ debugmsg("color constructor from ex,unsigned char", LOGLEVEL_CONSTRUCT);
+ tinfo_key = TINFO_color;
}
-color::color(color_types const t, ex const & i1, unsigned const rl)
- : indexed(i1), type(t), representation_label(rl)
+/** Construct object with one color index. This constructor is for internal
+ * use only. Use the color_T() function instead.
+ * @see color_T */
+color::color(const ex & b, const ex & i1, unsigned char rl) : inherited(b, i1), representation_label(rl)
{
- debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
+ debugmsg("color constructor from ex,ex,unsigned char", LOGLEVEL_CONSTRUCT);
+ tinfo_key = TINFO_color;
}
-color::color(color_types const t, ex const & i1, ex const & i2, unsigned const rl)
- : indexed(i1,i2), type(t), representation_label(rl)
+color::color(unsigned char rl, const exvector & v, bool discardable) : inherited(indexed::unknown, v, discardable), representation_label(rl)
{
- debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
+ debugmsg("color constructor from unsigned char,exvector", LOGLEVEL_CONSTRUCT);
+ tinfo_key = TINFO_color;
}
-color::color(color_types const t, ex const & i1, ex const & i2, ex const & i3,
- unsigned const rl) : indexed(i1,i2,i3), type(t), representation_label(rl)
+color::color(unsigned char rl, exvector * vp) : inherited(indexed::unknown, vp), representation_label(rl)
{
- debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
+ debugmsg("color constructor from unsigned char,exvector *", LOGLEVEL_CONSTRUCT);
+ tinfo_key = TINFO_color;
}
-color::color(color_types const t, exvector const & iv, unsigned const rl)
- : indexed(iv), type(t), representation_label(rl)
+//////////
+// archiving
+//////////
+
+color::color(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
{
- debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
+ debugmsg("color constructor from archive_node", LOGLEVEL_CONSTRUCT);
+ unsigned rl;
+ n.find_unsigned("label", rl);
+ representation_label = rl;
}
-color::color(color_types const t, exvector * ivp, unsigned const rl)
- : indexed(ivp), type(t), representation_label(rl)
+void color::archive(archive_node &n) const
{
- debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
+ inherited::archive(n);
+ n.add_unsigned("label", representation_label);
}
+DEFAULT_UNARCHIVE(color)
+DEFAULT_ARCHIVING(su3one)
+DEFAULT_ARCHIVING(su3t)
+DEFAULT_ARCHIVING(su3f)
+DEFAULT_ARCHIVING(su3d)
+
//////////
// functions overriding virtual functions from bases classes
//////////
-// public
-
-basic * color::duplicate() const
+int color::compare_same_type(const basic & other) const
{
- debugmsg("color duplicate",LOGLEVEL_DUPLICATE);
- return new color(*this);
-}
+ GINAC_ASSERT(other.tinfo() == TINFO_color);
+ const color &o = static_cast<const color &>(other);
-void color::printraw(ostream & os) const
-{
- debugmsg("color printraw",LOGLEVEL_PRINT);
- os << "color(type=" << (unsigned)type
- << ",representation_label=" << representation_label
- << ",indices=";
- printrawindices(os);
- os << ",hash=" << hashvalue << ",flags=" << flags << ")";
+ if (representation_label != o.representation_label) {
+ // different representation label
+ return representation_label < o.representation_label ? -1 : 1;
+ }
+
+ return inherited::compare_same_type(other);
}
-void color::printtree(ostream & os, unsigned indent) const
+DEFAULT_COMPARE(su3one)
+DEFAULT_COMPARE(su3t)
+DEFAULT_COMPARE(su3f)
+DEFAULT_COMPARE(su3d)
+
+DEFAULT_PRINT_LATEX(su3one, "ONE", "\\mathbb{1}")
+DEFAULT_PRINT(su3t, "T")
+DEFAULT_PRINT(su3f, "f")
+DEFAULT_PRINT(su3d, "d")
+
+/** Perform automatic simplification on noncommutative product of color
+ * objects. This removes superfluous ONEs. */
+ex color::simplify_ncmul(const exvector & v) const
{
- debugmsg("color printtree",LOGLEVEL_PRINT);
- os << string(indent,' ') << "color object: "
- << "type=" << (unsigned)type
- << ",representation_label=" << representation_label << ", ";
- os << seq.size() << " indices" << endl;
- printtreeindices(os,indent);
- os << string(indent,' ') << "hash=" << hashvalue
- << " (0x" << hex << hashvalue << dec << ")"
- << ", flags=" << flags << endl;
+ //!! TODO: sort by representation label
+ exvector s;
+ s.reserve(v.size());
+
+ exvector::const_iterator it = v.begin(), itend = v.end();
+ while (it != itend) {
+ if (!is_ex_of_type(it->op(0), su3one))
+ s.push_back(*it);
+ it++;
+ }
+
+ if (s.size() == 0)
+ return color(su3one());
+ else if (s.size() == v.size())
+ return simplified_ncmul(v);
+ else
+ return simplified_ncmul(s);
}
-void color::print(ostream & os, unsigned upper_precedence) const
+ex color::thisexprseq(const exvector & v) const
{
- debugmsg("color print",LOGLEVEL_PRINT);
- switch (type) {
- case color_T:
- os << "T";
- if (representation_label!=0) {
- os << "^(" << representation_label << ")";
- }
- break;
- case color_f:
- os << "f";
- break;
- case color_d:
- os << "d";
- break;
- case color_delta8:
- os << "delta8";
- break;
- case color_ONE:
- os << "color_ONE";
- break;
- case invalid:
- default:
- os << "INVALID_COLOR_OBJECT";
- break;
- }
- printindices(os);
+ return color(representation_label, v);
}
-void color::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+ex color::thisexprseq(exvector * vp) const
{
- debugmsg("color print csrc",LOGLEVEL_PRINT);
- print(os,upper_precedence);
+ return color(representation_label, vp);
}
-bool color::info(unsigned inf) const
+/** Given a vector iv3 of three indices and a vector iv2 of two indices that
+ * is a subset of iv3, return the (free) index that is in iv3 but not in
+ * iv2 and the sign introduced by permuting that index to the front.
+ *
+ * @param iv3 Vector of 3 indices
+ * @param iv2 Vector of 2 indices, must be a subset of iv3
+ * @param sig Returs sign introduced by index permutation
+ * @return the free index (the one that is in iv3 but not in iv2) */
+static ex permute_free_index_to_front(const exvector & iv3, const exvector & iv2, int & sig)
{
- return indexed::info(inf);
-}
+ GINAC_ASSERT(iv3.size() == 3);
+ GINAC_ASSERT(iv2.size() == 2);
-#define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
+ sig = 1;
-ex color::eval(int level) const
-{
- // canonicalize indices
-
- bool antisymmetric=false;
-
- switch (type) {
- case color_f:
- antisymmetric=true; // no break here!
- case color_d:
- case color_delta8:
- {
- exvector iv=seq;
- int sig=canonicalize_indices(iv,antisymmetric);
- if (sig!=INT_MAX) {
- // something has changed while sorting indices, more evaluations later
- if (sig==0) return exZERO();
- return ex(sig)*color(type,iv,representation_label);
- }
- }
- break;
- default:
- // nothing to canonicalize
- break;
- }
-
- switch (type) {
- case color_delta8:
- {
- ASSERT(seq.size()==2);
- coloridx const & idx1=ex_to_coloridx(seq[0]);
- coloridx const & idx2=ex_to_coloridx(seq[1]);
-
- // check for delta8_{a,a} where a is a symbolic index, replace by 8
- if ((idx1.is_symbolic())&&(idx1.is_equal_same_type(idx2))) {
- return ex(COLOR_EIGHT);
- }
-
- // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1
- if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
- if ((idx1.get_value()!=idx2.get_value())) {
- return exONE();
- } else {
- return exZERO();
- }
- }
+#define TEST_PERMUTATION(A,B,C,P) \
+ if (iv3[B].is_equal(iv2[0]) && iv3[C].is_equal(iv2[1])) { \
+ sig = P; \
+ return iv3[A]; \
+ }
+
+ TEST_PERMUTATION(0,1,2, 1);
+ TEST_PERMUTATION(0,2,1, -1);
+ TEST_PERMUTATION(1,0,2, -1);
+ TEST_PERMUTATION(1,2,0, 1);
+ TEST_PERMUTATION(2,0,1, 1);
+ TEST_PERMUTATION(2,1,0, -1);
+
+ throw(std::logic_error("permute_free_index_to_front(): no valid permutation found"));
+}
+
+/** Automatic symbolic evaluation of indexed symmetric structure constant. */
+ex su3d::eval_indexed(const basic & i) const
+{
+ GINAC_ASSERT(is_of_type(i, indexed));
+ GINAC_ASSERT(i.nops() == 4);
+ GINAC_ASSERT(is_ex_of_type(i.op(0), su3d));
+
+ // Convolutions are zero
+ if (static_cast<const indexed &>(i).get_dummy_indices().size() != 0)
+ return _ex0();
+
+ // Numeric evaluation
+ if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
+
+ // Sort indices
+ int v[3];
+ for (unsigned j=0; j<3; j++)
+ v[j] = ex_to_numeric(ex_to_idx(i.op(j + 1)).get_value()).to_int();
+ if (v[0] > v[1]) std::swap(v[0], v[1]);
+ if (v[0] > v[2]) std::swap(v[0], v[2]);
+ if (v[1] > v[2]) std::swap(v[1], v[2]);
+
+#define CMPINDICES(A,B,C) ((v[0] == (A)) && (v[1] == (B)) && (v[2] == (C)))
+
+ // Check for non-zero elements
+ if (CMPINDICES(1,4,6) || CMPINDICES(1,5,7) || CMPINDICES(2,5,6)
+ || CMPINDICES(3,4,4) || CMPINDICES(3,5,5))
+ return _ex1_2();
+ else if (CMPINDICES(2,4,7) || CMPINDICES(3,6,6) || CMPINDICES(3,7,7))
+ return _ex_1_2();
+ else if (CMPINDICES(1,1,8) || CMPINDICES(2,2,8) || CMPINDICES(3,3,8))
+ return sqrt(_ex3())/3;
+ else if (CMPINDICES(8,8,8))
+ return -sqrt(_ex3())/3;
+ else if (CMPINDICES(4,4,8) || CMPINDICES(5,5,8)
+ || CMPINDICES(6,6,8) || CMPINDICES(7,7,8))
+ return -sqrt(_ex3())/6;
+ else
+ return _ex0();
}
- break;
- case color_d:
- // check for d_{a,a,c} (=0) when a is symbolic
- {
- ASSERT(seq.size()==3);
- coloridx const & idx1=ex_to_coloridx(seq[0]);
- coloridx const & idx2=ex_to_coloridx(seq[1]);
- coloridx const & idx3=ex_to_coloridx(seq[2]);
-
- if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
- return exZERO();
- } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
- return exZERO();
- }
-
- // check for three numeric indices
- if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
- ASSERT(idx1.get_value()<=idx2.get_value());
- ASSERT(idx2.get_value()<=idx3.get_value());
- if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)||
- CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) {
- return exHALF();
- } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
- return -exHALF();
- } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) {
- return 1/sqrt(numeric(3));
- } else if (CMPINDICES(8,8,8)) {
- return -1/sqrt(numeric(3));
- } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
- return -1/(2*sqrt(numeric(3)));
- }
- return exZERO();
- }
- }
- break;
- case color_f:
- {
- ASSERT(seq.size()==3);
- coloridx const & idx1=ex_to_coloridx(seq[0]);
- coloridx const & idx2=ex_to_coloridx(seq[1]);
- coloridx const & idx3=ex_to_coloridx(seq[2]);
-
- // check for three numeric indices
- if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
- ASSERT(idx1.get_value()<=idx2.get_value());
- ASSERT(idx2.get_value()<=idx3.get_value());
- if (CMPINDICES(1,2,3)) {
- return exONE();
- } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
- CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
- return exHALF();
- } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
- return -exHALF();
- } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) {
- return sqrt(numeric(3))/2;
- } else if (CMPINDICES(8,8,8)) {
- return -1/sqrt(numeric(3));
- } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
- return -1/(2*sqrt(numeric(3)));
- }
- return exZERO();
- }
- break;
- }
- default:
- // nothing to evaluate
- break;
- }
-
- return this->hold();
-}
-
-// protected
-int color::compare_same_type(basic const & other) const
-{
- ASSERT(other.tinfo() == TINFO_color);
- const color *o = static_cast<const color *>(&other);
- if (type==o->type) {
- if (representation_label==o->representation_label) {
- return indexed::compare_same_type(other);
- }
- return representation_label < o->representation_label ? -1 : 1;
- }
- return type < o->type ? -1 : 1;
-}
+ // No further simplifications
+ return i.hold();
+}
+
+/** Automatic symbolic evaluation of indexed antisymmetric structure constant. */
+ex su3f::eval_indexed(const basic & i) const
+{
+ GINAC_ASSERT(is_of_type(i, indexed));
+ GINAC_ASSERT(i.nops() == 4);
+ GINAC_ASSERT(is_ex_of_type(i.op(0), su3f));
+
+ // Numeric evaluation
+ if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
+
+ // Sort indices, remember permutation sign
+ int v[3];
+ for (unsigned j=0; j<3; j++)
+ v[j] = ex_to_numeric(ex_to_idx(i.op(j + 1)).get_value()).to_int();
+ int sign = 1;
+ if (v[0] > v[1]) { std::swap(v[0], v[1]); sign = -sign; }
+ if (v[0] > v[2]) { std::swap(v[0], v[2]); sign = -sign; }
+ if (v[1] > v[2]) { std::swap(v[1], v[2]); sign = -sign; }
+
+ // Check for non-zero elements
+ if (CMPINDICES(1,2,3))
+ return sign;
+ else if (CMPINDICES(1,4,7) || CMPINDICES(2,4,6)
+ || CMPINDICES(2,5,7) || CMPINDICES(3,4,5))
+ return _ex1_2() * sign;
+ else if (CMPINDICES(1,5,6) || CMPINDICES(3,6,7))
+ return _ex_1_2() * sign;
+ else if (CMPINDICES(4,5,8) || CMPINDICES(6,7,8))
+ return sqrt(_ex3())/2 * sign;
+ else
+ return _ex0();
+ }
-bool color::is_equal_same_type(basic const & other) const
-{
- ASSERT(other.tinfo() == TINFO_color);
- const color *o = static_cast<const color *>(&other);
- if (type!=o->type) return false;
- if (representation_label!=o->representation_label) return false;
- return indexed::is_equal_same_type(other);
+ // No further simplifications
+ return i.hold();
}
-#include <iostream>
-ex color::simplify_ncmul(exvector const & v) const
+/** Contraction of an indexed symmetric structure constant with something else. */
+bool su3d::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
{
- // simplifications: contract delta8_{a,b} where possible
- // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...
- // remove superfluous ONEs
-
- // contract indices of delta8_{a,b} if they are different and symbolic
-
- exvector v_contracted=v;
- unsigned replacements;
- bool something_changed=false;
-
- exvector::iterator it=v_contracted.begin();
- while (it!=v_contracted.end()) {
- // process only delta8 objects
- if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) {
- color & d8=ex_to_nonconst_color(*it);
- ASSERT(d8.seq.size()==2);
- coloridx const & first_idx=ex_to_coloridx(d8.seq[0]);
- coloridx const & second_idx=ex_to_coloridx(d8.seq[1]);
- // delta8_{a,a} should have been contracted in color::eval()
- ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic()));
- ex saved_delta8=*it; // save to restore it later
-
- // try to contract first index
- replacements=1;
- if (first_idx.is_symbolic()) {
- replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx);
- if (replacements==1) {
- // not contracted except in itself, restore delta8 object
- *it=saved_delta8;
- } else {
- // a contracted index should occur exactly twice
- ASSERT(replacements==2);
- *it=exONE();
- something_changed=true;
- }
- }
-
- // try second index only if first was not contracted
- if ((replacements==1)&&(second_idx.is_symbolic())) {
- // first index not contracted, *it is guaranteed to be the original delta8 object
- replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx);
- if (replacements==1) {
- // not contracted except in itself, restore delta8 object
- *it=saved_delta8;
- } else {
- // a contracted index should occur exactly twice
- ASSERT(replacements==2);
- *it=exONE();
- something_changed=true;
- }
- }
- }
- ++it;
- }
-
- if (something_changed) {
- // do more simplifications later
- return nonsimplified_ncmul(v_contracted);
- }
-
- // there were no indices to contract
- // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown
- // (if there is at least one unknown object, all Ts will be unknown to not change the order)
-
- exvector delta8vec;
- exvector fvec;
- exvector dvec;
- vector<exvector> Tvecs;
- Tvecs.resize(MAX_REPRESENTATION_LABELS);
- vector<exvector> ONEvecs;
- ONEvecs.resize(MAX_REPRESENTATION_LABELS);
- exvector unknownvec;
-
- split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
-
- // d_{a,k,l} f_{b,k,l}=0 (includes case a=b)
- if ((dvec.size()>=1)&&(fvec.size()>=1)) {
- for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) {
- for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it1,color));
- ASSERT(is_ex_exactly_of_type(*it2,color));
- color const & col1=ex_to_color(*it1);
- color const & col2=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
- if (iv_intersect.size()>=2) return exZERO();
- }
- }
- }
-
- // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b)
- if (dvec.size()>=2) {
- for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) {
- for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it1,color));
- ASSERT(is_ex_exactly_of_type(*it2,color));
- color const & col1=ex_to_color(*it1);
- color const & col2=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
- if (iv_intersect.size()>=2) {
- if (iv_intersect.size()==3) {
- *it1=numeric(40)/numeric(3);
- *it2=exONE();
- } else {
- int sig1, sig2; // unimportant, since symmetric
- ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,false,&sig1);
- ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,false,&sig2);
- *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
- *it2=exONE();
- }
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
- }
- }
-
- // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b)
- if (fvec.size()>=2) {
- for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) {
- for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it1,color));
- ASSERT(is_ex_exactly_of_type(*it2,color));
- color const & col1=ex_to_color(*it1);
- color const & col2=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
- if (iv_intersect.size()>=2) {
- if (iv_intersect.size()==3) {
- *it1=numeric(24);
- *it2=exONE();
- } else {
- int sig1, sig2;
- ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,true,&sig1);
- ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,true,&sig2);
- *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
- *it2=exONE();
- }
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
- }
- }
-
- // d_{a,b,c} T_b T_c = 5/6 T_a
- // f_{a,b,c} T_b T_c = 3/2 I T_a
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- if ((Tvecs[rl].size()>=2)&&((dvec.size()>=1)||(fvec.size()>=1))) {
- for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) {
- exvector iv;
- ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
- ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T);
- iv.push_back(ex_to_color(*it1).seq[0]);
- iv.push_back(ex_to_color(*(it1+1)).seq[0]);
-
- // d_{a,b,c} T_b T_c = 5/6 T_a
- for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
- color const & dref=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(dref.seq,iv);
- if (iv_intersect.size()==2) {
- int sig; // unimportant, since symmetric
- ex free_idx=permute_free_index_to_front(dref.seq,iv,false,&sig);
- *it1=color(color_T,free_idx,rl);
- *(it1+1)=color(color_ONE,rl);
- *it2=numeric(5)/numeric(6);
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
-
- // f_{a,b,c} T_b T_c = 3/2 I T_a
- for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
- color const & fref=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(fref.seq,iv);
- if (iv_intersect.size()==2) {
- int sig;
- ex free_idx=permute_free_index_to_front(fref.seq,iv,true,&sig);
- *it1=color(color_T,free_idx,rl);
- *(it1+1)=color(color_ONE,rl);
- *it2=numeric(sig*3)/numeric(2)*I;
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
- }
- }
- }
-
- // clear all ONEs when there is at least one corresponding color_T
- // in this representation, retain one ONE otherwise
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- if (Tvecs[rl].size()!=0) {
- ONEvecs[rl].clear();
- } else if (ONEvecs[rl].size()!=0) {
- ONEvecs[rl].clear();
- ONEvecs[rl].push_back(color(color_ONE,rl));
- }
- }
-
- // return a sorted vector
- return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
- ONEvecs,unknownvec));
-}
+ GINAC_ASSERT(is_ex_of_type(*self, indexed));
+ GINAC_ASSERT(is_ex_of_type(*other, indexed));
+ GINAC_ASSERT(self->nops() == 4);
+ GINAC_ASSERT(is_ex_of_type(self->op(0), su3d));
-ex color::thisexprseq(exvector const & v) const
-{
- return color(type,v,representation_label);
-}
+ if (is_ex_exactly_of_type(other->op(0), su3d)) {
-ex color::thisexprseq(exvector * vp) const
-{
- return color(type,vp,representation_label);
-}
+ // Find the dummy indices of the contraction
+ exvector dummy_indices;
+ dummy_indices = ex_to_indexed(*self).get_dummy_indices(ex_to_indexed(*other));
-bool color::all_of_type_coloridx(void) const
-{
- // used only inside of ASSERTs
- for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
- if (!is_ex_of_type(*cit,coloridx)) return false;
- }
- return true;
-}
+ // d.abc*d.abc=40/3
+ if (dummy_indices.size() == 3) {
+ *self = numeric(40, 3);
+ *other = _ex1();
+ return true;
-//////////
-// virtual functions which can be overridden by derived classes
-//////////
+ // d.akl*d.bkl=5/3*delta.ab
+ } else if (dummy_indices.size() == 2) {
+ exvector a = index_set_difference(ex_to_indexed(*self).get_indices(), dummy_indices);
+ exvector b = index_set_difference(ex_to_indexed(*other).get_indices(), dummy_indices);
+ GINAC_ASSERT(a.size() > 0);
+ GINAC_ASSERT(b.size() > 0);
+ *self = numeric(5, 3) * delta_tensor(a[0], b[0]);
+ *other = _ex1();
+ return true;
+ }
+ }
-// none
+ return false;
+}
-//////////
-// non-virtual functions in this class
-//////////
+/** Contraction of an indexed antisymmetric structure constant with something else. */
+bool su3f::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+ GINAC_ASSERT(is_ex_of_type(*self, indexed));
+ GINAC_ASSERT(is_ex_of_type(*other, indexed));
+ GINAC_ASSERT(self->nops() == 4);
+ GINAC_ASSERT(is_ex_of_type(self->op(0), su3f));
-//////////
-// static member variables
-//////////
+ if (is_ex_exactly_of_type(other->op(0), su3f)) { // f*d is handled by su3d class
-// none
+ // Find the dummy indices of the contraction
+ exvector dummy_indices;
+ dummy_indices = ex_to_indexed(*self).get_dummy_indices(ex_to_indexed(*other));
-//////////
-// global constants
-//////////
+ // f.abc*f.abc=24
+ if (dummy_indices.size() == 3) {
+ *self = 24;
+ *other = _ex1();
+ return true;
-const color some_color;
-type_info const & typeid_color=typeid(some_color);
+ // f.akl*f.bkl=3*delta.ab
+ } else if (dummy_indices.size() == 2) {
+ int sign1, sign2;
+ ex a = permute_free_index_to_front(ex_to_indexed(*self).get_indices(), dummy_indices, sign1);
+ ex b = permute_free_index_to_front(ex_to_indexed(*other).get_indices(), dummy_indices, sign2);
+ *self = sign1 * sign2 * 3 * delta_tensor(a, b);
+ *other = _ex1();
+ return true;
+ }
+ }
+
+ return false;
+}
//////////
-// friend functions
+// global functions
//////////
-color color_ONE(unsigned const rl)
+ex color_ONE(unsigned char rl)
{
- return color(color::color_ONE,rl);
+ return color(su3one(), rl);
}
-color color_T(ex const & a, unsigned const rl)
+ex color_T(const ex & a, unsigned char rl)
{
- return color(color::color_T,a,rl);
-}
+ if (!is_ex_of_type(a, idx))
+ throw(std::invalid_argument("indices of color_T must be of type idx"));
+ if (!ex_to_idx(a).get_dim().is_equal(8))
+ throw(std::invalid_argument("index dimension for color_T must be 8"));
-color color_f(ex const & a, ex const & b, ex const & c)
-{
- return color(color::color_f,a,b,c);
+ return color(su3t(), a, rl);
}
-color color_d(ex const & a, ex const & b, ex const & c)
+ex color_f(const ex & a, const ex & b, const ex & c)
{
- return color(color::color_d,a,b,c);
-}
+ if (!is_ex_of_type(a, idx) || !is_ex_of_type(b, idx) || !is_ex_of_type(c, idx))
+ throw(std::invalid_argument("indices of color_f must be of type idx"));
+ if (!ex_to_idx(a).get_dim().is_equal(8) || !ex_to_idx(b).get_dim().is_equal(8) || !ex_to_idx(c).get_dim().is_equal(8))
+ throw(std::invalid_argument("index dimension for color_f must be 8"));
-ex color_h(ex const & a, ex const & b, ex const & c)
-{
- return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
+ return indexed(su3f(), indexed::antisymmetric, a, b, c);
}
-color color_delta8(ex const & a, ex const & b)
+ex color_d(const ex & a, const ex & b, const ex & c)
{
- return color(color::color_delta8,a,b);
-}
+ if (!is_ex_of_type(a, idx) || !is_ex_of_type(b, idx) || !is_ex_of_type(c, idx))
+ throw(std::invalid_argument("indices of color_d must be of type idx"));
+ if (!ex_to_idx(a).get_dim().is_equal(8) || !ex_to_idx(b).get_dim().is_equal(8) || !ex_to_idx(c).get_dim().is_equal(8))
+ throw(std::invalid_argument("index dimension for color_d must be 8"));
-void split_color_string_in_parts(exvector const & v, exvector & delta8vec,
- exvector & fvec, exvector & dvec,
- vector<exvector> & Tvecs,
- vector<exvector> & ONEvecs,
- exvector & unknownvec)
-{
- // if not all elements are of type color, put all Ts in unknownvec to
- // retain the ordering
- bool all_color=true;
- for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
- if (!is_ex_exactly_of_type(*cit,color)) {
- all_color=false;
- break;
- }
- }
-
- for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
- if (is_ex_exactly_of_type(*cit,color)) {
- switch (ex_to_color(*cit).type) {
- case color::color_delta8:
- delta8vec.push_back(*cit);
- break;
- case color::color_f:
- fvec.push_back(*cit);
- break;
- case color::color_d:
- dvec.push_back(*cit);
- break;
- case color::color_T:
- ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
- if (all_color) {
- Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
- } else {
- unknownvec.push_back(*cit);
- }
- break;
- case color::color_ONE:
- ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
- ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
- break;
- default:
- throw(std::logic_error("invalid type in split_color_string_in_parts()"));
- }
- } else {
- unknownvec.push_back(*cit);
- }
- }
-}
-
-exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
- exvector & dvec, vector<exvector> & Tvecs,
- vector<exvector> & ONEvecs, exvector & unknownvec)
-{
- unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- sz += Tvecs[rl].size();
- sz += ONEvecs[rl].size();
- }
- exvector v;
- v.reserve(sz);
-
- append_exvector_to_exvector(v,delta8vec);
- append_exvector_to_exvector(v,fvec);
- append_exvector_to_exvector(v,dvec);
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- append_exvector_to_exvector(v,Tvecs[rl]);
- append_exvector_to_exvector(v,ONEvecs[rl]);
- }
- append_exvector_to_exvector(v,unknownvec);
- return v;
+ return indexed(su3d(), indexed::symmetric, a, b, c);
}
-ex color_trace_of_one_representation_label(exvector const & v)
+ex color_h(const ex & a, const ex & b, const ex & c)
{
- if (v.size()==0) {
- return numeric(COLOR_THREE);
- } else if (v.size()==1) {
- ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
- return exZERO();
- }
- exvector v1=v;
- ex last_element=v1.back();
- ASSERT(is_ex_exactly_of_type(last_element,color));
- ASSERT(ex_to_color(last_element).type==color::color_T);
- v1.pop_back();
- ex next_to_last_element=v1.back();
- ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
- ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
- v1.pop_back();
- exvector v2=v1;
-
- ex const & last_index=ex_to_color(last_element).seq[0];
- ex const & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
- ex summation_index=coloridx();
-
- v2.push_back(color_T(summation_index)); // don't care about the representation_label
-
- // check this formula for SU(N) with N!=3 !!!!!!!!!
- return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
- % color_trace_of_one_representation_label(v1)
- +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
- % color_trace_of_one_representation_label(v2);
- /*
- ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
- % color_trace_of_one_representation_label(v1);
- cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
- ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
- % color_trace_of_one_representation_label(v2);
- cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
- return term1+term2;
- */
+ return color_d(a, b, c) + I * color_f(a, b, c);
}
-ex color_trace(exvector const & v, unsigned const rl)
+/** 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)
{
- ASSERT(rl<MAX_REPRESENTATION_LABELS);
-
- exvector v_rest;
- v_rest.reserve(v.size()+1); // max size if trace is empty
-
- exvector delta8vec;
- exvector fvec;
- exvector dvec;
- vector<exvector> Tvecs;
- Tvecs.resize(MAX_REPRESENTATION_LABELS);
- vector<exvector> ONEvecs;
- ONEvecs.resize(MAX_REPRESENTATION_LABELS);
- exvector unknownvec;
-
- split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
-
- if (unknownvec.size()!=0) {
- throw(std::invalid_argument("color_trace(): expression must be expanded"));
- }
-
- append_exvector_to_exvector(v_rest,delta8vec);
- append_exvector_to_exvector(v_rest,fvec);
- append_exvector_to_exvector(v_rest,dvec);
- for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
- if (i!=rl) {
- append_exvector_to_exvector(v_rest,Tvecs[i]);
- append_exvector_to_exvector(v_rest,ONEvecs[i]);
- } else {
- if (Tvecs[i].size()!=0) {
- v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
- } else if (ONEvecs[i].size()!=0) {
- v_rest.push_back(numeric(COLOR_THREE));
- } else {
- throw(std::logic_error("color_trace(): representation_label not in color string"));
- }
- }
- }
-
- return nonsimplified_ncmul(v_rest);
+ return ti == (TINFO_color + rl);
}
-ex simplify_pure_color_string(ex const & e)
-{
- ASSERT(is_ex_exactly_of_type(e,ncmul));
-
- exvector delta8vec;
- exvector fvec;
- exvector dvec;
- vector<exvector> Tvecs;
- Tvecs.resize(MAX_REPRESENTATION_LABELS);
- vector<exvector> ONEvecs;
- ONEvecs.resize(MAX_REPRESENTATION_LABELS);
- exvector unknownvec;
-
- split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
-
- // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- if (Tvecs[rl].size()>=2) {
- for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
- for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
- ex & t1=Tvecs[rl][i];
- ex & t2=Tvecs[rl][j];
- ASSERT(is_ex_exactly_of_type(t1,color)&&
- (ex_to_color(t1).type==color::color_T)&&
- (ex_to_color(t1).seq.size()==1));
- ASSERT(is_ex_exactly_of_type(t2,color)&&
- (ex_to_color(t2).type==color::color_T)&&
- (ex_to_color(t2).seq.size()==1));
- coloridx const & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
- coloridx const & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
-
- if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
- exvector S;
- for (unsigned k=i+1; k<j; ++k) {
- S.push_back(Tvecs[rl][k]);
- }
- t1=exONE();
- t2=exONE();
- ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- for (unsigned k=i+1; k<j; ++k) {
- S.push_back(exONE());
- }
- t1=color_trace_of_one_representation_label(S);
- ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- return simplify_color(term1+term2);
- }
- }
- }
- }
- }
-
- // TODO: higher contractions!!!!!!!!!!!!!
-
- return e;
-}
-
-ex simplify_color(ex const & e)
+ex color_trace(const ex & e, unsigned char rl)
{
- // all simplification is done on expanded objects
- ex e_expanded=e.expand();
-
- // simplification of sum=sum of simplifications
- if (is_ex_exactly_of_type(e_expanded,add)) {
- ex sum=exZERO();
- for (int i=0; i<e_expanded.nops(); ++i) {
- sum += simplify_color(e_expanded.op(i));
- }
- return sum;
- }
-
- // simplification of commutative product=commutative product of simplifications
- if (is_ex_exactly_of_type(e_expanded,mul)) {
- ex prod=exONE();
- for (int i=0; i<e_expanded.nops(); ++i) {
- prod *= simplify_color(e_expanded.op(i));
- }
- return prod;
- }
-
- // simplification of noncommutative product: test if everything is color
- if (is_ex_exactly_of_type(e_expanded,ncmul)) {
- bool all_color=true;
- for (int i=0; i<e_expanded.nops(); ++i) {
- if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
- all_color=false;
- break;
- }
- }
- if (all_color) {
- return simplify_pure_color_string(e_expanded);
- }
- }
-
- // cannot do anything
- return e_expanded;
-}
+ if (is_ex_of_type(e, color)) {
-ex brute_force_sum_color_indices(ex const & e)
-{
- exvector iv_all=e.get_indices();
- exvector iv_double;
-
- // find double symbolic indices
- if (iv_all.size()<2) return e;
- for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
- ASSERT(is_ex_of_type(*cit1,coloridx));
- for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
- ASSERT(is_ex_of_type(*cit2,coloridx));
- if (ex_to_coloridx(*cit1).is_symbolic() &&
- ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
- iv_double.push_back(*cit1);
- break;
- }
- }
- }
-
- vector<int> counter;
- counter.resize(iv_double.size());
- int l;
- for (l=0; unsigned(l)<iv_double.size(); ++l) {
- counter[l]=1;
- }
-
- ex sum=exZERO();
-
- while (1) {
- ex term=e;
- for (l=0; unsigned(l)<iv_double.size(); ++l) {
- term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
- //iv_double[l].print(cout);
- //cout << " " << counter[l] << " ";
- }
- //cout << endl;
- sum += term;
-
- // increment counter[]
- l=iv_double.size()-1;
- while ((l>=0)&&((++counter[l])>COLOR_EIGHT)) {
- counter[l]=1;
- l--;
- }
- if (l<2) { cout << counter[0] << counter[1] << endl; }
- if (l<0) break;
- }
-
- return sum;
-}
+ if (ex_to_color(e).get_representation_label() == rl
+ && is_ex_of_type(e.op(0), su3one))
+ return _ex3();
+ else
+ return _ex0();
-void append_exvector_to_exvector(exvector & dest, exvector const & source)
-{
- for (exvector::const_iterator cit=source.begin(); cit!=source.end(); ++cit) {
- dest.push_back(*cit);
- }
-}
+ } 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
git://www.ginac.de / ginac.git / blobdiff