ginac/simp_lor.cpp

   1 /** @file simp_lor.cpp
   2  *
   3  *  Implementation of GiNaC's simp_lor objects.
   4  *  No real implementation yet, to be done.     */
   5
   6 #include <string>
   7 #include <list>
   8 #include <algorithm>
   9 #include <iostream>
  10 #include <stdexcept>
  11 #include <map>
  12
  13 #include "ginac.h"
  14
  15 //////////
  16 // default constructor, destructor, copy constructor assignment operator and helpers
  17 //////////
  18
  19 // public
  20
  21 simp_lor::simp_lor() : type(invalid)
  22 {
  23     debugmsg("simp_lor default constructor",LOGLEVEL_CONSTRUCT);
  24     tinfo_key=TINFO_SIMP_LOR;
  25 }
  26
  27 simp_lor::~simp_lor()
  28 {
  29     debugmsg("simp_lor destructor",LOGLEVEL_DESTRUCT);
  30     destroy(0);
  31 }
  32
  33 simp_lor::simp_lor(simp_lor const & other)
  34 {
  35     debugmsg("simp_lor copy constructor",LOGLEVEL_CONSTRUCT);
  36     copy (other);
  37 }
  38
  39 simp_lor const & simp_lor::operator=(simp_lor const & other)
  40 {
  41     debugmsg("simp_lor operator=",LOGLEVEL_ASSIGNMENT);
  42     if (this != &other) {
  43         destroy(1);
  44         copy(other);
  45     }
  46     return *this;
  47 }
  48
  49 // protected
  50
  51 void simp_lor::copy(simp_lor const & other)
  52 {
  53     indexed::copy(other);
  54     type=other.type;
  55     name=other.name;
  56 }
  57
  58 void simp_lor::destroy(bool call_parent)
  59 {
  60     if (call_parent) {
  61         indexed::destroy(call_parent);
  62     }
  63 }
  64
  65 //////////
  66 // other constructors
  67 //////////
  68
  69 // protected
  70
  71 simp_lor::simp_lor(simp_lor_types const t) : type(t)
  72 {
  73     debugmsg("simp_lor constructor from simp_lor_types",LOGLEVEL_CONSTRUCT);
  74     tinfo_key=TINFO_SIMP_LOR;
  75 }
  76
  77 simp_lor::simp_lor(simp_lor_types const t, ex const & i1, ex const & i2) :
  78     indexed(i1,i2), type(t)
  79 {
  80     debugmsg("simp_lor constructor from simp_lor_types,ex,ex",LOGLEVEL_CONSTRUCT);
  81     tinfo_key=TINFO_SIMP_LOR;
  82     ASSERT(all_of_type_lorentzidx());
  83 }
  84
  85 simp_lor::simp_lor(simp_lor_types const t, string const & n, ex const & i1) :
  86     indexed(i1), type(t), name(n)
  87 {
  88     debugmsg("simp_lor constructor from simp_lor_types,string,ex",LOGLEVEL_CONSTRUCT);
  89     tinfo_key=TINFO_SIMP_LOR;
  90     ASSERT(all_of_type_lorentzidx());
  91 }
  92
  93 simp_lor::simp_lor(simp_lor_types const t, string const & n, exvector const & iv) :
  94     indexed(iv), type(t), name(n)
  95 {
  96     debugmsg("simp_lor constructor from simp_lor_types,string,exvector",LOGLEVEL_CONSTRUCT);
  97     tinfo_key=TINFO_SIMP_LOR;
  98     ASSERT(all_of_type_lorentzidx());
  99 }
 100
 101 simp_lor::simp_lor(simp_lor_types const t, string const & n, exvector * ivp) :
 102     indexed(ivp), type(t), name(n)
 103 {
 104     debugmsg("simp_lor constructor from simp_lor_types,string,exvector*",LOGLEVEL_CONSTRUCT);
 105     tinfo_key=TINFO_SIMP_LOR;
 106     ASSERT(all_of_type_lorentzidx());
 107 }
 108
 109 //////////
 110 // functions overriding virtual functions from bases classes
 111 //////////
 112
 113 // public
 114
 115 basic * simp_lor::duplicate() const
 116 {
 117     debugmsg("simp_lor duplicate",LOGLEVEL_DUPLICATE);
 118     return new simp_lor(*this);
 119 }
 120
 121 void simp_lor::printraw(ostream & os) const
 122 {
 123     debugmsg("simp_lor printraw",LOGLEVEL_PRINT);
 124     os << "simp_lor(type=" << (unsigned)type
 125        << ",name=" << name << ",indices=";
 126     printrawindices(os);
 127     os << ",hash=" << hashvalue << ",flags=" << flags << ")";
 128 }
 129
 130 void simp_lor::printtree(ostream & os, unsigned indent) const
 131 {
 132     debugmsg("simp_lor printtree",LOGLEVEL_PRINT);
 133     os << string(indent,' ') << "simp_lor object: "
 134        << "type=" << (unsigned)type
 135        << ", name=" << name << ", ";
 136     os << seq.size() << " indices" << endl;
 137     printtreeindices(os,indent);
 138     os << string(indent,' ') << "hash=" << hashvalue
 139        << " (0x" << hex << hashvalue << dec << ")"
 140        << ", flags=" << flags << endl;
 141 }
 142
 143 void simp_lor::print(ostream & os, unsigned upper_precedence) const
 144 {
 145     debugmsg("simp_lor print",LOGLEVEL_PRINT);
 146     switch (type) {
 147     case simp_lor_g:
 148         os << "g";
 149         break;
 150     case simp_lor_vec:
 151         os << name;
 152         break;
 153     case invalid:
 154     default:
 155         os << "INVALID_SIMP_LOR_OBJECT";
 156         break;
 157     }
 158     printindices(os);
 159 }
 160
 161 void simp_lor::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
 162 {
 163     debugmsg("simp_lor print csrc",LOGLEVEL_PRINT);
 164     print(os,upper_precedence);
 165 }
 166
 167 bool simp_lor::info(unsigned inf) const
 168 {
 169     return indexed::info(inf);
 170 }
 171
 172 ex simp_lor::eval(int level) const
 173 {
 174     if (type==simp_lor_g) {
 175         // canonicalize indices
 176         exvector iv=seq;
 177         int sig=canonicalize_indices(iv,false); // symmetric
 178         if (sig!=INT_MAX) {
 179             // something has changed while sorting indices, more evaluations later
 180             if (sig==0) return exZERO();
 181             return ex(sig)*simp_lor(type,name,iv);
 182         }
 183         lorentzidx const & idx1=ex_to_lorentzidx(seq[0]);
 184         lorentzidx const & idx2=ex_to_lorentzidx(seq[1]);
 185         if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
 186             // both indices are numeric
 187             if ((idx1.get_value()==idx2.get_value())) {
 188                 // both on diagonal
 189                 if (idx1.get_value()==0) {
 190                     // (0,0)
 191                     return exONE();
 192                 } else {
 193                     if (idx1.is_covariant()!=idx2.is_covariant()) {
 194                         // (_i,~i) or (~i,_i), i=1..3
 195                         return exONE();
 196                     } else {
 197                         // (_i,_i) or (~i,~i), i=1..3
 198                         return exMINUSONE();
 199                     }
 200                 }
 201             } else {
 202                 // at least one off-diagonal
 203                 return exZERO();
 204             }
 205         } else if (idx1.is_symbolic() &&
 206                    idx1.is_co_contra_pair(idx2)) {
 207             return Dim()-idx1.get_dim_parallel_space();
 208         }
 209     }
 210
 211     return this->hold();
 212 }
 213
 214 // protected
 215
 216 int simp_lor::compare_same_type(basic const & other) const
 217 {
 218     ASSERT(other.tinfo() == TINFO_SIMP_LOR);
 219     const simp_lor *o = static_cast<const simp_lor *>(&other);
 220     if (type==o->type) {
 221         if (name==o->name) {
 222             return indexed::compare_same_type(other);
 223         }
 224         return name.compare(o->name);
 225     }
 226     return type < o->type ? -1 : 1;
 227 }
 228
 229 bool simp_lor::is_equal_same_type(basic const & other) const
 230 {
 231     ASSERT(other.tinfo() == TINFO_SIMP_LOR);
 232     const simp_lor *o = static_cast<const simp_lor *>(&other);
 233     if (type!=o->type) return false;
 234     if (name!=o->name) return false;
 235     return indexed::is_equal_same_type(other);
 236 }
 237
 238 unsigned simp_lor::return_type(void) const
 239 {
 240     return return_types::commutative;
 241 }
 242
 243 unsigned simp_lor::return_type_tinfo(void) const
 244 {
 245     return tinfo_key;
 246 }
 247
 248 ex simp_lor::thisexprseq(exvector const & v) const
 249 {
 250     return simp_lor(type,name,v);
 251 }
 252
 253 ex simp_lor::thisexprseq(exvector * vp) const
 254 {
 255     return simp_lor(type,name,vp);
 256 }
 257
 258 //////////
 259 // virtual functions which can be overridden by derived classes
 260 //////////
 261
 262 // none
 263
 264 //////////
 265 // non-virtual functions in this class
 266 //////////
 267
 268 // protected
 269
 270 bool simp_lor::all_of_type_lorentzidx(void) const
 271 {
 272     // used only inside of ASSERTs
 273     for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 274         if (!is_ex_of_type(*cit,lorentzidx)) return false;
 275     }
 276     return true;
 277 }
 278
 279 //////////
 280 // static member variables
 281 //////////
 282
 283 // none
 284
 285 //////////
 286 // global constants
 287 //////////
 288
 289 const simp_lor some_simp_lor;
 290 type_info const & typeid_simp_lor=typeid(some_simp_lor);
 291
 292 //////////
 293 // friend functions
 294 //////////
 295
 296 simp_lor lor_g(ex const & mu, ex const & nu)
 297 {
 298     return simp_lor(simp_lor::simp_lor_g,mu,nu);
 299 }
 300
 301 simp_lor lor_vec(string const & n, ex const & mu)
 302 {
 303     return simp_lor(simp_lor::simp_lor_vec,n,mu);
 304 }
 305
 306 ex simplify_simp_lor_mul(ex const & m, scalar_products const & sp)
 307 {
 308     ASSERT(is_ex_exactly_of_type(m,mul));
 309     exvector v_contracted;
 310
 311     // collect factors in an exvector, store squares twice
 312     int n=m.nops();
 313     v_contracted.reserve(2*n);
 314     for (int i=0; i<n; ++i) {
 315         ex f=m.op(i);
 316         if (is_ex_exactly_of_type(f,power)&&f.op(1)==2) {
 317             v_contracted.push_back(f.op(0));
 318             v_contracted.push_back(f.op(0));
 319         } else {
 320             v_contracted.push_back(f);
 321         }
 322     }
 323
 324     unsigned replacements;
 325     bool something_changed=false;
 326
 327     exvector::iterator it=v_contracted.begin();
 328     while (it!=v_contracted.end()) {
 329         // process only lor_g objects
 330         if (is_ex_exactly_of_type(*it,simp_lor) &&
 331             (ex_to_simp_lor(*it).type==simp_lor::simp_lor_g)) {
 332             simp_lor const & g=ex_to_simp_lor(*it);
 333             ASSERT(g.seq.size()==2);
 334             idx const & first_idx=ex_to_lorentzidx(g.seq[0]);
 335             idx const & second_idx=ex_to_lorentzidx(g.seq[1]);
 336             // g_{mu,mu} should have been contracted in simp_lor::eval()
 337             ASSERT(!first_idx.is_equal(second_idx));
 338             ex saved_g=*it; // save to restore it later
 339
 340             // try to contract first index
 341             replacements=0;
 342             if (first_idx.is_symbolic()) {
 343                 replacements = subs_index_in_exvector(v_contracted,
 344                                    first_idx.toggle_covariant(),second_idx);
 345                 if (replacements==0) {
 346                     // not contracted, restore g object
 347                     *it=saved_g;
 348                 } else {
 349                     // a contracted index should occur exactly once
 350                     ASSERT(replacements==1);
 351                     *it=exONE();
 352                     something_changed=true;
 353                 }
 354             }
 355
 356             // try second index only if first was not contracted
 357             if ((replacements==0)&&(second_idx.is_symbolic())) {
 358                 // first index not contracted, *it is again the original g object
 359                 replacements = subs_index_in_exvector(v_contracted,
 360                                    second_idx.toggle_covariant(),first_idx);
 361                 if (replacements==0) {
 362                     // not contracted except in itself, restore g object
 363                     *it=saved_g;
 364                 } else {
 365                     // a contracted index should occur exactly once
 366                     ASSERT(replacements==1);
 367                     *it=exONE();
 368                     something_changed=true;
 369                 }
 370             }
 371         }
 372         ++it;
 373     }
 374
 375     // process only lor_vec objects
 376     bool jump_to_next=false;
 377     exvector::iterator it1=v_contracted.begin();
 378     while (it1!=v_contracted.end()-1) {
 379         if (is_ex_exactly_of_type(*it1,simp_lor) &&
 380             (ex_to_simp_lor(*it1).type==simp_lor::simp_lor_vec)) {
 381             exvector::iterator it2=it1+1;
 382             while ((it2!=v_contracted.end())&&!jump_to_next) {
 383                 if (is_ex_exactly_of_type(*it2,simp_lor) &&
 384                     (ex_to_simp_lor(*it2).type==simp_lor::simp_lor_vec)) {
 385                     simp_lor const & vec1=ex_to_simp_lor(*it1);
 386                     simp_lor const & vec2=ex_to_simp_lor(*it2);
 387                     ASSERT(vec1.seq.size()==1);
 388                     ASSERT(vec2.seq.size()==1);
 389                     lorentzidx const & idx1=ex_to_lorentzidx(vec1.seq[0]);
 390                     lorentzidx const & idx2=ex_to_lorentzidx(vec2.seq[0]);
 391                     if (idx1.is_symbolic() &&
 392                         idx1.is_co_contra_pair(idx2) &&
 393                         sp.is_defined(vec1,vec2)) {
 394                         *it1=sp.evaluate(vec1,vec2);
 395                         *it2=exONE();
 396                         something_changed=true;
 397                         jump_to_next=true;
 398                     }
 399                 }
 400                 ++it2;
 401             }
 402             jump_to_next=false;
 403         }
 404         ++it1;
 405     }
 406     if (something_changed) {
 407         return mul(v_contracted);
 408     }
 409     return m;
 410 }
 411
 412 ex simplify_simp_lor(ex const & e, scalar_products const & sp)
 413 {
 414     // all simplification is done on expanded objects
 415     ex e_expanded=e.expand();
 416
 417     // simplification of sum=sum of simplifications
 418     if (is_ex_exactly_of_type(e_expanded,add)) {
 419         ex sum=exZERO();
 420         for (int i=0; i<e_expanded.nops(); ++i) {
 421             sum += simplify_simp_lor(e_expanded.op(i),sp);
 422         }
 423         return sum;
 424     }
 425
 426     // simplification of commutative product=commutative product of simplifications
 427     if (is_ex_exactly_of_type(e_expanded,mul)) {
 428         return simplify_simp_lor_mul(e,sp);
 429     }
 430
 431     // cannot do anything
 432     return e_expanded;
 433 }
 434
 435 ex Dim(void)
 436 {
 437     static symbol * d=new symbol("dim");
 438     return *d;
 439 }
 440
 441 //////////
 442 // helper classes
 443 //////////
 444
 445 void scalar_products::reg(simp_lor const & v1, simp_lor const & v2,
 446                           ex const & sp)
 447 {
 448     if (v1.compare_same_type(v2)>0) {
 449         reg(v2,v1,sp);
 450         return;
 451     }
 452     spm[make_key(v1,v2)]=sp;
 453 }
 454
 455 bool scalar_products::is_defined(simp_lor const & v1, simp_lor const & v2) const
 456 {
 457     if (v1.compare_same_type(v2)>0) {
 458         return is_defined(v2,v1);
 459     }
 460     return spm.find(make_key(v1,v2))!=spm.end();
 461 }
 462
 463 ex scalar_products::evaluate(simp_lor const & v1, simp_lor const & v2) const
 464 {
 465     if (v1.compare_same_type(v2)>0) {
 466         return evaluate(v2,v1);
 467     }
 468     return spm.find(make_key(v1,v2))->second;
 469 }
 470
 471 void scalar_products::debugprint(void) const
 472 {
 473     cerr << "map size=" << spm.size() << endl;
 474     for (spmap::const_iterator cit=spm.begin(); cit!=spm.end(); ++cit) {
 475         spmapkey const & k=(*cit).first;
 476         cerr << "item key=((" << k.first.first
 477              << "," << k.first.second << "),";
 478         k.second.printraw(cerr);
 479         cerr << ") value=" << (*cit).second << endl;
 480     }
 481 }
 482
 483 spmapkey scalar_products::make_key(simp_lor const & v1, simp_lor const & v2)
 484 {
 485     ASSERT(v1.type==simp_lor::simp_lor_vec);
 486     ASSERT(v2.type==simp_lor::simp_lor_vec);
 487     lorentzidx anon=ex_to_lorentzidx(v1.seq[0]).create_anonymous_representative();
 488     ASSERT(anon.is_equal_same_type(ex_to_lorentzidx(v2.seq[0]).create_anonymous_representative()));
 489     return spmapkey(strstrpair(v1.name,v2.name),anon);
 490 }
 491
 492
 493