3 * Implementation of GiNaC's color objects.
4 * No real implementation yet, to be done. */
7 * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
35 #include "relational.h"
42 GINAC_IMPLEMENT_REGISTERED_CLASS(color, indexed)
45 // default constructor, destructor, copy constructor assignment operator and helpers
50 color::color() : inherited(TINFO_color), type(invalid), representation_label(0)
52 debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
57 void color::copy(const color & other)
59 inherited::copy(other);
61 representation_label=other.representation_label;
64 void color::destroy(bool call_parent)
67 inherited::destroy(call_parent);
77 /** Construct object without any color index. This constructor is for internal
78 * use only. Use the color_ONE() function instead.
80 color::color(color_types const t, unsigned rl) : type(t), representation_label(rl)
82 debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT);
83 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
84 tinfo_key=TINFO_color;
85 GINAC_ASSERT(all_of_type_coloridx());
88 /** Construct object with one color index. This constructor is for internal
89 * use only. Use the color_T() function instead.
91 color::color(color_types const t, const ex & i1, unsigned rl)
92 : inherited(i1), type(t), representation_label(rl)
94 debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
95 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
96 tinfo_key=TINFO_color;
97 GINAC_ASSERT(all_of_type_coloridx());
100 /** Construct object with two color indices. This constructor is for internal
101 * use only. Use the color_delta8() function instead.
102 * @see color_delta8 */
103 color::color(color_types const t, const ex & i1, const ex & i2, unsigned rl)
104 : inherited(i1,i2), type(t), representation_label(rl)
106 debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
107 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
108 tinfo_key=TINFO_color;
109 GINAC_ASSERT(all_of_type_coloridx());
112 /** Construct object with three color indices. This constructor is for internal
113 * use only. Use the color_f(), color_d() and color_h() functions instead.
117 color::color(color_types const t, const ex & i1, const ex & i2, const ex & i3,
118 unsigned rl) : inherited(i1,i2,i3), type(t), representation_label(rl)
120 debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
121 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
122 tinfo_key=TINFO_color;
123 GINAC_ASSERT(all_of_type_coloridx());
126 /** Construct object with arbitrary number of color indices. This
127 * constructor is for internal use only. */
128 color::color(color_types const t, const exvector & iv, unsigned rl)
129 : inherited(iv), type(t), representation_label(rl)
131 debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
132 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
133 tinfo_key=TINFO_color;
134 GINAC_ASSERT(all_of_type_coloridx());
137 color::color(color_types const t, exvector * ivp, unsigned rl)
138 : inherited(ivp), type(t), representation_label(rl)
140 debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
141 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
142 tinfo_key=TINFO_color;
143 GINAC_ASSERT(all_of_type_coloridx());
150 /** Construct object from archive_node. */
151 color::color(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
153 debugmsg("color constructor from archive_node", LOGLEVEL_CONSTRUCT);
155 if (!(n.find_unsigned("type", ty)))
156 throw (std::runtime_error("unknown color type in archive"));
157 type = (color_types)ty;
158 if (!(n.find_unsigned("representation", representation_label)))
159 throw (std::runtime_error("unknown color representation label in archive"));
162 /** Unarchive the object. */
163 ex color::unarchive(const archive_node &n, const lst &sym_lst)
165 return (new color(n, sym_lst))->setflag(status_flags::dynallocated);
168 /** Archive the object. */
169 void color::archive(archive_node &n) const
171 inherited::archive(n);
172 n.add_unsigned("type", type);
173 n.add_unsigned("representation", representation_label);
177 // functions overriding virtual functions from bases classes
182 void color::printraw(std::ostream & os) const
184 debugmsg("color printraw",LOGLEVEL_PRINT);
185 os << "color(type=" << (unsigned)type
186 << ",representation_label=" << representation_label
189 os << ",hash=" << hashvalue << ",flags=" << flags << ")";
192 void color::printtree(std::ostream & os, unsigned indent) const
194 debugmsg("color printtree",LOGLEVEL_PRINT);
195 os << std::string(indent,' ') << "color object: "
196 << "type=" << (unsigned)type
197 << ",representation_label=" << representation_label << ", ";
198 os << seq.size() << " indices" << std::endl;
199 printtreeindices(os,indent);
200 os << std::string(indent,' ') << "hash=" << hashvalue
201 << " (0x" << std::hex << hashvalue << std::dec << ")"
202 << ", flags=" << flags << std::endl;
205 void color::print(std::ostream & os, unsigned upper_precedence) const
207 debugmsg("color print",LOGLEVEL_PRINT);
211 if (representation_label!=0) {
212 os << "^(" << representation_label << ")";
229 os << "INVALID_COLOR_OBJECT";
235 bool color::info(unsigned inf) const
237 return inherited::info(inf);
240 #define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
242 ex color::eval(int level) const
244 // canonicalize indices
246 bool antisymmetric=false;
250 antisymmetric=true; // no break here!
255 int sig=canonicalize_indices(iv,antisymmetric);
257 // something has changed while sorting indices, more evaluations later
258 if (sig==0) return _ex0();
259 return ex(sig)*color(type,iv,representation_label);
264 // nothing to canonicalize
271 GINAC_ASSERT(seq.size()==2);
272 const coloridx & idx1=ex_to_coloridx(seq[0]);
273 const coloridx & idx2=ex_to_coloridx(seq[1]);
275 // check for delta8_{a,a} where a is a symbolic index, replace by 8
276 if ((idx1.is_symbolic())&&(idx1.is_equal_same_type(idx2))) {
277 return ex(COLOR_EIGHT);
280 // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1
281 if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
282 if ((idx1.get_value()!=idx2.get_value())) {
291 // check for d_{a,a,c} (=0) when a is symbolic
293 GINAC_ASSERT(seq.size()==3);
294 const coloridx & idx1=ex_to_coloridx(seq[0]);
295 const coloridx & idx2=ex_to_coloridx(seq[1]);
296 const coloridx & idx3=ex_to_coloridx(seq[2]);
298 if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
300 } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
304 // check for three numeric indices
305 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
306 GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
307 GINAC_ASSERT(idx2.get_value()<=idx3.get_value());
308 if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)||
309 CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) {
311 } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
313 } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) {
314 return 1/sqrt(numeric(3));
315 } else if (CMPINDICES(8,8,8)) {
316 return -1/sqrt(numeric(3));
317 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
318 return -1/(2*sqrt(numeric(3)));
326 GINAC_ASSERT(seq.size()==3);
327 const coloridx & idx1=ex_to_coloridx(seq[0]);
328 const coloridx & idx2=ex_to_coloridx(seq[1]);
329 const coloridx & idx3=ex_to_coloridx(seq[2]);
331 // check for three numeric indices
332 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
333 GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
334 GINAC_ASSERT(idx2.get_value()<=idx3.get_value());
335 if (CMPINDICES(1,2,3)) {
337 } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
338 CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
340 } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
342 } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) {
343 return sqrt(numeric(3))/2;
344 } else if (CMPINDICES(8,8,8)) {
345 return -1/sqrt(numeric(3));
346 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
347 return -1/(2*sqrt(numeric(3)));
354 // nothing to evaluate
363 int color::compare_same_type(const basic & other) const
365 GINAC_ASSERT(other.tinfo() == TINFO_color);
366 const color &o = static_cast<const color &>(other);
368 if (type != o.type) {
370 return type < o.type ? -1 : 1;
373 if (representation_label != o.representation_label) {
374 // different representation label
375 return representation_label < o.representation_label ? -1 : 1;
378 return inherited::compare_same_type(other);
381 bool color::is_equal_same_type(const basic & other) const
383 GINAC_ASSERT(other.tinfo() == TINFO_color);
384 const color &o = static_cast<const color &>(other);
386 if (type != o.type) return false;
387 if (representation_label != o.representation_label) return false;
388 return inherited::is_equal_same_type(other);
393 ex color::simplify_ncmul(const exvector & v) const
395 // simplifications: contract delta8_{a,b} where possible
396 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...
397 // remove superfluous ONEs
399 // contract indices of delta8_{a,b} if they are different and symbolic
401 exvector v_contracted=v;
402 unsigned replacements;
403 bool something_changed=false;
405 exvector::iterator it=v_contracted.begin();
406 while (it!=v_contracted.end()) {
407 // process only delta8 objects
408 if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) {
409 color & d8=ex_to_nonconst_color(*it);
410 GINAC_ASSERT(d8.seq.size()==2);
411 const coloridx & first_idx=ex_to_coloridx(d8.seq[0]);
412 const coloridx & second_idx=ex_to_coloridx(d8.seq[1]);
413 // delta8_{a,a} should have been contracted in color::eval()
414 GINAC_ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic()));
415 ex saved_delta8=*it; // save to restore it later
417 // try to contract first index
419 if (first_idx.is_symbolic()) {
420 replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx);
421 if (replacements==1) {
422 // not contracted except in itself, restore delta8 object
425 // a contracted index should occur exactly twice
426 GINAC_ASSERT(replacements==2);
428 something_changed=true;
432 // try second index only if first was not contracted
433 if ((replacements==1)&&(second_idx.is_symbolic())) {
434 // first index not contracted, *it is guaranteed to be the original delta8 object
435 replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx);
436 if (replacements==1) {
437 // not contracted except in itself, restore delta8 object
440 // a contracted index should occur exactly twice
441 GINAC_ASSERT(replacements==2);
443 something_changed=true;
450 if (something_changed) {
451 // do more simplifications later
452 return nonsimplified_ncmul(v_contracted);
455 // there were no indices to contract
456 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown
457 // (if there is at least one unknown object, all Ts will be unknown to not change the order)
462 exvectorvector Tvecs;
463 Tvecs.resize(MAX_REPRESENTATION_LABELS);
464 exvectorvector ONEvecs;
465 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
468 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
470 // d_{a,k,l} f_{b,k,l}=0 (includes case a=b)
471 if ((dvec.size()>=1)&&(fvec.size()>=1)) {
472 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) {
473 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
474 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
475 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
476 const color & col1=ex_to_color(*it1);
477 const color & col2=ex_to_color(*it2);
478 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
479 if (iv_intersect.size()>=2) return _ex0();
484 // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b)
485 if (dvec.size()>=2) {
486 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) {
487 for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) {
488 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
489 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
490 const color & col1=ex_to_color(*it1);
491 const color & col2=ex_to_color(*it2);
492 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
493 if (iv_intersect.size()>=2) {
494 if (iv_intersect.size()==3) {
495 *it1=numeric(40)/numeric(3);
498 int dummy; // sign unimportant, since symmetric
499 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&dummy);
500 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&dummy);
501 *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
504 return nonsimplified_ncmul(recombine_color_string(
505 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
511 // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b)
512 if (fvec.size()>=2) {
513 for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) {
514 for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) {
515 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
516 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
517 const color & col1=ex_to_color(*it1);
518 const color & col2=ex_to_color(*it2);
519 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
520 if (iv_intersect.size()>=2) {
521 if (iv_intersect.size()==3) {
526 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&sig1);
527 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&sig2);
528 *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
531 return nonsimplified_ncmul(recombine_color_string(
532 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
538 // d_{a,b,c} T_b T_c = 5/6 T_a
539 // f_{a,b,c} T_b T_c = 3/2 I T_a
540 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
541 if ((Tvecs[rl].size()>=2)&&((dvec.size()>=1)||(fvec.size()>=1))) {
542 for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) {
544 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
545 GINAC_ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T);
546 iv.push_back(ex_to_color(*it1).seq[0]);
547 iv.push_back(ex_to_color(*(it1+1)).seq[0]);
549 // d_{a,b,c} T_b T_c = 5/6 T_a
550 for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) {
551 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
552 const color & dref=ex_to_color(*it2);
553 exvector iv_intersect=idx_intersect(dref.seq,iv);
554 if (iv_intersect.size()==2) {
555 int dummy; // sign unimportant, since symmetric
556 ex free_idx=permute_free_index_to_front(dref.seq,iv,&dummy);
557 *it1=color(color_T,free_idx,rl);
558 *(it1+1)=color(color_ONE,rl);
559 *it2=numeric(5)/numeric(6);
560 return nonsimplified_ncmul(recombine_color_string(
561 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
565 // f_{a,b,c} T_b T_c = 3/2 I T_a
566 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
567 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
568 const color & fref=ex_to_color(*it2);
569 exvector iv_intersect=idx_intersect(fref.seq,iv);
570 if (iv_intersect.size()==2) {
572 ex free_idx=permute_free_index_to_front(fref.seq,iv,&sig);
573 *it1=color(color_T,free_idx,rl);
574 *(it1+1)=color(color_ONE,rl);
575 *it2=numeric(sig*3)/numeric(2)*I;
576 return nonsimplified_ncmul(recombine_color_string(
577 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
584 // clear all ONEs when there is at least one corresponding color_T
585 // in this representation, retain one ONE otherwise
586 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
587 if (Tvecs[rl].size()!=0) {
589 } else if (ONEvecs[rl].size()!=0) {
591 ONEvecs[rl].push_back(color(color_ONE,rl));
595 // return a sorted vector
596 return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
597 ONEvecs,unknownvec));
600 ex color::thisexprseq(const exvector & v) const
602 return color(type,v,representation_label);
605 ex color::thisexprseq(exvector * vp) const
607 return color(type,vp,representation_label);
610 /** Check whether all indices are of class coloridx or a subclass. This
611 * function is used internally to make sure that all constructed color
612 * objects really carry color indices and not some other classes. */
613 bool color::all_of_type_coloridx(void) const
615 for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
616 if (!is_ex_of_type(*cit,coloridx)) return false;
625 /** Construct an object representing the unity element of su(3).
627 * @param rl Representation label
628 * @return newly constructed object */
629 color color_ONE(unsigned rl)
631 return color(color::color_ONE,rl);
634 /** Construct an object representing the generators T_a of SU(3). The index
635 * must be of class coloridx.
638 * @param rl Representation label
639 * @return newly constructed object */
640 color color_T(const ex & a, unsigned rl)
642 return color(color::color_T,a,rl);
645 /** Construct an object representing the antisymmetric structure constants
646 * f_abc of SU(3). The indices must be of class coloridx.
648 * @param a First index
649 * @param b Second index
650 * @param c Third index
651 * @return newly constructed object */
652 color color_f(const ex & a, const ex & b, const ex & c)
654 return color(color::color_f,a,b,c);
657 /** Construct an object representing the symmetric structure constants d_abc
658 * of SU(3). The indices must be of class coloridx.
660 * @param a First index
661 * @param b Second index
662 * @param c Third index
663 * @return newly constructed object */
664 color color_d(const ex & a, const ex & b, const ex & c)
666 return color(color::color_d,a,b,c);
669 /** This returns the linear combination d_abc+I*f_abc.
671 * @param a First index
672 * @param b Second index
673 * @param c Third index
674 * @return newly constructed object */
675 ex color_h(const ex & a, const ex & b, const ex & c)
677 return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
680 /** Construct an object representing the unity matrix delta8_ab in su(3).
681 * The indices must be of class coloridx.
683 * @param a First index
684 * @param b Second index
685 * @return newly constructed object */
686 color color_delta8(const ex & a, const ex & b)
688 return color(color::color_delta8,a,b);
691 /** Given a vector of color (and possible other) objects, split it up
692 * according to the object type (structure constant, generator etc.) and
693 * representation label while preserving the order within each group. If
694 * there are non-color objetcs in the vector, the SU(3) generators T_a get
695 * sorted into the "unknown" group together with the non-color objects
696 * because we don't know whether these objects commute with the generators.
698 * @param v Source vector of expressions
699 * @param delta8vec Vector of unity matrices (returned)
700 * @param fvec Vector of antisymmetric structure constants (returned)
701 * @param dvec Vector of symmetric structure constants (returned)
702 * @param Tvecs Vectors of generators, one for each representation label (returned)
703 * @param ONEvecs Vectors of unity elements, one for each representation label (returned)
704 * @param unknownvec Vector of all non-color objects (returned)
706 * @see color::color_types
707 * @see recombine_color_string */
708 void split_color_string_in_parts(const exvector & v, exvector & delta8vec,
709 exvector & fvec, exvector & dvec,
710 exvectorvector & Tvecs,
711 exvectorvector & ONEvecs,
712 exvector & unknownvec)
714 // if not all elements are of type color, put all Ts in unknownvec to
715 // retain the ordering
717 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
718 if (!is_ex_exactly_of_type(*cit,color)) {
724 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
725 if (is_ex_exactly_of_type(*cit,color)) {
726 switch (ex_to_color(*cit).type) {
727 case color::color_delta8:
728 delta8vec.push_back(*cit);
731 fvec.push_back(*cit);
734 dvec.push_back(*cit);
737 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
739 Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
741 unknownvec.push_back(*cit);
744 case color::color_ONE:
745 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
746 ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
749 throw(std::logic_error("invalid type in split_color_string_in_parts()"));
752 unknownvec.push_back(*cit);
757 /** Merge vectors of color objects sorted by object type into one vector,
758 * retaining the order within each group. This is the inverse operation of
759 * split_color_string_in_parts().
761 * @param delta8vec Vector of unity matrices
762 * @param fvec Vector of antisymmetric structure constants
763 * @param dvec Vector of symmetric structure constants
764 * @param Tvecs Vectors of generators, one for each representation label
765 * @param ONEvecs Vectors of unity elements, one for each representation label
766 * @param unknownvec Vector of all non-color objects
767 * @return merged vector
769 * @see color::color_types
770 * @see split_color_string_in_parts */
771 exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
772 exvector & dvec, exvectorvector & Tvecs,
773 exvectorvector & ONEvecs, exvector & unknownvec)
775 unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
776 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
777 sz += Tvecs[rl].size();
778 sz += ONEvecs[rl].size();
783 append_exvector_to_exvector(v,delta8vec);
784 append_exvector_to_exvector(v,fvec);
785 append_exvector_to_exvector(v,dvec);
786 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
787 append_exvector_to_exvector(v,Tvecs[rl]);
788 append_exvector_to_exvector(v,ONEvecs[rl]);
790 append_exvector_to_exvector(v,unknownvec);
794 ex color_trace_of_one_representation_label(const exvector & v)
797 return numeric(COLOR_THREE);
798 } else if (v.size()==1) {
799 GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
803 ex last_element=v1.back();
804 GINAC_ASSERT(is_ex_exactly_of_type(last_element,color));
805 GINAC_ASSERT(ex_to_color(last_element).type==color::color_T);
807 ex next_to_last_element=v1.back();
808 GINAC_ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
809 GINAC_ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
813 const ex & last_index=ex_to_color(last_element).seq[0];
814 const ex & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
815 ex summation_index=coloridx();
817 v2.push_back(color_T(summation_index)); // don't care about the representation_label
819 // FIXME: check this formula for SU(N) with N!=3
820 return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
821 % color_trace_of_one_representation_label(v1)
822 +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
823 % color_trace_of_one_representation_label(v2);
825 ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
826 % color_trace_of_one_representation_label(v1);
827 cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
828 ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
829 % color_trace_of_one_representation_label(v2);
830 cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
835 /** Calculate the trace over the (hidden) indices of the su(3) Lie algebra
836 * elements (the SU(3) generators and the unity element) of a specified
837 * representation label in a string of color objects.
839 * @param v Vector of color objects
840 * @param rl Representation label
841 * @return value of the trace */
842 ex color_trace(const exvector & v, unsigned rl)
844 GINAC_ASSERT(rl<MAX_REPRESENTATION_LABELS);
847 v_rest.reserve(v.size()+1); // max size if trace is empty
852 exvectorvector Tvecs;
853 Tvecs.resize(MAX_REPRESENTATION_LABELS);
854 exvectorvector ONEvecs;
855 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
858 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
860 if (unknownvec.size()!=0) {
861 throw(std::invalid_argument("color_trace(): expression must be expanded"));
864 append_exvector_to_exvector(v_rest,delta8vec);
865 append_exvector_to_exvector(v_rest,fvec);
866 append_exvector_to_exvector(v_rest,dvec);
867 for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
869 append_exvector_to_exvector(v_rest,Tvecs[i]);
870 append_exvector_to_exvector(v_rest,ONEvecs[i]);
872 if (Tvecs[i].size()!=0) {
873 v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
874 } else if (ONEvecs[i].size()!=0) {
875 v_rest.push_back(numeric(COLOR_THREE));
877 throw(std::logic_error("color_trace(): representation_label not in color string"));
882 return nonsimplified_ncmul(v_rest);
885 ex simplify_pure_color_string(const ex & e)
887 GINAC_ASSERT(is_ex_exactly_of_type(e,ncmul));
892 exvectorvector Tvecs;
893 Tvecs.resize(MAX_REPRESENTATION_LABELS);
894 exvectorvector ONEvecs;
895 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
898 split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
900 // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
901 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
902 if (Tvecs[rl].size()>=2) {
903 for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
904 for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
905 ex & t1=Tvecs[rl][i];
906 ex & t2=Tvecs[rl][j];
907 GINAC_ASSERT(is_ex_exactly_of_type(t1,color)&&
908 (ex_to_color(t1).type==color::color_T)&&
909 (ex_to_color(t1).seq.size()==1));
910 GINAC_ASSERT(is_ex_exactly_of_type(t2,color)&&
911 (ex_to_color(t2).type==color::color_T)&&
912 (ex_to_color(t2).seq.size()==1));
913 const coloridx & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
914 const coloridx & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
916 if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
918 for (unsigned k=i+1; k<j; ++k) {
919 S.push_back(Tvecs[rl][k]);
923 ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
924 for (unsigned k=i+1; k<j; ++k) {
927 t1=color_trace_of_one_representation_label(S);
928 ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
929 return simplify_color(term1+term2);
936 // FIXME: higher contractions
941 /** Perform some simplifications on an expression containing color objects. */
942 ex simplify_color(const ex & e)
944 // all simplification is done on expanded objects
945 ex e_expanded=e.expand();
947 // simplification of sum=sum of simplifications
948 if (is_ex_exactly_of_type(e_expanded,add)) {
950 for (unsigned i=0; i<e_expanded.nops(); ++i)
951 sum += simplify_color(e_expanded.op(i));
956 // simplification of commutative product=commutative product of simplifications
957 if (is_ex_exactly_of_type(e_expanded,mul)) {
959 for (unsigned i=0; i<e_expanded.nops(); ++i)
960 prod *= simplify_color(e_expanded.op(i));
965 // simplification of noncommutative product: test if everything is color
966 if (is_ex_exactly_of_type(e_expanded,ncmul)) {
968 for (unsigned i=0; i<e_expanded.nops(); ++i) {
969 if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
975 return simplify_pure_color_string(e_expanded);
979 // cannot do anything
983 ex brute_force_sum_color_indices(const ex & e)
985 exvector iv_all=e.get_indices();
988 // find double symbolic indices
989 if (iv_all.size()<2) return e;
990 for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
991 GINAC_ASSERT(is_ex_of_type(*cit1,coloridx));
992 for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
993 GINAC_ASSERT(is_ex_of_type(*cit2,coloridx));
994 if (ex_to_coloridx(*cit1).is_symbolic() &&
995 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
996 iv_double.push_back(*cit1);
1002 std::vector<int> counter;
1003 counter.resize(iv_double.size());
1005 for (l=0; unsigned(l)<iv_double.size(); ++l) {
1013 for (l=0; unsigned(l)<iv_double.size(); ++l) {
1014 term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
1015 //iv_double[l].print(cout);
1016 //cout << " " << counter[l] << " ";
1021 // increment counter[]
1022 l = iv_double.size()-1;
1023 while ((l>=0)&&((++counter[l])>(int)COLOR_EIGHT)) {
1027 if (l<2) { std::cout << counter[0] << counter[1] << std::endl; }
1034 } // namespace GiNaC