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