Add support for Texinfo-5.0.
[ginac.git] / ginac / color.cpp
index 9e355cb..34ab4d8 100644 (file)
@@ -1,10 +1,9 @@
 /** @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-2001 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2011 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
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  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 "symmetry.h"
+#include "operators.h"
 #include "numeric.h"
-#include "relational.h"
+#include "mul.h"
+#include "power.h" // for sqrt()
+#include "symbol.h"
 #include "archive.h"
-#include "debugmsg.h"
 #include "utils.h"
 
+#include <iostream>
+#include <stdexcept>
+
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(color, indexed)
 
-//////////
-// default constructor, destructor, copy constructor assignment operator and helpers
-//////////
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3one, tensor,
+  print_func<print_dflt>(&su3one::do_print).
+  print_func<print_latex>(&su3one::do_print_latex))
 
-// public
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3t, tensor,
+  print_func<print_dflt>(&su3t::do_print).
+  print_func<print_latex>(&su3t::do_print))
 
-color::color() : inherited(TINFO_color), type(invalid), representation_label(0)
-{
-       debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
-}
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3f, tensor,
+  print_func<print_dflt>(&su3f::do_print).
+  print_func<print_latex>(&su3f::do_print))
 
-// protected
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3d, tensor,
+  print_func<print_dflt>(&su3d::do_print).
+  print_func<print_latex>(&su3d::do_print))
 
-void color::copy(const color & other)
-{
-       inherited::copy(other);
-       type=other.type;
-       representation_label=other.representation_label;
-}
+//////////
+// default constructors
+//////////
 
-void color::destroy(bool call_parent)
+color::color() : representation_label(0)
 {
-       if (call_parent) {
-               inherited::destroy(call_parent);
-       }
 }
 
+DEFAULT_CTOR(su3one)
+DEFAULT_CTOR(su3t)
+DEFAULT_CTOR(su3f)
+DEFAULT_CTOR(su3d)
+
 //////////
 // other constructors
 //////////
 
-// protected
-
-/** Construct object without any color index. This constructor is for internal
- *  use only. Use the color_ONE() function instead.
+/** Construct object without any color index. This constructor is for
+ *  internal use only. Use the color_ONE() function instead.
  *  @see color_ONE */
-color::color(color_types const t, unsigned rl) : type(t), representation_label(rl)
+color::color(const ex & b, unsigned char rl) : inherited(b), representation_label(rl)
 {
-       debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT);
-       GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
-       tinfo_key=TINFO_color;
-       GINAC_ASSERT(all_of_type_coloridx());
 }
 
 /** Construct object with one color index. This constructor is for internal
  *  use only. Use the color_T() function instead.
  *  @see color_T */
-color::color(color_types const t, const ex & i1, unsigned rl)
-  : inherited(i1), type(t), representation_label(rl)
-{
-       debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
-       GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
-       tinfo_key=TINFO_color;
-       GINAC_ASSERT(all_of_type_coloridx());
-}
-
-/** Construct object with two color indices. This constructor is for internal
- *  use only. Use the color_delta8() function instead.
- *  @see color_delta8 */
-color::color(color_types const t, const ex & i1, const ex & i2, unsigned rl)
-  : inherited(i1,i2), type(t), representation_label(rl)
+color::color(const ex & b, const ex & i1, unsigned char rl) : inherited(b, i1), representation_label(rl)
 {
-       debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
-       GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
-       tinfo_key=TINFO_color;
-       GINAC_ASSERT(all_of_type_coloridx());
 }
 
-/** Construct object with three color indices. This constructor is for internal
- *  use only. Use the color_f(), color_d() and color_h() functions instead.
- *  @see color_f
- *  @see color_d
- *  @see color_h */
-color::color(color_types const t, const ex & i1, const ex & i2, const ex & i3,
-             unsigned rl) : inherited(i1,i2,i3), type(t), representation_label(rl)
+color::color(unsigned char rl, const exvector & v, bool discardable) : inherited(not_symmetric(), v, discardable), representation_label(rl)
 {
-       debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
-       GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
-       tinfo_key=TINFO_color;
-       GINAC_ASSERT(all_of_type_coloridx());
 }
 
