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