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