-/** Construct object with arbitrary number of color indices. This
- *  constructor is for internal use only. */
-color::color(color_types const t, const exvector & iv, unsigned rl)
-  : inherited(iv), type(t), representation_label(rl)
+color::color(unsigned char rl, std::auto_ptr<exvector> vp) : inherited(not_symmetric(), vp), representation_label(rl)
 {
-       debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
-       GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
-       tinfo_key=TINFO_color;
-       GINAC_ASSERT(all_of_type_coloridx());
 }
 
-color::color(color_types const t, exvector * ivp, unsigned rl)
-  : inherited(ivp), type(t), representation_label(rl)
+return_type_t color::return_type_tinfo() const
 {
-       debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
-       GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
-       tinfo_key=TINFO_color;
-       GINAC_ASSERT(all_of_type_coloridx());
+       return make_return_type_t<color>(representation_label);
 }
 
 //////////
 // archiving
 //////////
 
-/** Construct object from archive_node. */
-color::color(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
-{
-       debugmsg("color constructor from archive_node", LOGLEVEL_CONSTRUCT);
-       unsigned int ty;
-       if (!(n.find_unsigned("type", ty)))
-               throw (std::runtime_error("unknown color type in archive"));
-       type = (color_types)ty;
-       if (!(n.find_unsigned("representation", representation_label)))
-               throw (std::runtime_error("unknown color representation label in archive"));
-}
-
-/** Unarchive the object. */
-ex color::unarchive(const archive_node &n, const lst &sym_lst)
+void color::read_archive(const archive_node& n, lst& sym_lst)
 {
-       return (new color(n, sym_lst))->setflag(status_flags::dynallocated);
+       inherited::read_archive(n, sym_lst);
+       unsigned rl;
+       n.find_unsigned("label", rl);
+       representation_label = rl;
 }
 
-/** Archive the object. */
 void color::archive(archive_node &n) const
 {
        inherited::archive(n);
-       n.add_unsigned("type", type);
-       n.add_unsigned("representation", representation_label);
+       n.add_unsigned("label", representation_label);
 }
 
+GINAC_BIND_UNARCHIVER(color);
+GINAC_BIND_UNARCHIVER(su3one);
+GINAC_BIND_UNARCHIVER(su3t);
+GINAC_BIND_UNARCHIVER(su3f);
+GINAC_BIND_UNARCHIVER(su3d);
+
 //////////
-// functions overriding virtual functions from bases classes
+// functions overriding virtual functions from base classes
 //////////
 
-// public
-
-void color::printraw(std::ostream & os) const
+int color::compare_same_type(const basic & other) const
 {
-       debugmsg("color printraw",LOGLEVEL_PRINT);
-       os << "color(type=" << (unsigned)type
-          << ",representation_label=" << representation_label
-          << ",indices=";
-       printrawindices(os);
-       os << ",hash=" << hashvalue << ",flags=" << flags << ")";
+       GINAC_ASSERT(is_a<color>(other));
+       const color &o = static_cast<const color &>(other);
+
+       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(std::ostream & os, unsigned indent) const
+bool color::match_same_type(const basic & other) const
 {
-       debugmsg("color printtree",LOGLEVEL_PRINT);
-       os << std::string(indent,' ') << "color object: "
-          << "type=" << (unsigned)type
-          << ",representation_label=" << representation_label << ", ";
-       os << seq.size() << " indices" << std::endl;
-       printtreeindices(os,indent);
-       os << std::string(indent,' ') << "hash=" << hashvalue
-          << " (0x" << std::hex << hashvalue << std::dec << ")"
-          << ", flags=" << flags << std::endl;
+       GINAC_ASSERT(is_a<color>(other));
+       const color &o = static_cast<const color &>(other);
+
+       return representation_label == o.representation_label;
 }
 
-void color::print(std::ostream & os, unsigned upper_precedence) 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::eval_ncmul(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;
+       exvector s;
+       s.reserve(v.size());
+
+       // Remove superfluous ONEs
+       exvector::const_iterator it = v.begin(), itend = v.end();
+       while (it != itend) {
+               if (!is_a<su3one>(it->op(0)))
+                       s.push_back(*it);
+               it++;
        }
-       printindices(os);
+
+       if (s.empty())
+               return color(su3one(), representation_label);
+       else
+               return hold_ncmul(s);
 }
 
-bool color::info(unsigned inf) const
+ex color::thiscontainer(const exvector & v) const
 {
-       return inherited::info(inf);
+       return color(representation_label, v);
 }
 
-#define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
+ex color::thiscontainer(std::auto_ptr<exvector> vp) const
+{
+       return color(representation_label, vp);
+}
 
-ex color::eval(int level) 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)
 {
-       // 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 _ex0();
-                               return ex(sig)*color(type,iv,representation_label);
-                       }
-               }
-               break;
-       default:
-               // nothing to canonicalize
-               break;
-       }
+       GINAC_ASSERT(iv3.size() == 3);
+       GINAC_ASSERT(iv2.size() == 2);
 
