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