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