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-2000 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 "symbol.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 //////////
  50 // default constructor, destructor, copy constructor assignment operator and helpers
  51 //////////
  52
  53 // public
  54
  55 lortensor::lortensor()
  56 {
  57         debugmsg("lortensor default constructor",LOGLEVEL_CONSTRUCT);
  58         serial=next_serial++;
  59         name=autoname_prefix()+ToString(serial);
  60         tinfo_key=TINFO_lortensor;
  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         indexed::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) {
  98                 indexed::destroy(call_parent);
  99         }
 100 }
 101
 102 //////////
 103 // other constructors
 104 //////////
 105
 106 // protected
 107
 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         serial=next_serial++;
 112         tinfo_key=TINFO_lortensor;
 113 }
 114
 115 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu) : indexed(mu), type(lt), name(n)
 116 {
 117         debugmsg("lortensor constructor from lortensor_types,string,ex",LOGLEVEL_CONSTRUCT);
 118         serial=next_serial++;
 119         GINAC_ASSERT(all_of_type_lorentzidx());
 120         tinfo_key=TINFO_lortensor;
 121 }
 122
 123 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu) : indexed(mu,nu), type(lt), name(n)
 124 {
 125         debugmsg("lortensor constructor from lortensor_types,string,ex,ex",LOGLEVEL_CONSTRUCT);
 126         serial=next_serial++;
 127         GINAC_ASSERT(all_of_type_lorentzidx());
 128         tinfo_key=TINFO_lortensor;
 129 }
 130
 131 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu, const ex & rho) : indexed(mu,nu,rho), type(lt), name(n)
 132 {
 133         debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex",LOGLEVEL_CONSTRUCT);
 134         serial=next_serial++;
 135         GINAC_ASSERT(all_of_type_lorentzidx());
 136         tinfo_key=TINFO_lortensor;
 137 }
 138
 139 lortensor::lortensor(lortensor_types const lt, const std::string & n, const ex & mu, const ex & nu, const ex & rho, const ex & sigma) : indexed(mu,nu,rho,sigma), type(lt), name(n)
 140 {
 141         debugmsg("lortensor constructor from lortensor_types,string,ex,ex,ex,ex",LOGLEVEL_CONSTRUCT);
 142         serial=next_serial++;
 143         GINAC_ASSERT(all_of_type_lorentzidx());
 144         tinfo_key=TINFO_lortensor;
 145 }
 146
 147 lortensor::lortensor(lortensor_types const lt, const std::string & n, const exvector & iv) : indexed(iv), type(lt), name(n)
 148 {
 149         debugmsg("lortensor constructor from lortensor_types,string,exvector",LOGLEVEL_CONSTRUCT);
 150         serial=next_serial++;
 151         GINAC_ASSERT(all_of_type_lorentzidx());
 152         tinfo_key=TINFO_lortensor;
 153 }
 154
 155 lortensor::lortensor(lortensor_types const lt, const std::string & n, unsigned s, const exvector & iv) : indexed(iv), type(lt), name(n), serial(s)
 156 {
 157         debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
 158         GINAC_ASSERT(all_of_type_lorentzidx());
 159         tinfo_key=TINFO_lortensor;
 160 }
 161
 162 lortensor::lortensor(lortensor_types const lt, const std::string & n, unsigned s, exvector *ivp) : indexed(ivp), type(lt), name(n), serial(s)
 163 {
 164         debugmsg("lortensor constructor from lortensor_types,string,unsigned,exvector",LOGLEVEL_CONSTRUCT);
 165         GINAC_ASSERT(all_of_type_lorentzidx());
 166         tinfo_key=TINFO_lortensor;
 167 }
 168
 169 //////////
 170 // functions overriding virtual functions from bases classes
 171 //////////
 172
 173 //public
 174
 175 basic * lortensor::duplicate() const
 176 {
 177         debugmsg("lortensor duplicate",LOGLEVEL_DUPLICATE);
 178         return new lortensor(*this);
 179 }
 180
 181 void lortensor::printraw(std::ostream & os) const
 182 {
 183         debugmsg("lortensor printraw",LOGLEVEL_PRINT);
 184         os << "lortensor(type=" << (unsigned)type
 185            << ",indices=";
 186         printrawindices(os);
 187         os << ",hash=" << hashvalue << ",flags=" << flags << ")";
 188 }
 189
 190 void lortensor::printtree(std::ostream & os, unsigned indent) const
 191 {
 192         debugmsg("lortensor printtree",LOGLEVEL_PRINT);
 193         os << std::string(indent,' ') <<"lortensor object: "
 194            << "type=" << (unsigned)type << ","
 195            << seq.size() << " indices" << std::endl;
 196         printtreeindices(os,indent);
 197         os << std::string(indent,' ') << "hash=" << hashvalue
 198            << " (0x" << std::hex << hashvalue << std::dec << ")"
 199            << ", flags=" << flags << std::endl;
 200 }
 201
 202 void lortensor::print(std::ostream & os, unsigned upper_precedence) const
 203 {
 204         debugmsg("lortensor print",LOGLEVEL_PRINT);
 205         switch (type) {
 206         case lortensor_g:
 207                 os << "g";
 208                 break;
 209         case lortensor_rankn:
 210                 os << name;
 211                 break;
 212         case lortensor_rank1:
 213                 os << name;
 214                 break;
 215         case lortensor_rank2:
 216                 os << name;
 217                 break;
 218         case lortensor_epsilon:
 219                 os << "epsilon";
 220                 break;
 221         case invalid:
 222         default:
 223                 os << "INVALID_LORTENSOR_OBJECT";
 224                 break;
 225         }
 226         printindices(os);
 227 }
 228
 229 void lortensor::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
 230 {
 231         debugmsg("lortensor print csrc",LOGLEVEL_PRINT);
 232         print(os,upper_precedence);
 233 }
 234
 235 bool lortensor::info(unsigned inf) const
 236 {
 237         return indexed::info(inf);
 238 }
 239
 240 ex lortensor::eval(int level) const
 241 {
 242         if (type==lortensor_g) {
 243                 // canonicalize indices
 244                 exvector iv=seq;
 245                 int sig=canonicalize_indices(iv,false); //symmetric
 246                 if (sig!=INT_MAX) {
 247                         //something has changed while sorting indices, more evaluations later
 248                         return ex(sig) *lortensor(type,name,iv);
 249                 }
 250                 const lorentzidx & idx1=ex_to_lorentzidx(seq[0]);
 251                 const lorentzidx & idx2=ex_to_lorentzidx(seq[1]);
 252                 if ((!idx1.is_symbolic()) && (!idx2.is_symbolic())) {
 253                         //both indices are numeric
 254                         if ((idx1.get_value()==idx2.get_value())) {
 255                                 //both on diagonal
 256                                 if (idx1.get_value()==0){
 257                                         // (0,0)
 258                                         return _ex1();
 259                                 } else {
 260                                         if (idx1.is_covariant() != idx2.is_covariant()) {
 261                                                 // (_i,~i) or (~i,_i), i = 1...3
 262                                                 return _ex1();
 263                                         } else {
 264                                                 // (_i,_i) or (~i,~i), i= 1...3
 265                                                 return _ex_1();
 266                                         }
 267                                 }
 268                         } else {
 269                                 // at least one off-diagonal
 270                                 return _ex0();
 271                         }
 272                 } else if (idx1.is_symbolic() && idx1.is_co_contra_pair(idx2)) {
 273                         return Dim()-idx1.get_dim_parallel_space();
 274                 }
 275         }
 276         return this -> hold();
 277 }
 278
 279 //protected
 280
 281 int lortensor::compare_same_type(const basic & other) const
 282 {
 283         GINAC_ASSERT(is_of_type(other,lortensor));
 284         const lortensor *o = static_cast <const lortensor *> (&other);
 285         if (type==o->type) {
 286                 if (type==lortensor_rankn) {
 287                         if (serial!=o->serial) {
 288                                 return serial < o->serial ? -1 : 1;
 289                         }
 290                 }
 291                 return indexed::compare_same_type(other);
 292         }
 293         return type < o->type ? -1 : 1;
 294 }
 295
 296 bool lortensor::is_equal_same_type(const basic & other) const
 297 {
 298         GINAC_ASSERT(is_of_type(other,lortensor));
 299         const lortensor *o=static_cast<const lortensor *> (&other);
 300         if (type!=o->type) return false;
 301         if (type==lortensor_rankn) {
 302                 if (serial!=o->serial) return false;
 303         }
 304         return indexed::is_equal_same_type(other);
 305 }
 306
 307 unsigned lortensor::return_type(void) const
 308 {
 309         return return_types::commutative;
 310 }
 311
 312 unsigned lortensor::return_type_tinfo(void) const
 313 {
 314         return tinfo_key;
 315 }
 316
 317 ex lortensor::thisexprseq(const exvector & v) const
 318 {
 319         return lortensor(type,name,serial,v);
 320 }
 321
 322 ex lortensor::thisexprseq(exvector *vp) const
 323 {
 324         return lortensor(type,name,serial,vp);
 325 }
 326
 327 //////////
 328 // non-virtual functions in this class
 329 //////////
 330
 331 // protected
 332
 333 void lortensor::setname(const std::string & n)
 334 {
 335         name = n;
 336 }
 337
 338 bool lortensor::all_of_type_lorentzidx(void) const
 339 {
 340         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++ cit) {
 341                 if (!is_ex_of_type(*cit,lorentzidx)) {
 342                         return false;
 343                 }
 344         }
 345         return true;
 346 }
 347
 348 // private
 349
 350 std::string & lortensor::autoname_prefix(void)
 351 {
 352         static std::string * s = new std::string("lortensor");
 353         return *s;
 354 }
 355
 356 //////////
 357 // static member variables
 358 //////////
 359
 360 // private
 361
 362 unsigned lortensor::next_serial=0;
 363
 364 //////////
 365 // friend functions
 366 //////////
 367
 368 lortensor lortensor_g(const ex & mu, const ex & nu)
 369 {
 370         return lortensor(lortensor::lortensor_g,"",mu,nu);
 371 }
 372
 373 lortensor lortensor_epsilon(const ex & mu, const ex & nu, const ex & rho, const ex & sigma)
 374 {
 375         return lortensor(lortensor::lortensor_epsilon,"",mu,nu,rho,sigma);
 376 }
 377
 378 lortensor lortensor_rank1(const std::string & n, const ex & mu)
 379 {
 380         return lortensor(lortensor::lortensor_rank1,n,mu);
 381 }
 382
 383 lortensor lortensor_rank2(const std::string & n, const ex & mu, const ex & nu)
 384 {
 385         return lortensor(lortensor::lortensor_rank2,n,mu,nu);
 386 }
 387
 388 ex simplify_lortensor_mul(const ex & m)
 389 {
 390         GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
 391         exvector v_contracted;
 392
 393         // collect factors in an exvector, store squares twice
 394         int n=m.nops();
 395         v_contracted.reserve(2*n);
 396         for (int i=0; i<n; ++i) {
 397                 ex f=m.op(i);
 398                 if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
 399                         v_contracted.push_back(f.op(0));
 400                         v_contracted.push_back(f.op(0));
 401                 } else {
 402                         v_contracted.push_back(f);
 403                 }
 404         }
 405
 406         unsigned replacements;
 407         bool something_changed=false;
 408
 409         exvector::iterator it=v_contracted.begin();
 410         while (it!=v_contracted.end()) {
 411                 // process only lor_g objects
 412                 if (is_ex_exactly_of_type(*it,lortensor) &&
 413                         (ex_to_lortensor(*it).type==lortensor::lortensor_g)) {
 414                         const lortensor & g=ex_to_lortensor(*it);
 415                         GINAC_ASSERT(g.seq.size()==2);
 416                         const idx & first_idx=ex_to_lorentzidx(g.seq[0]);
 417                         const idx & second_idx=ex_to_lorentzidx(g.seq[1]);
 418                         // g_{mu,mu} should have been contracted in lortensor::eval()
 419                         GINAC_ASSERT(!first_idx.is_equal(second_idx));
 420                         ex saved_g=*it; // save to restore it later
 421
 422                         // try to contract first index
 423                         replacements=0;
 424                         if (first_idx.is_symbolic()) {
 425                                 replacements = subs_index_in_exvector(v_contracted,
 426                                                                       first_idx.toggle_covariant(),second_idx);
 427                                 if (replacements==0) {
 428                                         // not contracted, restore g object
 429                                         *it=saved_g;
 430                                 } else {
 431                                         // a contracted index should occur exactly once
 432                                         GINAC_ASSERT(replacements==1);
 433                                         *it=_ex1();
 434                                         something_changed=true;
 435                                 }
 436                         }
 437
 438                         // try second index only if first was not contracted
 439                         if ((replacements==0)&&(second_idx.is_symbolic())) {
 440                                 // first index not contracted, *it is again the original g object
 441                                 replacements = subs_index_in_exvector(v_contracted,
 442                                                                       second_idx.toggle_covariant(),first_idx);
 443                                 if (replacements==0) {
 444                                         // not contracted except in itself, restore g object
 445                                         *it=saved_g;
 446                                 } else {
 447                                         // a contracted index should occur exactly once
 448                                         GINAC_ASSERT(replacements==1);
 449                                         *it=_ex1();
 450                                         something_changed=true;
 451                                 }
 452                         }
 453                 }
 454                 ++it;
 455         }
 456         if (something_changed) {
 457                 return mul(v_contracted);
 458         }
 459         return m;
 460 }
 461
 462 ex simplify_lortensor(const ex & e)
 463 {
 464         // all simplification is done on expanded objects
 465         ex e_expanded=e.expand();
 466
 467         // simplification of sum=sum of simplifications
 468         if (is_ex_exactly_of_type(e_expanded,add)) {
 469                 ex sum=_ex0();
 470                 for (unsigned i=0; i<e_expanded.nops(); ++i) {
 471                         sum += simplify_lortensor(e_expanded.op(i));
 472                 }
 473                 return sum;
 474         }
 475
 476         // simplification of commutative product=commutative product of simplifications
 477         if (is_ex_exactly_of_type(e_expanded,mul)) {
 478                 return simplify_lortensor_mul(e);
 479         }
 480
 481         // cannot do anything
 482         return e_expanded;
 483 }
 484
 485 ex Dim(void)
 486 {
 487         static symbol * d=new symbol("dim");
 488         return *d;
 489 }
 490
 491 //////////
 492 // global constants
 493 //////////
 494
 495 const lortensor some_lortensor;
 496 const std::type_info & typeid_lortensor = typeid(some_lortensor);
 497
 498 #ifndef NO_NAMESPACE_GINAC
 499 } // namespace GiNaC
 500 #endif // ndef NO_NAMESPACE_GINAC