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