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