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