- switched to automake build environment
[ginac.git] / ginac / color.cpp
1 /** @file color.cpp
2  *
3  *  Implementation of GiNaC's color objects.
4  *  No real implementation yet, to be done.     */
5
6 #include <string>
7 #include <list>
8 #include <algorithm>
9 #include <iostream>
10 #include <stdexcept>
11
12 #include "ginac.h"
13
14 //////////
15 // default constructor, destructor, copy constructor assignment operator and helpers
16 //////////
17
18 // public
19
20 color::color() : type(invalid), representation_label(0)
21 {
22     debugmsg("color default constructor",LOGLEVEL_CONSTRUCT);
23     tinfo_key=TINFO_COLOR;
24 }
25
26 color::~color()
27 {
28     debugmsg("color destructor",LOGLEVEL_DESTRUCT);
29     destroy(0);
30 }
31
32 color::color(color const & other)
33 {
34     debugmsg("color copy constructor",LOGLEVEL_CONSTRUCT);
35     copy (other);
36 }
37
38 color const & color::operator=(color const & other)
39 {
40     debugmsg("color operator=",LOGLEVEL_ASSIGNMENT);
41     if (this != &other) {
42         destroy(1);
43        copy(other);
44     }
45     return *this;
46 }
47
48 // protected
49
50 void color::copy(color const & other)
51 {
52     indexed::copy(other);
53     type=other.type;
54     representation_label=other.representation_label;
55 }
56
57 void color::destroy(bool call_parent)
58 {
59     if (call_parent) {
60         indexed::destroy(call_parent);
61     }
62 }
63
64 //////////
65 // other constructors
66 //////////
67
68 // protected
69
70 color::color(color_types const t, unsigned const rl) : type(t), representation_label(rl)
71 {
72     debugmsg("color constructor from color_types,unsigned",LOGLEVEL_CONSTRUCT);
73     ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
74     tinfo_key=TINFO_COLOR;
75     ASSERT(all_of_type_coloridx());
76 }
77
78 color::color(color_types const t, ex const & i1, unsigned const rl)
79     : indexed(i1), type(t), representation_label(rl)
80 {
81     debugmsg("color constructor from color_types,ex,unsigned",LOGLEVEL_CONSTRUCT);
82     ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
83     tinfo_key=TINFO_COLOR;
84     ASSERT(all_of_type_coloridx());
85 }
86
87 color::color(color_types const t, ex const & i1, ex const & i2, unsigned const rl)
88     : indexed(i1,i2), type(t), representation_label(rl)
89 {
90     debugmsg("color constructor from color_types,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
91     ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
92     tinfo_key=TINFO_COLOR;
93     ASSERT(all_of_type_coloridx());
94 }
95
96 color::color(color_types const t, ex const & i1, ex const & i2, ex const & i3,
97              unsigned const rl) : indexed(i1,i2,i3), type(t), representation_label(rl)
98 {
99     debugmsg("color constructor from color_types,ex,ex,ex,unsigned",LOGLEVEL_CONSTRUCT);
100     ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
101     tinfo_key=TINFO_COLOR;
102     ASSERT(all_of_type_coloridx());
103 }
104
105 color::color(color_types const t, exvector const & iv, unsigned const rl)
106     : indexed(iv), type(t), representation_label(rl)
107 {
108     debugmsg("color constructor from color_types,exvector,unsigned",LOGLEVEL_CONSTRUCT);
109     ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
110     tinfo_key=TINFO_COLOR;
111     ASSERT(all_of_type_coloridx());
112 }
113
114 color::color(color_types const t, exvector * ivp, unsigned const rl)
115     : indexed(ivp), type(t), representation_label(rl)
116 {
117     debugmsg("color constructor from color_types,exvector *,unsigned",LOGLEVEL_CONSTRUCT);
118     ASSERT(representation_label<MAX_REPRESENTATION_LABELS);
119     tinfo_key=TINFO_COLOR;
120     ASSERT(all_of_type_coloridx());
121 }
122
123 //////////
124 // functions overriding virtual functions from bases classes
125 //////////
126
127 // public
128
129 basic * color::duplicate() const
130 {
131     debugmsg("color duplicate",LOGLEVEL_DUPLICATE);
132     return new color(*this);
133 }
134
135 void color::printraw(ostream & os) const
136 {
137     debugmsg("color printraw",LOGLEVEL_PRINT);
138     os << "color(type=" << (unsigned)type
139        << ",representation_label=" << representation_label
140        << ",indices=";
141     printrawindices(os);
142     os << ",hash=" << hashvalue << ",flags=" << flags << ")";
143 }
144
145 void color::printtree(ostream & os, unsigned indent) const
146 {
147     debugmsg("color printtree",LOGLEVEL_PRINT);
148     os << string(indent,' ') << "color object: "
149        << "type=" << (unsigned)type
150        << ",representation_label=" << representation_label << ", ";
151     os << seq.size() << " indices" << endl;
152     printtreeindices(os,indent);
153     os << string(indent,' ') << "hash=" << hashvalue
154        << " (0x" << hex << hashvalue << dec << ")"
155        << ", flags=" << flags << endl;
156 }
157
158 void color::print(ostream & os, unsigned upper_precedence) const
159 {
160     debugmsg("color print",LOGLEVEL_PRINT);
161     switch (type) {
162     case color_T:
163         os << "T";
164         if (representation_label!=0) {
165             os << "^(" << representation_label << ")";
166         }
167         break;
168     case color_f:
169         os << "f";
170         break;
171     case color_d:
172         os << "d";
173         break;
174     case color_delta8:
175         os << "delta8";
176         break;
177     case color_ONE:
178         os << "color_ONE";
179         break;
180     case invalid:
181     default:
182         os << "INVALID_COLOR_OBJECT";
183         break;
184     }
185     printindices(os);
186 }
187
188 void color::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
189 {
190     debugmsg("color print csrc",LOGLEVEL_PRINT);
191     print(os,upper_precedence);
192 }
193
194 bool color::info(unsigned inf) const
195 {
196     return indexed::info(inf);
197 }
198
199 #define CMPINDICES(A,B,C) ((idx1.get_value()==(A))&&(idx2.get_value()==(B))&&(idx3.get_value()==(C)))
200
201 ex color::eval(int level) const
202 {
203     // canonicalize indices
204     
205     bool antisymmetric=false;
206     
207     switch (type) {
208     case color_f:
209         antisymmetric=true; // no break here!
210     case color_d:
211     case color_delta8:
212         {
213             exvector iv=seq;
214             int sig=canonicalize_indices(iv,antisymmetric);
215             if (sig!=INT_MAX) {
216                 // something has changed while sorting indices, more evaluations later
217                 if (sig==0) return exZERO();
218                 return ex(sig)*color(type,iv,representation_label);
219             }
220         }
221         break;
222     default:
223         // nothing to canonicalize
224         break;
225     }
226
227     switch (type) {
228     case color_delta8:
229         {
230             ASSERT(seq.size()==2);
231             coloridx const & idx1=ex_to_coloridx(seq[0]);
232             coloridx const & idx2=ex_to_coloridx(seq[1]);
233             
234             // check for delta8_{a,a} where a is a symbolic index, replace by 8
235             if ((idx1.is_symbolic())&&(idx1.is_equal_same_type(idx2))) {
236                 return ex(COLOR_EIGHT);
237             }
238
239             // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1
240             if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
241                 if ((idx1.get_value()!=idx2.get_value())) {
242                     return exONE();
243                 } else {
244                     return exZERO();
245                 }
246             }
247         }
248         break;
249     case color_d:
250         // check for d_{a,a,c} (=0) when a is symbolic
251         {
252             ASSERT(seq.size()==3);
253             coloridx const & idx1=ex_to_coloridx(seq[0]);
254             coloridx const & idx2=ex_to_coloridx(seq[1]);
255             coloridx const & idx3=ex_to_coloridx(seq[2]);
256             
257             if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) {
258                 return exZERO();
259             } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) {
260                 return exZERO();
261             }
262             
263             // check for three numeric indices
264             if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
265                 ASSERT(idx1.get_value()<=idx2.get_value());
266                 ASSERT(idx2.get_value()<=idx3.get_value());
267                 if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)||
268                     CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) {
269                     return exHALF();
270                 } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) {
271                     return -exHALF();
272                 } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) {
273                     return 1/sqrt(numeric(3));
274                 } else if (CMPINDICES(8,8,8)) {
275                     return -1/sqrt(numeric(3));
276                 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
277                     return -1/(2*sqrt(numeric(3)));
278                 }
279                 return exZERO();
280             }
281         }
282         break;
283     case color_f:
284         {
285             ASSERT(seq.size()==3);
286             coloridx const & idx1=ex_to_coloridx(seq[0]);
287             coloridx const & idx2=ex_to_coloridx(seq[1]);
288             coloridx const & idx3=ex_to_coloridx(seq[2]);
289             
290             // check for three numeric indices
291             if (!(idx1.is_symbolic()||idx2.is_symbolic()||idx3.is_symbolic())) {
292                 ASSERT(idx1.get_value()<=idx2.get_value());
293                 ASSERT(idx2.get_value()<=idx3.get_value());
294                 if (CMPINDICES(1,2,3)) {
295                     return exONE();
296                 } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)||
297                            CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) {
298                     return exHALF();
299                 } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) {
300                     return -exHALF();
301                 } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) {
302                     return sqrt(numeric(3))/2;
303                 } else if (CMPINDICES(8,8,8)) {
304                     return -1/sqrt(numeric(3));
305                 } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) {
306                     return -1/(2*sqrt(numeric(3)));
307                 }
308                 return exZERO();
309             }
310             break;
311         }
312     default:
313         // nothing to evaluate
314         break;
315     }
316     
317     return this->hold();
318 }
319     
320 // protected
321
322 int color::compare_same_type(basic const & other) const
323 {
324     ASSERT(other.tinfo() == TINFO_COLOR);
325     const color *o = static_cast<const color *>(&other);
326     if (type==o->type) {
327         if (representation_label==o->representation_label) {
328             return indexed::compare_same_type(other);
329         }
330         return representation_label < o->representation_label ? -1 : 1;
331     }
332     return type < o->type ? -1 : 1;
333 }
334
335 bool color::is_equal_same_type(basic const & other) const
336 {
337     ASSERT(other.tinfo() == TINFO_COLOR);
338     const color *o = static_cast<const color *>(&other);
339     if (type!=o->type) return false;
340     if (representation_label!=o->representation_label) return false;
341     return indexed::is_equal_same_type(other);
342 }
343
344 #include <iostream>
345
346 ex color::simplify_ncmul(exvector const & v) const
347 {
348     // simplifications: contract delta8_{a,b} where possible
349     //                  sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...
350     //                  remove superfluous ONEs
351     
352     // contract indices of delta8_{a,b} if they are different and symbolic
353
354     exvector v_contracted=v;
355     unsigned replacements;
356     bool something_changed=false;
357
358     exvector::iterator it=v_contracted.begin();
359     while (it!=v_contracted.end()) {
360         // process only delta8 objects
361         if (is_ex_exactly_of_type(*it,color) && (ex_to_color(*it).type==color_delta8)) {
362             color & d8=ex_to_nonconst_color(*it);
363             ASSERT(d8.seq.size()==2);
364             coloridx const & first_idx=ex_to_coloridx(d8.seq[0]);
365             coloridx const & second_idx=ex_to_coloridx(d8.seq[1]);
366             // delta8_{a,a} should have been contracted in color::eval()
367             ASSERT((!first_idx.is_equal(second_idx))||(!first_idx.is_symbolic()));
368             ex saved_delta8=*it; // save to restore it later
369
370             // try to contract first index
371             replacements=1;
372             if (first_idx.is_symbolic()) {
373                 replacements = subs_index_in_exvector(v_contracted,first_idx,second_idx);
374                 if (replacements==1) {
375                     // not contracted except in itself, restore delta8 object
376                     *it=saved_delta8;
377                 } else {
378                     // a contracted index should occur exactly twice
379                     ASSERT(replacements==2);
380                     *it=exONE();
381                     something_changed=true;
382                 }
383             }
384
385             // try second index only if first was not contracted
386             if ((replacements==1)&&(second_idx.is_symbolic())) {
387                 // first index not contracted, *it is guaranteed to be the original delta8 object
388                 replacements = subs_index_in_exvector(v_contracted,second_idx,first_idx);
389                 if (replacements==1) {
390                     // not contracted except in itself, restore delta8 object
391                     *it=saved_delta8;
392                 } else {
393                     // a contracted index should occur exactly twice
394                     ASSERT(replacements==2);
395                     *it=exONE();
396                     something_changed=true;
397                 }
398             }
399         }
400         ++it;
401     }
402
403     if (something_changed) {
404         // do more simplifications later
405         return nonsimplified_ncmul(v_contracted);
406     }
407
408     // there were no indices to contract
409     // sort delta8,f,d,T(rl=0),T(rl=1),...,ONE(rl=0),ONE(rl=1),...,unknown
410     // (if there is at least one unknown object, all Ts will be unknown to not change the order)
411     
412     exvector delta8vec;
413     exvector fvec;
414     exvector dvec;
415     vector<exvector> Tvecs;
416     Tvecs.resize(MAX_REPRESENTATION_LABELS);
417     vector<exvector> ONEvecs;
418     ONEvecs.resize(MAX_REPRESENTATION_LABELS);
419     exvector unknownvec;
420     
421     split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
422
423     // d_{a,k,l} f_{b,k,l}=0 (includes case a=b)
424     if ((dvec.size()>=1)&&(fvec.size()>=1)) {
425         for (exvector::iterator it1=dvec.begin(); it1!=dvec.end(); ++it1) {
426             for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
427                 ASSERT(is_ex_exactly_of_type(*it1,color));
428                 ASSERT(is_ex_exactly_of_type(*it2,color));
429                 color const & col1=ex_to_color(*it1);
430                 color const & col2=ex_to_color(*it2);
431                 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
432                 if (iv_intersect.size()>=2) return exZERO();
433             }
434         }
435     }
436     
437     // d_{a,k,l} d_{b,k,l}=5/3 delta8_{a,b} (includes case a=b)
438     if (dvec.size()>=2) {
439         for (exvector::iterator it1=dvec.begin(); it1!=dvec.end()-1; ++it1) {
440             for (exvector::iterator it2=it1+1; it2!=dvec.end(); ++it2) {
441                 ASSERT(is_ex_exactly_of_type(*it1,color));
442                 ASSERT(is_ex_exactly_of_type(*it2,color));
443                 color const & col1=ex_to_color(*it1);
444                 color const & col2=ex_to_color(*it2);
445                 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
446                 if (iv_intersect.size()>=2) {
447                     if (iv_intersect.size()==3) {
448                         *it1=numeric(40)/numeric(3);
449                         *it2=exONE();
450                     } else {
451                         int sig1, sig2; // unimportant, since symmetric
452                         ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,false,&sig1);
453                         ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,false,&sig2);
454                         *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2);
455                         *it2=exONE();
456                     }
457                     return nonsimplified_ncmul(recombine_color_string(
458                            delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
459                 }
460             }
461         }
462     }
463
464     // f_{a,k,l} f_{b,k,l}=3 delta8_{a,b} (includes case a=b)
465     if (fvec.size()>=2) {
466         for (exvector::iterator it1=fvec.begin(); it1!=fvec.end()-1; ++it1) {
467             for (exvector::iterator it2=it1+1; it2!=fvec.end(); ++it2) {
468                 ASSERT(is_ex_exactly_of_type(*it1,color));
469                 ASSERT(is_ex_exactly_of_type(*it2,color));
470                 color const & col1=ex_to_color(*it1);
471                 color const & col2=ex_to_color(*it2);
472                 exvector iv_intersect=idx_intersect(col1.seq,col2.seq);
473                 if (iv_intersect.size()>=2) {
474                     if (iv_intersect.size()==3) {
475                         *it1=numeric(24);
476                         *it2=exONE();
477                     } else {
478                         int sig1, sig2;
479                         ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,true,&sig1);
480                         ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,true,&sig2);
481                         *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2);
482                         *it2=exONE();
483                     }
484                     return nonsimplified_ncmul(recombine_color_string(
485                            delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
486                 }
487             }
488         }
489     }
490
491     // d_{a,b,c} T_b T_c = 5/6 T_a
492     // f_{a,b,c} T_b T_c = 3/2 I T_a
493     for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
494         if ((Tvecs[rl].size()>=2)&&((dvec.size()>=1)||(fvec.size()>=1))) {
495             for (exvector::iterator it1=Tvecs[rl].begin(); it1!=Tvecs[rl].end()-1; ++it1) {
496                 exvector iv;
497                 ASSERT(is_ex_exactly_of_type(*it1,color)&&ex_to_color(*it1).type==color_T);
498                 ASSERT(is_ex_exactly_of_type(*(it1+1),color)&&ex_to_color(*(it1+1)).type==color_T);
499                 iv.push_back(ex_to_color(*it1).seq[0]);
500                 iv.push_back(ex_to_color(*(it1+1)).seq[0]);
501                 
502                 // d_{a,b,c} T_b T_c = 5/6 T_a
503                 for (exvector::iterator it2=dvec.begin(); it2!=dvec.end(); ++it2) {
504                     ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_d);
505                     color const & dref=ex_to_color(*it2);
506                     exvector iv_intersect=idx_intersect(dref.seq,iv);
507                     if (iv_intersect.size()==2) {
508                         int sig; // unimportant, since symmetric
509                         ex free_idx=permute_free_index_to_front(dref.seq,iv,false,&sig);
510                         *it1=color(color_T,free_idx,rl);
511                         *(it1+1)=color(color_ONE,rl);
512                         *it2=numeric(5)/numeric(6);
513                         return nonsimplified_ncmul(recombine_color_string(
514                                delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
515                     }
516                 }
517
518                 // f_{a,b,c} T_b T_c = 3/2 I T_a
519                 for (exvector::iterator it2=fvec.begin(); it2!=fvec.end(); ++it2) {
520                     ASSERT(is_ex_exactly_of_type(*it2,color)&&ex_to_color(*it2).type==color_f);
521                     color const & fref=ex_to_color(*it2);
522                     exvector iv_intersect=idx_intersect(fref.seq,iv);
523                     if (iv_intersect.size()==2) {
524                         int sig;
525                         ex free_idx=permute_free_index_to_front(fref.seq,iv,true,&sig);
526                         *it1=color(color_T,free_idx,rl);
527                         *(it1+1)=color(color_ONE,rl);
528                         *it2=numeric(sig*3)/numeric(2)*I;
529                         return nonsimplified_ncmul(recombine_color_string(
530                                delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
531                     }
532                 }
533             }
534         }
535     }
536     
537     // clear all ONEs when there is at least one corresponding color_T
538     // in this representation, retain one ONE otherwise
539     for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
540         if (Tvecs[rl].size()!=0) {
541             ONEvecs[rl].clear();
542         } else if (ONEvecs[rl].size()!=0) {
543             ONEvecs[rl].clear();
544             ONEvecs[rl].push_back(color(color_ONE,rl));
545         }
546     }
547
548     // return a sorted vector
549     return simplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,
550                                                    ONEvecs,unknownvec));
551 }
552
553 ex color::thisexprseq(exvector const & v) const
554 {
555     return color(type,v,representation_label);
556 }
557
558 ex color::thisexprseq(exvector * vp) const
559 {
560     return color(type,vp,representation_label);
561 }
562
563 bool color::all_of_type_coloridx(void) const
564 {
565     // used only inside of ASSERTs
566     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
567         if (!is_ex_of_type(*cit,coloridx)) return false;
568     }
569     return true;
570 }
571
572 //////////
573 // virtual functions which can be overridden by derived classes
574 //////////
575
576 // none
577
578 //////////
579 // non-virtual functions in this class
580 //////////
581
582 //////////
583 // static member variables
584 //////////
585
586 // none
587
588 //////////
589 // global constants
590 //////////
591
592 const color some_color;
593 type_info const & typeid_color=typeid(some_color);
594
595 //////////
596 // friend functions
597 //////////
598
599 color color_ONE(unsigned const rl)
600 {
601     return color(color::color_ONE,rl);
602 }
603
604 color color_T(ex const & a, unsigned const rl)
605 {
606     return color(color::color_T,a,rl);
607 }
608
609 color color_f(ex const & a, ex const & b, ex const & c)
610 {
611     return color(color::color_f,a,b,c);
612 }
613
614 color color_d(ex const & a, ex const & b, ex const & c)
615 {
616     return color(color::color_d,a,b,c);
617 }
618
619 ex color_h(ex const & a, ex const & b, ex const & c)
620 {
621     return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
622 }
623
624 color color_delta8(ex const & a, ex const & b)
625 {
626     return color(color::color_delta8,a,b);
627 }
628
629 void split_color_string_in_parts(exvector const & v, exvector & delta8vec,
630                                  exvector & fvec, exvector & dvec,
631                                  vector<exvector> & Tvecs,
632                                  vector<exvector> & ONEvecs,
633                                  exvector & unknownvec)
634 {
635     // if not all elements are of type color, put all Ts in unknownvec to
636     // retain the ordering
637     bool all_color=true;
638     for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
639         if (!is_ex_exactly_of_type(*cit,color)) {
640             all_color=false;
641             break;
642         }
643     }
644     
645     for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
646         if (is_ex_exactly_of_type(*cit,color)) {
647             switch (ex_to_color(*cit).type) {
648             case color::color_delta8:
649                 delta8vec.push_back(*cit);
650                 break;
651             case color::color_f:
652                 fvec.push_back(*cit);
653                 break;
654             case color::color_d:
655                 dvec.push_back(*cit);
656                 break;
657             case color::color_T:
658                 ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
659                 if (all_color) {
660                     Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
661                 } else {
662                     unknownvec.push_back(*cit);
663                 }
664                 break;
665             case color::color_ONE:
666                 ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
667                 ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
668                 break;
669             default:
670                 throw(std::logic_error("invalid type in split_color_string_in_parts()"));
671             }
672         } else {
673             unknownvec.push_back(*cit);
674         }
675     }
676 }    
677
678 exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
679                                 exvector & dvec, vector<exvector> & Tvecs,
680                                 vector<exvector> & ONEvecs, exvector & unknownvec)
681 {
682     unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
683     for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
684         sz += Tvecs[rl].size();
685         sz += ONEvecs[rl].size();
686     }
687     exvector v;
688     v.reserve(sz);
689     
690     append_exvector_to_exvector(v,delta8vec);
691     append_exvector_to_exvector(v,fvec);
692     append_exvector_to_exvector(v,dvec);
693     for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
694         append_exvector_to_exvector(v,Tvecs[rl]);
695         append_exvector_to_exvector(v,ONEvecs[rl]);
696     }
697     append_exvector_to_exvector(v,unknownvec);
698     return v;
699 }
700
701 ex color_trace_of_one_representation_label(exvector const & v)
702 {
703     if (v.size()==0) {
704         return numeric(COLOR_THREE);
705     } else if (v.size()==1) {
706         ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
707         return exZERO();
708     }
709     exvector v1=v;
710     ex last_element=v1.back();
711     ASSERT(is_ex_exactly_of_type(last_element,color));
712     ASSERT(ex_to_color(last_element).type==color::color_T);
713     v1.pop_back();
714     ex next_to_last_element=v1.back();
715     ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
716     ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
717     v1.pop_back();
718     exvector v2=v1;
719
720     ex const & last_index=ex_to_color(last_element).seq[0];
721     ex const & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
722     ex summation_index=coloridx();
723
724     v2.push_back(color_T(summation_index)); // don't care about the representation_label
725     
726     // check this formula for SU(N) with N!=3 !!!!!!!!!
727     return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
728            % color_trace_of_one_representation_label(v1)
729           +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
730            % color_trace_of_one_representation_label(v2);
731     /*
732     ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
733            % color_trace_of_one_representation_label(v1);
734     cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
735     ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
736            % color_trace_of_one_representation_label(v2);
737     cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
738     return term1+term2;
739     */
740 }
741
742 ex color_trace(exvector const & v, unsigned const rl)
743 {
744     ASSERT(rl<MAX_REPRESENTATION_LABELS);
745     
746     exvector v_rest;
747     v_rest.reserve(v.size()+1); // max size if trace is empty
748     
749     exvector delta8vec;
750     exvector fvec;
751     exvector dvec;
752     vector<exvector> Tvecs;
753     Tvecs.resize(MAX_REPRESENTATION_LABELS);
754     vector<exvector> ONEvecs;
755     ONEvecs.resize(MAX_REPRESENTATION_LABELS);
756     exvector unknownvec;
757
758     split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
759
760     if (unknownvec.size()!=0) {
761         throw(std::invalid_argument("color_trace(): expression must be expanded"));
762     }
763
764     append_exvector_to_exvector(v_rest,delta8vec);
765     append_exvector_to_exvector(v_rest,fvec);
766     append_exvector_to_exvector(v_rest,dvec);
767     for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
768         if (i!=rl) {
769             append_exvector_to_exvector(v_rest,Tvecs[i]);
770             append_exvector_to_exvector(v_rest,ONEvecs[i]);
771         } else {
772             if (Tvecs[i].size()!=0) {
773                 v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
774             } else if (ONEvecs[i].size()!=0) {
775                 v_rest.push_back(numeric(COLOR_THREE));
776             } else {
777                 throw(std::logic_error("color_trace(): representation_label not in color string"));
778             }
779         }
780     }
781
782     return nonsimplified_ncmul(v_rest);
783 }
784
785 ex simplify_pure_color_string(ex const & e)
786 {
787     ASSERT(is_ex_exactly_of_type(e,ncmul));
788
789     exvector delta8vec;
790     exvector fvec;
791     exvector dvec;
792     vector<exvector> Tvecs;
793     Tvecs.resize(MAX_REPRESENTATION_LABELS);
794     vector<exvector> ONEvecs;
795     ONEvecs.resize(MAX_REPRESENTATION_LABELS);
796     exvector unknownvec;
797
798     split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
799
800     // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
801     for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
802         if (Tvecs[rl].size()>=2) {
803             for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
804                 for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
805                     ex & t1=Tvecs[rl][i];
806                     ex & t2=Tvecs[rl][j];
807                     ASSERT(is_ex_exactly_of_type(t1,color)&&
808                            (ex_to_color(t1).type==color::color_T)&&
809                            (ex_to_color(t1).seq.size()==1));
810                     ASSERT(is_ex_exactly_of_type(t2,color)&&
811                            (ex_to_color(t2).type==color::color_T)&&
812                            (ex_to_color(t2).seq.size()==1));
813                     coloridx const & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
814                     coloridx const & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
815                     
816                     if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
817                         exvector S;
818                         for (unsigned k=i+1; k<j; ++k) {
819                             S.push_back(Tvecs[rl][k]);
820                         }
821                         t1=exONE();
822                         t2=exONE();
823                         ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(
824                                  delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
825                         for (unsigned k=i+1; k<j; ++k) {
826                             S.push_back(exONE());
827                         }
828                         t1=color_trace_of_one_representation_label(S);
829                         ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(
830                                  delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
831                         return simplify_color(term1+term2);
832                     }
833                 }
834             }
835         }
836     }
837     
838     // TODO: higher contractions!!!!!!!!!!!!!
839     
840     return e;
841 }
842     
843 ex simplify_color(ex const & e)
844 {
845     // all simplification is done on expanded objects
846     ex e_expanded=e.expand();
847
848     // simplification of sum=sum of simplifications
849     if (is_ex_exactly_of_type(e_expanded,add)) {
850         ex sum=exZERO();
851         for (int i=0; i<e_expanded.nops(); ++i) {
852             sum += simplify_color(e_expanded.op(i));
853         }
854         return sum;
855     }
856
857     // simplification of commutative product=commutative product of simplifications
858     if (is_ex_exactly_of_type(e_expanded,mul)) {
859         ex prod=exONE();
860         for (int i=0; i<e_expanded.nops(); ++i) {
861             prod *= simplify_color(e_expanded.op(i));
862         }
863         return prod;
864     }
865
866     // simplification of noncommutative product: test if everything is color
867     if (is_ex_exactly_of_type(e_expanded,ncmul)) {
868         bool all_color=true;
869         for (int i=0; i<e_expanded.nops(); ++i) {
870             if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
871                 all_color=false;
872                 break;
873             }
874         }
875         if (all_color) {
876             return simplify_pure_color_string(e_expanded);
877         }
878     }
879
880     // cannot do anything
881     return e_expanded;
882 }
883
884 ex brute_force_sum_color_indices(ex const & e)
885 {
886     exvector iv_all=e.get_indices();
887     exvector iv_double;
888     
889     // find double symbolic indices
890     if (iv_all.size()<2) return e;
891     for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
892         ASSERT(is_ex_of_type(*cit1,coloridx));
893         for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
894             ASSERT(is_ex_of_type(*cit2,coloridx));
895             if (ex_to_coloridx(*cit1).is_symbolic() && 
896                 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
897                 iv_double.push_back(*cit1);
898                 break;
899             }
900         }
901     }
902
903     vector<int> counter;
904     counter.resize(iv_double.size());
905     int l;
906     for (l=0; unsigned(l)<iv_double.size(); ++l) {
907         counter[l]=1;
908     }
909
910     ex sum=exZERO();
911     
912     while (1) {
913         ex term=e;
914         for (l=0; unsigned(l)<iv_double.size(); ++l) {
915             term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
916             //iv_double[l].print(cout);
917             //cout << " " << counter[l] << " ";
918         }
919         //cout << endl;
920         sum += term;
921         
922         // increment counter[]
923         l=iv_double.size()-1;
924         while ((l>=0)&&((++counter[l])>COLOR_EIGHT)) {
925             counter[l]=1;    
926             l--;
927         }
928         if (l<2) { cout << counter[0] << counter[1] << endl; }
929         if (l<0) break;
930     }
931     
932     return sum;
933 }
934
935 void append_exvector_to_exvector(exvector & dest, exvector const & source)
936 {
937     for (exvector::const_iterator cit=source.begin(); cit!=source.end(); ++cit) {
938         dest.push_back(*cit);
939     }
940 }
941
942
943
944
945
946
947
948
949
950
951
952
953