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"
40 #ifndef NO_NAMESPACE_GINAC
42 #endif // ndef NO_NAMESPACE_GINAC
44 GINAC_IMPLEMENT_REGISTERED_CLASS(color, indexed)
47 // default constructor, destructor, copy constructor assignment operator and helpers
52 color::color() : inherited(TINFO_color), type(invalid), representation_label(0)
54 debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
59 debugmsg("color destructor",LOGLEVEL_DESTRUCT);
63 color::color(const color & other)
65 debugmsg("color copy constructor",LOGLEVEL_CONSTRUCT);
69 const color & color::operator=(const color & other)
71 debugmsg("color operator=",LOGLEVEL_ASSIGNMENT);
81 void color::copy(const color & other)
83 inherited::copy(other);
85 representation_label=other.representation_label;
88 void color::destroy(bool call_parent)
91 inherited::destroy(call_parent);
101 /** Construct object without any color index. This constructor is for internal
102 * use only. Use the color_ONE() function instead.
104 color::color(color_types const t, unsigned rl) : type(t), representation_label(rl)
106 debugmsg("color constructor from color_types,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 one color index. This constructor is for internal
113 * use only. Use the color_T() function instead.
115 color::color(color_types const t, const ex & i1, unsigned rl)
116 : inherited(i1), type(t), representation_label(rl)
118 debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
119 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
120 tinfo_key=TINFO_color;
121 GINAC_ASSERT(all_of_type_coloridx());
124 /** Construct object with two color indices. This constructor is for internal
125 * use only. Use the color_delta8() function instead.
126 * @see color_delta8 */
127 color::color(color_types const t, const ex & i1, const ex & i2, unsigned rl)
128 : inherited(i1,i2), type(t), representation_label(rl)
130 debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
131 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
132 tinfo_key=TINFO_color;
133 GINAC_ASSERT(all_of_type_coloridx());
136 /** Construct object with three color indices. This constructor is for internal
137 * use only. Use the color_f(), color_d() and color_h() functions instead.
141 color::color(color_types const t, const ex & i1, const ex & i2, const ex & i3,
142 unsigned rl) : inherited(i1,i2,i3), type(t), representation_label(rl)
144 debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
145 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
146 tinfo_key=TINFO_color;
147 GINAC_ASSERT(all_of_type_coloridx());
150 /** Construct object with arbitrary number of color indices. This
151 * constructor is for internal use only. */
152 color::color(color_types const t, const exvector & iv, unsigned rl)
153 : inherited(iv), type(t), representation_label(rl)
155 debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
156 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
157 tinfo_key=TINFO_color;
158 GINAC_ASSERT(all_of_type_coloridx());
161 color::color(color_types const t, exvector * ivp, unsigned rl)
162 : inherited(ivp), type(t), representation_label(rl)
164 debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
165 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
166 tinfo_key=TINFO_color;
167 GINAC_ASSERT(all_of_type_coloridx());
174 /** Construct object from archive_node. */
175 color::color(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
177 debugmsg("color constructor from archive_node", LOGLEVEL_CONSTRUCT);
179 if (!(n.find_unsigned("type", ty)))
180 throw (std::runtime_error("unknown color type in archive"));
181 type = (color_types)ty;
182 if (!(n.find_unsigned("representation", representation_label)))
183 throw (std::runtime_error("unknown color representation label in archive"));
186 /** Unarchive the object. */
187 ex color::unarchive(const archive_node &n, const lst &sym_lst)
189 return (new color(n, sym_lst))->setflag(status_flags::dynallocated);
192 /** Archive the object. */
193 void color::archive(archive_node &n) const
195 inherited::archive(n);
196 n.add_unsigned("type", type);
197 n.add_unsigned("representation", representation_label);
201 // functions overriding virtual functions from bases classes
206 basic * color::duplicate() const
208 debugmsg("color duplicate",LOGLEVEL_DUPLICATE);
209 return new color(*this);
212 void color::printraw(std::ostream & os) const
214 debugmsg("color printraw",LOGLEVEL_PRINT);
215 os << "color(type=" << (unsigned)type
216 << ",representation_label=" << representation_label
219 os << ",hash=" << hashvalue << ",flags=" << flags << ")";
222 void color::printtree(std::ostream & os, unsigned indent) const
224 debugmsg("color printtree",LOGLEVEL_PRINT);
225 os << std::string(indent,' ') << "color object: "
226 << "type=" << (unsigned)type
227 << ",representation_label=" << representation_label << ", ";
228 os << seq.size() << " indices" << std::endl;
229 printtreeindices(os,indent);
230 os << std::string(indent,' ') << "hash=" << hashvalue
231 << " (0x" << std::hex << hashvalue << std::dec << ")"
232 << ", flags=" << flags << std::endl;
235 void color::print(std::ostream & os, unsigned upper_precedence) const
237 debugmsg("color print",LOGLEVEL_PRINT);
241 if (representation_label!=0) {
242 os << "^(" << representation_label << ")";
259 os << "INVALID_COLOR_OBJECT";
265 bool color::info(unsigned inf) const
267 return inherited::info(inf);
270 #define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
272 ex color::eval(int level) const
274 // canonicalize indices
276 bool antisymmetric=false;
280 antisymmetric=true; // no break here!
285 int sig=canonicalize_indices(iv,antisymmetric);
287 // something has changed while sorting indices, more evaluations later
288 if (sig==0) return _ex0();
289 return ex(sig)*color(type,iv,representation_label);
294 // nothing to canonicalize
301 GINAC_ASSERT(seq.size()==2);
302 const coloridx & idx1=ex_to_coloridx(seq[0]);
303 const coloridx & idx2=ex_to_coloridx(seq[1]);
305 // check for delta8_{a,a} where a is a symbolic index, replace by 8
306 if ((idx1.is_symbolic())&&(idx1.is_equal_same_type(idx2))) {
307 return ex(COLOR_EIGHT);
310 // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1
311 if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
312 if ((idx1.get_value()!=idx2.get_value())) {
321 // check for d_{a,a,c} (=0) when a is symbolic
323 GINAC_ASSERT(seq.size()==3);
324 const coloridx & idx1=ex_to_coloridx(seq[0]);
325 const coloridx & idx2=ex_to_coloridx(seq[1]);
326 const coloridx & idx3=ex_to_coloridx(seq[2]);
328 if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
330 } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
334 // check for three numeric indices
335 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
336 GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
337 GINAC_ASSERT(idx2.get_value()<=idx3.get_value());
338 if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)||
339 CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) {
341 } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
343 } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) {
344 return 1/sqrt(numeric(3));
345 } else if (CMPINDICES(8,8,8)) {
346 return -1/sqrt(numeric(3));
347 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
348 return -1/(2*sqrt(numeric(3)));
356 GINAC_ASSERT(seq.size()==3);
357 const coloridx & idx1=ex_to_coloridx(seq[0]);
358 const coloridx & idx2=ex_to_coloridx(seq[1]);
359 const coloridx & idx3=ex_to_coloridx(seq[2]);
361 // check for three numeric indices
362 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
363 GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
364 GINAC_ASSERT(idx2.get_value()<=idx3.get_value());
365 if (CMPINDICES(1,2,3)) {
367 } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
368 CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
370 } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
372 } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) {
373 return sqrt(numeric(3))/2;
374 } else if (CMPINDICES(8,8,8)) {
375 return -1/sqrt(numeric(3));
376 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
377 return -1/(2*sqrt(numeric(3)));
384 // nothing to evaluate
393 int color::compare_same_type(const basic & other) const
395 GINAC_ASSERT(other.tinfo() == TINFO_color);
396 const color &o = static_cast<const color &>(other);
398 if (type != o.type) {
400 return type < o.type ? -1 : 1;
403 if (representation_label != o.representation_label) {
404 // different representation label
405 return representation_label < o.representation_label ? -1 : 1;
408 return inherited::compare_same_type(other);
411 bool color::is_equal_same_type(const basic & other) const
413 GINAC_ASSERT(other.tinfo() == TINFO_color);
414 const color &o = static_cast<const color &>(other);
416 if (type != o.type) return false;
417 if (representation_label != o.representation_label) return false;
418 return inherited::is_equal_same_type(other);
423 ex color::simplify_ncmul(const exvector & v) const
425 // simplifications: contract delta8_{a,b} where possible
426 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...
427 // remove superfluous ONEs
429 // contract indices of delta8_{a,b} if they are different and symbolic
431 exvector v_contracted=v;
432 unsigned replacements;
433 bool something_changed=false;
435 exvector::iterator it=v_contracted.begin();
436 while (it!=v_contracted.end()) {
437 // process only delta8 objects
438 if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) {
439 color & d8=ex_to_nonconst_color(*it);
440 GINAC_ASSERT(d8.seq.size()==2);
441 const coloridx & first_idx=ex_to_coloridx(d8.seq[0]);
442 const coloridx & second_idx=ex_to_coloridx(d8.seq[1]);
443 // delta8_{a,a} should have been contracted in color::eval()
444 GINAC_ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic()));
445 ex saved_delta8=*it; // save to restore it later
447 // try to contract first index
449 if (first_idx.is_symbolic()) {
450 replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx);
451 if (replacements==1) {
452 // not contracted except in itself, restore delta8 object
455 // a contracted index should occur exactly twice
456 GINAC_ASSERT(replacements==2);
458 something_changed=true;
462 // try second index only if first was not contracted
463 if ((replacements==1)&&(second_idx.is_symbolic())) {
464 // first index not contracted, *it is guaranteed to be the original delta8 object
465 replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx);
466 if (replacements==1) {
467 // not contracted except in itself, restore delta8 object
470 // a contracted index should occur exactly twice
471 GINAC_ASSERT(replacements==2);
473 something_changed=true;
480 if (something_changed) {
481 // do more simplifications later
482 return nonsimplified_ncmul(v_contracted);
485 // there were no indices to contract
486 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown
487 // (if there is at least one unknown object, all Ts will be unknown to not change the order)
492 exvectorvector Tvecs;
493 Tvecs.resize(MAX_REPRESENTATION_LABELS);
494 exvectorvector ONEvecs;
495 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
498 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
500 // d_{a,k,l} f_{b,k,l}=0 (includes case a=b)
501 if ((dvec.size()>=1)&&(fvec.size()>=1)) {
502 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) {
503 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
504 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
505 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
506 const color & col1=ex_to_color(*it1);
507 const color & col2=ex_to_color(*it2);
508 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
509 if (iv_intersect.size()>=2) return _ex0();
514 // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b)
515 if (dvec.size()>=2) {
516 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) {
517 for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) {
518 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
519 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
520 const color & col1=ex_to_color(*it1);
521 const color & col2=ex_to_color(*it2);
522 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
523 if (iv_intersect.size()>=2) {
524 if (iv_intersect.size()==3) {
525 *it1=numeric(40)/numeric(3);
528 int dummy; // sign unimportant, since symmetric
529 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&dummy);
530 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&dummy);
531 *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
534 return nonsimplified_ncmul(recombine_color_string(
535 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
541 // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b)
542 if (fvec.size()>=2) {
543 for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) {
544 for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) {
545 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
546 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
547 const color & col1=ex_to_color(*it1);
548 const color & col2=ex_to_color(*it2);
549 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
550 if (iv_intersect.size()>=2) {
551 if (iv_intersect.size()==3) {
556 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&sig1);
557 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&sig2);
558 *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
561 return nonsimplified_ncmul(recombine_color_string(
562 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
568 // d_{a,b,c} T_b T_c = 5/6 T_a
569 // f_{a,b,c} T_b T_c = 3/2 I T_a
570 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
571 if ((Tvecs[rl].size()>=2)&&((dvec.size()>=1)||(fvec.size()>=1))) {
572 for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) {
574 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
575 GINAC_ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T);
576 iv.push_back(ex_to_color(*it1).seq[0]);
577 iv.push_back(ex_to_color(*(it1+1)).seq[0]);
579 // d_{a,b,c} T_b T_c = 5/6 T_a
580 for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) {
581 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
582 const color & dref=ex_to_color(*it2);
583 exvector iv_intersect=idx_intersect(dref.seq,iv);
584 if (iv_intersect.size()==2) {
585 int dummy; // sign unimportant, since symmetric
586 ex free_idx=permute_free_index_to_front(dref.seq,iv,&dummy);
587 *it1=color(color_T,free_idx,rl);
588 *(it1+1)=color(color_ONE,rl);
589 *it2=numeric(5)/numeric(6);
590 return nonsimplified_ncmul(recombine_color_string(
591 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
595 // f_{a,b,c} T_b T_c = 3/2 I T_a
596 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
597 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
598 const color & fref=ex_to_color(*it2);
599 exvector iv_intersect=idx_intersect(fref.seq,iv);
600 if (iv_intersect.size()==2) {
602 ex free_idx=permute_free_index_to_front(fref.seq,iv,&sig);
603 *it1=color(color_T,free_idx,rl);
604 *(it1+1)=color(color_ONE,rl);
605 *it2=numeric(sig*3)/numeric(2)*I;
606 return nonsimplified_ncmul(recombine_color_string(
607 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
614 // clear all ONEs when there is at least one corresponding color_T
615 // in this representation, retain one ONE otherwise
616 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
617 if (Tvecs[rl].size()!=0) {
619 } else if (ONEvecs[rl].size()!=0) {
621 ONEvecs[rl].push_back(color(color_ONE,rl));
625 // return a sorted vector
626 return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
627 ONEvecs,unknownvec));
630 ex color::thisexprseq(const exvector & v) const
632 return color(type,v,representation_label);
635 ex color::thisexprseq(exvector * vp) const
637 return color(type,vp,representation_label);
640 /** Check whether all indices are of class coloridx or a subclass. This
641 * function is used internally to make sure that all constructed color
642 * objects really carry color indices and not some other classes. */
643 bool color::all_of_type_coloridx(void) const
645 for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
646 if (!is_ex_of_type(*cit,coloridx)) return false;
652 // virtual functions which can be overridden by derived classes
658 // non-virtual functions in this class
662 // static member variables
671 const color some_color;
672 const std::type_info & typeid_color = typeid(some_color);
678 /** Construct an object representing the unity element of su(3).
680 * @param rl Representation label
681 * @return newly constructed object */
682 color color_ONE(unsigned rl)
684 return color(color::color_ONE,rl);
687 /** Construct an object representing the generators T_a of SU(3). The index
688 * must be of class coloridx.
691 * @param rl Representation label
692 * @return newly constructed object */
693 color color_T(const ex & a, unsigned rl)
695 return color(color::color_T,a,rl);
698 /** Construct an object representing the antisymmetric structure constants
699 * f_abc of SU(3). The indices must be of class coloridx.
701 * @param a First index
702 * @param b Second index
703 * @param c Third index
704 * @return newly constructed object */
705 color color_f(const ex & a, const ex & b, const ex & c)
707 return color(color::color_f,a,b,c);
710 /** Construct an object representing the symmetric structure constants d_abc
711 * of SU(3). The indices must be of class coloridx.
713 * @param a First index
714 * @param b Second index
715 * @param c Third index
716 * @return newly constructed object */
717 color color_d(const ex & a, const ex & b, const ex & c)
719 return color(color::color_d,a,b,c);
722 /** This returns the linear combination d_abc+I*f_abc.
724 * @param a First index
725 * @param b Second index
726 * @param c Third index
727 * @return newly constructed object */
728 ex color_h(const ex & a, const ex & b, const ex & c)
730 return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
733 /** Construct an object representing the unity matrix delta8_ab in su(3).
734 * The indices must be of class coloridx.
736 * @param a First index
737 * @param b Second index
738 * @return newly constructed object */
739 color color_delta8(const ex & a, const ex & b)
741 return color(color::color_delta8,a,b);
744 /** Given a vector of color (and possible other) objects, split it up
745 * according to the object type (structure constant, generator etc.) and
746 * representation label while preserving the order within each group. If
747 * there are non-color objetcs in the vector, the SU(3) generators T_a get
748 * sorted into the "unknown" group together with the non-color objects
749 * because we don't know whether these objects commute with the generators.
751 * @param v Source vector of expressions
752 * @param delta8vec Vector of unity matrices (returned)
753 * @param fvec Vector of antisymmetric structure constants (returned)
754 * @param dvec Vector of symmetric structure constants (returned)
755 * @param Tvecs Vectors of generators, one for each representation label (returned)
756 * @param ONEvecs Vectors of unity elements, one for each representation label (returned)
757 * @param unknownvec Vector of all non-color objects (returned)
759 * @see color::color_types
760 * @see recombine_color_string */
761 void split_color_string_in_parts(const exvector & v, exvector & delta8vec,
762 exvector & fvec, exvector & dvec,
763 exvectorvector & Tvecs,
764 exvectorvector & ONEvecs,
765 exvector & unknownvec)
767 // if not all elements are of type color, put all Ts in unknownvec to
768 // retain the ordering
770 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
771 if (!is_ex_exactly_of_type(*cit,color)) {
777 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
778 if (is_ex_exactly_of_type(*cit,color)) {
779 switch (ex_to_color(*cit).type) {
780 case color::color_delta8:
781 delta8vec.push_back(*cit);
784 fvec.push_back(*cit);
787 dvec.push_back(*cit);
790 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
792 Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
794 unknownvec.push_back(*cit);
797 case color::color_ONE:
798 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
799 ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
802 throw(std::logic_error("invalid type in split_color_string_in_parts()"));
805 unknownvec.push_back(*cit);
810 /** Merge vectors of color objects sorted by object type into one vector,
811 * retaining the order within each group. This is the inverse operation of
812 * split_color_string_in_parts().
814 * @param delta8vec Vector of unity matrices
815 * @param fvec Vector of antisymmetric structure constants
816 * @param dvec Vector of symmetric structure constants
817 * @param Tvecs Vectors of generators, one for each representation label
818 * @param ONEvecs Vectors of unity elements, one for each representation label
819 * @param unknownvec Vector of all non-color objects
820 * @return merged vector
822 * @see color::color_types
823 * @see split_color_string_in_parts */
824 exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
825 exvector & dvec, exvectorvector & Tvecs,
826 exvectorvector & ONEvecs, exvector & unknownvec)
828 unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
829 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
830 sz += Tvecs[rl].size();
831 sz += ONEvecs[rl].size();
836 append_exvector_to_exvector(v,delta8vec);
837 append_exvector_to_exvector(v,fvec);
838 append_exvector_to_exvector(v,dvec);
839 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
840 append_exvector_to_exvector(v,Tvecs[rl]);
841 append_exvector_to_exvector(v,ONEvecs[rl]);
843 append_exvector_to_exvector(v,unknownvec);
847 ex color_trace_of_one_representation_label(const exvector & v)
850 return numeric(COLOR_THREE);
851 } else if (v.size()==1) {
852 GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
856 ex last_element=v1.back();
857 GINAC_ASSERT(is_ex_exactly_of_type(last_element,color));
858 GINAC_ASSERT(ex_to_color(last_element).type==color::color_T);
860 ex next_to_last_element=v1.back();
861 GINAC_ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
862 GINAC_ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
866 const ex & last_index=ex_to_color(last_element).seq[0];
867 const ex & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
868 ex summation_index=coloridx();
870 v2.push_back(color_T(summation_index)); // don't care about the representation_label
872 // FIXME: check this formula for SU(N) with N!=3
873 return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
874 % color_trace_of_one_representation_label(v1)
875 +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
876 % color_trace_of_one_representation_label(v2);
878 ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
879 % color_trace_of_one_representation_label(v1);
880 cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
881 ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
882 % color_trace_of_one_representation_label(v2);
883 cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
888 /** Calculate the trace over the (hidden) indices of the su(3) Lie algebra
889 * elements (the SU(3) generators and the unity element) of a specified
890 * representation label in a string of color objects.
892 * @param v Vector of color objects
893 * @param rl Representation label
894 * @return value of the trace */
895 ex color_trace(const exvector & v, unsigned rl)
897 GINAC_ASSERT(rl<MAX_REPRESENTATION_LABELS);
900 v_rest.reserve(v.size()+1); // max size if trace is empty
905 exvectorvector Tvecs;
906 Tvecs.resize(MAX_REPRESENTATION_LABELS);
907 exvectorvector ONEvecs;
908 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
911 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
913 if (unknownvec.size()!=0) {
914 throw(std::invalid_argument("color_trace(): expression must be expanded"));
917 append_exvector_to_exvector(v_rest,delta8vec);
918 append_exvector_to_exvector(v_rest,fvec);
919 append_exvector_to_exvector(v_rest,dvec);
920 for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
922 append_exvector_to_exvector(v_rest,Tvecs[i]);
923 append_exvector_to_exvector(v_rest,ONEvecs[i]);
925 if (Tvecs[i].size()!=0) {
926 v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
927 } else if (ONEvecs[i].size()!=0) {
928 v_rest.push_back(numeric(COLOR_THREE));
930 throw(std::logic_error("color_trace(): representation_label not in color string"));
935 return nonsimplified_ncmul(v_rest);
938 ex simplify_pure_color_string(const ex & e)
940 GINAC_ASSERT(is_ex_exactly_of_type(e,ncmul));
945 exvectorvector Tvecs;
946 Tvecs.resize(MAX_REPRESENTATION_LABELS);
947 exvectorvector ONEvecs;
948 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
951 split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
953 // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
954 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
955 if (Tvecs[rl].size()>=2) {
956 for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
957 for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
958 ex & t1=Tvecs[rl][i];
959 ex & t2=Tvecs[rl][j];
960 GINAC_ASSERT(is_ex_exactly_of_type(t1,color)&&
961 (ex_to_color(t1).type==color::color_T)&&
962 (ex_to_color(t1).seq.size()==1));
963 GINAC_ASSERT(is_ex_exactly_of_type(t2,color)&&
964 (ex_to_color(t2).type==color::color_T)&&
965 (ex_to_color(t2).seq.size()==1));
966 const coloridx & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
967 const coloridx & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
969 if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
971 for (unsigned k=i+1; k<j; ++k) {
972 S.push_back(Tvecs[rl][k]);
976 ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
977 for (unsigned k=i+1; k<j; ++k) {
980 t1=color_trace_of_one_representation_label(S);
981 ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
982 return simplify_color(term1+term2);
989 // FIXME: higher contractions
994 /** Perform some simplifications on an expression containing color objects. */
995 ex simplify_color(const ex & e)
997 // all simplification is done on expanded objects
998 ex e_expanded=e.expand();
1000 // simplification of sum=sum of simplifications
1001 if (is_ex_exactly_of_type(e_expanded,add)) {
1003 for (unsigned i=0; i<e_expanded.nops(); ++i)
1004 sum += simplify_color(e_expanded.op(i));
1009 // simplification of commutative product=commutative product of simplifications
1010 if (is_ex_exactly_of_type(e_expanded,mul)) {
1012 for (unsigned i=0; i<e_expanded.nops(); ++i)
1013 prod *= simplify_color(e_expanded.op(i));
1018 // simplification of noncommutative product: test if everything is color
1019 if (is_ex_exactly_of_type(e_expanded,ncmul)) {
1020 bool all_color=true;
1021 for (unsigned i=0; i<e_expanded.nops(); ++i) {
1022 if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
1028 return simplify_pure_color_string(e_expanded);
1032 // cannot do anything
1036 ex brute_force_sum_color_indices(const ex & e)
1038 exvector iv_all=e.get_indices();
1041 // find double symbolic indices
1042 if (iv_all.size()<2) return e;
1043 for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
1044 GINAC_ASSERT(is_ex_of_type(*cit1,coloridx));
1045 for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
1046 GINAC_ASSERT(is_ex_of_type(*cit2,coloridx));
1047 if (ex_to_coloridx(*cit1).is_symbolic() &&
1048 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
1049 iv_double.push_back(*cit1);
1055 std::vector<int> counter;
1056 counter.resize(iv_double.size());
1058 for (l=0; unsigned(l)<iv_double.size(); ++l) {
1066 for (l=0; unsigned(l)<iv_double.size(); ++l) {
1067 term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
1068 //iv_double[l].print(cout);
1069 //cout << " " << counter[l] << " ";
1074 // increment counter[]
1075 l = iv_double.size()-1;
1076 while ((l>=0)&&((++counter[l])>(int)COLOR_EIGHT)) {
1080 if (l<2) { std::cout << counter[0] << counter[1] << std::endl; }
1087 #ifndef NO_NAMESPACE_GINAC
1088 } // namespace GiNaC
1089 #endif // ndef NO_NAMESPACE_GINAC