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