ginac/simp_lor.cpp

   1 /** @file simp_lor.cpp
   2  *
   3  *  Implementation of GiNaC's simp_lor objects.
   4  *  No real implementation yet, to 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 "simp_lor.h"
  32 #include "ex.h"
  33 #include "mul.h"
  34 #include "symbol.h"
  35 #include "debugmsg.h"
  36 #include "utils.h"
  37
  38 #ifndef NO_NAMESPACE_GINAC
  39 namespace GiNaC {
  40 #endif // ndef NO_NAMESPACE_GINAC
  41
  42 //////////
  43 // default constructor, destructor, copy constructor assignment operator and helpers
  44 //////////
  45
  46 // public
  47
  48 simp_lor::simp_lor() : type(invalid)
  49 {
  50         debugmsg("simp_lor default constructor",LOGLEVEL_CONSTRUCT);
  51         tinfo_key=TINFO_simp_lor;
  52 }
  53
  54 simp_lor::~simp_lor()
  55 {
  56         debugmsg("simp_lor destructor",LOGLEVEL_DESTRUCT);
  57         destroy(false);
  58 }
  59
  60 simp_lor::simp_lor(const simp_lor & other)
  61 {
  62         debugmsg("simp_lor copy constructor",LOGLEVEL_CONSTRUCT);
  63         copy (other);
  64 }
  65
  66 const simp_lor & simp_lor::operator=(const simp_lor & other)
  67 {
  68         debugmsg("simp_lor operator=",LOGLEVEL_ASSIGNMENT);
  69         if (this != &other) {
  70                 destroy(true);
  71                 copy(other);
  72         }
  73         return *this;
  74 }
  75
  76 // protected
  77
  78 void simp_lor::copy(const simp_lor & other)
  79 {
  80         indexed::copy(other);
  81         type=other.type;
  82         name=other.name;
  83 }
  84
  85 void simp_lor::destroy(bool call_parent)
  86 {
  87         if (call_parent) {
  88                 indexed::destroy(call_parent);
  89         }
  90 }
  91
  92 //////////
  93 // other constructors
  94 //////////
  95
  96 // protected
  97
  98 simp_lor::simp_lor(simp_lor_types const t) : type(t)
  99 {
 100         debugmsg("simp_lor constructor from simp_lor_types",LOGLEVEL_CONSTRUCT);
 101         tinfo_key=TINFO_simp_lor;
 102 }
 103
 104 simp_lor::simp_lor(simp_lor_types const t, const ex & i1, const ex & i2)
 105   : indexed(i1,i2), type(t)
 106 {
 107         debugmsg("simp_lor constructor from simp_lor_types,ex,ex",LOGLEVEL_CONSTRUCT);
 108         tinfo_key=TINFO_simp_lor;
 109         GINAC_ASSERT(all_of_type_lorentzidx());
 110 }
 111
 112 simp_lor::simp_lor(simp_lor_types const t, const std::string & n, const ex & i1)
 113   : indexed(i1), type(t), name(n)
 114 {
 115         debugmsg("simp_lor constructor from simp_lor_types,string,ex",LOGLEVEL_CONSTRUCT);
 116         tinfo_key=TINFO_simp_lor;
 117         GINAC_ASSERT(all_of_type_lorentzidx());
 118 }
 119
 120 simp_lor::simp_lor(simp_lor_types const t, const std::string & n, const exvector & iv)
 121   : indexed(iv), type(t), name(n)
 122 {
 123         debugmsg("simp_lor constructor from simp_lor_types,string,exvector",LOGLEVEL_CONSTRUCT);
 124         tinfo_key=TINFO_simp_lor;
 125         GINAC_ASSERT(all_of_type_lorentzidx());
 126 }
 127
 128 simp_lor::simp_lor(simp_lor_types const t, const std::string & n, exvector * ivp)
 129   : indexed(ivp), type(t), name(n)
 130 {
 131         debugmsg("simp_lor constructor from simp_lor_types,string,exvector*",LOGLEVEL_CONSTRUCT);
 132         tinfo_key=TINFO_simp_lor;
 133         GINAC_ASSERT(all_of_type_lorentzidx());
 134 }
 135
 136 //////////
 137 // functions overriding virtual functions from bases classes
 138 //////////
 139
 140 // public
 141
 142 basic * simp_lor::duplicate() const
 143 {
 144         debugmsg("simp_lor duplicate",LOGLEVEL_DUPLICATE);
 145         return new simp_lor(*this);
 146 }
 147
 148 void simp_lor::printraw(std::ostream & os) const
 149 {
 150         debugmsg("simp_lor printraw",LOGLEVEL_PRINT);
 151         os << "simp_lor(type=" << (unsigned)type
 152            << ",name=" << name << ",indices=";
 153         printrawindices(os);
 154         os << ",hash=" << hashvalue << ",flags=" << flags << ")";
 155 }
 156
 157 void simp_lor::printtree(std::ostream & os, unsigned indent) const
 158 {
 159         debugmsg("simp_lor printtree",LOGLEVEL_PRINT);
 160         os << std::string(indent,' ') << "simp_lor object: "
 161            << "type=" << (unsigned)type
 162            << ", name=" << name << ", ";
 163         os << seq.size() << " indices" << std::endl;
 164         printtreeindices(os,indent);
 165         os << std::string(indent,' ') << "hash=" << hashvalue
 166            << " (0x" << std::hex << hashvalue << std::dec << ")"
 167            << ", flags=" << flags << std::endl;
 168 }
 169
 170 void simp_lor::print(std::ostream & os, unsigned upper_precedence) const
 171 {
 172         debugmsg("simp_lor print",LOGLEVEL_PRINT);
 173         switch (type) {
 174         case simp_lor_g:
 175                 os << "g";
 176                 break;
 177         case simp_lor_vec:
 178                 os << name;
 179                 break;
 180         case invalid:
 181         default:
 182                 os << "INVALID_SIMP_LOR_OBJECT";
 183                 break;
 184         }
 185         printindices(os);
 186 }
 187
 188 void simp_lor::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
 189 {
 190         debugmsg("simp_lor print csrc",LOGLEVEL_PRINT);
 191         print(os,upper_precedence);
 192 }
 193
 194 bool simp_lor::info(unsigned inf) const
 195 {
 196         return indexed::info(inf);
 197 }
 198
 199 ex simp_lor::eval(int level) const
 200 {
 201         if (type==simp_lor_g) {
 202                 // canonicalize indices
 203                 exvector iv=seq;
 204                 int sig=canonicalize_indices(iv,false); // symmetric
 205                 if (sig!=INT_MAX) {
 206                         // something has changed while sorting indices, more evaluations later
 207                         if (sig==0) return _ex0();
 208                         return ex(sig)*simp_lor(type,name,iv);
 209                 }
 210                 const lorentzidx & idx1=ex_to_lorentzidx(seq[0]);
 211                 const lorentzidx & idx2=ex_to_lorentzidx(seq[1]);
 212                 if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
 213                         // both indices are numeric
 214                         if ((idx1.get_value()==idx2.get_value())) {
 215                                 // both on diagonal
 216                                 if (idx1.get_value()==0) {
 217                                         // (0,0)
 218                                         return _ex1();
 219                                 } else {
 220                                         if (idx1.is_covariant()!=idx2.is_covariant()) {
 221                                                 // (_i,~i) or (~i,_i), i=1..3
 222                                                 return _ex1();
 223                                         } else {
 224                                                 // (_i,_i) or (~i,~i), i=1..3
 225                                                 return _ex_1();
 226                                         }
 227                                 }
 228                         } else {
 229                                 // at least one off-diagonal
 230                                 return _ex0();
 231                         }
 232                 } else if (idx1.is_symbolic() &&
 233                            idx1.is_co_contra_pair(idx2)) {
 234                         return Dim()-idx1.get_dim_parallel_space();
 235                 }
 236         }
 237
 238         return this->hold();
 239 }
 240
 241 // protected
 242
 243 int simp_lor::compare_same_type(const basic & other) const
 244 {
 245         GINAC_ASSERT(other.tinfo() == TINFO_simp_lor);
 246         const simp_lor *o = static_cast<const simp_lor *>(&other);
 247         if (type==o->type) {
 248                 if (name==o->name) {
 249                         return indexed::compare_same_type(other);
 250                 }
 251                 return name.compare(o->name);
 252         }
 253         return type < o->type ? -1 : 1;
 254 }
 255
 256 bool simp_lor::is_equal_same_type(const basic & other) const
 257 {
 258         GINAC_ASSERT(other.tinfo() == TINFO_simp_lor);
 259         const simp_lor *o = static_cast<const simp_lor *>(&other);
 260         if (type!=o->type) return false;
 261         if (name!=o->name) return false;
 262         return indexed::is_equal_same_type(other);
 263 }
 264
 265 unsigned simp_lor::return_type(void) const
 266 {
 267         return return_types::commutative;
 268 }
 269
 270 unsigned simp_lor::return_type_tinfo(void) const
 271 {
 272         return tinfo_key;
 273 }
 274
 275 ex simp_lor::thisexprseq(const exvector & v) const
 276 {
 277         return simp_lor(type,name,v);
 278 }
 279
 280 ex simp_lor::thisexprseq(exvector * vp) const
 281 {
 282         return simp_lor(type,name,vp);
 283 }
 284
 285 //////////
 286 // virtual functions which can be overridden by derived classes
 287 //////////
 288
 289 // none
 290
 291 //////////
 292 // non-virtual functions in this class
 293 //////////
 294
 295 // protected
 296
 297 bool simp_lor::all_of_type_lorentzidx(void) const
 298 {
 299         // used only inside of ASSERTs
 300         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 301                 if (!is_ex_of_type(*cit,lorentzidx)) return false;
 302         }
 303         return true;
 304 }
 305
 306 //////////
 307 // static member variables
 308 //////////
 309
 310 // none
 311
 312 //////////
 313 // global constants
 314 //////////
 315
 316 const simp_lor some_simp_lor;
 317 const std::type_info & typeid_simp_lor = typeid(some_simp_lor);
 318
 319 //////////
 320 // friend functions
 321 //////////
 322
 323 simp_lor lor_g(const ex & mu, const ex & nu)
 324 {
 325         return simp_lor(simp_lor::simp_lor_g,mu,nu);
 326 }
 327
 328 simp_lor lor_vec(const std::string & n, const ex & mu)
 329 {
 330         return simp_lor(simp_lor::simp_lor_vec,n,mu);
 331 }
 332
 333 ex simplify_simp_lor_mul(const ex & m, const scalar_products & sp)
 334 {
 335         GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
 336         exvector v_contracted;
 337
 338         // collect factors in an exvector, store squares twice
 339         unsigned n=m.nops();
 340         v_contracted.reserve(2*n);
 341         for (unsigned i=0; i<n; ++i) {
 342                 ex f=m.op(i);
 343                 if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
 344                         v_contracted.push_back(f.op(0));
 345                         v_contracted.push_back(f.op(0));
 346                 } else {
 347                         v_contracted.push_back(f);
 348                 }
 349         }
 350
 351         unsigned replacements;
 352         bool something_changed=false;
 353
 354         exvector::iterator it=v_contracted.begin();
 355         while (it!=v_contracted.end()) {
 356                 // process only lor_g objects
 357                 if (is_ex_exactly_of_type(*it,simp_lor) &&
 358                         (ex_to_simp_lor(*it).type==simp_lor::simp_lor_g)) {
 359                         const simp_lor & g=ex_to_simp_lor(*it);
 360                         GINAC_ASSERT(g.seq.size()==2);
 361                         const idx & first_idx=ex_to_lorentzidx(g.seq[0]);
 362                         const idx & second_idx=ex_to_lorentzidx(g.seq[1]);
 363                         // g_{mu,mu} should have been contracted in simp_lor::eval()
 364                         GINAC_ASSERT(!first_idx.is_equal(second_idx));
 365                         ex saved_g=*it; // save to restore it later
 366
 367                         // try to contract first index
 368                         replacements=0;
 369                         if (first_idx.is_symbolic()) {
 370                                 replacements = subs_index_in_exvector(v_contracted, first_idx.toggle_covariant(),second_idx);
 371                                 if (replacements==0) {
 372                                         // not contracted, restore g object
 373                                         *it=saved_g;
 374                                 } else {
 375                                         // a contracted index should occur exactly once
 376                                         GINAC_ASSERT(replacements==1);
 377                                         *it=_ex1();
 378                                         something_changed=true;
 379                                 }
 380                         }
 381
 382                         // try second index only if first was not contracted
 383                         if ((replacements==0)&&(second_idx.is_symbolic())) {
 384                                 // first index not contracted, *it is again the original g object
 385                                 replacements = subs_index_in_exvector(v_contracted, second_idx.toggle_covariant(),first_idx);
 386                                 if (replacements==0) {
 387                                         // not contracted except in itself, restore g object
 388                                         *it=saved_g;
 389                                 } else {
 390                                         // a contracted index should occur exactly once
 391                                         GINAC_ASSERT(replacements==1);
 392                                         *it=_ex1();
 393                                         something_changed=true;
 394                                 }
 395                         }
 396                 }
 397                 ++it;
 398         }
 399
 400         // process only lor_vec objects
 401         bool jump_to_next=false;
 402         exvector::iterator it1=v_contracted.begin();
 403         while (it1!=v_contracted.end()-1) {
 404                 if (is_ex_exactly_of_type(*it1,simp_lor) &&
 405                         (ex_to_simp_lor(*it1).type==simp_lor::simp_lor_vec)) {
 406                         exvector::iterator it2=it1+1;
 407                         while ((it2!=v_contracted.end())&&!jump_to_next) {
 408                                 if (is_ex_exactly_of_type(*it2,simp_lor) &&
 409                                         (ex_to_simp_lor(*it2).type==simp_lor::simp_lor_vec)) {
 410                                         const simp_lor & vec1=ex_to_simp_lor(*it1);
 411                                         const simp_lor & vec2=ex_to_simp_lor(*it2);
 412                                         GINAC_ASSERT(vec1.seq.size()==1);
 413                                         GINAC_ASSERT(vec2.seq.size()==1);
 414                                         const lorentzidx & idx1=ex_to_lorentzidx(vec1.seq[0]);
 415                                         const lorentzidx & idx2=ex_to_lorentzidx(vec2.seq[0]);
 416                                         if (idx1.is_symbolic() &&
 417                                                 idx1.is_co_contra_pair(idx2) &&
 418                                                 sp.is_defined(vec1,vec2)) {
 419                                                 *it1=sp.evaluate(vec1,vec2);
 420                                                 *it2=_ex1();
 421                                                 something_changed=true;
 422                                                 jump_to_next=true;
 423                                         }
 424                                 }
 425                                 ++it2;
 426                         }
 427                         jump_to_next=false;
 428                 }
 429                 ++it1;
 430         }
 431         if (something_changed) {
 432                 return mul(v_contracted);
 433         }
 434         return m;
 435 }
 436
 437 ex simplify_simp_lor(const ex & e, const scalar_products & sp)
 438 {
 439         // all simplification is done on expanded objects
 440         ex e_expanded=e.expand();
 441
 442         // simplification of sum=sum of simplifications
 443         if (is_ex_exactly_of_type(e_expanded,add)) {
 444                 ex sum=_ex0();
 445                 for (unsigned i=0; i<e_expanded.nops(); ++i)
 446                         sum += simplify_simp_lor(e_expanded.op(i),sp);
 447
 448                 return sum;
 449         }
 450
 451         // simplification of commutative product=commutative product of simplifications
 452         if (is_ex_exactly_of_type(e_expanded,mul)) {
 453                 return simplify_simp_lor_mul(e,sp);
 454         }
 455
 456         // cannot do anything
 457         return e_expanded;
 458 }
 459
 460 //ex Dim(void)   // FIXME: what's going on here?
 461 //{
 462 //    static symbol * d=new symbol("dim");
 463 //    return *d;
 464 //}
 465
 466 //////////
 467 // helper classes
 468 //////////
 469
 470 void scalar_products::reg(const simp_lor & v1, const simp_lor & v2,
 471                           const ex & sp)
 472 {
 473         if (v1.compare_same_type(v2)>0) {
 474                 reg(v2,v1,sp);
 475                 return;
 476         }
 477         spm[make_key(v1,v2)]=sp;
 478 }
 479
 480 bool scalar_products::is_defined(const simp_lor & v1, const simp_lor & v2) const
 481 {
 482         if (v1.compare_same_type(v2)>0) {
 483                 return is_defined(v2,v1);
 484         }
 485         return spm.find(make_key(v1,v2))!=spm.end();
 486 }
 487
 488 ex scalar_products::evaluate(const simp_lor & v1, const simp_lor & v2) const
 489 {
 490         if (v1.compare_same_type(v2)>0)
 491                 return evaluate(v2, v1);
 492
 493         return (*spm.find(make_key(v1,v2))).second;
 494 }
 495
 496 void scalar_products::debugprint(void) const
 497 {
 498         std::cerr << "map size=" << spm.size() << std::endl;
 499         for (spmap::const_iterator cit=spm.begin(); cit!=spm.end(); ++cit) {
 500                 const spmapkey & k=(*cit).first;
 501                 std::cerr << "item key=((" << k.first.first
 502                           << "," << k.first.second << "),";
 503                 k.second.printraw(std::cerr);
 504                 std::cerr << ") value=" << (*cit).second << std::endl;
 505         }
 506 }
 507
 508 spmapkey scalar_products::make_key(const simp_lor & v1, const simp_lor & v2)
 509 {
 510         GINAC_ASSERT(v1.type==simp_lor::simp_lor_vec);
 511         GINAC_ASSERT(v2.type==simp_lor::simp_lor_vec);
 512         lorentzidx anon=ex_to_lorentzidx(v1.seq[0]).create_anonymous_representative();
 513         GINAC_ASSERT(anon.is_equal_same_type(ex_to_lorentzidx(v2.seq[0]).create_anonymous_representative()));
 514         return spmapkey(strstrpair(v1.name,v2.name),anon);
 515 }
 516
 517 #ifndef NO_NAMESPACE_GINAC
 518 } // namespace GiNaC
 519 #endif // ndef NO_NAMESPACE_GINAC