X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fcolor.cpp;h=9e355cb01a4118184d994bd1d865a7962a9dfc4d;hp=468645ddadfaa5e87912c16eb13b9e71469e0de7;hb=dad107ff48f68d45e72469a8716df375ae145cf3;hpb=afdd7fa8c6c0a587f7c80789198551383e8beb7b

diff --git a/ginac/color.cpp b/ginac/color.cpp
index 468645dd..9e355cb0 100644
--- a/ginac/color.cpp
+++ b/ginac/color.cpp
@@ -4,7 +4,7 @@
  *  No real implementation yet, to be done.     */
 
 /*
- *  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
@@ -33,11 +33,13 @@
 #include "ncmul.h"
 #include "numeric.h"
 #include "relational.h"
+#include "archive.h"
 #include "debugmsg.h"
+#include "utils.h"
 
-#ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
-#endif // ndef NO_GINAC_NAMESPACE
+
+GINAC_IMPLEMENT_REGISTERED_CLASS(color, indexed)
 
 //////////
 // default constructor, destructor, copy constructor assignment operator and helpers
@@ -45,48 +47,25 @@ namespace GiNaC {
 
 // public
 
-color::color() : type(invalid), representation_label(0)
-{
-    debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
-    tinfo_key=TINFO_color;
-}
-
-color::~color()
+color::color() : inherited(TINFO_color), type(invalid), representation_label(0)
 {
-    debugmsg("color destructor",LOGLEVEL_DESTRUCT);
-    destroy(0);
-}
-
-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;
+	debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
 }
 
 // protected
 
-void color::copy(color const & other)
+void color::copy(const color & other)
 {
-    indexed::copy(other);
-    type=other.type;
-    representation_label=other.representation_label;
+	inherited::copy(other);
+	type=other.type;
+	representation_label=other.representation_label;
 }
 
 void color::destroy(bool call_parent)
 {
-    if (call_parent) {
-        indexed::destroy(call_parent);
-    }
+	if (call_parent) {
+		inherited::destroy(call_parent);
+	}
 }
 
 //////////
@@ -95,879 +74,961 @@ void color::destroy(bool call_parent)
 
 // 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(color_types const t, unsigned rl) : type(t), 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());
+	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());
 }
 
-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(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());
+	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());
 }
 
-color::color(color_types const t, ex const & i1, ex const & i2, unsigned const rl)
-    : indexed(i1,i2), type(t), representation_label(rl)
+/** 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)
 {
-    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());
+	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());
 }
 
-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)
+/** 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)
 {
-    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());
+	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());
 }
 
-color::color(color_types const t, exvector const & iv, unsigned const rl)
-    : indexed(iv), type(t), representation_label(rl)
+/** 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)
 {
-    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());
+	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 const rl)
-    : indexed(ivp), type(t), representation_label(rl)
+color::color(color_types const t, exvector * ivp, unsigned rl)
+  : inherited(ivp), type(t), 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());
+	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());
 }
 
 //////////
-// functions overriding virtual functions from bases classes
+// archiving
 //////////
 
-// public
+/** 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"));
+}
 
-basic * color::duplicate() const
+/** Unarchive the object. */
+ex color::unarchive(const archive_node &n, const lst &sym_lst)
 {
-    debugmsg("color duplicate",LOGLEVEL_DUPLICATE);
-    return new color(*this);
+	return (new color(n, sym_lst))->setflag(status_flags::dynallocated);
 }
 
-void color::printraw(ostream & os) const
+/** Archive the object. */
+void color::archive(archive_node &n) const
 {
-    debugmsg("color printraw",LOGLEVEL_PRINT);
-    os << "color(type=" << (unsigned)type
-       << ",representation_label=" << representation_label
-       << ",indices=";
-    printrawindices(os);
-    os << ",hash=" << hashvalue << ",flags=" << flags << ")";
+	inherited::archive(n);
+	n.add_unsigned("type", type);
+	n.add_unsigned("representation", representation_label);
 }
 
-void color::printtree(ostream & os, unsigned indent) const
+//////////
+// functions overriding virtual functions from bases classes
+//////////
+
+// public
+
+void color::printraw(std::ostream & os) 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;
+	debugmsg("color printraw",LOGLEVEL_PRINT);
+	os << "color(type=" << (unsigned)type
+	   << ",representation_label=" << representation_label
+	   << ",indices=";
+	printrawindices(os);
+	os << ",hash=" << hashvalue << ",flags=" << flags << ")";
 }
 
-void color::print(ostream & os, unsigned upper_precedence) const
+void color::printtree(std::ostream & os, unsigned indent) 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);
+	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;
 }
 
-void color::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+void color::print(std::ostream & os, unsigned upper_precedence) const
 {
-    debugmsg("color print csrc",LOGLEVEL_PRINT);
-    print(os,upper_precedence);
+	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);
 }
 
 bool color::info(unsigned inf) const
 {
-    return indexed::info(inf);
+	return inherited::info(inf);
 }
 
 #define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
 
 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:
-        {
-            GINAC_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();
-                }
-            }
+	// 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;
 	}
-        break;
-    case color_d:
-        // check for d_{a,a,c} (=0) when a is symbolic
-        {
-            GINAC_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())) {
-                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 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:
-        {
-            GINAC_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())) {
-                GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
-                GINAC_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
-{
-    GINAC_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;
+	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);
+			}
+
+			// 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;
+	}
+	
+	return this->hold();
 }
+	
+// protected
 
-bool color::is_equal_same_type(basic const & other) const
+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) return false;
-    if (representation_label!=o->representation_label) return false;
-    return indexed::is_equal_same_type(other);
-}
+	GINAC_ASSERT(other.tinfo() == TINFO_color);
+	const color &o = static_cast<const color &>(other);
 
