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