ginac/color.cpp

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