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(false);
  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(true);
  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(std::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(std::ostream & os, unsigned indent) const
 208 {
 209         debugmsg("color printtree",LOGLEVEL_PRINT);
 210         os << std::string(indent,' ') << "color object: "
 211            << "type=" << (unsigned)type
 212            << ",representation_label=" << representation_label << ", ";
 213         os << seq.size() << " indices" << std::endl;
 214         printtreeindices(os,indent);
 215         os << std::string(indent,' ') << "hash=" << hashvalue
 216            << " (0x" << std::hex << hashvalue << std::dec << ")"
 217            << ", flags=" << flags << std::endl;
 218 }
 219
 220 void color::print(std::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(std::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 std::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(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
 886                                                 for (unsigned k=i+1; k<j; ++k) {
 887                                                         S.push_back(_ex1());
 888                                                 }
 889                                                 t1=color_trace_of_one_representation_label(S);
 890                                                 ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
 891                                                 return simplify_color(term1+term2);
 892                                         }
 893                                 }
 894                         }
 895                 }
 896         }
 897
 898         // FIXME: higher contractions
 899
 900         return e;
 901 }
 902
 903 ex simplify_color(const ex & e)
 904 {
 905         // all simplification is done on expanded objects
 906         ex e_expanded=e.expand();
 907
 908         // simplification of sum=sum of simplifications
 909         if (is_ex_exactly_of_type(e_expanded,add)) {
 910                 ex sum=_ex0();
 911                 for (unsigned i=0; i<e_expanded.nops(); ++i)
 912                         sum += simplify_color(e_expanded.op(i));
 913
 914                 return sum;
 915         }
 916
 917         // simplification of commutative product=commutative product of simplifications
 918         if (is_ex_exactly_of_type(e_expanded,mul)) {
 919                 ex prod=_ex1();
 920                 for (unsigned i=0; i<e_expanded.nops(); ++i)
 921                         prod *= simplify_color(e_expanded.op(i));
 922
 923                 return prod;
 924         }
 925
 926         // simplification of noncommutative product: test if everything is color
 927         if (is_ex_exactly_of_type(e_expanded,ncmul)) {
 928                 bool all_color=true;
 929                 for (unsigned i=0; i<e_expanded.nops(); ++i) {
 930                         if (!is_ex_exactly_of_type(e_expanded.op(i),color)) {
 931                                 all_color=false;
 932                                 break;
 933                         }
 934                 }
 935                 if (all_color) {
 936                         return simplify_pure_color_string(e_expanded);
 937                 }
 938         }
 939
 940         // cannot do anything
 941         return e_expanded;
 942 }
 943
 944 ex brute_force_sum_color_indices(const ex & e)
 945 {
 946         exvector iv_all=e.get_indices();
 947         exvector iv_double;
 948
 949         // find double symbolic indices
 950         if (iv_all.size()<2) return e;
 951         for (exvector::const_iterator cit1=iv_all.begin(); cit1!=iv_all.end()-1; ++cit1) {
 952                 GINAC_ASSERT(is_ex_of_type(*cit1,coloridx));
 953                 for (exvector::const_iterator cit2=cit1+1; cit2!=iv_all.end(); ++cit2) {
 954                         GINAC_ASSERT(is_ex_of_type(*cit2,coloridx));
 955                         if (ex_to_coloridx(*cit1).is_symbolic() &&
 956                                 ex_to_coloridx(*cit1).is_equal(ex_to_coloridx(*cit2))) {
 957                                 iv_double.push_back(*cit1);
 958                                 break;
 959                         }
 960                 }
 961         }
 962
 963         std::vector<int> counter;
 964         counter.resize(iv_double.size());
 965         int l;
 966         for (l=0; unsigned(l)<iv_double.size(); ++l) {
 967                 counter[l]=1;
 968         }
 969
 970         ex sum;
 971
 972         while (1) {
 973                 ex term = e;
 974                 for (l=0; unsigned(l)<iv_double.size(); ++l) {
 975                         term=term.subs(iv_double[l]==coloridx((unsigned)(counter[l])));
 976                         //iv_double[l].print(cout);
 977                         //cout << " " << counter[l] << " ";
 978                 }
 979                 //cout << endl;
 980                 sum += term;
 981
 982                 // increment counter[]
 983                 l = iv_double.size()-1;
 984                 while ((l>=0)&&((++counter[l])>(int)COLOR_EIGHT)) {
 985                         counter[l]=1;
 986                         l--;
 987                 }
 988                 if (l<2) { std::cout << counter[0] << counter[1] << std::endl; }
 989                 if (l<0) break;
 990         }
 991
 992         return sum;
 993 }
 994
 995 void append_exvector_to_exvector(exvector & dest, const exvector & source)
 996 {
 997         for (exvector::const_iterator cit=source.begin(); cit!=source.end(); ++cit) {
 998                 dest.push_back(*cit);
 999         }
1000 }
1001
1002 #ifndef NO_NAMESPACE_GINAC
1003 } // namespace GiNaC
1004 #endif // ndef NO_NAMESPACE_GINAC
1005