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