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