-#include <iostream>
+	if (type != o.type) {
+		// different type
+		return type < o.type ? -1 : 1;
+	}
+
+	if (representation_label != o.representation_label) {
+		// different representation label
+		return representation_label < o.representation_label ? -1 : 1;
+	}
 
-ex color::simplify_ncmul(exvector const & 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);
-            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()
-            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=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
-                    GINAC_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) {
-                GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
-                GINAC_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) {
-                GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
-                GINAC_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) {
-                GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
-                GINAC_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;
-                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);
-                    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) {
-                    GINAC_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));
-}
-
-ex color::thisexprseq(exvector const & v) const
-{
-    return color(type,v,representation_label);
+	return inherited::compare_same_type(other);
 }
 
-ex color::thisexprseq(exvector * vp) const
+bool color::is_equal_same_type(const basic & other) const
 {
-    return color(type,vp,representation_label);
+	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);
 }
 
-bool color::all_of_type_coloridx(void) const
+#include <iostream>
+
+ex color::simplify_ncmul(const exvector & v) 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;
-}
+	// 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;
+				}
+			}
+
+			// 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;
+	}
 
-//////////
-// virtual functions which can be overridden by derived classes
-//////////
+	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;
+	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));
+				}
+			}
+		}
+	}
 
-// none
+	// 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));
+				}
+			}
+		}
+	}
 
-//////////
-// non-virtual functions in this class
-//////////
+	// 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));
+					}
+				}
+			}
+		}
+	}
+	
+	// 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));
+		}
+	}
 
-//////////
-// static member variables
-//////////
+	// return a sorted vector
+	return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
+												   ONEvecs,unknownvec));
+}
 
-// none
+ex color::thisexprseq(const exvector & v) const
+{
+	return color(type,v,representation_label);
+}
 
-//////////
-// global constants
-//////////
+ex color::thisexprseq(exvector * vp) const
+{
+	return color(type,vp,representation_label);
+}
 
-const color some_color;
-type_info const & typeid_color=typeid(some_color);
+/** 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;
+	}
+	return true;
+}
 
 //////////
 // friend functions
 //////////
 
-color color_ONE(unsigned const rl)
+/** Construct an object representing the unity element of su(3).
+ *
+ *  @param rl Representation label
+ *  @return newly constructed object */
+color color_ONE(unsigned rl)
 {
-    return color(color::color_ONE,rl);
+	return color(color::color_ONE,rl);
 }
 
-color color_T(ex const & a, unsigned const 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)
 {
-    return color(color::color_T,a,rl);
+	return color(color::color_T,a,rl);
 }
 
-color color_f(ex const & a, ex const & b, ex const & c)
+/** 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)
 {
-    return color(color::color_f,a,b,c);
+	return color(color::color_f,a,b,c);
 }
 
-color color_d(ex const & a, ex const & b, ex const & 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)
 {
-    return color(color::color_d,a,b,c);
+	return color(color::color_d,a,b,c);
 }
 
-ex color_h(ex const & a, ex const & b, ex const & 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(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
 }
 
-color color_delta8(ex const & a, ex const & b)
+/** 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)
 {
-    return color(color::color_delta8,a,b);
+	return color(color::color_delta8,a,b);
 }
 
-void split_color_string_in_parts(exvector const & v, exvector & delta8vec,
+/** 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,
-                                 vector<exvector> & Tvecs,
-                                 vector<exvector> & ONEvecs,
+                                 exvectorvector & Tvecs,
+                                 exvectorvector & 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:
-                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);
-        }
-    }
+	// 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, 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;
-}
-
-ex color_trace_of_one_representation_label(exvector const & v)
-{
-    if (v.size()==0) {
-        return numeric(COLOR_THREE);
-    } else if (v.size()==1) {
-        GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
-        return exZERO();
-    }
-    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;
-
-    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
-    
-    // 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;
-    */
-}
-
-ex color_trace(exvector const & v, unsigned const 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;
-    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);
-}
-
-ex simplify_pure_color_string(ex const & e)
-{
-    GINAC_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];
-                    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));
-                    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);
-                    }
-                }
-            }
-        }
-    }
-    
-    // FIXME: higher contractions
-    
-    return e;
-}
-    
-ex simplify_color(ex const & 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=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;
-}
-
-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) {
-        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;
-            }
-        }
-    }
-
-    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])>(int)COLOR_EIGHT)) {
-            counter[l]=1;    
-            l--;
-        }
-        if (l<2) { cout << counter[0] << counter[1] << endl; }
-        if (l<0) break;
-    }
-    
-    return sum;
-}
-
-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);
-    }
-}
-
-#ifndef NO_GINAC_NAMESPACE
-} // namespace GiNaC
-#endif // ndef NO_GINAC_NAMESPACE
+                                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;
+}
 
+ex color_trace_of_one_representation_label(const exvector & v)
+{
+	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;
+	*/
+}
+
+/** 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"));
+	}
+
+	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);
+}
+
+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;
+	}
+
+	// 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));
+		
+		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);
+		}
+	}
+
+	// cannot do anything
+	return e_expanded;
+}
+
+ex brute_force_sum_color_indices(const ex & 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) {
+		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;
+			}
+		}
+	}
+
+	std::vector<int> counter;
+	counter.resize(iv_double.size());
+	int l;
+	for (l=0; unsigned(l)<iv_double.size(); ++l) {
+		counter[l]=1;
+	}
+
+	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;
+}
+
+} // namespace GiNaC