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