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