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