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