-       switch (type) {
-       case color_delta8:
-               {
-                       GINAC_ASSERT(seq.size()==2);
-                       const coloridx & idx1=ex_to_coloridx(seq[0]);
-                       const coloridx & 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);
-                       }
+       sig = 1;
 
-                       // 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 _ex1();
-                               } else {
-                                       return _ex0();
-                               }
-                       }
-       }
-               break;
-       case color_d:
-               // check for d_{a,a,c} (=0) when a is symbolic
-               {
-                       GINAC_ASSERT(seq.size()==3);
-                       const coloridx & idx1=ex_to_coloridx(seq[0]);
-                       const coloridx & idx2=ex_to_coloridx(seq[1]);
-                       const coloridx & idx3=ex_to_coloridx(seq[2]);
-                       
-                       if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
-                               return _ex0();
-                       } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
-                               return _ex0();
-                       }
-                       
-                       // check for three numeric indices
-                       if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
-                               GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
-                               GINAC_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 _ex1_2();
-                               } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
-                                       return -_ex1_2();
-                               } 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 _ex0();
-                       }
-               }
-               break;
-       case color_f:
-               {
-                       GINAC_ASSERT(seq.size()==3);
-                       const coloridx & idx1=ex_to_coloridx(seq[0]);
-                       const coloridx & idx2=ex_to_coloridx(seq[1]);
-                       const coloridx & idx3=ex_to_coloridx(seq[2]);
-                       
-                       // check for three numeric indices
-                       if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
-                               GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
-                               GINAC_ASSERT(idx2.get_value()<=idx3.get_value());
-                               if (CMPINDICES(1,2,3)) {
-                                       return _ex1();
-                               } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
-                                          CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
-                                       return _ex1_2();
-                               } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
-                                       return -_ex1_2();
-                               } 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 _ex0();
-                       }
-                       break;
-               }
-       default:
-               // nothing to evaluate
-               break;
+#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]; \
        }
        
-       return this->hold();
-}
-       
-// protected
-
-int color::compare_same_type(const basic & other) const
-{
-       GINAC_ASSERT(other.tinfo() == TINFO_color);
-       const color &o = static_cast<const color &>(other);
-
-       if (type != o.type) {
-               // different type
-               return type < o.type ? -1 : 1;
+       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_a<indexed>(i));
+       GINAC_ASSERT(i.nops() == 4);
+       GINAC_ASSERT(is_a<su3d>(i.op(0)));
+
+       // Convolutions are zero
+       if (!(static_cast<const indexed &>(i).get_dummy_indices().empty()))
+               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)*_ex1_3;
+               else if (CMPINDICES(8,8,8))
+                       return sqrt(_ex3)*_ex_1_3;
+               else if (CMPINDICES(4,4,8) || CMPINDICES(5,5,8)
+                     || CMPINDICES(6,6,8) || CMPINDICES(7,7,8))
+                       return sqrt(_ex3)/_ex_6;
+               else
+                       return _ex0;
        }
 
