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