ginac/color.cpp

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