-       if (representation_label != o.representation_label) {
-               // different representation label
-               return representation_label < o.representation_label ? -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_a<indexed>(i));
+       GINAC_ASSERT(i.nops() == 4);
+       GINAC_ASSERT(is_a<su3f>(i.op(0)));
+
+       // 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;
        }
 
-       return inherited::compare_same_type(other);
+       // No further simplifications
+       return i.hold();
 }
 
-bool color::is_equal_same_type(const basic & other) const
-{
-       GINAC_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 inherited::is_equal_same_type(other);
-}
+/** Contraction of generator with something else. */
+bool su3t::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+       GINAC_ASSERT(is_a<indexed>(*self));
+       GINAC_ASSERT(is_a<indexed>(*other));
+       GINAC_ASSERT(self->nops() == 2);
+       GINAC_ASSERT(is_a<su3t>(self->op(0)));
+       unsigned char rl = ex_to<color>(*self).get_representation_label();
 
-#include <iostream>
+       if (is_exactly_a<su3t>(other->op(0))) {
 
-ex color::simplify_ncmul(const 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);
-                       GINAC_ASSERT(d8.seq.size()==2);
-                       const coloridx & first_idx=ex_to_coloridx(d8.seq[0]);
-                       const coloridx & second_idx=ex_to_coloridx(d8.seq[1]);
-                       // delta8_{a,a} should have been contracted in color::eval()
-                       GINAC_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
-                                       GINAC_ASSERT(replacements==2);
-                                       *it=_ex1();
-                                       something_changed=true;
-                               }
-                       }
+               // Contraction only makes sense if the represenation labels are equal
+               GINAC_ASSERT(is_a<color>(*other));
+               if (ex_to<color>(*other).get_representation_label() != rl)
+                       return false;
 
-                       // 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
-                                       GINAC_ASSERT(replacements==2);
-                                       *it=_ex1();
-                                       something_changed=true;
-                               }
-                       }
-               }
-               ++it;
-       }
+               // T.a T.a = 4/3 ONE
+               if (other - self == 1) {
+                       *self = numeric(4, 3);
+                       *other = color_ONE(rl);
+                       return true;
 
-       if (something_changed) {
-               // do more simplifications later
-               return nonsimplified_ncmul(v_contracted);
-       }
+               // T.a T.b T.a = -1/6 T.b
+               } else if (other - self == 2
+                       && is_a<color>(self[1])) {
+                       *self = numeric(-1, 6);
+                       *other = _ex1;
+                       return true;
 
-       // 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;
-       exvectorvector Tvecs;
-       Tvecs.resize(MAX_REPRESENTATION_LABELS);
-       exvectorvector 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) {
-                               GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
-                               GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
-                               const color & col1=ex_to_color(*it1);
-                               const color & col2=ex_to_color(*it2);
-                               exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
-                               if (iv_intersect.size()>=2) return _ex0();
-                       }
-               }
-       }
-       
-       // 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) {
-                               GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
-                               GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
-                               const color & col1=ex_to_color(*it1);
-                               const color & 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=_ex1();
-                                       } else {
-                                               int dummy; // sign unimportant, since symmetric
-                                               ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&dummy);
-                                               ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&dummy);
-                                               *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
-                                               *it2=_ex1();
-                                       }
-                                       return nonsimplified_ncmul(recombine_color_string(
-                                                  delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
+               // T.a S T.a = 1/2 Tr(S) - 1/6 S
+               } else {
+                       exvector::iterator it = self + 1;
+                       while (it != other) {
+                               if (!is_a<color>(*it)) {
+                                       return false;
                                }
+                               it++;
                        }
-               }
-       }
 
