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(0);
  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(1);
  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 unsigned lortensor::return_type_tinfo(void) const
 312 {
 313         return tinfo_key;
 314 }
 315 ex lortensor::thisexprseq(const exvector & v) const
 316 {
 317         return lortensor(type,name,serial,v);
 318 }
 319 ex lortensor::thisexprseq(exvector *vp) const
 320 {
 321         return lortensor(type,name,serial,vp);
 322 }
 323
 324 //////////
 325 // non-virtual functions in this class
 326 //////////
 327
 328 // protected
 329
 330 void lortensor::setname(const std::string & n)
 331 {
 332         name = n;
 333 }
 334
 335 bool lortensor::all_of_type_lorentzidx(void) const
 336 {
 337         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++ cit) {
 338                 if (!is_ex_of_type(*cit,lorentzidx)) {
 339                         return false;
 340                 }
 341         }
 342         return true;
 343 }
 344
 345 // private
 346
 347 std::string & lortensor::autoname_prefix(void)
 348 {
 349         static std::string * s = new std::string("lortensor");
 350         return *s;
 351 }
 352
 353 //////////
 354 // static member variables
 355 //////////
 356
 357 // private
 358
 359 unsigned lortensor::next_serial=0;
 360
 361 //////////
 362 // friend functions
 363 //////////
 364
 365 lortensor lortensor_g(const ex & mu, const ex & nu)
 366 {
 367         return lortensor(lortensor::lortensor_g,"",mu,nu);
 368 }
 369
 370 lortensor lortensor_epsilon(const ex & mu, const ex & nu, const ex & rho, const ex & sigma)
 371 {
 372         return lortensor(lortensor::lortensor_epsilon,"",mu,nu,rho,sigma);
 373 }
 374
 375 lortensor lortensor_rank1(const std::string & n, const ex & mu)
 376 {
 377         return lortensor(lortensor::lortensor_rank1,n,mu);
 378 }
 379
 380 lortensor lortensor_rank2(const std::string & n, const ex & mu, const ex & nu)
 381 {
 382         return lortensor(lortensor::lortensor_rank2,n,mu,nu);
 383 }
 384
 385 ex simplify_lortensor_mul(const ex & m)
 386 {
 387         GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
 388         exvector v_contracted;
 389
 390         // collect factors in an exvector, store squares twice
 391         int n=m.nops();
 392         v_contracted.reserve(2*n);
 393         for (int i=0; i<n; ++i) {
 394                 ex f=m.op(i);
 395                 if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
 396                         v_contracted.push_back(f.op(0));
 397                         v_contracted.push_back(f.op(0));
 398                 } else {
 399                         v_contracted.push_back(f);
 400                 }
 401         }
 402
 403         unsigned replacements;
 404         bool something_changed=false;
 405
 406         exvector::iterator it=v_contracted.begin();
 407         while (it!=v_contracted.end()) {
 408                 // process only lor_g objects
 409                 if (is_ex_exactly_of_type(*it,lortensor) &&
 410                         (ex_to_lortensor(*it).type==lortensor::lortensor_g)) {
 411                         const lortensor & g=ex_to_lortensor(*it);
 412                         GINAC_ASSERT(g.seq.size()==2);
 413                         const idx & first_idx=ex_to_lorentzidx(g.seq[0]);
 414                         const idx & second_idx=ex_to_lorentzidx(g.seq[1]);
 415                         // g_{mu,mu} should have been contracted in lortensor::eval()
 416                         GINAC_ASSERT(!first_idx.is_equal(second_idx));
 417                         ex saved_g=*it; // save to restore it later
 418
 419                         // try to contract first index
 420                         replacements=0;
 421                         if (first_idx.is_symbolic()) {
 422                                 replacements = subs_index_in_exvector(v_contracted,
 423                                                                    first_idx.toggle_covariant(),second_idx);
 424                                 if (replacements==0) {
 425                                         // not contracted, restore g object
 426                                         *it=saved_g;
 427                                 } else {
 428                                         // a contracted index should occur exactly once
 429                                         GINAC_ASSERT(replacements==1);
 430                                         *it=_ex1();
 431                                         something_changed=true;
 432                                 }
 433                         }
 434
 435                         // try second index only if first was not contracted
 436                         if ((replacements==0)&&(second_idx.is_symbolic())) {
 437                                 // first index not contracted, *it is again the original g object
 438                                 replacements = subs_index_in_exvector(v_contracted,
 439                                                                    second_idx.toggle_covariant(),first_idx);
 440                                 if (replacements==0) {
 441                                         // not contracted except in itself, restore g object
 442                                         *it=saved_g;
 443                                 } else {
 444                                         // a contracted index should occur exactly once
 445                                         GINAC_ASSERT(replacements==1);
 446                                         *it=_ex1();
 447                                         something_changed=true;
 448                                 }
 449                         }
 450                 }
 451                 ++it;
 452         }
 453         if (something_changed) {
 454                 return mul(v_contracted);
 455         }
 456         return m;
 457 }
 458
 459 ex simplify_lortensor(const ex & e)
 460 {
 461         // all simplification is done on expanded objects
 462         ex e_expanded=e.expand();
 463
 464         // simplification of sum=sum of simplifications
 465         if (is_ex_exactly_of_type(e_expanded,add)) {
 466                 ex sum=_ex0();
 467                 for (unsigned i=0; i<e_expanded.nops(); ++i) {
 468                         sum += simplify_lortensor(e_expanded.op(i));
 469                 }
 470                 return sum;
 471         }
 472
 473         // simplification of commutative product=commutative product of simplifications
 474         if (is_ex_exactly_of_type(e_expanded,mul)) {
 475                 return simplify_lortensor_mul(e);
 476         }
 477
 478         // cannot do anything
 479         return e_expanded;
 480 }
 481
 482 ex Dim(void)
 483 {
 484         static symbol * d=new symbol("dim");
 485         return *d;
 486 }
 487
 488 //////////
 489 // global constants
 490 //////////
 491
 492 const lortensor some_lortensor;
 493 const type_info & typeid_lortensor=typeid(some_lortensor);
 494
 495 #ifndef NO_NAMESPACE_GINAC
 496 } // namespace GiNaC
 497 #endif // ndef NO_NAMESPACE_GINAC