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