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 void color::copy(const color & other)
61 inherited::copy(other);
63 representation_label=other.representation_label;
66 void color::destroy(bool call_parent)
69 inherited::destroy(call_parent);
79 /** Construct object without any color index. This constructor is for internal
80 * use only. Use the color_ONE() function instead.
82 color::color(color_types const t, unsigned rl) : type(t), representation_label(rl)
84 debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT);
85 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
86 tinfo_key=TINFO_color;
87 GINAC_ASSERT(all_of_type_coloridx());
90 /** Construct object with one color index. This constructor is for internal
91 * use only. Use the color_T() function instead.
93 color::color(color_types const t, const ex & i1, unsigned rl)
94 : inherited(i1), type(t), representation_label(rl)
96 debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
97 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
98 tinfo_key=TINFO_color;
99 GINAC_ASSERT(all_of_type_coloridx());
102 /** Construct object with two color indices. This constructor is for internal
103 * use only. Use the color_delta8() function instead.
104 * @see color_delta8 */
105 color::color(color_types const t, const ex & i1, const ex & i2, unsigned rl)
106 : inherited(i1,i2), type(t), representation_label(rl)
108 debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
109 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
110 tinfo_key=TINFO_color;
111 GINAC_ASSERT(all_of_type_coloridx());
114 /** Construct object with three color indices. This constructor is for internal
115 * use only. Use the color_f(), color_d() and color_h() functions instead.
119 color::color(color_types const t, const ex & i1, const ex & i2, const ex & i3,
120 unsigned rl) : inherited(i1,i2,i3), type(t), representation_label(rl)
122 debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
123 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
124 tinfo_key=TINFO_color;
125 GINAC_ASSERT(all_of_type_coloridx());
128 /** Construct object with arbitrary number of color indices. This
129 * constructor is for internal use only. */
130 color::color(color_types const t, const exvector & iv, unsigned rl)
131 : inherited(iv), type(t), representation_label(rl)
133 debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
134 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
135 tinfo_key=TINFO_color;
136 GINAC_ASSERT(all_of_type_coloridx());
139 color::color(color_types const t, exvector * ivp, unsigned rl)
140 : inherited(ivp), type(t), representation_label(rl)
142 debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
143 GINAC_ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
144 tinfo_key=TINFO_color;
145 GINAC_ASSERT(all_of_type_coloridx());
152 /** Construct object from archive_node. */
153 color::color(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
155 debugmsg("color constructor from archive_node", LOGLEVEL_CONSTRUCT);
157 if (!(n.find_unsigned("type", ty)))
158 throw (std::runtime_error("unknown color type in archive"));
159 type = (color_types)ty;
160 if (!(n.find_unsigned("representation", representation_label)))
161 throw (std::runtime_error("unknown color representation label in archive"));
164 /** Unarchive the object. */
165 ex color::unarchive(const archive_node &n, const lst &sym_lst)
167 return (new color(n, sym_lst))->setflag(status_flags::dynallocated);
170 /** Archive the object. */
171 void color::archive(archive_node &n) const
173 inherited::archive(n);
174 n.add_unsigned("type", type);
175 n.add_unsigned("representation", representation_label);
179 // functions overriding virtual functions from bases classes
184 void color::printraw(std::ostream & os) const
186 debugmsg("color printraw",LOGLEVEL_PRINT);
187 os << "color(type=" << (unsigned)type
188 << ",representation_label=" << representation_label
191 os << ",hash=" << hashvalue << ",flags=" << flags << ")";
194 void color::printtree(std::ostream & os, unsigned indent) const
196 debugmsg("color printtree",LOGLEVEL_PRINT);
197 os << std::string(indent,' ') << "color object: "
198 << "type=" << (unsigned)type
199 << ",representation_label=" << representation_label << ", ";
200 os << seq.size() << " indices" << std::endl;
201 printtreeindices(os,indent);
202 os << std::string(indent,' ') << "hash=" << hashvalue
203 << " (0x" << std::hex << hashvalue << std::dec << ")"
204 << ", flags=" << flags << std::endl;
207 void color::print(std::ostream & os, unsigned upper_precedence) const
209 debugmsg("color print",LOGLEVEL_PRINT);
213 if (representation_label!=0) {
214 os << "^(" << representation_label << ")";
231 os << "INVALID_COLOR_OBJECT";
237 bool color::info(unsigned inf) const
239 return inherited::info(inf);
242 #define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
244 ex color::eval(int level) const
246 // canonicalize indices
248 bool antisymmetric=false;
252 antisymmetric=true; // no break here!
257 int sig=canonicalize_indices(iv,antisymmetric);
259 // something has changed while sorting indices, more evaluations later
260 if (sig==0) return _ex0();
261 return ex(sig)*color(type,iv,representation_label);
266 // nothing to canonicalize
273 GINAC_ASSERT(seq.size()==2);
274 const coloridx & idx1=ex_to_coloridx(seq[0]);
275 const coloridx & idx2=ex_to_coloridx(seq[1]);
277 // check for delta8_{a,a} where a is a symbolic index, replace by 8
278 if ((idx1.is_symbolic())&&(idx1.is_equal_same_type(idx2))) {
279 return ex(COLOR_EIGHT);
282 // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1
283 if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
284 if ((idx1.get_value()!=idx2.get_value())) {
293 // check for d_{a,a,c} (=0) when a is symbolic
295 GINAC_ASSERT(seq.size()==3);
296 const coloridx & idx1=ex_to_coloridx(seq[0]);
297 const coloridx & idx2=ex_to_coloridx(seq[1]);
298 const coloridx & idx3=ex_to_coloridx(seq[2]);
300 if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
302 } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
306 // check for three numeric indices
307 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
308 GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
309 GINAC_ASSERT(idx2.get_value()<=idx3.get_value());
310 if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)||
311 CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) {
313 } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
315 } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) {
316 return 1/sqrt(numeric(3));
317 } else if (CMPINDICES(8,8,8)) {
318 return -1/sqrt(numeric(3));
319 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
320 return -1/(2*sqrt(numeric(3)));
328 GINAC_ASSERT(seq.size()==3);
329 const coloridx & idx1=ex_to_coloridx(seq[0]);
330 const coloridx & idx2=ex_to_coloridx(seq[1]);
331 const coloridx & idx3=ex_to_coloridx(seq[2]);
333 // check for three numeric indices
334 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
335 GINAC_ASSERT(idx1.get_value()<=idx2.get_value());
336 GINAC_ASSERT(idx2.get_value()<=idx3.get_value());
337 if (CMPINDICES(1,2,3)) {
339 } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
340 CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
342 } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
344 } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) {
345 return sqrt(numeric(3))/2;
346 } else if (CMPINDICES(8,8,8)) {
347 return -1/sqrt(numeric(3));
348 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
349 return -1/(2*sqrt(numeric(3)));
356 // nothing to evaluate
365 int color::compare_same_type(const basic & other) const
367 GINAC_ASSERT(other.tinfo() == TINFO_color);
368 const color &o = static_cast<const color &>(other);
370 if (type != o.type) {
372 return type < o.type ? -1 : 1;
375 if (representation_label != o.representation_label) {
376 // different representation label
377 return representation_label < o.representation_label ? -1 : 1;
380 return inherited::compare_same_type(other);
383 bool color::is_equal_same_type(const basic & other) const
385 GINAC_ASSERT(other.tinfo() == TINFO_color);
386 const color &o = static_cast<const color &>(other);
388 if (type != o.type) return false;
389 if (representation_label != o.representation_label) return false;
390 return inherited::is_equal_same_type(other);
395 ex color::simplify_ncmul(const exvector & v) const
397 // simplifications: contract delta8_{a,b} where possible
398 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...
399 // remove superfluous ONEs
401 // contract indices of delta8_{a,b} if they are different and symbolic
403 exvector v_contracted=v;
404 unsigned replacements;
405 bool something_changed=false;
407 exvector::iterator it=v_contracted.begin();
408 while (it!=v_contracted.end()) {
409 // process only delta8 objects
410 if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) {
411 color & d8=ex_to_nonconst_color(*it);
412 GINAC_ASSERT(d8.seq.size()==2);
413 const coloridx & first_idx=ex_to_coloridx(d8.seq[0]);
414 const coloridx & second_idx=ex_to_coloridx(d8.seq[1]);
415 // delta8_{a,a} should have been contracted in color::eval()
416 GINAC_ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic()));
417 ex saved_delta8=*it; // save to restore it later
419 // try to contract first index
421 if (first_idx.is_symbolic()) {
422 replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx);
423 if (replacements==1) {
424 // not contracted except in itself, restore delta8 object
427 // a contracted index should occur exactly twice
428 GINAC_ASSERT(replacements==2);
430 something_changed=true;
434 // try second index only if first was not contracted
435 if ((replacements==1)&&(second_idx.is_symbolic())) {
436 // first index not contracted, *it is guaranteed to be the original delta8 object
437 replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx);
438 if (replacements==1) {
439 // not contracted except in itself, restore delta8 object
442 // a contracted index should occur exactly twice
443 GINAC_ASSERT(replacements==2);
445 something_changed=true;
452 if (something_changed) {
453 // do more simplifications later
454 return nonsimplified_ncmul(v_contracted);
457 // there were no indices to contract
458 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown
459 // (if there is at least one unknown object, all Ts will be unknown to not change the order)
464 exvectorvector Tvecs;
465 Tvecs.resize(MAX_REPRESENTATION_LABELS);
466 exvectorvector ONEvecs;
467 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
470 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
472 // d_{a,k,l} f_{b,k,l}=0 (includes case a=b)
473 if ((dvec.size()>=1)&&(fvec.size()>=1)) {
474 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) {
475 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
476 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
477 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
478 const color & col1=ex_to_color(*it1);
479 const color & col2=ex_to_color(*it2);
480 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
481 if (iv_intersect.size()>=2) return _ex0();
486 // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b)
487 if (dvec.size()>=2) {
488 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) {
489 for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) {
490 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
491 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
492 const color & col1=ex_to_color(*it1);
493 const color & col2=ex_to_color(*it2);
494 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
495 if (iv_intersect.size()>=2) {
496 if (iv_intersect.size()==3) {
497 *it1=numeric(40)/numeric(3);
500 int dummy; // sign unimportant, since symmetric
501 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&dummy);
502 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&dummy);
503 *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
506 return nonsimplified_ncmul(recombine_color_string(
507 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
513 // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b)
514 if (fvec.size()>=2) {
515 for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) {
516 for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) {
517 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color));
518 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color));
519 const color & col1=ex_to_color(*it1);
520 const color & col2=ex_to_color(*it2);
521 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
522 if (iv_intersect.size()>=2) {
523 if (iv_intersect.size()==3) {
528 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,&sig1);
529 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,&sig2);
530 *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
533 return nonsimplified_ncmul(recombine_color_string(
534 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
540 // d_{a,b,c} T_b T_c = 5/6 T_a
541 // f_{a,b,c} T_b T_c = 3/2 I T_a
542 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
543 if ((Tvecs[rl].size()>=2)&&((dvec.size()>=1)||(fvec.size()>=1))) {
544 for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) {
546 GINAC_ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
547 GINAC_ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T);
548 iv.push_back(ex_to_color(*it1).seq[0]);
549 iv.push_back(ex_to_color(*(it1+1)).seq[0]);
551 // d_{a,b,c} T_b T_c = 5/6 T_a
552 for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) {
553 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
554 const color & dref=ex_to_color(*it2);
555 exvector iv_intersect=idx_intersect(dref.seq,iv);
556 if (iv_intersect.size()==2) {
557 int dummy; // sign unimportant, since symmetric
558 ex free_idx=permute_free_index_to_front(dref.seq,iv,&dummy);
559 *it1=color(color_T,free_idx,rl);
560 *(it1+1)=color(color_ONE,rl);
561 *it2=numeric(5)/numeric(6);
562 return nonsimplified_ncmul(recombine_color_string(
563 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
567 // f_{a,b,c} T_b T_c = 3/2 I T_a
568 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
569 GINAC_ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
570 const color & fref=ex_to_color(*it2);
571 exvector iv_intersect=idx_intersect(fref.seq,iv);
572 if (iv_intersect.size()==2) {
574 ex free_idx=permute_free_index_to_front(fref.seq,iv,&sig);
575 *it1=color(color_T,free_idx,rl);
576 *(it1+1)=color(color_ONE,rl);
577 *it2=numeric(sig*3)/numeric(2)*I;
578 return nonsimplified_ncmul(recombine_color_string(
579 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
586 // clear all ONEs when there is at least one corresponding color_T
587 // in this representation, retain one ONE otherwise
588 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
589 if (Tvecs[rl].size()!=0) {
591 } else if (ONEvecs[rl].size()!=0) {
593 ONEvecs[rl].push_back(color(color_ONE,rl));
597 // return a sorted vector
598 return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
599 ONEvecs,unknownvec));
602 ex color::thisexprseq(const exvector & v) const
604 return color(type,v,representation_label);
607 ex color::thisexprseq(exvector * vp) const
609 return color(type,vp,representation_label);
612 /** Check whether all indices are of class coloridx or a subclass. This
613 * function is used internally to make sure that all constructed color
614 * objects really carry color indices and not some other classes. */
615 bool color::all_of_type_coloridx(void) const
617 for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
618 if (!is_ex_of_type(*cit,coloridx)) return false;
627 /** Construct an object representing the unity element of su(3).
629 * @param rl Representation label
630 * @return newly constructed object */
631 color color_ONE(unsigned rl)
633 return color(color::color_ONE,rl);
636 /** Construct an object representing the generators T_a of SU(3). The index
637 * must be of class coloridx.
640 * @param rl Representation label
641 * @return newly constructed object */
642 color color_T(const ex & a, unsigned rl)
644 return color(color::color_T,a,rl);
647 /** Construct an object representing the antisymmetric structure constants
648 * f_abc of SU(3). The indices must be of class coloridx.
650 * @param a First index
651 * @param b Second index
652 * @param c Third index
653 * @return newly constructed object */
654 color color_f(const ex & a, const ex & b, const ex & c)
656 return color(color::color_f,a,b,c);
659 /** Construct an object representing the symmetric structure constants d_abc
660 * of SU(3). The indices must be of class coloridx.
662 * @param a First index
663 * @param b Second index
664 * @param c Third index
665 * @return newly constructed object */
666 color color_d(const ex & a, const ex & b, const ex & c)
668 return color(color::color_d,a,b,c);
671 /** This returns the linear combination d_abc+I*f_abc.
673 * @param a First index
674 * @param b Second index
675 * @param c Third index
676 * @return newly constructed object */
677 ex color_h(const ex & a, const ex & b, const ex & c)
679 return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
682 /** Construct an object representing the unity matrix delta8_ab in su(3).
683 * The indices must be of class coloridx.
685 * @param a First index
686 * @param b Second index
687 * @return newly constructed object */
688 color color_delta8(const ex & a, const ex & b)
690 return color(color::color_delta8,a,b);
693 /** Given a vector of color (and possible other) objects, split it up
694 * according to the object type (structure constant, generator etc.) and
695 * representation label while preserving the order within each group. If
696 * there are non-color objetcs in the vector, the SU(3) generators T_a get
697 * sorted into the "unknown" group together with the non-color objects
698 * because we don't know whether these objects commute with the generators.
700 * @param v Source vector of expressions
701 * @param delta8vec Vector of unity matrices (returned)
702 * @param fvec Vector of antisymmetric structure constants (returned)
703 * @param dvec Vector of symmetric structure constants (returned)
704 * @param Tvecs Vectors of generators, one for each representation label (returned)
705 * @param ONEvecs Vectors of unity elements, one for each representation label (returned)
706 * @param unknownvec Vector of all non-color objects (returned)
708 * @see color::color_types
709 * @see recombine_color_string */
710 void split_color_string_in_parts(const exvector & v, exvector & delta8vec,
711 exvector & fvec, exvector & dvec,
712 exvectorvector & Tvecs,
713 exvectorvector & ONEvecs,
714 exvector & unknownvec)
716 // if not all elements are of type color, put all Ts in unknownvec to
717 // retain the ordering
719 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
720 if (!is_ex_exactly_of_type(*cit,color)) {
726 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
727 if (is_ex_exactly_of_type(*cit,color)) {
728 switch (ex_to_color(*cit).type) {
729 case color::color_delta8:
730 delta8vec.push_back(*cit);
733 fvec.push_back(*cit);
736 dvec.push_back(*cit);
739 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
741 Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
743 unknownvec.push_back(*cit);
746 case color::color_ONE:
747 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
748 ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
751 throw(std::logic_error("invalid type in split_color_string_in_parts()"));
754 unknownvec.push_back(*cit);
759 /** Merge vectors of color objects sorted by object type into one vector,
760 * retaining the order within each group. This is the inverse operation of
761 * split_color_string_in_parts().
763 * @param delta8vec Vector of unity matrices
764 * @param fvec Vector of antisymmetric structure constants
765 * @param dvec Vector of symmetric structure constants
766 * @param Tvecs Vectors of generators, one for each representation label
767 * @param ONEvecs Vectors of unity elements, one for each representation label
768 * @param unknownvec Vector of all non-color objects
769 * @return merged vector
771 * @see color::color_types
772 * @see split_color_string_in_parts */
773 exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
774 exvector & dvec, exvectorvector & Tvecs,
775 exvectorvector & ONEvecs, exvector & unknownvec)
777 unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
778 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
779 sz += Tvecs[rl].size();
780 sz += ONEvecs[rl].size();
785 append_exvector_to_exvector(v,delta8vec);
786 append_exvector_to_exvector(v,fvec);
787 append_exvector_to_exvector(v,dvec);
788 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
789 append_exvector_to_exvector(v,Tvecs[rl]);
790 append_exvector_to_exvector(v,ONEvecs[rl]);
792 append_exvector_to_exvector(v,unknownvec);
796 ex color_trace_of_one_representation_label(const exvector & v)
799 return numeric(COLOR_THREE);
800 } else if (v.size()==1) {
801 GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
805 ex last_element=v1.back();
806 GINAC_ASSERT(is_ex_exactly_of_type(last_element,color));
807 GINAC_ASSERT(ex_to_color(last_element).type==color::color_T);
809 ex next_to_last_element=v1.back();
810 GINAC_ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
811 GINAC_ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
815 const ex & last_index=ex_to_color(last_element).seq[0];
816 const ex & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
817 ex summation_index=coloridx();
819 v2.push_back(color_T(summation_index)); // don't care about the representation_label
821 // FIXME: check this formula for SU(N) with N!=3
822 return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
823 % color_trace_of_one_representation_label(v1)
824 +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
825 % color_trace_of_one_representation_label(v2);
827 ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
828 % color_trace_of_one_representation_label(v1);
829 cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
830 ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
831 % color_trace_of_one_representation_label(v2);
832 cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
837 /** Calculate the trace over the (hidden) indices of the su(3) Lie algebra
838 * elements (the SU(3) generators and the unity element) of a specified
839 * representation label in a string of color objects.
841 * @param v Vector of color objects
842 * @param rl Representation label
843 * @return value of the trace */
844 ex color_trace(const exvector & v, unsigned rl)
846 GINAC_ASSERT(rl<MAX_REPRESENTATION_LABELS);
849 v_rest.reserve(v.size()+1); // max size if trace is empty
854 exvectorvector Tvecs;
855 Tvecs.resize(MAX_REPRESENTATION_LABELS);
856 exvectorvector ONEvecs;
857 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
860 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
862 if (unknownvec.size()!=0) {
863 throw(std::invalid_argument("color_trace(): expression must be expanded"));
866 append_exvector_to_exvector(v_rest,delta8vec);
867 append_exvector_to_exvector(v_rest,fvec);
868 append_exvector_to_exvector(v_rest,dvec);
869 for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
871 append_exvector_to_exvector(v_rest,Tvecs[i]);
872 append_exvector_to_exvector(v_rest,ONEvecs[i]);
874 if (Tvecs[i].size()!=0) {
875 v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
876 } else if (ONEvecs[i].size()!=0) {
877 v_rest.push_back(numeric(COLOR_THREE));
879 throw(std::logic_error("color_trace(): representation_label not in color string"));
884 return nonsimplified_ncmul(v_rest);
887 ex simplify_pure_color_string(const ex & e)
889 GINAC_ASSERT(is_ex_exactly_of_type(e,ncmul));
894 exvectorvector Tvecs;
895 Tvecs.resize(MAX_REPRESENTATION_LABELS);
896 exvectorvector ONEvecs;
897 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
900 split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
902 // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
903 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
904 if (Tvecs[rl].size()>=2) {
905 for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
906 for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
907 ex & t1=Tvecs[rl][i];
908 ex & t2=Tvecs[rl][j];
909 GINAC_ASSERT(is_ex_exactly_of_type(t1,color)&&
910 (ex_to_color(t1).type==color::color_T)&&
911 (ex_to_color(t1).seq.size()==1));
912 GINAC_ASSERT(is_ex_exactly_of_type(t2,color)&&
913 (ex_to_color(t2).type==color::color_T)&&
914 (ex_to_color(t2).seq.size()==1));
915 const coloridx & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
916 const coloridx & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
918 if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
920 for (unsigned k=i+1; k<j; ++k) {
921 S.push_back(Tvecs[rl][k]);
925 ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
926 for (unsigned k=i+1; k<j; ++k) {
929 t1=color_trace_of_one_representation_label(S);
930 ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
931 return simplify_color(term1+term2);
938 // FIXME: higher contractions
943 /** Perform some simplifications on an expression containing color objects. */
944 ex simplify_color(const ex & e)
946 // all simplification is done on expanded objects
947 ex e_expanded=e.expand();
949 // simplification of sum=sum of simplifications
950 if (is_ex_exactly_of_type(e_expanded,add)) {
952 for (unsigned i=0; i<e_expanded.nops(); ++i)
953 sum += simplify_color(e_expanded.op(i));
958 // simplification of commutative product=commutative product of simplifications
959 if (is_ex_exactly_of_type(e_expanded,mul)) {
961 for (unsigned i=0; i<e_expanded.nops(); ++i)
962 prod *= simplify_color(e_expanded.op(i));
967 // simplification of noncommutative product: test if everything is color
968 if (is_ex_exactly_of_type(e_expanded,ncmul)) {
970 for (unsigned i=0; i<e_expanded.nops(); ++i) {
971 if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
977 return simplify_pure_color_string(e_expanded);
981 // cannot do anything
985 ex brute_force_sum_color_indices(const ex & e)
987 exvector iv_all=e.get_indices();
990 // find double symbolic indices
991 if (iv_all.size()<2) return e;
992 for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
993 GINAC_ASSERT(is_ex_of_type(*cit1,coloridx));
994 for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
995 GINAC_ASSERT(is_ex_of_type(*cit2,coloridx));
996 if (ex_to_coloridx(*cit1).is_symbolic() &&
997 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
998 iv_double.push_back(*cit1);
1004 std::vector<int> counter;
1005 counter.resize(iv_double.size());
1007 for (l=0; unsigned(l)<iv_double.size(); ++l) {
1015 for (l=0; unsigned(l)<iv_double.size(); ++l) {
1016 term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
1017 //iv_double[l].print(cout);
1018 //cout << " " << counter[l] << " ";
1023 // increment counter[]
1024 l = iv_double.size()-1;
1025 while ((l>=0)&&((++counter[l])>(int)COLOR_EIGHT)) {
1029 if (l<2) { std::cout << counter[0] << counter[1] << std::endl; }
1036 #ifndef NO_NAMESPACE_GINAC
1037 } // namespace GiNaC
1038 #endif // ndef NO_NAMESPACE_GINAC