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