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