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