ginac/lortensor.cpp

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