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