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