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