/** @file color.cpp
*
- * Implementation of GiNaC's color objects.
- * No real implementation yet, to be done.
- *
- * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
+ * Implementation of GiNaC's color (SU(3) Lie algebra) objects. */
+
+/*
+ * 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 "utils.h"
-//////////
-// default constructor, destructor, copy constructor assignment operator and helpers
-//////////
+#include <iostream>
+#include <stdexcept>
-// public
+namespace GiNaC {
-color::color() : type(invalid), representation_label(0)
-{
- debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
- tinfo_key=TINFO_color;
-}
+GINAC_IMPLEMENT_REGISTERED_CLASS(color, indexed)
-color::~color()
-{
- debugmsg("color destructor",LOGLEVEL_DESTRUCT);
- destroy(0);
-}
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3one, tensor,
+ print_func<print_dflt>(&su3one::do_print).
+ print_func<print_latex>(&su3one::do_print_latex))
-color::color(color const & other)
-{
- debugmsg("color copy constructor",LOGLEVEL_CONSTRUCT);
- copy (other);
-}
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3t, tensor,
+ print_func<print_dflt>(&su3t::do_print).
+ print_func<print_latex>(&su3t::do_print))
-color const & color::operator=(color const & other)
-{
- debugmsg("color operator=",LOGLEVEL_ASSIGNMENT);
- if (this != &other) {
- destroy(1);
- copy(other);
- }
- return *this;
-}
+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(color const & other)
-{
- indexed::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) {
- indexed::destroy(call_parent);
- }
}
+DEFAULT_CTOR(su3one)
+DEFAULT_CTOR(su3t)
+DEFAULT_CTOR(su3f)
+DEFAULT_CTOR(su3d)
+
//////////
// other constructors
//////////
-// protected
-
-color::color(color_types const t, unsigned const rl) : type(t), representation_label(rl)
+/** Construct object without any color index. This constructor is for
+ * internal use only. Use the color_ONE() function instead.
+ * @see color_ONE */
+color::color(const ex & b, unsigned char rl) : inherited(b), representation_label(rl)
{
- debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
}
-color::color(color_types const t, ex const & i1, unsigned const rl)
- : indexed(i1), type(t), representation_label(rl)
+/** Construct object with one color index. This constructor is for internal
+ * use only. Use the color_T() function instead.
+ * @see color_T */
+color::color(const ex & b, const ex & i1, unsigned char rl) : inherited(b, i1), representation_label(rl)
{
- debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
}
-color::color(color_types const t, ex const & i1, ex const & i2, unsigned const rl)
- : indexed(i1,i2), type(t), representation_label(rl)
+color::color(unsigned char rl, const exvector & v, bool discardable) : inherited(not_symmetric(), v, discardable), representation_label(rl)
{
- debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
}
-color::color(color_types const t, ex const & i1, ex const & i2, ex const & i3,
- unsigned const rl) : indexed(i1,i2,i3), type(t), representation_label(rl)
+color::color(unsigned char rl, std::auto_ptr<exvector> vp) : inherited(not_symmetric(), vp), representation_label(rl)
{
- debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
}
-color::color(color_types const t, exvector const & iv, unsigned const rl)
- : indexed(iv), type(t), representation_label(rl)
+return_type_t color::return_type_tinfo() const
{
- debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
-}
-
-color::color(color_types const t, exvector * ivp, unsigned const rl)
- : indexed(ivp), type(t), representation_label(rl)
-{
- debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
- ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
- tinfo_key=TINFO_color;
- ASSERT(all_of_type_coloridx());
+ return make_return_type_t<color>(representation_label);
}
//////////
-// functions overriding virtual functions from bases classes
+// archiving
//////////
-// public
-
-basic * color::duplicate() const
-{
- debugmsg("color duplicate",LOGLEVEL_DUPLICATE);
- return new color(*this);
-}
-
-void color::printraw(ostream & os) const
+void color::read_archive(const archive_node& n, lst& sym_lst)
{
- debugmsg("color printraw",LOGLEVEL_PRINT);
- os << "color(type=" << (unsigned)type
- << ",representation_label=" << representation_label
- << ",indices=";
- printrawindices(os);
- os << ",hash=" << hashvalue << ",flags=" << flags << ")";
+ inherited::read_archive(n, sym_lst);
+ unsigned rl;
+ n.find_unsigned("label", rl);
+ representation_label = rl;
}
-void color::printtree(ostream & os, unsigned indent) const
+void color::archive(archive_node &n) 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;
+ inherited::archive(n);
+ n.add_unsigned("label", representation_label);
}
-void color::print(ostream & os, unsigned upper_precedence) 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);
-}
+GINAC_BIND_UNARCHIVER(color);
+GINAC_BIND_UNARCHIVER(su3one);
+GINAC_BIND_UNARCHIVER(su3t);
+GINAC_BIND_UNARCHIVER(su3f);
+GINAC_BIND_UNARCHIVER(su3d);
-void color::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("color print csrc",LOGLEVEL_PRINT);
- print(os,upper_precedence);
-}
+//////////
+// functions overriding virtual functions from base classes
+//////////
-bool color::info(unsigned inf) const
+int color::compare_same_type(const basic & other) const
{
- return indexed::info(inf);
-}
+ GINAC_ASSERT(is_a<color>(other));
+ const color &o = static_cast<const color &>(other);
-#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:
- {
- 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();
- }
- }
+ if (representation_label != o.representation_label) {
+ // different representation label
+ return representation_label < o.representation_label ? -1 : 1;
}
- break;
- case color_d:
- // check for d_{a,a,c} (=0) when a is symbolic
- {
- ASSERT(seq.size()==3);
- coloridx const & idx1=ex_to_coloridx(seq[0]);
- coloridx const & idx2=ex_to_coloridx(seq[1]);
- coloridx const & idx3=ex_to_coloridx(seq[2]);
-
- if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
- return exZERO();
- } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
- return exZERO();
- }
-
- // check for three numeric indices
- if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
- ASSERT(idx1.get_value()<=idx2.get_value());
- ASSERT(idx2.get_value()<=idx3.get_value());
- if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)||
- CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) {
- return exHALF();
- } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
- return -exHALF();
- } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) {
- return 1/sqrt(numeric(3));
- } else if (CMPINDICES(8,8,8)) {
- return -1/sqrt(numeric(3));
- } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
- return -1/(2*sqrt(numeric(3)));
- }
- return exZERO();
- }
- }
- break;
- case color_f:
- {
- ASSERT(seq.size()==3);
- coloridx const & idx1=ex_to_coloridx(seq[0]);
- coloridx const & idx2=ex_to_coloridx(seq[1]);
- coloridx const & idx3=ex_to_coloridx(seq[2]);
-
- // check for three numeric indices
- if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
- ASSERT(idx1.get_value()<=idx2.get_value());
- ASSERT(idx2.get_value()<=idx3.get_value());
- if (CMPINDICES(1,2,3)) {
- return exONE();
- } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
- CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
- return exHALF();
- } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
- return -exHALF();
- } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) {
- return sqrt(numeric(3))/2;
- } else if (CMPINDICES(8,8,8)) {
- return -1/sqrt(numeric(3));
- } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
- return -1/(2*sqrt(numeric(3)));
- }
- return exZERO();
- }
- break;
- }
- default:
- // nothing to evaluate
- break;
- }
-
- return this->hold();
-}
-
-// protected
-int color::compare_same_type(basic const & other) const
-{
- ASSERT(other.tinfo() == TINFO_color);
- const color *o = static_cast<const color *>(&other);
- if (type==o->type) {
- if (representation_label==o->representation_label) {
- return indexed::compare_same_type(other);
- }
- return representation_label < o->representation_label ? -1 : 1;
- }
- return type < o->type ? -1 : 1;
+ return inherited::compare_same_type(other);
}
-bool color::is_equal_same_type(basic const & other) const
+bool color::match_same_type(const basic & other) const
{
- ASSERT(other.tinfo() == TINFO_color);
- const color *o = static_cast<const color *>(&other);
- if (type!=o->type) return false;
- if (representation_label!=o->representation_label) return false;
- return indexed::is_equal_same_type(other);
+ GINAC_ASSERT(is_a<color>(other));
+ const color &o = static_cast<const color &>(other);
+
+ return representation_label == o.representation_label;
}
-#include <iostream>
+DEFAULT_COMPARE(su3one)
+DEFAULT_COMPARE(su3t)
+DEFAULT_COMPARE(su3f)
+DEFAULT_COMPARE(su3d)
-ex color::simplify_ncmul(exvector const & v) const
+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
{
- // simplifications: contract delta8_{a,b} where possible
- // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...
- // remove superfluous ONEs
-
- // contract indices of delta8_{a,b} if they are different and symbolic
-
- exvector v_contracted=v;
- unsigned replacements;
- bool something_changed=false;
-
- exvector::iterator it=v_contracted.begin();
- while (it!=v_contracted.end()) {
- // process only delta8 objects
- if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) {
- color & d8=ex_to_nonconst_color(*it);
- ASSERT(d8.seq.size()==2);
- coloridx const & first_idx=ex_to_coloridx(d8.seq[0]);
- coloridx const & second_idx=ex_to_coloridx(d8.seq[1]);
- // delta8_{a,a} should have been contracted in color::eval()
- ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic()));
- ex saved_delta8=*it; // save to restore it later
-
- // try to contract first index
- replacements=1;
- if (first_idx.is_symbolic()) {
- replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx);
- if (replacements==1) {
- // not contracted except in itself, restore delta8 object
- *it=saved_delta8;
- } else {
- // a contracted index should occur exactly twice
- ASSERT(replacements==2);
- *it=exONE();
- something_changed=true;
- }
- }
-
- // try second index only if first was not contracted
- if ((replacements==1)&&(second_idx.is_symbolic())) {
- // first index not contracted, *it is guaranteed to be the original delta8 object
- replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx);
- if (replacements==1) {
- // not contracted except in itself, restore delta8 object
- *it=saved_delta8;
- } else {
- // a contracted index should occur exactly twice
- ASSERT(replacements==2);
- *it=exONE();
- something_changed=true;
- }
- }
- }
- ++it;
- }
-
- if (something_changed) {
- // do more simplifications later
- return nonsimplified_ncmul(v_contracted);
- }
-
- // there were no indices to contract
- // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown
- // (if there is at least one unknown object, all Ts will be unknown to not change the order)
-
- exvector delta8vec;
- exvector fvec;
- exvector dvec;
- vector<exvector> Tvecs;
- Tvecs.resize(MAX_REPRESENTATION_LABELS);
- vector<exvector> ONEvecs;
- ONEvecs.resize(MAX_REPRESENTATION_LABELS);
- exvector unknownvec;
-
- split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
-
- // d_{a,k,l} f_{b,k,l}=0 (includes case a=b)
- if ((dvec.size()>=1)&&(fvec.size()>=1)) {
- for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) {
- for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it1,color));
- ASSERT(is_ex_exactly_of_type(*it2,color));
- color const & col1=ex_to_color(*it1);
- color const & col2=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
- if (iv_intersect.size()>=2) return exZERO();
- }
- }
- }
-
- // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b)
- if (dvec.size()>=2) {
- for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) {
- for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it1,color));
- ASSERT(is_ex_exactly_of_type(*it2,color));
- color const & col1=ex_to_color(*it1);
- color const & col2=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
- if (iv_intersect.size()>=2) {
- if (iv_intersect.size()==3) {
- *it1=numeric(40)/numeric(3);
- *it2=exONE();
- } else {
- int sig1, sig2; // unimportant, since symmetric
- ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,false,&sig1);
- ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,false,&sig2);
- *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
- *it2=exONE();
- }
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
- }
- }
-
- // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b)
- if (fvec.size()>=2) {
- for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) {
- for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it1,color));
- ASSERT(is_ex_exactly_of_type(*it2,color));
- color const & col1=ex_to_color(*it1);
- color const & col2=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
- if (iv_intersect.size()>=2) {
- if (iv_intersect.size()==3) {
- *it1=numeric(24);
- *it2=exONE();
- } else {
- int sig1, sig2;
- ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,true,&sig1);
- ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,true,&sig2);
- *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
- *it2=exONE();
- }
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
- }
- }
-
- // d_{a,b,c} T_b T_c = 5/6 T_a
- // f_{a,b,c} T_b T_c = 3/2 I T_a
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- if ((Tvecs[rl].size()>=2)&&((dvec.size()>=1)||(fvec.size()>=1))) {
- for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) {
- exvector iv;
- ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
- ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T);
- iv.push_back(ex_to_color(*it1).seq[0]);
- iv.push_back(ex_to_color(*(it1+1)).seq[0]);
-
- // d_{a,b,c} T_b T_c = 5/6 T_a
- for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
- color const & dref=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(dref.seq,iv);
- if (iv_intersect.size()==2) {
- int sig; // unimportant, since symmetric
- ex free_idx=permute_free_index_to_front(dref.seq,iv,false,&sig);
- *it1=color(color_T,free_idx,rl);
- *(it1+1)=color(color_ONE,rl);
- *it2=numeric(5)/numeric(6);
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
-
- // f_{a,b,c} T_b T_c = 3/2 I T_a
- for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
- ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
- color const & fref=ex_to_color(*it2);
- exvector iv_intersect=idx_intersect(fref.seq,iv);
- if (iv_intersect.size()==2) {
- int sig;
- ex free_idx=permute_free_index_to_front(fref.seq,iv,true,&sig);
- *it1=color(color_T,free_idx,rl);
- *(it1+1)=color(color_ONE,rl);
- *it2=numeric(sig*3)/numeric(2)*I;
- return nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- }
- }
- }
- }
- }
-
- // clear all ONEs when there is at least one corresponding color_T
- // in this representation, retain one ONE otherwise
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- if (Tvecs[rl].size()!=0) {
- ONEvecs[rl].clear();
- } else if (ONEvecs[rl].size()!=0) {
- ONEvecs[rl].clear();
- ONEvecs[rl].push_back(color(color_ONE,rl));
- }
- }
-
- // return a sorted vector
- return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
- ONEvecs,unknownvec));
+ 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++;
+ }
+
+ if (s.empty())
+ return color(su3one(), representation_label);
+ else
+ return hold_ncmul(s);
}
-ex color::thisexprseq(exvector const & v) const
+ex color::thiscontainer(const exvector & v) const
{
- return color(type,v,representation_label);
+ return color(representation_label, v);
}
-ex color::thisexprseq(exvector * vp) const
+ex color::thiscontainer(std::auto_ptr<exvector> vp) const
{
- return color(type,vp,representation_label);
+ return color(representation_label, vp);
}
-bool color::all_of_type_coloridx(void) 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)
{
- // 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;
-}
+ GINAC_ASSERT(iv3.size() == 3);
+ GINAC_ASSERT(iv2.size() == 2);
-//////////
-// virtual functions which can be overridden by derived classes
-//////////
+ sig = 1;
-// none
+#define TEST_PERMUTATION(A,B,C,P) \
+ if (iv3[B].is_equal(iv2[0]) && iv3[C].is_equal(iv2[1])) { \
+ sig = P; \
+ return iv3[A]; \
+ }
+
+ TEST_PERMUTATION(0,1,2, 1);
+ TEST_PERMUTATION(0,2,1, -1);
+ TEST_PERMUTATION(1,0,2, -1);
+ TEST_PERMUTATION(1,2,0, 1);
+ TEST_PERMUTATION(2,0,1, 1);
+ TEST_PERMUTATION(2,1,0, -1);
+
+ throw(std::logic_error("permute_free_index_to_front(): no valid permutation found"));
+}
+
+/** Automatic symbolic evaluation of indexed symmetric structure constant. */
+ex su3d::eval_indexed(const basic & i) const
+{
+ GINAC_ASSERT(is_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;
+ }
-//////////
-// non-virtual functions in this class
-//////////
+ // 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;
+ }
-//////////
-// static member variables
-//////////
+ // No further simplifications
+ return i.hold();
+}
+
+
+/** 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();
+
+ if (is_exactly_a<su3t>(other->op(0))) {
+
+ // 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;
+
+ // T.a T.a = 4/3 ONE
+ if (other - self == 1) {
+ *self = numeric(4, 3);
+ *other = color_ONE(rl);
+ return true;
+
+ // 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;
+
+ // 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++;
+ }
+
+ 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;
+ }
+ }
-// none
+ 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;
+ }
+
+ } 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;
+ }
+ }
-//////////
-// global constants
-//////////
+ return false;
+}
+
+/** Contraction of an indexed antisymmetric structure constant with something else. */
+bool su3f::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+ GINAC_ASSERT(is_a<indexed>(*self));
+ GINAC_ASSERT(is_a<indexed>(*other));
+ GINAC_ASSERT(self->nops() == 4);
+ GINAC_ASSERT(is_a<su3f>(self->op(0)));
+
+ if (is_exactly_a<su3f>(other->op(0))) { // f*d is handled by su3d class
+
+ // 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;
+ }
+ }
-const color some_color;
-type_info const & typeid_color=typeid(some_color);
+ return false;
+}
//////////
-// friend functions
+// global functions
//////////
-color color_ONE(unsigned const 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);
}
-color color_T(ex const & a, unsigned const rl)
+ex color_T(const ex & a, unsigned char rl)
{
- return color(color::color_T,a,rl);
-}
+ static ex t = (new su3t)->setflag(status_flags::dynallocated);
-color color_f(ex const & a, ex const & b, ex const & c)
-{
- return color(color::color_f,a,b,c);
-}
+ 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"));
-color color_d(ex const & a, ex const & b, ex const & c)
-{
- return color(color::color_d,a,b,c);
+ return color(t, a, rl);
}
-ex color_h(ex const & a, ex const & b, ex const & c)
+ex color_f(const ex & a, const ex & b, const ex & c)
{
- return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
-}
+ static ex f = (new su3f)->setflag(status_flags::dynallocated);
-color color_delta8(ex const & a, ex const & b)
-{
- return color(color::color_delta8,a,b);
-}
+ 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"));
-void split_color_string_in_parts(exvector const & v, exvector & delta8vec,
- exvector & fvec, exvector & dvec,
- vector<exvector> & Tvecs,
- vector<exvector> & ONEvecs,
- exvector & unknownvec)
-{
- // if not all elements are of type color, put all Ts in unknownvec to
- // retain the ordering
- bool all_color=true;
- for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
- if (!is_ex_exactly_of_type(*cit,color)) {
- all_color=false;
- break;
- }
- }
-
- for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
- if (is_ex_exactly_of_type(*cit,color)) {
- switch (ex_to_color(*cit).type) {
- case color::color_delta8:
- delta8vec.push_back(*cit);
- break;
- case color::color_f:
- fvec.push_back(*cit);
- break;
- case color::color_d:
- dvec.push_back(*cit);
- break;
- case color::color_T:
- ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
- if (all_color) {
- Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
- } else {
- unknownvec.push_back(*cit);
- }
- break;
- case color::color_ONE:
- ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
- ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
- break;
- default:
- throw(std::logic_error("invalid type in split_color_string_in_parts()"));
- }
- } else {
- unknownvec.push_back(*cit);
- }
- }
-}
-
-exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
- exvector & dvec, vector<exvector> & Tvecs,
- vector<exvector> & ONEvecs, exvector & unknownvec)
-{
- unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- sz += Tvecs[rl].size();
- sz += ONEvecs[rl].size();
- }
- exvector v;
- v.reserve(sz);
-
- append_exvector_to_exvector(v,delta8vec);
- append_exvector_to_exvector(v,fvec);
- append_exvector_to_exvector(v,dvec);
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- append_exvector_to_exvector(v,Tvecs[rl]);
- append_exvector_to_exvector(v,ONEvecs[rl]);
- }
- append_exvector_to_exvector(v,unknownvec);
- return v;
+ return indexed(f, antisymmetric3(), a, b, c);
}
-ex color_trace_of_one_representation_label(exvector const & v)
+ex color_d(const ex & a, const ex & b, const ex & c)
{
- if (v.size()==0) {
- return numeric(COLOR_THREE);
- } else if (v.size()==1) {
- ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
- return exZERO();
- }
- exvector v1=v;
- ex last_element=v1.back();
- ASSERT(is_ex_exactly_of_type(last_element,color));
- ASSERT(ex_to_color(last_element).type==color::color_T);
- v1.pop_back();
- ex next_to_last_element=v1.back();
- ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
- ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
- v1.pop_back();
- exvector v2=v1;
-
- ex const & last_index=ex_to_color(last_element).seq[0];
- ex const & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
- ex summation_index=coloridx();
-
- v2.push_back(color_T(summation_index)); // don't care about the representation_label
-
- // check this formula for SU(N) with N!=3 !!!!!!!!!
- return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
- % color_trace_of_one_representation_label(v1)
- +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
- % color_trace_of_one_representation_label(v2);
- /*
- ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
- % color_trace_of_one_representation_label(v1);
- cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
- ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
- % color_trace_of_one_representation_label(v2);
- cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
- return term1+term2;
- */
-}
+ static ex d = (new su3d)->setflag(status_flags::dynallocated);
-ex color_trace(exvector const & v, unsigned const rl)
-{
- ASSERT(rl<MAX_REPRESENTATION_LABELS);
-
- exvector v_rest;
- v_rest.reserve(v.size()+1); // max size if trace is empty
-
- exvector delta8vec;
- exvector fvec;
- exvector dvec;
- vector<exvector> Tvecs;
- Tvecs.resize(MAX_REPRESENTATION_LABELS);
- vector<exvector> ONEvecs;
- ONEvecs.resize(MAX_REPRESENTATION_LABELS);
- exvector unknownvec;
-
- split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
-
- if (unknownvec.size()!=0) {
- throw(std::invalid_argument("color_trace(): expression must be expanded"));
- }
-
- append_exvector_to_exvector(v_rest,delta8vec);
- append_exvector_to_exvector(v_rest,fvec);
- append_exvector_to_exvector(v_rest,dvec);
- for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
- if (i!=rl) {
- append_exvector_to_exvector(v_rest,Tvecs[i]);
- append_exvector_to_exvector(v_rest,ONEvecs[i]);
- } else {
- if (Tvecs[i].size()!=0) {
- v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
- } else if (ONEvecs[i].size()!=0) {
- v_rest.push_back(numeric(COLOR_THREE));
- } else {
- throw(std::logic_error("color_trace(): representation_label not in color string"));
- }
- }
- }
-
- return nonsimplified_ncmul(v_rest);
+ 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);
}
-ex simplify_pure_color_string(ex const & e)
+ex color_h(const ex & a, const ex & b, const ex & c)
{
- ASSERT(is_ex_exactly_of_type(e,ncmul));
-
- exvector delta8vec;
- exvector fvec;
- exvector dvec;
- vector<exvector> Tvecs;
- Tvecs.resize(MAX_REPRESENTATION_LABELS);
- vector<exvector> ONEvecs;
- ONEvecs.resize(MAX_REPRESENTATION_LABELS);
- exvector unknownvec;
-
- split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
-
- // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
- for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
- if (Tvecs[rl].size()>=2) {
- for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
- for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
- ex & t1=Tvecs[rl][i];
- ex & t2=Tvecs[rl][j];
- ASSERT(is_ex_exactly_of_type(t1,color)&&
- (ex_to_color(t1).type==color::color_T)&&
- (ex_to_color(t1).seq.size()==1));
- ASSERT(is_ex_exactly_of_type(t2,color)&&
- (ex_to_color(t2).type==color::color_T)&&
- (ex_to_color(t2).seq.size()==1));
- coloridx const & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
- coloridx const & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
-
- if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
- exvector S;
- for (unsigned k=i+1; k<j; ++k) {
- S.push_back(Tvecs[rl][k]);
- }
- t1=exONE();
- t2=exONE();
- ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- for (unsigned k=i+1; k<j; ++k) {
- S.push_back(exONE());
- }
- t1=color_trace_of_one_representation_label(S);
- ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(
- delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
- return simplify_color(term1+term2);
- }
- }
- }
- }
- }
-
- // TODO: higher contractions!!!!!!!!!!!!!
-
- return e;
+ return color_d(a, b, c) + I * color_f(a, b, c);
}
-
-ex simplify_color(ex const & e)
+
+/** 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)
{
- // 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;
+ return *(ti.tinfo) == typeid(color);
}
-ex brute_force_sum_color_indices(ex const & e)
+/** Extract representation label from tinfo key (as returned by
+ * return_type_tinfo()). */
+static unsigned char get_representation_label(const return_type_t& ti)
{
- exvector iv_all=e.get_indices();
- exvector iv_double;
-
- // find double symbolic indices
- if (iv_all.size()<2) return e;
- for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
- ASSERT(is_ex_of_type(*cit1,coloridx));
- for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
- ASSERT(is_ex_of_type(*cit2,coloridx));
- if (ex_to_coloridx(*cit1).is_symbolic() &&
- ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
- iv_double.push_back(*cit1);
- break;
- }
- }
- }
-
- vector<int> counter;
- counter.resize(iv_double.size());
- int l;
- for (l=0; unsigned(l)<iv_double.size(); ++l) {
- counter[l]=1;
- }
-
- ex sum=exZERO();
-
- while (1) {
- ex term=e;
- for (l=0; unsigned(l)<iv_double.size(); ++l) {
- term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
- //iv_double[l].print(cout);
- //cout << " " << counter[l] << " ";
- }
- //cout << endl;
- sum += term;
-
- // increment counter[]
- l=iv_double.size()-1;
- while ((l>=0)&&((++counter[l])>COLOR_EIGHT)) {
- counter[l]=1;
- l--;
- }
- if (l<2) { cout << counter[0] << counter[1] << endl; }
- if (l<0) break;
- }
-
- return sum;
+ return (unsigned char)ti.rl;
}
-void append_exvector_to_exvector(exvector & dest, exvector const & source)
+ex color_trace(const ex & e, const std::set<unsigned char> & rls)
{
- for (exvector::const_iterator cit=source.begin(); cit!=source.end(); ++cit) {
- dest.push_back(*cit);
- }
-}
+ if (is_a<color>(e)) {
+
+ unsigned char rl = ex_to<color>(e).get_representation_label();
+
+ // Are we taking the trace over this object's representation label?
+ if (rls.find(rl) == rls.end())
+ return e;
+
+ // Yes, all generators are traceless, except for color_ONE
+ if (is_a<su3one>(e.op(0)))
+ return _ex3;
+ else
+ return _ex0;
+
+ } else if (is_exactly_a<mul>(e)) {
+
+ // Trace of product: pull out non-color factors
+ ex prod = _ex1;
+ for (size_t i=0; i<e.nops(); i++) {
+ const ex &o = e.op(i);
+ if (is_color_tinfo(o.return_type_tinfo()))
+ prod *= color_trace(o, rls);
+ else
+ prod *= o;
+ }
+ return prod;
+
+ } else if (is_exactly_a<ncmul>(e)) {
+
+ unsigned char rl = get_representation_label(e.return_type_tinfo());
+
+ // Are we taking the trace over this string's representation label?
+ if (rls.find(rl) == rls.end())
+ return e;
+ // Yes, expand product if necessary
+ ex e_expanded = e.expand();
+ if (!is_a<ncmul>(e_expanded))
+ return color_trace(e_expanded, rls);
+ size_t num = e.nops();
+ if (num == 2) {
+ // Tr T_a T_b = 1/2 delta_a_b
+ return delta_tensor(e.op(0).op(1), e.op(1).op(1)) / 2;
+ } else if (num == 3) {
+ // Tr T_a T_b T_c = 1/4 h_a_b_c
+ return color_h(e.op(0).op(1), e.op(1).op(1), e.op(2).op(1)) / 4;
+ } else {
+ // Traces of 4 or more generators are computed recursively:
+ // Tr T_a1 .. T_an =
+ // 1/6 delta_a(n-1)_an Tr T_a1 .. T_a(n-2)
+ // + 1/2 h_a(n-1)_an_k Tr T_a1 .. T_a(n-2) T_k
+ const ex &last_index = e.op(num - 1).op(1);
+ const ex &next_to_last_index = e.op(num - 2).op(1);
+ idx summation_index((new symbol)->setflag(status_flags::dynallocated), 8);
+ exvector v1;
+ v1.reserve(num - 2);
+ for (size_t i=0; i<num-2; i++)
+ v1.push_back(e.op(i));
+ exvector v2 = v1;
+ v2.push_back(color_T(summation_index, rl));
+ return delta_tensor(next_to_last_index, last_index) * color_trace(ncmul(v1), rl) / 6
+ + color_h(next_to_last_index, last_index, summation_index) * color_trace(ncmul(v2), rl) / 2;
+ }
+ } else if (e.nops() > 0) {
+
+ // Trace maps to all other container classes (this includes sums)
+ pointer_to_map_function_1arg<const std::set<unsigned char> &> fcn(color_trace, rls);
+ return e.map(fcn);
+
+ } else
+ return _ex0;
+}
+
+ex color_trace(const ex & e, const lst & rll)
+{
+ // Convert list to set
+ std::set<unsigned char> rls;
+ for (lst::const_iterator i = rll.begin(); i != rll.end(); ++i) {
+ if (i->info(info_flags::nonnegint))
+ rls.insert(ex_to<numeric>(*i).to_int());
+ }
+
+ return color_trace(e, rls);
+}
+
+ex color_trace(const ex & e, unsigned char rl)
+{
+ // Convert label to set
+ std::set<unsigned char> rls;
+ rls.insert(rl);
+
+ return color_trace(e, rls);
+}
+} // namespace GiNaC
git://www.ginac.de / ginac.git / blobdiff