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