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