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