-       // 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) {
-                               GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
-                               GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
-                               const color & col1=ex_to_color(*it1);
-                               const color & 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=_ex1();
-                                       } else {
-                                               int sig1, sig2;
-                                               ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&sig1);
-                                               ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&sig2);
-                                               *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
-                                               *it2=_ex1();
-                                       }
-                                       return nonsimplified_ncmul(recombine_color_string(
-                                                  delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
-                               }
+                       it = self + 1;
+                       ex S = _ex1;
+                       while (it != other) {
+                               S *= *it;
+                               *it++ = _ex1;
                        }
+
+                       *self = color_trace(S, rl) * color_ONE(rl) / 2 - S / 6;
+                       *other = _ex1;
+                       return true;
                }
        }
 
-       // 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;
-                               GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
-                               GINAC_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) {
-                                       GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
-                                       const color & dref=ex_to_color(*it2);
-                                       exvector iv_intersect=idx_intersect(dref.seq,iv);
-                                       if (iv_intersect.size()==2) {
-                                               int dummy; // sign unimportant, since symmetric
-                                               ex free_idx=permute_free_index_to_front(dref.seq,iv,&dummy);
-                                               *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) {
-                                       GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
-                                       const color & 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,&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));
-                                       }
-                               }
-                       }
+       return false;
+}
+
+/** Contraction of an indexed symmetric structure constant with something else. */
+bool su3d::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+       GINAC_ASSERT(is_a<indexed>(*self));
+       GINAC_ASSERT(is_a<indexed>(*other));
+       GINAC_ASSERT(self->nops() == 4);
+       GINAC_ASSERT(is_a<su3d>(self->op(0)));
+
+       if (is_exactly_a<su3d>(other->op(0))) {
+
+               // Find the dummy indices of the contraction
+               exvector self_indices = ex_to<indexed>(*self).get_indices();
+               exvector other_indices = ex_to<indexed>(*other).get_indices();
+               exvector all_indices = self_indices;
+               all_indices.insert(all_indices.end(), other_indices.begin(), other_indices.end());
+               exvector free_indices, dummy_indices;
+               find_free_and_dummy(all_indices, free_indices, dummy_indices);
+
+               // d.abc d.abc = 40/3
+               if (dummy_indices.size() == 3) {
+                       *self = numeric(40, 3);
+                       *other = _ex1;
+                       return true;
+
+               // d.akl d.bkl = 5/3 delta.ab
+               } else if (dummy_indices.size() == 2) {
+                       exvector a;
+                       std::back_insert_iterator<exvector> ita(a);
+                       ita = set_difference(self_indices.begin(), self_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less());
+                       ita = set_difference(other_indices.begin(), other_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less());
+                       GINAC_ASSERT(a.size() == 2);
+                       *self = numeric(5, 3) * delta_tensor(a[0], a[1]);
+                       *other = _ex1;
+                       return true;
                }
-       }
-       
-       // 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));
+
+       } else if (is_exactly_a<su3t>(other->op(0))) {
+
+               // d.abc T.b T.c = 5/6 T.a
+               if (other+1 != v.end()
+                && is_exactly_a<su3t>(other[1].op(0))
+                && ex_to<indexed>(*self).has_dummy_index_for(other[1].op(1))) {
+
+                       exvector self_indices = ex_to<indexed>(*self).get_indices();
+                       exvector dummy_indices;
+                       dummy_indices.push_back(other[0].op(1));
+                       dummy_indices.push_back(other[1].op(1));
+                       int sig;
+                       ex a = permute_free_index_to_front(self_indices, dummy_indices, sig);
+                       *self = numeric(5, 6);
+                       other[0] = color_T(a, ex_to<color>(other[0]).get_representation_label());
+                       other[1] = _ex1;
+                       return true;
                }
        }
 
-       // return a sorted vector
-       return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
-                                                                                                  ONEvecs,unknownvec));
+       return false;
 }
 
-ex color::thisexprseq(const exvector & v) const
+/** Contraction of an indexed antisymmetric structure constant with something else. */
+bool su3f::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 {
-       return color(type,v,representation_label);
-}
+       GINAC_ASSERT(is_a<indexed>(*self));
+       GINAC_ASSERT(is_a<indexed>(*other));
+       GINAC_ASSERT(self->nops() == 4);
+       GINAC_ASSERT(is_a<su3f>(self->op(0)));
 
