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