3 * Implementation of GiNaC's color objects.
4 * No real implementation yet, to be done.
6 * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
34 #include "relational.h"
37 // default constructor, destructor, copy constructor assignment operator and helpers
42 color::color() : type(invalid), representation_label(0)
44 debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
45 tinfo_key=TINFO_color;
50 debugmsg("color destructor",LOGLEVEL_DESTRUCT);
54 color::color(color const & other)
56 debugmsg("color copy constructor",LOGLEVEL_CONSTRUCT);
60 color const & color::operator=(color const & other)
62 debugmsg("color operator=",LOGLEVEL_ASSIGNMENT);
72 void color::copy(color const & other)
76 representation_label=other.representation_label;
79 void color::destroy(bool call_parent)
82 indexed::destroy(call_parent);
92 color::color(color_types const t, unsigned const rl) : type(t), representation_label(rl)
94 debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT);
95 ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
96 tinfo_key=TINFO_color;
97 ASSERT(all_of_type_coloridx());
100 color::color(color_types const t, ex const & i1, unsigned const rl)
101 : indexed(i1), type(t), representation_label(rl)
103 debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
104 ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
105 tinfo_key=TINFO_color;
106 ASSERT(all_of_type_coloridx());
109 color::color(color_types const t, ex const & i1, ex const & i2, unsigned const rl)
110 : indexed(i1,i2), type(t), representation_label(rl)
112 debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
113 ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
114 tinfo_key=TINFO_color;
115 ASSERT(all_of_type_coloridx());
118 color::color(color_types const t, ex const & i1, ex const & i2, ex const & i3,
119 unsigned const rl) : indexed(i1,i2,i3), type(t), representation_label(rl)
121 debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
122 ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
123 tinfo_key=TINFO_color;
124 ASSERT(all_of_type_coloridx());
127 color::color(color_types const t, exvector const & iv, unsigned const rl)
128 : indexed(iv), type(t), representation_label(rl)
130 debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
131 ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
132 tinfo_key=TINFO_color;
133 ASSERT(all_of_type_coloridx());
136 color::color(color_types const t, exvector * ivp, unsigned const rl)
137 : indexed(ivp), type(t), representation_label(rl)
139 debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
140 ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
141 tinfo_key=TINFO_color;
142 ASSERT(all_of_type_coloridx());
146 // functions overriding virtual functions from bases classes
151 basic * color::duplicate() const
153 debugmsg("color duplicate",LOGLEVEL_DUPLICATE);
154 return new color(*this);
157 void color::printraw(ostream & os) const
159 debugmsg("color printraw",LOGLEVEL_PRINT);
160 os << "color(type=" << (unsigned)type
161 << ",representation_label=" << representation_label
164 os << ",hash=" << hashvalue << ",flags=" << flags << ")";
167 void color::printtree(ostream & os, unsigned indent) const
169 debugmsg("color printtree",LOGLEVEL_PRINT);
170 os << string(indent,' ') << "color object: "
171 << "type=" << (unsigned)type
172 << ",representation_label=" << representation_label << ", ";
173 os << seq.size() << " indices" << endl;
174 printtreeindices(os,indent);
175 os << string(indent,' ') << "hash=" << hashvalue
176 << " (0x" << hex << hashvalue << dec << ")"
177 << ", flags=" << flags << endl;
180 void color::print(ostream & os, unsigned upper_precedence) const
182 debugmsg("color print",LOGLEVEL_PRINT);
186 if (representation_label!=0) {
187 os << "^(" << representation_label << ")";
204 os << "INVALID_COLOR_OBJECT";
210 void color::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
212 debugmsg("color print csrc",LOGLEVEL_PRINT);
213 print(os,upper_precedence);
216 bool color::info(unsigned inf) const
218 return indexed::info(inf);
221 #define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
223 ex color::eval(int level) const
225 // canonicalize indices
227 bool antisymmetric=false;
231 antisymmetric=true; // no break here!
236 int sig=canonicalize_indices(iv,antisymmetric);
238 // something has changed while sorting indices, more evaluations later
239 if (sig==0) return exZERO();
240 return ex(sig)*color(type,iv,representation_label);
245 // nothing to canonicalize
252 ASSERT(seq.size()==2);
253 coloridx const & idx1=ex_to_coloridx(seq[0]);
254 coloridx const & idx2=ex_to_coloridx(seq[1]);
256 // check for delta8_{a,a} where a is a symbolic index, replace by 8
257 if ((idx1.is_symbolic())&&(idx1.is_equal_same_type(idx2))) {
258 return ex(COLOR_EIGHT);
261 // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1
262 if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
263 if ((idx1.get_value()!=idx2.get_value())) {
272 // check for d_{a,a,c} (=0) when a is symbolic
274 ASSERT(seq.size()==3);
275 coloridx const & idx1=ex_to_coloridx(seq[0]);
276 coloridx const & idx2=ex_to_coloridx(seq[1]);
277 coloridx const & idx3=ex_to_coloridx(seq[2]);
279 if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
281 } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
285 // check for three numeric indices
286 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
287 ASSERT(idx1.get_value()<=idx2.get_value());
288 ASSERT(idx2.get_value()<=idx3.get_value());
289 if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)||
290 CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) {
292 } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
294 } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) {
295 return 1/sqrt(numeric(3));
296 } else if (CMPINDICES(8,8,8)) {
297 return -1/sqrt(numeric(3));
298 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
299 return -1/(2*sqrt(numeric(3)));
307 ASSERT(seq.size()==3);
308 coloridx const & idx1=ex_to_coloridx(seq[0]);
309 coloridx const & idx2=ex_to_coloridx(seq[1]);
310 coloridx const & idx3=ex_to_coloridx(seq[2]);
312 // check for three numeric indices
313 if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
314 ASSERT(idx1.get_value()<=idx2.get_value());
315 ASSERT(idx2.get_value()<=idx3.get_value());
316 if (CMPINDICES(1,2,3)) {
318 } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
319 CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
321 } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
323 } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) {
324 return sqrt(numeric(3))/2;
325 } else if (CMPINDICES(8,8,8)) {
326 return -1/sqrt(numeric(3));
327 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
328 return -1/(2*sqrt(numeric(3)));
335 // nothing to evaluate
344 int color::compare_same_type(basic const & other) const
346 ASSERT(other.tinfo() == TINFO_color);
347 const color *o = static_cast<const color *>(&other);
349 if (representation_label==o->representation_label) {
350 return indexed::compare_same_type(other);
352 return representation_label < o->representation_label ? -1 : 1;
354 return type < o->type ? -1 : 1;
357 bool color::is_equal_same_type(basic const & other) const
359 ASSERT(other.tinfo() == TINFO_color);
360 const color *o = static_cast<const color *>(&other);
361 if (type!=o->type) return false;
362 if (representation_label!=o->representation_label) return false;
363 return indexed::is_equal_same_type(other);
368 ex color::simplify_ncmul(exvector const & v) const
370 // simplifications: contract delta8_{a,b} where possible
371 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...
372 // remove superfluous ONEs
374 // contract indices of delta8_{a,b} if they are different and symbolic
376 exvector v_contracted=v;
377 unsigned replacements;
378 bool something_changed=false;
380 exvector::iterator it=v_contracted.begin();
381 while (it!=v_contracted.end()) {
382 // process only delta8 objects
383 if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) {
384 color & d8=ex_to_nonconst_color(*it);
385 ASSERT(d8.seq.size()==2);
386 coloridx const & first_idx=ex_to_coloridx(d8.seq[0]);
387 coloridx const & second_idx=ex_to_coloridx(d8.seq[1]);
388 // delta8_{a,a} should have been contracted in color::eval()
389 ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic()));
390 ex saved_delta8=*it; // save to restore it later
392 // try to contract first index
394 if (first_idx.is_symbolic()) {
395 replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx);
396 if (replacements==1) {
397 // not contracted except in itself, restore delta8 object
400 // a contracted index should occur exactly twice
401 ASSERT(replacements==2);
403 something_changed=true;
407 // try second index only if first was not contracted
408 if ((replacements==1)&&(second_idx.is_symbolic())) {
409 // first index not contracted, *it is guaranteed to be the original delta8 object
410 replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx);
411 if (replacements==1) {
412 // not contracted except in itself, restore delta8 object
415 // a contracted index should occur exactly twice
416 ASSERT(replacements==2);
418 something_changed=true;
425 if (something_changed) {
426 // do more simplifications later
427 return nonsimplified_ncmul(v_contracted);
430 // there were no indices to contract
431 // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown
432 // (if there is at least one unknown object, all Ts will be unknown to not change the order)
437 vector<exvector> Tvecs;
438 Tvecs.resize(MAX_REPRESENTATION_LABELS);
439 vector<exvector> ONEvecs;
440 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
443 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
445 // d_{a,k,l} f_{b,k,l}=0 (includes case a=b)
446 if ((dvec.size()>=1)&&(fvec.size()>=1)) {
447 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) {
448 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
449 ASSERT(is_ex_exactly_of_type(*it1,color));
450 ASSERT(is_ex_exactly_of_type(*it2,color));
451 color const & col1=ex_to_color(*it1);
452 color const & col2=ex_to_color(*it2);
453 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
454 if (iv_intersect.size()>=2) return exZERO();
459 // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b)
460 if (dvec.size()>=2) {
461 for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) {
462 for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) {
463 ASSERT(is_ex_exactly_of_type(*it1,color));
464 ASSERT(is_ex_exactly_of_type(*it2,color));
465 color const & col1=ex_to_color(*it1);
466 color const & col2=ex_to_color(*it2);
467 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
468 if (iv_intersect.size()>=2) {
469 if (iv_intersect.size()==3) {
470 *it1=numeric(40)/numeric(3);
473 int sig1, sig2; // unimportant, since symmetric
474 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,false,&sig1);
475 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,false,&sig2);
476 *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
479 return nonsimplified_ncmul(recombine_color_string(
480 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
486 // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b)
487 if (fvec.size()>=2) {
488 for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) {
489 for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) {
490 ASSERT(is_ex_exactly_of_type(*it1,color));
491 ASSERT(is_ex_exactly_of_type(*it2,color));
492 color const & col1=ex_to_color(*it1);
493 color const & 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) {
501 ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,true,&sig1);
502 ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,true,&sig2);
503 *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
506 return nonsimplified_ncmul(recombine_color_string(
507 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
513 // d_{a,b,c} T_b T_c = 5/6 T_a
514 // f_{a,b,c} T_b T_c = 3/2 I T_a
515 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
516 if ((Tvecs[rl].size()>=2)&&((dvec.size()>=1)||(fvec.size()>=1))) {
517 for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) {
519 ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
520 ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T);
521 iv.push_back(ex_to_color(*it1).seq[0]);
522 iv.push_back(ex_to_color(*(it1+1)).seq[0]);
524 // d_{a,b,c} T_b T_c = 5/6 T_a
525 for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) {
526 ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
527 color const & dref=ex_to_color(*it2);
528 exvector iv_intersect=idx_intersect(dref.seq,iv);
529 if (iv_intersect.size()==2) {
530 int sig; // unimportant, since symmetric
531 ex free_idx=permute_free_index_to_front(dref.seq,iv,false,&sig);
532 *it1=color(color_T,free_idx,rl);
533 *(it1+1)=color(color_ONE,rl);
534 *it2=numeric(5)/numeric(6);
535 return nonsimplified_ncmul(recombine_color_string(
536 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
540 // f_{a,b,c} T_b T_c = 3/2 I T_a
541 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
542 ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
543 color const & fref=ex_to_color(*it2);
544 exvector iv_intersect=idx_intersect(fref.seq,iv);
545 if (iv_intersect.size()==2) {
547 ex free_idx=permute_free_index_to_front(fref.seq,iv,true,&sig);
548 *it1=color(color_T,free_idx,rl);
549 *(it1+1)=color(color_ONE,rl);
550 *it2=numeric(sig*3)/numeric(2)*I;
551 return nonsimplified_ncmul(recombine_color_string(
552 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
559 // clear all ONEs when there is at least one corresponding color_T
560 // in this representation, retain one ONE otherwise
561 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
562 if (Tvecs[rl].size()!=0) {
564 } else if (ONEvecs[rl].size()!=0) {
566 ONEvecs[rl].push_back(color(color_ONE,rl));
570 // return a sorted vector
571 return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
572 ONEvecs,unknownvec));
575 ex color::thisexprseq(exvector const & v) const
577 return color(type,v,representation_label);
580 ex color::thisexprseq(exvector * vp) const
582 return color(type,vp,representation_label);
585 bool color::all_of_type_coloridx(void) const
587 // used only inside of ASSERTs
588 for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
589 if (!is_ex_of_type(*cit,coloridx)) return false;
595 // virtual functions which can be overridden by derived classes
601 // non-virtual functions in this class
605 // static member variables
614 const color some_color;
615 type_info const & typeid_color=typeid(some_color);
621 color color_ONE(unsigned const rl)
623 return color(color::color_ONE,rl);
626 color color_T(ex const & a, unsigned const rl)
628 return color(color::color_T,a,rl);
631 color color_f(ex const & a, ex const & b, ex const & c)
633 return color(color::color_f,a,b,c);
636 color color_d(ex const & a, ex const & b, ex const & c)
638 return color(color::color_d,a,b,c);
641 ex color_h(ex const & a, ex const & b, ex const & c)
643 return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
646 color color_delta8(ex const & a, ex const & b)
648 return color(color::color_delta8,a,b);
651 void split_color_string_in_parts(exvector const & v, exvector & delta8vec,
652 exvector & fvec, exvector & dvec,
653 vector<exvector> & Tvecs,
654 vector<exvector> & ONEvecs,
655 exvector & unknownvec)
657 // if not all elements are of type color, put all Ts in unknownvec to
658 // retain the ordering
660 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
661 if (!is_ex_exactly_of_type(*cit,color)) {
667 for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
668 if (is_ex_exactly_of_type(*cit,color)) {
669 switch (ex_to_color(*cit).type) {
670 case color::color_delta8:
671 delta8vec.push_back(*cit);
674 fvec.push_back(*cit);
677 dvec.push_back(*cit);
680 ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
682 Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
684 unknownvec.push_back(*cit);
687 case color::color_ONE:
688 ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
689 ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
692 throw(std::logic_error("invalid type in split_color_string_in_parts()"));
695 unknownvec.push_back(*cit);
700 exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
701 exvector & dvec, vector<exvector> & Tvecs,
702 vector<exvector> & ONEvecs, exvector & unknownvec)
704 unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
705 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
706 sz += Tvecs[rl].size();
707 sz += ONEvecs[rl].size();
712 append_exvector_to_exvector(v,delta8vec);
713 append_exvector_to_exvector(v,fvec);
714 append_exvector_to_exvector(v,dvec);
715 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
716 append_exvector_to_exvector(v,Tvecs[rl]);
717 append_exvector_to_exvector(v,ONEvecs[rl]);
719 append_exvector_to_exvector(v,unknownvec);
723 ex color_trace_of_one_representation_label(exvector const & v)
726 return numeric(COLOR_THREE);
727 } else if (v.size()==1) {
728 ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
732 ex last_element=v1.back();
733 ASSERT(is_ex_exactly_of_type(last_element,color));
734 ASSERT(ex_to_color(last_element).type==color::color_T);
736 ex next_to_last_element=v1.back();
737 ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
738 ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
742 ex const & last_index=ex_to_color(last_element).seq[0];
743 ex const & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
744 ex summation_index=coloridx();
746 v2.push_back(color_T(summation_index)); // don't care about the representation_label
748 // check this formula for SU(N) with N!=3 !!!!!!!!!
749 return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
750 % color_trace_of_one_representation_label(v1)
751 +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
752 % color_trace_of_one_representation_label(v2);
754 ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
755 % color_trace_of_one_representation_label(v1);
756 cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
757 ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
758 % color_trace_of_one_representation_label(v2);
759 cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
764 ex color_trace(exvector const & v, unsigned const rl)
766 ASSERT(rl<MAX_REPRESENTATION_LABELS);
769 v_rest.reserve(v.size()+1); // max size if trace is empty
774 vector<exvector> Tvecs;
775 Tvecs.resize(MAX_REPRESENTATION_LABELS);
776 vector<exvector> ONEvecs;
777 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
780 split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
782 if (unknownvec.size()!=0) {
783 throw(std::invalid_argument("color_trace(): expression must be expanded"));
786 append_exvector_to_exvector(v_rest,delta8vec);
787 append_exvector_to_exvector(v_rest,fvec);
788 append_exvector_to_exvector(v_rest,dvec);
789 for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
791 append_exvector_to_exvector(v_rest,Tvecs[i]);
792 append_exvector_to_exvector(v_rest,ONEvecs[i]);
794 if (Tvecs[i].size()!=0) {
795 v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
796 } else if (ONEvecs[i].size()!=0) {
797 v_rest.push_back(numeric(COLOR_THREE));
799 throw(std::logic_error("color_trace(): representation_label not in color string"));
804 return nonsimplified_ncmul(v_rest);
807 ex simplify_pure_color_string(ex const & e)
809 ASSERT(is_ex_exactly_of_type(e,ncmul));
814 vector<exvector> Tvecs;
815 Tvecs.resize(MAX_REPRESENTATION_LABELS);
816 vector<exvector> ONEvecs;
817 ONEvecs.resize(MAX_REPRESENTATION_LABELS);
820 split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
822 // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
823 for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
824 if (Tvecs[rl].size()>=2) {
825 for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
826 for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
827 ex & t1=Tvecs[rl][i];
828 ex & t2=Tvecs[rl][j];
829 ASSERT(is_ex_exactly_of_type(t1,color)&&
830 (ex_to_color(t1).type==color::color_T)&&
831 (ex_to_color(t1).seq.size()==1));
832 ASSERT(is_ex_exactly_of_type(t2,color)&&
833 (ex_to_color(t2).type==color::color_T)&&
834 (ex_to_color(t2).seq.size()==1));
835 coloridx const & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
836 coloridx const & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
838 if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
840 for (unsigned k=i+1; k<j; ++k) {
841 S.push_back(Tvecs[rl][k]);
845 ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(
846 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
847 for (unsigned k=i+1; k<j; ++k) {
848 S.push_back(exONE());
850 t1=color_trace_of_one_representation_label(S);
851 ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(
852 delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
853 return simplify_color(term1+term2);
860 // TODO: higher contractions!!!!!!!!!!!!!
865 ex simplify_color(ex const & e)
867 // all simplification is done on expanded objects
868 ex e_expanded=e.expand();
870 // simplification of sum=sum of simplifications
871 if (is_ex_exactly_of_type(e_expanded,add)) {
873 for (int i=0; i<e_expanded.nops(); ++i) {
874 sum += simplify_color(e_expanded.op(i));
879 // simplification of commutative product=commutative product of simplifications
880 if (is_ex_exactly_of_type(e_expanded,mul)) {
882 for (int i=0; i<e_expanded.nops(); ++i) {
883 prod *= simplify_color(e_expanded.op(i));
888 // simplification of noncommutative product: test if everything is color
889 if (is_ex_exactly_of_type(e_expanded,ncmul)) {
891 for (int i=0; i<e_expanded.nops(); ++i) {
892 if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
898 return simplify_pure_color_string(e_expanded);
902 // cannot do anything
906 ex brute_force_sum_color_indices(ex const & e)
908 exvector iv_all=e.get_indices();
911 // find double symbolic indices
912 if (iv_all.size()<2) return e;
913 for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
914 ASSERT(is_ex_of_type(*cit1,coloridx));
915 for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
916 ASSERT(is_ex_of_type(*cit2,coloridx));
917 if (ex_to_coloridx(*cit1).is_symbolic() &&
918 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
919 iv_double.push_back(*cit1);
926 counter.resize(iv_double.size());
928 for (l=0; unsigned(l)<iv_double.size(); ++l) {
936 for (l=0; unsigned(l)<iv_double.size(); ++l) {
937 term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
938 //iv_double[l].print(cout);
939 //cout << " " << counter[l] << " ";
944 // increment counter[]
945 l=iv_double.size()-1;
946 while ((l>=0)&&((++counter[l])>COLOR_EIGHT)) {
950 if (l<2) { cout << counter[0] << counter[1] << endl; }
957 void append_exvector_to_exvector(exvector & dest, exvector const & source)
959 for (exvector::const_iterator cit=source.begin(); cit!=source.end(); ++cit) {
960 dest.push_back(*cit);