-ex color::thisexprseq(exvector * vp) const
-{
-       return color(type,vp,representation_label);
-}
+       if (is_exactly_a<su3f>(other->op(0))) { // f*d is handled by su3d class
 
-/** Check whether all indices are of class coloridx or a subclass. This
- *  function is used internally to make sure that all constructed color
- *  objects really carry color indices and not some other classes. */
-bool color::all_of_type_coloridx(void) const
-{
-       for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
-               if (!is_ex_of_type(*cit,coloridx)) return false;
+               // Find the dummy indices of the contraction
+               exvector dummy_indices;
+               dummy_indices = ex_to<indexed>(*self).get_dummy_indices(ex_to<indexed>(*other));
+
+               // f.abc f.abc = 24
+               if (dummy_indices.size() == 3) {
+                       *self = 24;
+                       *other = _ex1;
+                       return true;
+
+               // 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;
+               }
+
+       } else if (is_exactly_a<su3t>(other->op(0))) {
+
+               // f.abc T.b T.c = 3/2 I T.a
+               if (other+1 != v.end()
+                && is_exactly_a<su3t>(other[1].op(0))
+                && ex_to<indexed>(*self).has_dummy_index_for(other[1].op(1))) {
+
+                       exvector self_indices = ex_to<indexed>(*self).get_indices();
+                       exvector dummy_indices;
+                       dummy_indices.push_back(other[0].op(1));
+                       dummy_indices.push_back(other[1].op(1));
+                       int sig;
+                       ex a = permute_free_index_to_front(self_indices, dummy_indices, sig);
+                       *self = numeric(3, 2) * sig * I;
+                       other[0] = color_T(a, ex_to<color>(other[0]).get_representation_label());
+                       other[1] = _ex1;
+                       return true;
+               }
        }
-       return true;
+
+       return false;
 }
 
 //////////
-// friend functions
+// global functions
 //////////
 
-/** Construct an object representing the unity element of su(3).
- *
- *  @param rl Representation label
- *  @return newly constructed object */
-color color_ONE(unsigned rl)
+ex color_ONE(unsigned char rl)
 {
-       return color(color::color_ONE,rl);
+       static ex ONE = (new su3one)->setflag(status_flags::dynallocated);
+       return color(ONE, rl);
 }
 
-/** Construct an object representing the generators T_a of SU(3). The index
- *  must be of class coloridx.
- *
- *  @param a Index
- *  @param rl Representation label
- *  @return newly constructed object */
-color color_T(const ex & a, unsigned rl)
+ex color_T(const ex & a, unsigned char rl)
 {
-       return color(color::color_T,a,rl);
+       static ex t = (new su3t)->setflag(status_flags::dynallocated);
+
+       if (!is_a<idx>(a))
+               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"));
+
+       return color(t, a, rl);
 }
 
-/** Construct an object representing the antisymmetric structure constants
- *  f_abc of SU(3). The indices must be of class coloridx.
- *
- *  @param a First index
- *  @param b Second index
- *  @param c Third index
- *  @return newly constructed object */
-color color_f(const ex & a, const ex & b, const ex & c)
+ex color_f(const ex & a, const ex & b, const ex & c)
 {
-       return color(color::color_f,a,b,c);
+       static ex f = (new su3f)->setflag(status_flags::dynallocated);
+
+       if (!is_a<idx>(a) || !is_a<idx>(b) || !is_a<idx>(c))
+               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"));
+
+       return indexed(f, antisymmetric3(), a, b, c);
 }
 
-/** Construct an object representing the symmetric structure constants d_abc
- *  of SU(3). The indices must be of class coloridx.
- *
- *  @param a First index
- *  @param b Second index
- *  @param c Third index
- *  @return newly constructed object */
-color color_d(const ex & a, const ex & b, const ex & c)
+ex color_d(const ex & a, const ex & b, const ex & c)
 {
-       return color(color::color_d,a,b,c);
+       static ex d = (new su3d)->setflag(status_flags::dynallocated);
+
+       if (!is_a<idx>(a) || !is_a<idx>(b) || !is_a<idx>(c))
+               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"));
+
+       return indexed(d, symmetric3(), a, b, c);
 }
 
