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