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