-/** This returns the linear combination d_abc+I*f_abc.
- *
- *  @param a First index
- *  @param b Second index
- *  @param c Third index
- *  @return newly constructed object */
 ex color_h(const ex & a, const ex & b, const ex & c)
 {
-       return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
+       return color_d(a, b, c) + I * color_f(a, b, c);
 }
 
-/** Construct an object representing the unity matrix delta8_ab in su(3).
- *  The indices must be of class coloridx.
- *
- *  @param a First index
- *  @param b Second index
- *  @return newly constructed object */
-color color_delta8(const ex & a, const ex & b)
+/** 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(const return_type_t& ti)
 {
-       return color(color::color_delta8,a,b);
+       return *(ti.tinfo) == typeid(color);
 }
 
-/** Given a vector of color (and possible other) objects, split it up
- *  according to the object type (structure constant, generator etc.) and
- *  representation label while preserving the order within each group. If
- *  there are non-color objetcs in the vector, the SU(3) generators T_a get
- *  sorted into the "unknown" group together with the non-color objects
- *  because we don't know whether these objects commute with the generators.
- *
- *  @param v Source vector of expressions
- *  @param delta8vec Vector of unity matrices (returned)
- *  @param fvec Vector of antisymmetric structure constants (returned)
- *  @param dvec Vector of symmetric structure constants (returned)
- *  @param Tvecs Vectors of generators, one for each representation label (returned)
- *  @param ONEvecs Vectors of unity elements, one for each representation label (returned)
- *  @param unknownvec Vector of all non-color objects (returned)
- *
- *  @see color::color_types
- *  @see recombine_color_string */
-void split_color_string_in_parts(const exvector & v, exvector & delta8vec,
-                                 exvector & fvec, exvector & dvec,
-                                 exvectorvector & Tvecs,
-                                 exvectorvector & ONEvecs,
-                                 exvector & unknownvec)
+/** Extract representation label from tinfo key (as returned by
+ *  return_type_tinfo()). */
+static unsigned char get_representation_label(const return_type_t& ti)
 {
-       // 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:
-                               GINAC_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:
-                               GINAC_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);
-               }
-       }
-}    
-
-/** Merge vectors of color objects sorted by object type into one vector,
- *  retaining the order within each group. This is the inverse operation of
- *  split_color_string_in_parts().
- *
- *  @param delta8vec Vector of unity matrices
- *  @param fvec Vector of antisymmetric structure constants
- *  @param dvec Vector of symmetric structure constants
- *  @param Tvecs Vectors of generators, one for each representation label
- *  @param ONEvecs Vectors of unity elements, one for each representation label
- *  @param unknownvec Vector of all non-color objects
- *  @return merged vector
- *
- *  @see color::color_types
- *  @see split_color_string_in_parts */
-exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
-                                exvector & dvec, exvectorvector & Tvecs,
-                                exvectorvector & 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 (unsigned char)ti.rl;
 }
 
-ex color_trace_of_one_representation_label(const exvector & v)
+ex color_trace(const ex & e, const std::set<unsigned char> & rls)
 {
-       if (v.size()==0) {
-               return numeric(COLOR_THREE);
-       } else if (v.size()==1) {
-               GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
-               return _ex0();
-       }
-       exvector v1=v;
-       ex last_element=v1.back();
-       GINAC_ASSERT(is_ex_exactly_of_type(last_element,color));
-       GINAC_ASSERT(ex_to_color(last_element).type==color::color_T);
-       v1.pop_back();
-       ex next_to_last_element=v1.back();
-       GINAC_ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
-       GINAC_ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
-       v1.pop_back();
-       exvector v2=v1;
-
-       const ex & last_index=ex_to_color(last_element).seq[0];
-       const ex & 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
-       
-       // FIXME: 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;
-       */
-}
+       if (is_a<color>(e)) {
 
-/** Calculate the trace over the (hidden) indices of the su(3) Lie algebra
- *  elements (the SU(3) generators and the unity element) of a specified
- *  representation label in a string of color objects.
- *
- *  @param v Vector of color objects
- *  @param rl Representation label
- *  @return value of the trace */
-ex color_trace(const exvector & v, unsigned rl)
-{
-       GINAC_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;
-       exvectorvector Tvecs;
-       Tvecs.resize(MAX_REPRESENTATION_LABELS);
-       exvectorvector 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"));
-       }
+               unsigned char rl = ex_to<color>(e).get_representation_label();
 
-       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"));
-                       }
-               }
-       }
+               // Are we taking the trace over this object's representation label?
+               if (rls.find(rl) == rls.end())
+                       return e;
 
