- added documentation to the indexed, color, lortensor and clifford classes
[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 // virtual functions which can be overridden by derived classes
653 //////////
654
655 // none
656
657 //////////
658 // non-virtual functions in this class
659 //////////
660
661 //////////
662 // static member variables
663 //////////
664
665 // none
666
667 //////////
668 // global constants
669 //////////
670
671 const color some_color;
672 const std::type_info & typeid_color = typeid(some_color);
673
674 //////////
675 // friend functions
676 //////////
677
678 /** Construct an object representing the unity element of su(3).
679  *
680  *  @param rl Representation label
681  *  @return newly constructed object */
682 color color_ONE(unsigned rl)
683 {
684         return color(color::color_ONE,rl);
685 }
686
687 /** Construct an object representing the generators T_a of SU(3). The index
688  *  must be of class coloridx.
689  *
690  *  @param a Index
691  *  @param rl Representation label
692  *  @return newly constructed object */
693 color color_T(const ex & a, unsigned rl)
694 {
695         return color(color::color_T,a,rl);
696 }
697
698 /** Construct an object representing the antisymmetric structure constants
699  *  f_abc of SU(3). The indices must be of class coloridx.
700  *
701  *  @param a First index
702  *  @param b Second index
703  *  @param c Third index
704  *  @return newly constructed object */
705 color color_f(const ex & a, const ex & b, const ex & c)
706 {
707         return color(color::color_f,a,b,c);
708 }
709
710 /** Construct an object representing the symmetric structure constants d_abc
711  *  of SU(3). The indices must be of class coloridx.
712  *
713  *  @param a First index
714  *  @param b Second index
715  *  @param c Third index
716  *  @return newly constructed object */
717 color color_d(const ex & a, const ex & b, const ex & c)
718 {
719         return color(color::color_d,a,b,c);
720 }
721
722 /** This returns the linear combination d_abc+I*f_abc.
723  *
724  *  @param a First index
725  *  @param b Second index
726  *  @param c Third index
727  *  @return newly constructed object */
728 ex color_h(const ex & a, const ex & b, const ex & c)
729 {
730         return color(color::color_d,a,b,c)+I*color(color::color_f,a,b,c);
731 }
732
733 /** Construct an object representing the unity matrix delta8_ab in su(3).
734  *  The indices must be of class coloridx.
735  *
736  *  @param a First index
737  *  @param b Second index
738  *  @return newly constructed object */
739 color color_delta8(const ex & a, const ex & b)
740 {
741         return color(color::color_delta8,a,b);
742 }
743
744 /** Given a vector of color (and possible other) objects, split it up
745  *  according to the object type (structure constant, generator etc.) and
746  *  representation label while preserving the order within each group. If
747  *  there are non-color objetcs in the vector, the SU(3) generators T_a get
748  *  sorted into the "unknown" group together with the non-color objects
749  *  because we don't know whether these objects commute with the generators.
750  *
751  *  @param v Source vector of expressions
752  *  @param delta8vec Vector of unity matrices (returned)
753  *  @param fvec Vector of antisymmetric structure constants (returned)
754  *  @param dvec Vector of symmetric structure constants (returned)
755  *  @param Tvecs Vectors of generators, one for each representation label (returned)
756  *  @param ONEvecs Vectors of unity elements, one for each representation label (returned)
757  *  @param unknownvec Vector of all non-color objects (returned)
758  *
759  *  @see color::color_types
760  *  @see recombine_color_string */
761 void split_color_string_in_parts(const exvector & v, exvector & delta8vec,
762                                  exvector & fvec, exvector & dvec,
763                                  exvectorvector & Tvecs,
764                                  exvectorvector & ONEvecs,
765                                  exvector & unknownvec)
766 {
767         // if not all elements are of type color, put all Ts in unknownvec to
768         // retain the ordering
769         bool all_color=true;
770         for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
771                 if (!is_ex_exactly_of_type(*cit,color)) {
772                         all_color=false;
773                         break;
774                 }
775         }
776         
777         for (exvector::const_iterator cit=v.begin(); cit!=v.end(); ++cit) {
778                 if (is_ex_exactly_of_type(*cit,color)) {
779                         switch (ex_to_color(*cit).type) {
780                         case color::color_delta8:
781                                 delta8vec.push_back(*cit);
782                                 break;
783                         case color::color_f:
784                                 fvec.push_back(*cit);
785                                 break;
786                         case color::color_d:
787                                 dvec.push_back(*cit);
788                                 break;
789                         case color::color_T:
790                                 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
791                                 if (all_color) {
792                                         Tvecs[ex_to_color(*cit).representation_label].push_back(*cit);
793                                 } else {
794                                         unknownvec.push_back(*cit);
795                                 }
796                                 break;
797                         case color::color_ONE:
798                                 GINAC_ASSERT(ex_to_color(*cit).representation_label<MAX_REPRESENTATION_LABELS);
799                                 ONEvecs[ex_to_color(*cit).representation_label].push_back(*cit);
800                                 break;
801                         default:
802                                 throw(std::logic_error("invalid type in split_color_string_in_parts()"));
803                         }
804                 } else {
805                         unknownvec.push_back(*cit);
806                 }
807         }
808 }    
809
810 /** Merge vectors of color objects sorted by object type into one vector,
811  *  retaining the order within each group. This is the inverse operation of
812  *  split_color_string_in_parts().
813  *
814  *  @param delta8vec Vector of unity matrices
815  *  @param fvec Vector of antisymmetric structure constants
816  *  @param dvec Vector of symmetric structure constants
817  *  @param Tvecs Vectors of generators, one for each representation label
818  *  @param ONEvecs Vectors of unity elements, one for each representation label
819  *  @param unknownvec Vector of all non-color objects
820  *  @return merged vector
821  *
822  *  @see color::color_types
823  *  @see split_color_string_in_parts */
824 exvector recombine_color_string(exvector & delta8vec, exvector & fvec,
825                                 exvector & dvec, exvectorvector & Tvecs,
826                                 exvectorvector & ONEvecs, exvector & unknownvec)
827 {
828         unsigned sz=delta8vec.size()+fvec.size()+dvec.size()+unknownvec.size();
829         for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
830                 sz += Tvecs[rl].size();
831                 sz += ONEvecs[rl].size();
832         }
833         exvector v;
834         v.reserve(sz);
835         
836         append_exvector_to_exvector(v,delta8vec);
837         append_exvector_to_exvector(v,fvec);
838         append_exvector_to_exvector(v,dvec);
839         for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
840                 append_exvector_to_exvector(v,Tvecs[rl]);
841                 append_exvector_to_exvector(v,ONEvecs[rl]);
842         }
843         append_exvector_to_exvector(v,unknownvec);
844         return v;
845 }
846
847 ex color_trace_of_one_representation_label(const exvector & v)
848 {
849         if (v.size()==0) {
850                 return numeric(COLOR_THREE);
851         } else if (v.size()==1) {
852                 GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color));
853                 return _ex0();
854         }
855         exvector v1=v;
856         ex last_element=v1.back();
857         GINAC_ASSERT(is_ex_exactly_of_type(last_element,color));
858         GINAC_ASSERT(ex_to_color(last_element).type==color::color_T);
859         v1.pop_back();
860         ex next_to_last_element=v1.back();
861         GINAC_ASSERT(is_ex_exactly_of_type(next_to_last_element,color));
862         GINAC_ASSERT(ex_to_color(next_to_last_element).type==color::color_T);
863         v1.pop_back();
864         exvector v2=v1;
865
866         const ex & last_index=ex_to_color(last_element).seq[0];
867         const ex & next_to_last_index=ex_to_color(next_to_last_element).seq[0];
868         ex summation_index=coloridx();
869
870         v2.push_back(color_T(summation_index)); // don't care about the representation_label
871         
872         // FIXME: check this formula for SU(N) with N!=3
873         return numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
874                % color_trace_of_one_representation_label(v1)
875                +numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
876                % color_trace_of_one_representation_label(v2);
877         /*
878         ex term1=numeric(1)/numeric(2*COLOR_THREE)*color_delta8(next_to_last_index,last_index)
879                    % color_trace_of_one_representation_label(v1);
880         cout << "term 1 of trace of " << v.size() << " ts=" << term1 << endl;
881         ex term2=numeric(1)/numeric(2)*color_h(next_to_last_index,last_index,summation_index)
882                    % color_trace_of_one_representation_label(v2);
883         cout << "term 2 of trace of " << v.size() << " ts=" << term2 << endl;
884         return term1+term2;
885         */
886 }
887
888 /** Calculate the trace over the (hidden) indices of the su(3) Lie algebra
889  *  elements (the SU(3) generators and the unity element) of a specified
890  *  representation label in a string of color objects.
891  *
892  *  @param v Vector of color objects
893  *  @param rl Representation label
894  *  @return value of the trace */
895 ex color_trace(const exvector & v, unsigned rl)
896 {
897         GINAC_ASSERT(rl<MAX_REPRESENTATION_LABELS);
898         
899         exvector v_rest;
900         v_rest.reserve(v.size()+1); // max size if trace is empty
901         
902         exvector delta8vec;
903         exvector fvec;
904         exvector dvec;
905         exvectorvector Tvecs;
906         Tvecs.resize(MAX_REPRESENTATION_LABELS);
907         exvectorvector ONEvecs;
908         ONEvecs.resize(MAX_REPRESENTATION_LABELS);
909         exvector unknownvec;
910
911         split_color_string_in_parts(v,delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
912
913         if (unknownvec.size()!=0) {
914                 throw(std::invalid_argument("color_trace(): expression must be expanded"));
915         }
916
917         append_exvector_to_exvector(v_rest,delta8vec);
918         append_exvector_to_exvector(v_rest,fvec);
919         append_exvector_to_exvector(v_rest,dvec);
920         for (unsigned i=0; i<MAX_REPRESENTATION_LABELS; ++i) {
921                 if (i!=rl) {
922                         append_exvector_to_exvector(v_rest,Tvecs[i]);
923                         append_exvector_to_exvector(v_rest,ONEvecs[i]);
924                 } else {
925                         if (Tvecs[i].size()!=0) {
926                                 v_rest.push_back(color_trace_of_one_representation_label(Tvecs[i]));
927                         } else if (ONEvecs[i].size()!=0) {
928                                 v_rest.push_back(numeric(COLOR_THREE));
929                         } else {
930                                 throw(std::logic_error("color_trace(): representation_label not in color string"));
931                         }
932                 }
933         }
934
935         return nonsimplified_ncmul(v_rest);
936 }
937
938 ex simplify_pure_color_string(const ex & e)
939 {
940         GINAC_ASSERT(is_ex_exactly_of_type(e,ncmul));
941
942         exvector delta8vec;
943         exvector fvec;
944         exvector dvec;
945         exvectorvector Tvecs;
946         Tvecs.resize(MAX_REPRESENTATION_LABELS);
947         exvectorvector ONEvecs;
948         ONEvecs.resize(MAX_REPRESENTATION_LABELS);
949         exvector unknownvec;
950
951         split_color_string_in_parts(ex_to_ncmul(e).get_factors(),delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec);
952
953         // search for T_k S T_k (=1/2 Tr(S) - 1/6 S)
954         for (unsigned rl=0; rl<MAX_REPRESENTATION_LABELS; ++rl) {
955                 if (Tvecs[rl].size()>=2) {
956                         for (unsigned i=0; i<Tvecs[rl].size()-1; ++i) {
957                                 for (unsigned j=i+1; j<Tvecs[rl].size(); ++j) {
958                                         ex & t1=Tvecs[rl][i];
959                                         ex & t2=Tvecs[rl][j];
960                                         GINAC_ASSERT(is_ex_exactly_of_type(t1,color)&&
961                                                    (ex_to_color(t1).type==color::color_T)&&
962                                                    (ex_to_color(t1).seq.size()==1));
963                                         GINAC_ASSERT(is_ex_exactly_of_type(t2,color)&&
964                                                    (ex_to_color(t2).type==color::color_T)&&
965                                                    (ex_to_color(t2).seq.size()==1));
966                                         const coloridx & idx1=ex_to_coloridx(ex_to_color(t1).seq[0]);
967                                         const coloridx & idx2=ex_to_coloridx(ex_to_color(t2).seq[0]);
968                                         
969                                         if (idx1.is_equal(idx2) && idx1.is_symbolic()) {
970                                                 exvector S;
971                                                 for (unsigned k=i+1; k<j; ++k) {
972                                                         S.push_back(Tvecs[rl][k]);
973                                                 }
974                                                 t1=_ex1();
975                                                 t2=_ex1();
976                                                 ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
977                                                 for (unsigned k=i+1; k<j; ++k) {
978                                                         S.push_back(_ex1());
979                                                 }
980                                                 t1=color_trace_of_one_representation_label(S);
981                                                 ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
982                                                 return simplify_color(term1+term2);
983                                         }
984                                 }
985                         }
986                 }
987         }
988         
989         // FIXME: higher contractions
990         
991         return e;
992 }
993         
994 /** Perform some simplifications on an expression containing color objects. */
995 ex simplify_color(const ex & e)
996 {
997         // all simplification is done on expanded objects
998         ex e_expanded=e.expand();
999
1000         // simplification of sum=sum of simplifications
1001         if (is_ex_exactly_of_type(e_expanded,add)) {
1002                 ex sum=_ex0();
1003                 for (unsigned i=0; i<e_expanded.nops(); ++i)
1004                         sum += simplify_color(e_expanded.op(i));
1005                 
1006                 return sum;
1007         }
1008
1009         // simplification of commutative product=commutative product of simplifications
1010         if (is_ex_exactly_of_type(e_expanded,mul)) {
1011                 ex prod=_ex1();
1012                 for (unsigned i=0; i<e_expanded.nops(); ++i)
1013                         prod *= simplify_color(e_expanded.op(i));
1014                 
1015                 return prod;
1016         }
1017
1018         // simplification of noncommutative product: test if everything is color
1019         if (is_ex_exactly_of_type(e_expanded,ncmul)) {
1020                 bool all_color=true;
1021                 for (unsigned i=0; i<e_expanded.nops(); ++i) {
1022                         if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
1023                                 all_color=false;
1024                                 break;
1025                         }
1026                 }
1027                 if (all_color) {
1028                         return simplify_pure_color_string(e_expanded);
1029                 }
1030         }
1031
1032         // cannot do anything
1033         return e_expanded;
1034 }
1035
1036 ex brute_force_sum_color_indices(const ex & e)
1037 {
1038         exvector iv_all=e.get_indices();
1039         exvector iv_double;
1040         
1041         // find double symbolic indices
1042         if (iv_all.size()<2) return e;
1043         for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
1044                 GINAC_ASSERT(is_ex_of_type(*cit1,coloridx));
1045                 for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
1046                         GINAC_ASSERT(is_ex_of_type(*cit2,coloridx));
1047                         if (ex_to_coloridx(*cit1).is_symbolic() && 
1048                                 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
1049                                 iv_double.push_back(*cit1);
1050                                 break;
1051                         }
1052                 }
1053         }
1054
1055         std::vector<int> counter;
1056         counter.resize(iv_double.size());
1057         int l;
1058         for (l=0; unsigned(l)<iv_double.size(); ++l) {
1059                 counter[l]=1;
1060         }
1061
1062         ex sum;
1063         
1064         while (1) {
1065                 ex term = e;
1066                 for (l=0; unsigned(l)<iv_double.size(); ++l) {
1067                         term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
1068                         //iv_double[l].print(cout);
1069                         //cout << " " << counter[l] << " ";
1070                 }
1071                 //cout << endl;
1072                 sum += term;
1073                 
1074                 // increment counter[]
1075                 l = iv_double.size()-1;
1076                 while ((l>=0)&&((++counter[l])>(int)COLOR_EIGHT)) {
1077                         counter[l]=1;    
1078                         l--;
1079                 }
1080                 if (l<2) { std::cout << counter[0] << counter[1] << std::endl; }
1081                 if (l<0) break;
1082         }
1083         
1084         return sum;
1085 }
1086
1087 #ifndef NO_NAMESPACE_GINAC
1088 } // namespace GiNaC
1089 #endif // ndef NO_NAMESPACE_GINAC