- more indentation fixes
[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(false);
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(true);
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(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
886                                                 for (unsigned k=i+1; k<j; ++k) {
887                                                         S.push_back(_ex1());
888                                                 }
889                                                 t1=color_trace_of_one_representation_label(S);
890                                                 ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
891                                                 return simplify_color(term1+term2);
892                                         }
893                                 }
894                         }
895                 }
896         }
897         
898         // FIXME: higher contractions
899         
900         return e;
901 }
902         
903 ex simplify_color(const ex & e)
904 {
905         // all simplification is done on expanded objects
906         ex e_expanded=e.expand();
907
908         // simplification of sum=sum of simplifications
909         if (is_ex_exactly_of_type(e_expanded,add)) {
910                 ex sum=_ex0();
911                 for (unsigned i=0; i<e_expanded.nops(); ++i)
912                         sum += simplify_color(e_expanded.op(i));
913                 
914                 return sum;
915         }
916
917         // simplification of commutative product=commutative product of simplifications
918         if (is_ex_exactly_of_type(e_expanded,mul)) {
919                 ex prod=_ex1();
920                 for (unsigned i=0; i<e_expanded.nops(); ++i)
921                         prod *= simplify_color(e_expanded.op(i));
922                 
923                 return prod;
924         }
925
926         // simplification of noncommutative product: test if everything is color
927         if (is_ex_exactly_of_type(e_expanded,ncmul)) {
928                 bool all_color=true;
929                 for (unsigned i=0; i<e_expanded.nops(); ++i) {
930                         if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
931                                 all_color=false;
932                                 break;
933                         }
934                 }
935                 if (all_color) {
936                         return simplify_pure_color_string(e_expanded);
937                 }
938         }
939
940         // cannot do anything
941         return e_expanded;
942 }
943
944 ex brute_force_sum_color_indices(const ex & e)
945 {
946         exvector iv_all=e.get_indices();
947         exvector iv_double;
948         
949         // find double symbolic indices
950         if (iv_all.size()<2) return e;
951         for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
952                 GINAC_ASSERT(is_ex_of_type(*cit1,coloridx));
953                 for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
954                         GINAC_ASSERT(is_ex_of_type(*cit2,coloridx));
955                         if (ex_to_coloridx(*cit1).is_symbolic() && 
956                                 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
957                                 iv_double.push_back(*cit1);
958                                 break;
959                         }
960                 }
961         }
962
963         std::vector<int> counter;
964         counter.resize(iv_double.size());
965         int l;
966         for (l=0; unsigned(l)<iv_double.size(); ++l) {
967                 counter[l]=1;
968         }
969
970         ex sum;
971         
972         while (1) {
973                 ex term = e;
974                 for (l=0; unsigned(l)<iv_double.size(); ++l) {
975                         term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
976                         //iv_double[l].print(cout);
977                         //cout << " " << counter[l] << " ";
978                 }
979                 //cout << endl;
980                 sum += term;
981                 
982                 // increment counter[]
983                 l = iv_double.size()-1;
984                 while ((l>=0)&&((++counter[l])>(int)COLOR_EIGHT)) {
985                         counter[l]=1;    
986                         l--;
987                 }
988                 if (l<2) { std::cout << counter[0] << counter[1] << std::endl; }
989                 if (l<0) break;
990         }
991         
992         return sum;
993 }
994
995 void append_exvector_to_exvector(exvector & dest, const exvector & source)
996 {
997         for (exvector::const_iterator cit=source.begin(); cit!=source.end(); ++cit) {
998                 dest.push_back(*cit);
999         }
1000 }
1001
1002 #ifndef NO_NAMESPACE_GINAC
1003 } // namespace GiNaC
1004 #endif // ndef NO_NAMESPACE_GINAC
1005