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