-       return nonsimplified_ncmul(v_rest);
-}
+               // Yes, all generators are traceless, except for color_ONE
+               if (is_a<su3one>(e.op(0)))
+                       return _ex3;
+               else
+                       return _ex0;
 
-ex simplify_pure_color_string(const ex & e)
-{
-       GINAC_ASSERT(is_ex_exactly_of_type(e,ncmul));
-
-       exvector delta8vec;
-       exvector fvec;
-       exvector dvec;
-       exvectorvector Tvecs;
-       Tvecs.resize(MAX_REPRESENTATION_LABELS);
-       exvectorvector 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];
-                                       GINAC_ASSERT(is_ex_exactly_of_type(t1,color)&&
-                                                  (ex_to_color(t1).type==color::color_T)&&
-                                                  (ex_to_color(t1).seq.size()==1));
-                                       GINAC_ASSERT(is_ex_exactly_of_type(t2,color)&&
-                                                  (ex_to_color(t2).type==color::color_T)&&
-                                                  (ex_to_color(t2).seq.size()==1));
-                                       const coloridx & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
-                                       const coloridx & 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=_ex1();
-                                               t2=_ex1();
-                                               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(_ex1());
-                                               }
-                                               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);
-                                       }
-                               }
-                       }
-               }
-       }
-       
-       // FIXME: higher contractions
-       
-       return e;
-}
-       
-/** Perform some simplifications on an expression containing color objects. */
-ex simplify_color(const ex & e)
-{
-       // 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=_ex0();
-               for (unsigned i=0; i<e_expanded.nops(); ++i)
-                       sum += simplify_color(e_expanded.op(i));
-               
-               return sum;
-       }
+       } else if (is_exactly_a<mul>(e)) {
 
-       // simplification of commutative product=commutative product of simplifications
-       if (is_ex_exactly_of_type(e_expanded,mul)) {
-               ex prod=_ex1();
-               for (unsigned i=0; i<e_expanded.nops(); ++i)
-                       prod *= simplify_color(e_expanded.op(i));
-               
+               // 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;
-       }
 
-       // simplification of noncommutative product: test if everything is color
-       if (is_ex_exactly_of_type(e_expanded,ncmul)) {
-               bool all_color=true;
-               for (unsigned 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);
+       } 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;
                }
-       }
 
-       // cannot do anything
-       return e_expanded;
+       } 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 brute_force_sum_color_indices(const ex & e)
+ex color_trace(const ex & e, const lst & rll)
 {
-       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) {
-               GINAC_ASSERT(is_ex_of_type(*cit1,coloridx));
-               for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
-                       GINAC_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;
-                       }
-               }
+       // 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());
        }
 
-       std::vector<int> counter;
-       counter.resize(iv_double.size());
-       int l;
-       for (l=0; unsigned(l)<iv_double.size(); ++l) {
-               counter[l]=1;
-       }
+       return color_trace(e, rls);
+}
 
-       ex sum;
-       
-       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])>(int)COLOR_EIGHT)) {
-                       counter[l]=1;    
-                       l--;
-               }
-               if (l<2) { std::cout << counter[0] << counter[1] << std::endl; }
-               if (l<0) break;
-       }
-       
-       return sum;
+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);
 }
 
 } // namespace GiNaC