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