ginac/lortensor.cpp

   1 /** @file lortensor.cpp
   2  *
   3  *  Implementation of GiNaC's Lorentz tensors. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2001 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 "basic.h"
  31 #include "add.h"
  32 #include "mul.h"
  33 #include "debugmsg.h"
  34 #include "lst.h"
  35 #include "lortensor.h"
  36 #include "operators.h"
  37 #include "tinfos.h"
  38 #include "power.h"
  39 #include "archive.h"
  40 #include "utils.h"
  41 #include "config.h"
  42
  43 #ifndef NO_NAMESPACE_GINAC
  44 namespace GiNaC {
  45 #endif // ndef NO_NAMESPACE_GINAC
  46
  47 GINAC_IMPLEMENT_REGISTERED_CLASS(lortensor, indexed)
  48
  49 //////////
  50 // default constructor, destructor, copy constructor assignment operator and helpers
  51 //////////
  52
  53 // public
  54
  55 lortensor::lortensor() : inherited(TINFO_lortensor), type(invalid)
  56 {
  57         debugmsg("lortensor default constructor",LOGLEVEL_CONSTRUCT);
  58         serial=next_serial++;
  59         name=autoname_prefix()+ToString(serial);
  60 }
  61
  62 //protected
  63
  64 void lortensor::copy(const lortensor & other)
  65 {
  66         inherited::copy(other);
  67         type=other.type;
  68         name=other.name;
  69         serial=other.serial;
  70 }
  71
  72 void lortensor::destroy(bool call_parent)
  73 {
  74         if (call_parent) inherited::destroy(call_parent);
  75 }
  76
  77 //////////
  78 // other constructors
  79 //////////
  80
  81 // protected
  82
  83 /** Construct object without any Lorentz index. This constructor is for
  84  *  internal use only. */
  85 lortensor::lortensor(lortensor_types const lt, const std::string & n) : type(lt), name(n)
  86 {
  87         debugmsg("lortensor constructor from lortensor_types,string",LOGLEVEL_CONSTRUCT);
  88         if (lt == lortensor_symbolic)
  89                 serial = next_serial++;
  90         else
  91                 serial = 0;
  92         tinfo_key = TINFO_lortensor;
  93 }
  94
  95 /** Construct object with one Lorentz index. This constructor is for
  96  *  internal use only. Use the lortensor_vector() or lortensor_symbolic()
  97  *  functions instead.
  98  *  @see lortensor_vector
  99  *  @see lortensor_symbolic */
 100 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu) : inherited(mu), type(lt), name(n)
 101 {
 102         debugmsg("lortensor constructor from lortensor_types,string,ex",LOGLEVEL_CONSTRUCT);
 103         GINAC_ASSERT(all_of_type_lorentzidx());
 104         if (lt == lortensor_symbolic)
 105                 serial = next_serial++;
 106         else
 107                 serial = 0;
 108         tinfo_key=TINFO_lortensor;
 109 }
 110
 111 /** Construct object with two Lorentz indices. This constructor is for
 112  *  internal use only. Use the lortensor_g(), lortensor_delta() or
 113  *  lortensor_symbolic() functions instead.
 114  *  @see lortensor_g
 115  *  @see lortensor_delta
 116  *  @see lortensor_symbolic */
 117 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu) : inherited(mu,nu), type(lt), name(n)
 118 {
 119         debugmsg("lortensor constructor from lortensor_types,string,ex,ex",LOGLEVEL_CONSTRUCT);
 120         GINAC_ASSERT(all_of_type_lorentzidx());
 121         if (lt == lortensor_symbolic)
 122                 serial = next_serial++;
 123         else
 124                 serial = 0;
 125         tinfo_key=TINFO_lortensor;
 126 }
 127
 128 /** Construct object with three Lorentz indices. This constructor is for
 129  *  internal use only. Use the lortensor_symbolic() function instead.
 130  *  @see lortensor_symbolic */
 131 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu, const ex & rho) : inherited(mu,nu,rho), type(lt), name(n)
 132 {
 133         debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex",LOGLEVEL_CONSTRUCT);
 134         GINAC_ASSERT(all_of_type_lorentzidx());
 135         if (lt == lortensor_symbolic)
 136                 serial = next_serial++;
 137         else
 138                 serial = 0;
 139         tinfo_key=TINFO_lortensor;
 140 }
 141
 142 /** Construct object with four Lorentz indices. This constructor is for
 143  *  internal use only. Use the lortensor_epsilon() or lortensor_symbolic()
 144  *  functions instead.
 145  *  @see lortensor_epsilon
 146  *  @see lortensor_symbolic */
 147 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu, const ex & rho, const ex & sigma) : inherited(mu,nu,rho,sigma), type(lt), name(n)
 148 {
 149         debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex,ex",LOGLEVEL_CONSTRUCT);
 150         GINAC_ASSERT(all_of_type_lorentzidx());
 151         if (lt == lortensor_symbolic)
 152                 serial = next_serial++;
 153         else
 154                 serial = 0;
 155         tinfo_key=TINFO_lortensor;
 156 }
 157
 158 /** Construct object with arbitrary number of Lorentz indices. This
 159  *  constructor is for internal use only. Use the lortensor_symbolic()
 160  *  function instead.
 161  *
 162  *  @see lortensor_symbolic */
 163 lortensor::lortensor(lortensor_types const lt, const std::string & n, const exvector & iv) : inherited(iv), type(lt), name(n)
 164 {
 165         debugmsg("lortensor constructor from lortensor_types,string,exvector",LOGLEVEL_CONSTRUCT);
 166         GINAC_ASSERT(all_of_type_lorentzidx());
 167         if (lt == lortensor_symbolic)
 168                 serial = next_serial++;
 169         else
 170                 serial = 0;
 171         tinfo_key=TINFO_lortensor;
 172 }
 173
 174 lortensor::lortensor(lortensor_types const lt, const std::string & n, unsigned s, const exvector & iv) : indexed(iv), type(lt), name(n), serial(s)
 175 {
 176         debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
 177         GINAC_ASSERT(all_of_type_lorentzidx());
 178         tinfo_key=TINFO_lortensor;
 179 }
 180
 181 lortensor::lortensor(lortensor_types const lt, const std::string & n, unsigned s, exvector *ivp) : indexed(ivp), type(lt), name(n), serial(s)
 182 {
 183         debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
 184         GINAC_ASSERT(all_of_type_lorentzidx());
 185         tinfo_key=TINFO_lortensor;
 186 }
 187
 188
 189 //////////
 190 // archiving
 191 //////////
 192
 193 /** Construct object from archive_node. */
 194 lortensor::lortensor(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 195 {
 196         debugmsg("lortensor constructor from archive_node", LOGLEVEL_CONSTRUCT);
 197         unsigned int ty;
 198         if (!(n.find_unsigned("type", ty)))
 199                 throw (std::runtime_error("unknown lortensor type in archive"));
 200         type = (lortensor_types)ty;
 201         if (type == lortensor_symbolic) {
 202                 serial = next_serial++;
 203                 if (!(n.find_string("name", name)))
 204                         name = autoname_prefix() + ToString(serial);
 205         } else
 206                 serial = 0;
 207 }
 208
 209 /** Unarchive the object. */
 210 ex lortensor::unarchive(const archive_node &n, const lst &sym_lst)
 211 {
 212         ex s = (new lortensor(n, sym_lst))->setflag(status_flags::dynallocated);
 213
 214         if (ex_to_lortensor(s).type == lortensor_symbolic) {
 215                 // If lortensor is in sym_lst, return the existing lortensor
 216                 for (unsigned i=0; i<sym_lst.nops(); i++) {
 217                         if (is_ex_of_type(sym_lst.op(i), lortensor) && (ex_to_lortensor(sym_lst.op(i)).name == ex_to_lortensor(s).name))
 218                                 return sym_lst.op(i);
 219                 }
 220         }
 221         return s;
 222 }
 223
 224 /** Archive the object. */
 225 void lortensor::archive(archive_node &n) const
 226 {
 227         inherited::archive(n);
 228         n.add_unsigned("type", type);
 229         if (type == lortensor_symbolic)
 230                 n.add_string("name", name);
 231 }
 232
 233
 234 //////////
 235 // functions overriding virtual functions from bases classes
 236 //////////
 237
 238 //public
 239
 240 void lortensor::printraw(std::ostream & os) const
 241 {
 242         debugmsg("lortensor printraw",LOGLEVEL_PRINT);
 243         os << "lortensor(type=" << (unsigned)type
 244            << ",indices=";
 245         printrawindices(os);
 246         os << ",serial=" << serial;
 247         os << ",hash=" << hashvalue << ",flags=" << flags << ")";
 248 }
 249
 250 void lortensor::printtree(std::ostream & os, unsigned indent) const
 251 {
 252         debugmsg("lortensor printtree",LOGLEVEL_PRINT);
 253         os << std::string(indent,' ') <<"lortensor object: "
 254            << "type=" << (unsigned)type << ","
 255            << seq.size() << " indices" << std::endl;
 256         printtreeindices(os,indent);
 257         os << std::string(indent,' ') << "hash=" << hashvalue
 258            << " (0x" << std::hex << hashvalue << std::dec << ")"
 259            << ", flags=" << flags << std::endl;
 260 }
 261
 262 void lortensor::print(std::ostream & os, unsigned upper_precedence) const
 263 {
 264         debugmsg("lortensor print",LOGLEVEL_PRINT);
 265         switch (type) {
 266         case lortensor_g:
 267                 os << "g";
 268                 break;
 269         case lortensor_delta:
 270                 os << "delta";
 271                 break;
 272         case lortensor_epsilon:
 273                 os << "epsilon";
 274                 break;
 275         case lortensor_symbolic:
 276                 os << name;
 277                 break;
 278         case invalid:
 279         default:
 280                 os << "INVALID_LORTENSOR_OBJECT";
 281                 break;
 282         }
 283         printindices(os);
 284 }
 285
 286 bool lortensor::info(unsigned inf) const
 287 {
 288         return inherited::info(inf);
 289 }
 290
 291 ex lortensor::eval(int level) const
 292 {
 293         if (type==lortensor_g) {
 294                 // canonicalize indices
 295                 exvector iv=seq;
 296                 int sig=canonicalize_indices(iv,false); //symmetric
 297                 if (sig!=INT_MAX) {
 298                         //something has changed while sorting indices, more evaluations later
 299                         return ex(sig) *lortensor(type,name,iv);
 300                 }
 301                 const lorentzidx & idx1=ex_to_lorentzidx(seq[0]);
 302                 const lorentzidx & idx2=ex_to_lorentzidx(seq[1]);
 303                 if ((!idx1.is_symbolic()) && (!idx2.is_symbolic())) {
 304                         //both indices are numeric
 305                         if ((idx1.get_value()==idx2.get_value())) {
 306                                 //both on diagonal
 307                                 if (idx1.get_value()==0){
 308                                         // (0,0)
 309                                         return _ex1();
 310                                 } else {
 311                                         if (idx1.is_covariant() != idx2.is_covariant()) {
 312                                                 // (_i,~i) or (~i,_i), i = 1...3
 313                                                 return _ex1();
 314                                         } else {
 315                                                 // (_i,_i) or (~i,~i), i= 1...3
 316                                                 return _ex_1();
 317                                         }
 318                                 }
 319                         } else {
 320                                 // at least one off-diagonal
 321                                 return _ex0();
 322                         }
 323                 } else if (idx1.is_symbolic() && idx1.is_co_contra_pair(idx2)) {
 324                         if (idx1.is_orthogonal_only())
 325                                 return Dim() - idx1.get_dim_parallel_space();
 326                         else
 327                                 return Dim();
 328                 }
 329         }
 330         return this -> hold();
 331 }
 332
 333 //protected
 334
 335 int lortensor::compare_same_type(const basic & other) const
 336 {
 337         GINAC_ASSERT(is_of_type(other,lortensor));
 338         const lortensor &o = static_cast<const lortensor &>(other);
 339
 340         if (type!=o.type) {
 341                 // different type
 342                 return type < o.type ? -1 : 1;
 343         }
 344
 345         if (type == lortensor_symbolic) {
 346                 // symbolic, compare serials
 347                 if (serial != o.serial) {
 348                         return serial < o.serial ? -1 : 1;
 349                 }
 350         }
 351
 352         return inherited::compare_same_type(other);
 353 }
 354
 355 bool lortensor::is_equal_same_type(const basic & other) const
 356 {
 357         GINAC_ASSERT(is_of_type(other,lortensor));
 358         const lortensor &o = static_cast<const lortensor &>(other);
 359
 360         if (type != o.type) return false;
 361         if (type == lortensor_symbolic && serial != o.serial) return false;
 362         return inherited::is_equal_same_type(other);
 363 }
 364
 365 unsigned lortensor::return_type(void) const
 366 {
 367         return return_types::commutative;
 368 }
 369
 370 unsigned lortensor::return_type_tinfo(void) const
 371 {
 372         return tinfo_key;
 373 }
 374
 375 ex lortensor::thisexprseq(const exvector & v) const
 376 {
 377         return lortensor(type,name,serial,v);
 378 }
 379
 380 ex lortensor::thisexprseq(exvector *vp) const
 381 {
 382         return lortensor(type,name,serial,vp);
 383 }
 384
 385 //////////
 386 // non-virtual functions in this class
 387 //////////
 388
 389 // protected
 390
 391 /** Check whether all indices are of class lorentzidx or a subclass. This
 392  *  function is used internally to make sure that all constructed Lorentz
 393  *  tensors really carry Lorentz indices and not some other classes. */
 394 bool lortensor::all_of_type_lorentzidx(void) const
 395 {
 396         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++ cit) {
 397                 if (!is_ex_of_type(*cit,lorentzidx)) return false;
 398         }
 399         return true;
 400 }
 401
 402 // private
 403
 404 std::string & lortensor::autoname_prefix(void)
 405 {
 406         static std::string * s = new std::string("lortensor");
 407         return *s;
 408 }
 409
 410 //////////
 411 // static member variables
 412 //////////
 413
 414 // private
 415
 416 unsigned lortensor::next_serial=0;
 417
 418 //////////
 419 // friend functions
 420 //////////
 421
 422 /** Construct an object representing the metric tensor g. The indices must
 423  *  be of class lorentzidx.
 424  *
 425  *  @param mu First index
 426  *  @param nu Second index
 427  *  @return newly constructed object */
 428 lortensor lortensor_g(const ex & mu, const ex & nu)
 429 {
 430         return lortensor(lortensor::lortensor_g,"",mu,nu);
 431 }
 432
 433 /** Construct an object representing the unity matrix delta. The indices
 434  *  must be of class lorentzidx.
 435  *
 436  *  @param mu First index
 437  *  @param nu Second index
 438  *  @return newly constructed object */
 439 lortensor lortensor_delta(const ex & mu, const ex & nu)
 440 {
 441         return lortensor(lortensor::lortensor_delta,"",mu,nu);
 442 }
 443
 444 /** Construct an object representing the four-dimensional totally
 445  *  antisymmetric tensor epsilon. The indices must be of class lorentzidx.
 446  *
 447  *  @param mu First index
 448  *  @param nu Second index
 449  *  @param rho Third index
 450  *  @param sigma Fourth index
 451  *  @return newly constructed object */
 452 lortensor lortensor_epsilon(const ex & mu, const ex & nu, const ex & rho, const ex & sigma)
 453 {
 454         return lortensor(lortensor::lortensor_epsilon,"",mu,nu,rho,sigma);
 455 }
 456
 457 /** Construct an object representing a symbolic Lorentz vector. The index
 458  *  must be of class lorentzidx.
 459  *
 460  *  @param n Symbolic name
 461  *  @param mu Index
 462  *  @return newly constructed object */
 463 lortensor lortensor_vector(const std::string & n, const ex & mu)
 464 {
 465         return lortensor(lortensor::lortensor_symbolic,n,mu);
 466 }
 467
 468 /** Construct an object representing a symbolic Lorentz tensor of arbitrary
 469  *  rank. The indices must be of class lorentzidx.
 470  *
 471  *  @param n Symbolic name
 472  *  @param iv Vector of indices
 473  *  @return newly constructed object */
 474 lortensor lortensor_symbolic(const std::string & n, const exvector & iv)
 475 {
 476         return lortensor(lortensor::lortensor_symbolic,n,iv);
 477 }
 478
 479 ex simplify_lortensor_mul(const ex & m)
 480 {
 481         GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
 482         exvector v_contracted;
 483
 484         // collect factors in an exvector, store squares twice
 485         int n=m.nops();
 486         v_contracted.reserve(2*n);
 487         for (int i=0; i<n; ++i) {
 488                 ex f=m.op(i);
 489                 if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
 490                         v_contracted.push_back(f.op(0));
 491                         v_contracted.push_back(f.op(0));
 492                 } else {
 493                         v_contracted.push_back(f);
 494                 }
 495         }
 496
 497         unsigned replacements;
 498         bool something_changed=false;
 499
 500         exvector::iterator it=v_contracted.begin();
 501         while (it!=v_contracted.end()) {
 502                 // process only lor_g objects
 503                 if (is_ex_exactly_of_type(*it,lortensor) &&
 504                         (ex_to_lortensor(*it).type==lortensor::lortensor_g)) {
 505                         const lortensor & g=ex_to_lortensor(*it);
 506                         GINAC_ASSERT(g.seq.size()==2);
 507                         const idx & first_idx=ex_to_lorentzidx(g.seq[0]);
 508                         const idx & second_idx=ex_to_lorentzidx(g.seq[1]);
 509                         // g_{mu,mu} should have been contracted in lortensor::eval()
 510                         GINAC_ASSERT(!first_idx.is_equal(second_idx));
 511                         ex saved_g=*it; // save to restore it later
 512
 513                         // try to contract first index
 514                         replacements=0;
 515                         if (first_idx.is_symbolic()) {
 516                                 replacements = subs_index_in_exvector(v_contracted,
 517                                                                       first_idx.toggle_covariant(),second_idx);
 518                                 if (replacements==0) {
 519                                         // not contracted, restore g object
 520                                         *it=saved_g;
 521                                 } else {
 522                                         // a contracted index should occur exactly once
 523                                         GINAC_ASSERT(replacements==1);
 524                                         *it=_ex1();
 525                                         something_changed=true;
 526                                 }
 527                         }
 528
 529                         // try second index only if first was not contracted
 530                         if ((replacements==0)&&(second_idx.is_symbolic())) {
 531                                 // first index not contracted, *it is again the original g object
 532                                 replacements = subs_index_in_exvector(v_contracted,
 533                                                                       second_idx.toggle_covariant(),first_idx);
 534                                 if (replacements==0) {
 535                                         // not contracted except in itself, restore g object
 536                                         *it=saved_g;
 537                                 } else {
 538                                         // a contracted index should occur exactly once
 539                                         GINAC_ASSERT(replacements==1);
 540                                         *it=_ex1();
 541                                         something_changed=true;
 542                                 }
 543                         }
 544                 }
 545                 ++it;
 546         }
 547         if (something_changed) {
 548                 return mul(v_contracted);
 549         }
 550         return m;
 551 }
 552
 553 /** Perform some simplifications on an expression containing Lorentz tensors. */
 554 ex simplify_lortensor(const ex & e)
 555 {
 556         // all simplification is done on expanded objects
 557         ex e_expanded=e.expand();
 558
 559         // simplification of sum=sum of simplifications
 560         if (is_ex_exactly_of_type(e_expanded,add)) {
 561                 ex sum=_ex0();
 562                 for (unsigned i=0; i<e_expanded.nops(); ++i) {
 563                         sum += simplify_lortensor(e_expanded.op(i));
 564                 }
 565                 return sum;
 566         }
 567
 568         // simplification of (commutative) product
 569         if (is_ex_exactly_of_type(e_expanded,mul)) {
 570                 return simplify_lortensor_mul(e);
 571         }
 572
 573         // cannot do anything
 574         return e_expanded;
 575 }
 576
 577 #ifndef NO_NAMESPACE_GINAC
 578 } // namespace GiNaC
 579 #endif // ndef NO_NAMESPACE_GINAC