X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fcolor.cpp;h=34ab4d82d7ef055a80b8b6734621ceed147c000c;hp=9e355cb01a4118184d994bd1d865a7962a9dfc4d;hb=073bf40a73e419a3dbcb6dfa190947ce2cc3bdce;hpb=4e3a4ac2bcb0837611ea31bc8fc05d84a20c33ac diff --git a/ginac/color.cpp b/ginac/color.cpp index 9e355cb0..34ab4d82 100644 --- a/ginac/color.cpp +++ b/ginac/color.cpp @@ -1,10 +1,9 @@ /** @file color.cpp * - * Implementation of GiNaC's color objects. - * No real implementation yet, to be done. */ + * Implementation of GiNaC's color (SU(3) Lie algebra) objects. */ /* - * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -18,1017 +17,629 @@ * * 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 -#include -#include -#include -#include - #include "color.h" -#include "ex.h" -#include "coloridx.h" +#include "idx.h" #include "ncmul.h" +#include "symmetry.h" +#include "operators.h" #include "numeric.h" -#include "relational.h" +#include "mul.h" +#include "power.h" // for sqrt() +#include "symbol.h" #include "archive.h" -#include "debugmsg.h" #include "utils.h" +#include +#include + namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS(color, indexed) -////////// -// default constructor, destructor, copy constructor assignment operator and helpers -////////// +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3one, tensor, + print_func(&su3one::do_print). + print_func(&su3one::do_print_latex)) -// public +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3t, tensor, + print_func(&su3t::do_print). + print_func(&su3t::do_print)) -color::color() : inherited(TINFO_color), type(invalid), representation_label(0) -{ - debugmsg("color default constructor",LOGLEVEL_CONSTRUCT); -} +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3f, tensor, + print_func(&su3f::do_print). + print_func(&su3f::do_print)) -// protected +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(su3d, tensor, + print_func(&su3d::do_print). + print_func(&su3d::do_print)) -void color::copy(const color & other) -{ - inherited::copy(other); - type=other.type; - representation_label=other.representation_label; -} +////////// +// default constructors +////////// -void color::destroy(bool call_parent) +color::color() : representation_label(0) { - if (call_parent) { - inherited::destroy(call_parent); - } } +DEFAULT_CTOR(su3one) +DEFAULT_CTOR(su3t) +DEFAULT_CTOR(su3f) +DEFAULT_CTOR(su3d) + ////////// // other constructors ////////// -// protected - -/** Construct object without any color index. This constructor is for internal - * use only. Use the color_ONE() function instead. +/** Construct object without any color index. This constructor is for + * internal use only. Use the color_ONE() function instead. * @see color_ONE */ -color::color(color_types const t, unsigned rl) : type(t), representation_label(rl) +color::color(const ex & b, unsigned char rl) : inherited(b), representation_label(rl) { - debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT); - GINAC_ASSERT(representation_label vp) : inherited(not_symmetric(), vp), representation_label(rl) { - debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT); - GINAC_ASSERT(representation_label(representation_label); } ////////// // archiving ////////// -/** Construct object from archive_node. */ -color::color(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) -{ - debugmsg("color constructor from archive_node", LOGLEVEL_CONSTRUCT); - unsigned int ty; - if (!(n.find_unsigned("type", ty))) - throw (std::runtime_error("unknown color type in archive")); - type = (color_types)ty; - if (!(n.find_unsigned("representation", representation_label))) - throw (std::runtime_error("unknown color representation label in archive")); -} - -/** Unarchive the object. */ -ex color::unarchive(const archive_node &n, const lst &sym_lst) +void color::read_archive(const archive_node& n, lst& sym_lst) { - return (new color(n, sym_lst))->setflag(status_flags::dynallocated); + inherited::read_archive(n, sym_lst); + unsigned rl; + n.find_unsigned("label", rl); + representation_label = rl; } -/** Archive the object. */ void color::archive(archive_node &n) const { inherited::archive(n); - n.add_unsigned("type", type); - n.add_unsigned("representation", representation_label); + n.add_unsigned("label", representation_label); } +GINAC_BIND_UNARCHIVER(color); +GINAC_BIND_UNARCHIVER(su3one); +GINAC_BIND_UNARCHIVER(su3t); +GINAC_BIND_UNARCHIVER(su3f); +GINAC_BIND_UNARCHIVER(su3d); + ////////// -// functions overriding virtual functions from bases classes +// functions overriding virtual functions from base classes ////////// -// public - -void color::printraw(std::ostream & os) const +int color::compare_same_type(const basic & other) const { - debugmsg("color printraw",LOGLEVEL_PRINT); - os << "color(type=" << (unsigned)type - << ",representation_label=" << representation_label - << ",indices="; - printrawindices(os); - os << ",hash=" << hashvalue << ",flags=" << flags << ")"; + GINAC_ASSERT(is_a(other)); + const color &o = static_cast(other); + + if (representation_label != o.representation_label) { + // different representation label + return representation_label < o.representation_label ? -1 : 1; + } + + return inherited::compare_same_type(other); } -void color::printtree(std::ostream & os, unsigned indent) const +bool color::match_same_type(const basic & other) const { - debugmsg("color printtree",LOGLEVEL_PRINT); - os << std::string(indent,' ') << "color object: " - << "type=" << (unsigned)type - << ",representation_label=" << representation_label << ", "; - os << seq.size() << " indices" << std::endl; - printtreeindices(os,indent); - os << std::string(indent,' ') << "hash=" << hashvalue - << " (0x" << std::hex << hashvalue << std::dec << ")" - << ", flags=" << flags << std::endl; + GINAC_ASSERT(is_a(other)); + const color &o = static_cast(other); + + return representation_label == o.representation_label; } -void color::print(std::ostream & os, unsigned upper_precedence) const +DEFAULT_COMPARE(su3one) +DEFAULT_COMPARE(su3t) +DEFAULT_COMPARE(su3f) +DEFAULT_COMPARE(su3d) + +DEFAULT_PRINT_LATEX(su3one, "ONE", "\\mathbb{1}") +DEFAULT_PRINT(su3t, "T") +DEFAULT_PRINT(su3f, "f") +DEFAULT_PRINT(su3d, "d") + +/** Perform automatic simplification on noncommutative product of color + * objects. This removes superfluous ONEs. */ +ex color::eval_ncmul(const exvector & v) const { - debugmsg("color print",LOGLEVEL_PRINT); - switch (type) { - case color_T: - os << "T"; - if (representation_label!=0) { - os << "^(" << representation_label << ")"; - } - break; - case color_f: - os << "f"; - break; - case color_d: - os << "d"; - break; - case color_delta8: - os << "delta8"; - break; - case color_ONE: - os << "color_ONE"; - break; - case invalid: - default: - os << "INVALID_COLOR_OBJECT"; - break; + exvector s; + s.reserve(v.size()); + + // Remove superfluous ONEs + exvector::const_iterator it = v.begin(), itend = v.end(); + while (it != itend) { + if (!is_a(it->op(0))) + s.push_back(*it); + it++; } - printindices(os); + + if (s.empty()) + return color(su3one(), representation_label); + else + return hold_ncmul(s); } -bool color::info(unsigned inf) const +ex color::thiscontainer(const exvector & v) const { - return inherited::info(inf); + return color(representation_label, v); } -#define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C))) +ex color::thiscontainer(std::auto_ptr vp) const +{ + return color(representation_label, vp); +} -ex color::eval(int level) const +/** Given a vector iv3 of three indices and a vector iv2 of two indices that + * is a subset of iv3, return the (free) index that is in iv3 but not in + * iv2 and the sign introduced by permuting that index to the front. + * + * @param iv3 Vector of 3 indices + * @param iv2 Vector of 2 indices, must be a subset of iv3 + * @param sig Returs sign introduced by index permutation + * @return the free index (the one that is in iv3 but not in iv2) */ +static ex permute_free_index_to_front(const exvector & iv3, const exvector & iv2, int & sig) { - // canonicalize indices - - bool antisymmetric=false; - - switch (type) { - case color_f: - antisymmetric=true; // no break here! - case color_d: - case color_delta8: - { - exvector iv=seq; - int sig=canonicalize_indices(iv,antisymmetric); - if (sig!=INT_MAX) { - // something has changed while sorting indices, more evaluations later - if (sig==0) return _ex0(); - return ex(sig)*color(type,iv,representation_label); - } - } - break; - default: - // nothing to canonicalize - break; - } + GINAC_ASSERT(iv3.size() == 3); + GINAC_ASSERT(iv2.size() == 2); - switch (type) { - case color_delta8: - { - GINAC_ASSERT(seq.size()==2); - const coloridx & idx1=ex_to_coloridx(seq[0]); - const coloridx & idx2=ex_to_coloridx(seq[1]); - - // check for delta8_{a,a} where a is a symbolic index, replace by 8 - if ((idx1.is_symbolic())&&(idx1.is_equal_same_type(idx2))) { - return ex(COLOR_EIGHT); - } + sig = 1; - // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1 - if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) { - if ((idx1.get_value()!=idx2.get_value())) { - return _ex1(); - } else { - return _ex0(); - } - } - } - break; - case color_d: - // check for d_{a,a,c} (=0) when a is symbolic - { - GINAC_ASSERT(seq.size()==3); - const coloridx & idx1=ex_to_coloridx(seq[0]); - const coloridx & idx2=ex_to_coloridx(seq[1]); - const coloridx & idx3=ex_to_coloridx(seq[2]); - - if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) { - return _ex0(); - } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) { - return _ex0(); - } - - // check for three numeric indices - if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) { - GINAC_ASSERT(idx1.get_value()<=idx2.get_value()); - GINAC_ASSERT(idx2.get_value()<=idx3.get_value()); - if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)|| - CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) { - return _ex1_2(); - } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) { - return -_ex1_2(); - } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) { - return 1/sqrt(numeric(3)); - } else if (CMPINDICES(8,8,8)) { - return -1/sqrt(numeric(3)); - } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) { - return -1/(2*sqrt(numeric(3))); - } - return _ex0(); - } - } - break; - case color_f: - { - GINAC_ASSERT(seq.size()==3); - const coloridx & idx1=ex_to_coloridx(seq[0]); - const coloridx & idx2=ex_to_coloridx(seq[1]); - const coloridx & idx3=ex_to_coloridx(seq[2]); - - // check for three numeric indices - if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) { - GINAC_ASSERT(idx1.get_value()<=idx2.get_value()); - GINAC_ASSERT(idx2.get_value()<=idx3.get_value()); - if (CMPINDICES(1,2,3)) { - return _ex1(); - } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)|| - CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) { - return _ex1_2(); - } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) { - return -_ex1_2(); - } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) { - return sqrt(numeric(3))/2; - } else if (CMPINDICES(8,8,8)) { - return -1/sqrt(numeric(3)); - } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) { - return -1/(2*sqrt(numeric(3))); - } - return _ex0(); - } - break; - } - default: - // nothing to evaluate - break; +#define TEST_PERMUTATION(A,B,C,P) \ + if (iv3[B].is_equal(iv2[0]) && iv3[C].is_equal(iv2[1])) { \ + sig = P; \ + return iv3[A]; \ } - return this->hold(); -} - -// protected - -int color::compare_same_type(const basic & other) const -{ - GINAC_ASSERT(other.tinfo() == TINFO_color); - const color &o = static_cast(other); - - if (type != o.type) { - // different type - return type < o.type ? -1 : 1; + TEST_PERMUTATION(0,1,2, 1); + TEST_PERMUTATION(0,2,1, -1); + TEST_PERMUTATION(1,0,2, -1); + TEST_PERMUTATION(1,2,0, 1); + TEST_PERMUTATION(2,0,1, 1); + TEST_PERMUTATION(2,1,0, -1); + + throw(std::logic_error("permute_free_index_to_front(): no valid permutation found")); +} + +/** Automatic symbolic evaluation of indexed symmetric structure constant. */ +ex su3d::eval_indexed(const basic & i) const +{ + GINAC_ASSERT(is_a(i)); + GINAC_ASSERT(i.nops() == 4); + GINAC_ASSERT(is_a(i.op(0))); + + // Convolutions are zero + if (!(static_cast(i).get_dummy_indices().empty())) + return _ex0; + + // Numeric evaluation + if (static_cast(i).all_index_values_are(info_flags::nonnegint)) { + + // Sort indices + int v[3]; + for (unsigned j=0; j<3; j++) + v[j] = ex_to(ex_to(i.op(j + 1)).get_value()).to_int(); + if (v[0] > v[1]) std::swap(v[0], v[1]); + if (v[0] > v[2]) std::swap(v[0], v[2]); + if (v[1] > v[2]) std::swap(v[1], v[2]); + +#define CMPINDICES(A,B,C) ((v[0] == (A)) && (v[1] == (B)) && (v[2] == (C))) + + // Check for non-zero elements + if (CMPINDICES(1,4,6) || CMPINDICES(1,5,7) || CMPINDICES(2,5,6) + || CMPINDICES(3,4,4) || CMPINDICES(3,5,5)) + return _ex1_2; + else if (CMPINDICES(2,4,7) || CMPINDICES(3,6,6) || CMPINDICES(3,7,7)) + return _ex_1_2; + else if (CMPINDICES(1,1,8) || CMPINDICES(2,2,8) || CMPINDICES(3,3,8)) + return sqrt(_ex3)*_ex1_3; + else if (CMPINDICES(8,8,8)) + return sqrt(_ex3)*_ex_1_3; + else if (CMPINDICES(4,4,8) || CMPINDICES(5,5,8) + || CMPINDICES(6,6,8) || CMPINDICES(7,7,8)) + return sqrt(_ex3)/_ex_6; + else + return _ex0; } - if (representation_label != o.representation_label) { - // different representation label - return representation_label < o.representation_label ? -1 : 1; + // No further simplifications + return i.hold(); +} + +/** Automatic symbolic evaluation of indexed antisymmetric structure constant. */ +ex su3f::eval_indexed(const basic & i) const +{ + GINAC_ASSERT(is_a(i)); + GINAC_ASSERT(i.nops() == 4); + GINAC_ASSERT(is_a(i.op(0))); + + // Numeric evaluation + if (static_cast(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(ex_to(i.op(j + 1)).get_value()).to_int(); + int sign = 1; + if (v[0] > v[1]) { std::swap(v[0], v[1]); sign = -sign; } + if (v[0] > v[2]) { std::swap(v[0], v[2]); sign = -sign; } + if (v[1] > v[2]) { std::swap(v[1], v[2]); sign = -sign; } + + // Check for non-zero elements + if (CMPINDICES(1,2,3)) + return sign; + else if (CMPINDICES(1,4,7) || CMPINDICES(2,4,6) + || CMPINDICES(2,5,7) || CMPINDICES(3,4,5)) + return _ex1_2 * sign; + else if (CMPINDICES(1,5,6) || CMPINDICES(3,6,7)) + return _ex_1_2 * sign; + else if (CMPINDICES(4,5,8) || CMPINDICES(6,7,8)) + return sqrt(_ex3)/2 * sign; + else + return _ex0; } - return inherited::compare_same_type(other); + // No further simplifications + return i.hold(); } -bool color::is_equal_same_type(const basic & other) const -{ - GINAC_ASSERT(other.tinfo() == TINFO_color); - const color &o = static_cast(other); - if (type != o.type) return false; - if (representation_label != o.representation_label) return false; - return inherited::is_equal_same_type(other); -} +/** Contraction of generator with something else. */ +bool su3t::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const +{ + GINAC_ASSERT(is_a(*self)); + GINAC_ASSERT(is_a(*other)); + GINAC_ASSERT(self->nops() == 2); + GINAC_ASSERT(is_a(self->op(0))); + unsigned char rl = ex_to(*self).get_representation_label(); -#include + if (is_exactly_a(other->op(0))) { -ex color::simplify_ncmul(const exvector & v) const -{ - // simplifications: contract delta8_{a,b} where possible - // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),... - // remove superfluous ONEs - - // contract indices of delta8_{a,b} if they are different and symbolic - - exvector v_contracted=v; - unsigned replacements; - bool something_changed=false; - - exvector::iterator it=v_contracted.begin(); - while (it!=v_contracted.end()) { - // process only delta8 objects - if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) { - color & d8=ex_to_nonconst_color(*it); - GINAC_ASSERT(d8.seq.size()==2); - const coloridx & first_idx=ex_to_coloridx(d8.seq[0]); - const coloridx & second_idx=ex_to_coloridx(d8.seq[1]); - // delta8_{a,a} should have been contracted in color::eval() - GINAC_ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic())); - ex saved_delta8=*it; // save to restore it later - - // try to contract first index - replacements=1; - if (first_idx.is_symbolic()) { - replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx); - if (replacements==1) { - // not contracted except in itself, restore delta8 object - *it=saved_delta8; - } else { - // a contracted index should occur exactly twice - GINAC_ASSERT(replacements==2); - *it=_ex1(); - something_changed=true; - } - } + // Contraction only makes sense if the represenation labels are equal + GINAC_ASSERT(is_a(*other)); + if (ex_to(*other).get_representation_label() != rl) + return false; - // try second index only if first was not contracted - if ((replacements==1)&&(second_idx.is_symbolic())) { - // first index not contracted, *it is guaranteed to be the original delta8 object - replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx); - if (replacements==1) { - // not contracted except in itself, restore delta8 object - *it=saved_delta8; - } else { - // a contracted index should occur exactly twice - GINAC_ASSERT(replacements==2); - *it=_ex1(); - something_changed=true; - } - } - } - ++it; - } + // T.a T.a = 4/3 ONE + if (other - self == 1) { + *self = numeric(4, 3); + *other = color_ONE(rl); + return true; - if (something_changed) { - // do more simplifications later - return nonsimplified_ncmul(v_contracted); - } + // T.a T.b T.a = -1/6 T.b + } else if (other - self == 2 + && is_a(self[1])) { + *self = numeric(-1, 6); + *other = _ex1; + return true; - // there were no indices to contract - // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown - // (if there is at least one unknown object, all Ts will be unknown to not change the order) - - exvector delta8vec; - exvector fvec; - exvector dvec; - exvectorvector Tvecs; - Tvecs.resize(MAX_REPRESENTATION_LABELS); - exvectorvector ONEvecs; - ONEvecs.resize(MAX_REPRESENTATION_LABELS); - exvector unknownvec; - - split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec); - - // d_{a,k,l} f_{b,k,l}=0 (includes case a=b) - if ((dvec.size()>=1)&&(fvec.size()>=1)) { - for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) { - for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) { - GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)); - GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)); - const color & col1=ex_to_color(*it1); - const color & col2=ex_to_color(*it2); - exvector iv_intersect=idx_intersect(col1.seq,col2.seq); - if (iv_intersect.size()>=2) return _ex0(); - } - } - } - - // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b) - if (dvec.size()>=2) { - for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) { - for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) { - GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)); - GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)); - const color & col1=ex_to_color(*it1); - const color & col2=ex_to_color(*it2); - exvector iv_intersect=idx_intersect(col1.seq,col2.seq); - if (iv_intersect.size()>=2) { - if (iv_intersect.size()==3) { - *it1=numeric(40)/numeric(3); - *it2=_ex1(); - } else { - int dummy; // sign unimportant, since symmetric - ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&dummy); - ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&dummy); - *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2); - *it2=_ex1(); - } - return nonsimplified_ncmul(recombine_color_string( - delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec)); + // T.a S T.a = 1/2 Tr(S) - 1/6 S + } else { + exvector::iterator it = self + 1; + while (it != other) { + if (!is_a(*it)) { + return false; } + it++; } - } - } - // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b) - if (fvec.size()>=2) { - for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) { - for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) { - GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)); - GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)); - const color & col1=ex_to_color(*it1); - const color & col2=ex_to_color(*it2); - exvector iv_intersect=idx_intersect(col1.seq,col2.seq); - if (iv_intersect.size()>=2) { - if (iv_intersect.size()==3) { - *it1=numeric(24); - *it2=_ex1(); - } else { - int sig1, sig2; - ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&sig1); - ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&sig2); - *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2); - *it2=_ex1(); - } - return nonsimplified_ncmul(recombine_color_string( - delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec)); - } + it = self + 1; + ex S = _ex1; + while (it != other) { + S *= *it; + *it++ = _ex1; } + + *self = color_trace(S, rl) * color_ONE(rl) / 2 - S / 6; + *other = _ex1; + return true; } } - // d_{a,b,c} T_b T_c = 5/6 T_a - // f_{a,b,c} T_b T_c = 3/2 I T_a - for (unsigned rl=0; rl=2)&&((dvec.size()>=1)||(fvec.size()>=1))) { - for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) { - exvector iv; - GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T); - GINAC_ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T); - iv.push_back(ex_to_color(*it1).seq[0]); - iv.push_back(ex_to_color(*(it1+1)).seq[0]); - - // d_{a,b,c} T_b T_c = 5/6 T_a - for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) { - GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d); - const color & dref=ex_to_color(*it2); - exvector iv_intersect=idx_intersect(dref.seq,iv); - if (iv_intersect.size()==2) { - int dummy; // sign unimportant, since symmetric - ex free_idx=permute_free_index_to_front(dref.seq,iv,&dummy); - *it1=color(color_T,free_idx,rl); - *(it1+1)=color(color_ONE,rl); - *it2=numeric(5)/numeric(6); - return nonsimplified_ncmul(recombine_color_string( - delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec)); - } - } - - // f_{a,b,c} T_b T_c = 3/2 I T_a - for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) { - GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f); - const color & fref=ex_to_color(*it2); - exvector iv_intersect=idx_intersect(fref.seq,iv); - if (iv_intersect.size()==2) { - int sig; - ex free_idx=permute_free_index_to_front(fref.seq,iv,&sig); - *it1=color(color_T,free_idx,rl); - *(it1+1)=color(color_ONE,rl); - *it2=numeric(sig*3)/numeric(2)*I; - return nonsimplified_ncmul(recombine_color_string( - delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec)); - } - } - } + return false; +} + +/** Contraction of an indexed symmetric structure constant with something else. */ +bool su3d::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const +{ + GINAC_ASSERT(is_a(*self)); + GINAC_ASSERT(is_a(*other)); + GINAC_ASSERT(self->nops() == 4); + GINAC_ASSERT(is_a(self->op(0))); + + if (is_exactly_a(other->op(0))) { + + // Find the dummy indices of the contraction + exvector self_indices = ex_to(*self).get_indices(); + exvector other_indices = ex_to(*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 ita(a); + ita = set_difference(self_indices.begin(), self_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less()); + ita = set_difference(other_indices.begin(), other_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less()); + GINAC_ASSERT(a.size() == 2); + *self = numeric(5, 3) * delta_tensor(a[0], a[1]); + *other = _ex1; + return true; } - } - - // clear all ONEs when there is at least one corresponding color_T - // in this representation, retain one ONE otherwise - for (unsigned rl=0; rl(other->op(0))) { + + // d.abc T.b T.c = 5/6 T.a + if (other+1 != v.end() + && is_exactly_a(other[1].op(0)) + && ex_to(*self).has_dummy_index_for(other[1].op(1))) { + + exvector self_indices = ex_to(*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(other[0]).get_representation_label()); + other[1] = _ex1; + return true; } } - // return a sorted vector - return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs, - ONEvecs,unknownvec)); + return false; } -ex color::thisexprseq(const exvector & v) const +/** Contraction of an indexed antisymmetric structure constant with something else. */ +bool su3f::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { - return color(type,v,representation_label); -} + GINAC_ASSERT(is_a(*self)); + GINAC_ASSERT(is_a(*other)); + GINAC_ASSERT(self->nops() == 4); + GINAC_ASSERT(is_a(self->op(0))); -ex color::thisexprseq(exvector * vp) const -{ - return color(type,vp,representation_label); -} + if (is_exactly_a(other->op(0))) { // f*d is handled by su3d class -/** Check whether all indices are of class coloridx or a subclass. This - * function is used internally to make sure that all constructed color - * objects really carry color indices and not some other classes. */ -bool color::all_of_type_coloridx(void) const -{ - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - if (!is_ex_of_type(*cit,coloridx)) return false; + // Find the dummy indices of the contraction + exvector dummy_indices; + dummy_indices = ex_to(*self).get_dummy_indices(ex_to(*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(*self).get_indices(), dummy_indices, sign1); + ex b = permute_free_index_to_front(ex_to(*other).get_indices(), dummy_indices, sign2); + *self = sign1 * sign2 * 3 * delta_tensor(a, b); + *other = _ex1; + return true; + } + + } else if (is_exactly_a(other->op(0))) { + + // f.abc T.b T.c = 3/2 I T.a + if (other+1 != v.end() + && is_exactly_a(other[1].op(0)) + && ex_to(*self).has_dummy_index_for(other[1].op(1))) { + + exvector self_indices = ex_to(*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(other[0]).get_representation_label()); + other[1] = _ex1; + return true; + } } - return true; + + return false; } ////////// -// friend functions +// global functions ////////// -/** Construct an object representing the unity element of su(3). - * - * @param rl Representation label - * @return newly constructed object */ -color color_ONE(unsigned rl) +ex color_ONE(unsigned char rl) { - return color(color::color_ONE,rl); + static ex ONE = (new su3one)->setflag(status_flags::dynallocated); + return color(ONE, rl); } -/** Construct an object representing the generators T_a of SU(3). The index - * must be of class coloridx. - * - * @param a Index - * @param rl Representation label - * @return newly constructed object */ -color color_T(const ex & a, unsigned rl) +ex color_T(const ex & a, unsigned char rl) { - return color(color::color_T,a,rl); + static ex t = (new su3t)->setflag(status_flags::dynallocated); + + if (!is_a(a)) + throw(std::invalid_argument("indices of color_T must be of type idx")); + if (!ex_to(a).get_dim().is_equal(8)) + throw(std::invalid_argument("index dimension for color_T must be 8")); + + return color(t, a, rl); } -/** Construct an object representing the antisymmetric structure constants - * f_abc of SU(3). The indices must be of class coloridx. - * - * @param a First index - * @param b Second index - * @param c Third index - * @return newly constructed object */ -color color_f(const ex & a, const ex & b, const ex & c) +ex color_f(const ex & a, const ex & b, const ex & c) { - return color(color::color_f,a,b,c); + static ex f = (new su3f)->setflag(status_flags::dynallocated); + + if (!is_a(a) || !is_a(b) || !is_a(c)) + throw(std::invalid_argument("indices of color_f must be of type idx")); + if (!ex_to(a).get_dim().is_equal(8) || !ex_to(b).get_dim().is_equal(8) || !ex_to(c).get_dim().is_equal(8)) + throw(std::invalid_argument("index dimension for color_f must be 8")); + + return indexed(f, antisymmetric3(), a, b, c); } -/** Construct an object representing the symmetric structure constants d_abc - * of SU(3). The indices must be of class coloridx. - * - * @param a First index - * @param b Second index - * @param c Third index - * @return newly constructed object */ -color color_d(const ex & a, const ex & b, const ex & c) +ex color_d(const ex & a, const ex & b, const ex & c) { - return color(color::color_d,a,b,c); + static ex d = (new su3d)->setflag(status_flags::dynallocated); + + if (!is_a(a) || !is_a(b) || !is_a(c)) + throw(std::invalid_argument("indices of color_d must be of type idx")); + if (!ex_to(a).get_dim().is_equal(8) || !ex_to(b).get_dim().is_equal(8) || !ex_to(c).get_dim().is_equal(8)) + throw(std::invalid_argument("index dimension for color_d must be 8")); + + return indexed(d, symmetric3(), a, b, c); } -/** This returns the linear combination d_abc+I*f_abc. - * - * @param a First index - * @param b Second index - * @param c Third index - * @return newly constructed object */ ex color_h(const ex & a, const ex & b, const ex & c) { - return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c); + return color_d(a, b, c) + I * color_f(a, b, c); } -/** Construct an object representing the unity matrix delta8_ab in su(3). - * The indices must be of class coloridx. - * - * @param a First index - * @param b Second index - * @return newly constructed object */ -color color_delta8(const ex & a, const ex & b) +/** Check whether a given tinfo key (as returned by return_type_tinfo() + * is that of a color object (with an arbitrary representation label). */ +static bool is_color_tinfo(const return_type_t& ti) { - return color(color::color_delta8,a,b); + return *(ti.tinfo) == typeid(color); } -/** Given a vector of color (and possible other) objects, split it up - * according to the object type (structure constant, generator etc.) and - * representation label while preserving the order within each group. If - * there are non-color objetcs in the vector, the SU(3) generators T_a get - * sorted into the "unknown" group together with the non-color objects - * because we don't know whether these objects commute with the generators. - * - * @param v Source vector of expressions - * @param delta8vec Vector of unity matrices (returned) - * @param fvec Vector of antisymmetric structure constants (returned) - * @param dvec Vector of symmetric structure constants (returned) - * @param Tvecs Vectors of generators, one for each representation label (returned) - * @param ONEvecs Vectors of unity elements, one for each representation label (returned) - * @param unknownvec Vector of all non-color objects (returned) - * - * @see color::color_types - * @see recombine_color_string */ -void split_color_string_in_parts(const exvector & v, exvector & delta8vec, - exvector & fvec, exvector & dvec, - exvectorvector & Tvecs, - exvectorvector & ONEvecs, - exvector & unknownvec) +/** Extract representation label from tinfo key (as returned by + * return_type_tinfo()). */ +static unsigned char get_representation_label(const return_type_t& ti) { - // if not all elements are of type color, put all Ts in unknownvec to - // retain the ordering - bool all_color=true; - for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) { - if (!is_ex_exactly_of_type(*cit,color)) { - all_color=false; - break; - } - } - - for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) { - if (is_ex_exactly_of_type(*cit,color)) { - switch (ex_to_color(*cit).type) { - case color::color_delta8: - delta8vec.push_back(*cit); - break; - case color::color_f: - fvec.push_back(*cit); - break; - case color::color_d: - dvec.push_back(*cit); - break; - case color::color_T: - GINAC_ASSERT(ex_to_color(*cit).representation_label & rls) { - if (v.size()==0) { - return numeric(COLOR_THREE); - } else if (v.size()==1) { - GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color)); - return _ex0(); - } - exvector v1=v; - ex last_element=v1.back(); - GINAC_ASSERT(is_ex_exactly_of_type(last_element,color)); - GINAC_ASSERT(ex_to_color(last_element).type==color::color_T); - v1.pop_back(); - ex next_to_last_element=v1.back(); - GINAC_ASSERT(is_ex_exactly_of_type(next_to_last_element,color)); - GINAC_ASSERT(ex_to_color(next_to_last_element).type==color::color_T); - v1.pop_back(); - exvector v2=v1; - - const ex & last_index=ex_to_color(last_element).seq[0]; - const ex & next_to_last_index=ex_to_color(next_to_last_element).seq[0]; - ex summation_index=coloridx(); - - v2.push_back(color_T(summation_index)); // don't care about the representation_label - - // FIXME: check this formula for SU(N) with N!=3 - return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index) - % color_trace_of_one_representation_label(v1) - +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index) - % color_trace_of_one_representation_label(v2); - /* - ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index) - % color_trace_of_one_representation_label(v1); - cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl; - ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index) - % color_trace_of_one_representation_label(v2); - cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl; - return term1+term2; - */ -} + if (is_a(e)) { -/** Calculate the trace over the (hidden) indices of the su(3) Lie algebra - * elements (the SU(3) generators and the unity element) of a specified - * representation label in a string of color objects. - * - * @param v Vector of color objects - * @param rl Representation label - * @return value of the trace */ -ex color_trace(const exvector & v, unsigned rl) -{ - GINAC_ASSERT(rl(e).get_representation_label(); - append_exvector_to_exvector(v_rest,delta8vec); - append_exvector_to_exvector(v_rest,fvec); - append_exvector_to_exvector(v_rest,dvec); - for (unsigned i=0; i(e.op(0))) + return _ex3; + else + return _ex0; -ex simplify_pure_color_string(const ex & e) -{ - GINAC_ASSERT(is_ex_exactly_of_type(e,ncmul)); - - exvector delta8vec; - exvector fvec; - exvector dvec; - exvectorvector Tvecs; - Tvecs.resize(MAX_REPRESENTATION_LABELS); - exvectorvector ONEvecs; - ONEvecs.resize(MAX_REPRESENTATION_LABELS); - exvector unknownvec; - - split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec); - - // search for T_k S T_k (=1/2 Tr(S) - 1/6 S) - for (unsigned rl=0; rl=2) { - for (unsigned i=0; i(e)) { - // simplification of commutative product=commutative product of simplifications - if (is_ex_exactly_of_type(e_expanded,mul)) { - ex prod=_ex1(); - for (unsigned i=0; i(e)) { + + 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(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 0) { + + // Trace maps to all other container classes (this includes sums) + pointer_to_map_function_1arg &> fcn(color_trace, rls); + return e.map(fcn); + + } else + return _ex0; } -ex brute_force_sum_color_indices(const ex & e) +ex color_trace(const ex & e, const lst & rll) { - exvector iv_all=e.get_indices(); - exvector iv_double; - - // find double symbolic indices - if (iv_all.size()<2) return e; - for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) { - GINAC_ASSERT(is_ex_of_type(*cit1,coloridx)); - for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) { - GINAC_ASSERT(is_ex_of_type(*cit2,coloridx)); - if (ex_to_coloridx(*cit1).is_symbolic() && - ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) { - iv_double.push_back(*cit1); - break; - } - } + // Convert list to set + std::set rls; + for (lst::const_iterator i = rll.begin(); i != rll.end(); ++i) { + if (i->info(info_flags::nonnegint)) + rls.insert(ex_to(*i).to_int()); } - std::vector counter; - counter.resize(iv_double.size()); - int l; - for (l=0; unsigned(l)=0)&&((++counter[l])>(int)COLOR_EIGHT)) { - counter[l]=1; - l--; - } - if (l<2) { std::cout << counter[0] << counter[1] << std::endl; } - if (l<0) break; - } - - return sum; +ex color_trace(const ex & e, unsigned char rl) +{ + // Convert label to set + std::set rls; + rls.insert(rl); + + return color_trace(e, rls); } } // namespace GiNaC