ginac/idx.cpp

   1 /** @file idx.cpp
   2  *
   3  *  Implementation of GiNaC's indices. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  21  */
  22
  23 #include <iostream>
  24 #include <sstream>
  25 #include <stdexcept>
  26
  27 #include "idx.h"
  28 #include "symbol.h"
  29 #include "lst.h"
  30 #include "relational.h"
  31 #include "operators.h"
  32 #include "archive.h"
  33 #include "utils.h"
  34
  35 namespace GiNaC {
  36
  37 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(idx, basic,
  38   print_func<print_context>(&idx::do_print).
  39   print_func<print_latex>(&idx::do_print_latex).
  40   print_func<print_csrc>(&idx::do_print_csrc).
  41   print_func<print_tree>(&idx::do_print_tree))
  42
  43 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(varidx, idx,
  44   print_func<print_context>(&varidx::do_print).
  45   print_func<print_latex>(&varidx::do_print_latex).
  46   print_func<print_tree>(&varidx::do_print_tree))
  47
  48 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(spinidx, varidx,
  49   print_func<print_context>(&spinidx::do_print).
  50   print_func<print_latex>(&spinidx::do_print_latex).
  51   print_func<print_tree>(&spinidx::do_print_tree))
  52
  53 //////////
  54 // default constructor
  55 //////////
  56
  57 idx::idx() : inherited(&idx::tinfo_static) {}
  58
  59 varidx::varidx() : covariant(false)
  60 {
  61         tinfo_key = &varidx::tinfo_static;
  62 }
  63
  64 spinidx::spinidx() : dotted(false)
  65 {
  66         tinfo_key = &spinidx::tinfo_static;
  67 }
  68
  69 //////////
  70 // other constructors
  71 //////////
  72
  73 idx::idx(const ex & v, const ex & d) : inherited(&idx::tinfo_static), value(v), dim(d)
  74 {
  75         if (is_dim_numeric())
  76                 if (!dim.info(info_flags::posint))
  77                         throw(std::invalid_argument("dimension of space must be a positive integer"));
  78 }
  79
  80 varidx::varidx(const ex & v, const ex & d, bool cov) : inherited(v, d), covariant(cov)
  81 {
  82         tinfo_key = &varidx::tinfo_static;
  83 }
  84
  85 spinidx::spinidx(const ex & v, const ex & d, bool cov, bool dot) : inherited(v, d, cov), dotted(dot)
  86 {
  87         tinfo_key = &spinidx::tinfo_static;
  88 }
  89
  90 //////////
  91 // archiving
  92 //////////
  93
  94 idx::idx(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
  95 {
  96         n.find_ex("value", value, sym_lst);
  97         n.find_ex("dim", dim, sym_lst);
  98 }
  99
 100 varidx::varidx(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
 101 {
 102         n.find_bool("covariant", covariant);
 103 }
 104
 105 spinidx::spinidx(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
 106 {
 107         n.find_bool("dotted", dotted);
 108 }
 109
 110 void idx::archive(archive_node &n) const
 111 {
 112         inherited::archive(n);
 113         n.add_ex("value", value);
 114         n.add_ex("dim", dim);
 115 }
 116
 117 void varidx::archive(archive_node &n) const
 118 {
 119         inherited::archive(n);
 120         n.add_bool("covariant", covariant);
 121 }
 122
 123 void spinidx::archive(archive_node &n) const
 124 {
 125         inherited::archive(n);
 126         n.add_bool("dotted", dotted);
 127 }
 128
 129 DEFAULT_UNARCHIVE(idx)
 130 DEFAULT_UNARCHIVE(varidx)
 131 DEFAULT_UNARCHIVE(spinidx)
 132
 133 //////////
 134 // functions overriding virtual functions from base classes
 135 //////////
 136
 137 void idx::print_index(const print_context & c, unsigned level) const
 138 {
 139         bool need_parens = !(is_exactly_a<numeric>(value) || is_a<symbol>(value));
 140         if (need_parens)
 141                 c.s << "(";
 142         value.print(c);
 143         if (need_parens)
 144                 c.s << ")";
 145         if (c.options & print_options::print_index_dimensions) {
 146                 c.s << "[";
 147                 dim.print(c);
 148                 c.s << "]";
 149         }
 150 }
 151
 152 void idx::do_print(const print_context & c, unsigned level) const
 153 {
 154         c.s << ".";
 155         print_index(c, level);
 156 }
 157
 158 void idx::do_print_latex(const print_latex & c, unsigned level) const
 159 {
 160         c.s << "{";
 161         print_index(c, level);
 162         c.s << "}";
 163 }
 164
 165 void idx::do_print_csrc(const print_csrc & c, unsigned level) const
 166 {
 167         c.s << "[";
 168         if (value.info(info_flags::integer))
 169                 c.s << ex_to<numeric>(value).to_int();
 170         else
 171                 value.print(c);
 172         c.s << "]";
 173 }
 174
 175 void idx::do_print_tree(const print_tree & c, unsigned level) const
 176 {
 177         c.s << std::string(level, ' ') << class_name() << " @" << this
 178             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 179             << std::endl;
 180         value.print(c, level +  c.delta_indent);
 181         dim.print(c, level + c.delta_indent);
 182 }
 183
 184 void varidx::do_print(const print_context & c, unsigned level) const
 185 {
 186         if (covariant)
 187                 c.s << ".";
 188         else
 189                 c.s << "~";
 190         print_index(c, level);
 191 }
 192
 193 void varidx::do_print_tree(const print_tree & c, unsigned level) const
 194 {
 195         c.s << std::string(level, ' ') << class_name() << " @" << this
 196             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 197             << (covariant ? ", covariant" : ", contravariant")
 198             << std::endl;
 199         value.print(c, level + c.delta_indent);
 200         dim.print(c, level + c.delta_indent);
 201 }
 202
 203 void spinidx::do_print(const print_context & c, unsigned level) const
 204 {
 205         if (covariant)
 206                 c.s << ".";
 207         else
 208                 c.s << "~";
 209         if (dotted)
 210                 c.s << "*";
 211         print_index(c, level);
 212 }
 213
 214 void spinidx::do_print_latex(const print_latex & c, unsigned level) const
 215 {
 216         if (dotted)
 217                 c.s << "\\dot{";
 218         else
 219                 c.s << "{";
 220         print_index(c, level);
 221         c.s << "}";
 222 }
 223
 224 void spinidx::do_print_tree(const print_tree & c, unsigned level) const
 225 {
 226         c.s << std::string(level, ' ') << class_name() << " @" << this
 227             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 228             << (covariant ? ", covariant" : ", contravariant")
 229             << (dotted ? ", dotted" : ", undotted")
 230             << std::endl;
 231         value.print(c, level + c.delta_indent);
 232         dim.print(c, level + c.delta_indent);
 233 }
 234
 235 bool idx::info(unsigned inf) const
 236 {
 237         switch(inf) {
 238                 case info_flags::idx:
 239                 case info_flags::has_indices:
 240                         return true;
 241         }
 242         return inherited::info(inf);
 243 }
 244
 245 size_t idx::nops() const
 246 {
 247         // don't count the dimension as that is not really a sub-expression
 248         return 1;
 249 }
 250
 251 ex idx::op(size_t i) const
 252 {
 253         GINAC_ASSERT(i == 0);
 254         return value;
 255 }
 256
 257 ex idx::map(map_function & f) const
 258 {
 259         const ex &mapped_value = f(value);
 260         if (are_ex_trivially_equal(value, mapped_value))
 261                 return *this;
 262         else {
 263                 idx *copy = duplicate();
 264                 copy->setflag(status_flags::dynallocated);
 265                 copy->clearflag(status_flags::hash_calculated);
 266                 copy->value = mapped_value;
 267                 return *copy;
 268         }
 269 }
 270
 271 /** Returns order relation between two indices of the same type. The order
 272  *  must be such that dummy indices lie next to each other. */
 273 int idx::compare_same_type(const basic & other) const
 274 {
 275         GINAC_ASSERT(is_a<idx>(other));
 276         const idx &o = static_cast<const idx &>(other);
 277
 278         int cmpval = value.compare(o.value);
 279         if (cmpval)
 280                 return cmpval;
 281         return dim.compare(o.dim);
 282 }
 283
 284 bool idx::match_same_type(const basic & other) const
 285 {
 286         GINAC_ASSERT(is_a<idx>(other));
 287         const idx &o = static_cast<const idx &>(other);
 288
 289         return dim.is_equal(o.dim);
 290 }
 291
 292 int varidx::compare_same_type(const basic & other) const
 293 {
 294         GINAC_ASSERT(is_a<varidx>(other));
 295         const varidx &o = static_cast<const varidx &>(other);
 296
 297         int cmpval = inherited::compare_same_type(other);
 298         if (cmpval)
 299                 return cmpval;
 300
 301         // Check variance last so dummy indices will end up next to each other
 302         if (covariant != o.covariant)
 303                 return covariant ? -1 : 1;
 304
 305         return 0;
 306 }
 307
 308 bool varidx::match_same_type(const basic & other) const
 309 {
 310         GINAC_ASSERT(is_a<varidx>(other));
 311         const varidx &o = static_cast<const varidx &>(other);
 312
 313         if (covariant != o.covariant)
 314                 return false;
 315
 316         return inherited::match_same_type(other);
 317 }
 318
 319 int spinidx::compare_same_type(const basic & other) const
 320 {
 321         GINAC_ASSERT(is_a<spinidx>(other));
 322         const spinidx &o = static_cast<const spinidx &>(other);
 323
 324         // Check dottedness first so dummy indices will end up next to each other
 325         if (dotted != o.dotted)
 326                 return dotted ? -1 : 1;
 327
 328         int cmpval = inherited::compare_same_type(other);
 329         if (cmpval)
 330                 return cmpval;
 331
 332         return 0;
 333 }
 334
 335 bool spinidx::match_same_type(const basic & other) const
 336 {
 337         GINAC_ASSERT(is_a<spinidx>(other));
 338         const spinidx &o = static_cast<const spinidx &>(other);
 339
 340         if (dotted != o.dotted)
 341                 return false;
 342         return inherited::match_same_type(other);
 343 }
 344
 345 unsigned idx::calchash() const
 346 {
 347         // NOTE: The code in simplify_indexed() assumes that canonically
 348         // ordered sequences of indices have the two members of dummy index
 349         // pairs lying next to each other. The hash values for indices must
 350         // be devised accordingly. The easiest (only?) way to guarantee the
 351         // desired ordering is to make indices with the same value have equal
 352         // hash keys. That is, the hash values must not depend on the index
 353         // dimensions or other attributes (variance etc.).
 354         // The compare_same_type() methods will take care of the rest.
 355         unsigned v = golden_ratio_hash((p_int)tinfo());
 356         v = rotate_left(v);
 357         v ^= value.gethash();
 358
 359         // Store calculated hash value only if object is already evaluated
 360         if (flags & status_flags::evaluated) {
 361                 setflag(status_flags::hash_calculated);
 362                 hashvalue = v;
 363         }
 364
 365         return v;
 366 }
 367
 368 /** By default, basic::evalf would evaluate the index value but we don't want
 369  *  a.1 to become a.(1.0). */
 370 ex idx::evalf(int level) const
 371 {
 372         return *this;
 373 }
 374
 375 ex idx::subs(const exmap & m, unsigned options) const
 376 {
 377         // First look for index substitutions
 378         exmap::const_iterator it = m.find(*this);
 379         if (it != m.end()) {
 380
 381                 // Substitution index->index
 382                 if (is_a<idx>(it->second) || (options & subs_options::really_subs_idx))
 383                         return it->second;
 384
 385                 // Otherwise substitute value
 386                 idx *i_copy = duplicate();
 387                 i_copy->value = it->second;
 388                 i_copy->clearflag(status_flags::hash_calculated);
 389                 return i_copy->setflag(status_flags::dynallocated);
 390         }
 391
 392         // None, substitute objects in value (not in dimension)
 393         const ex &subsed_value = value.subs(m, options);
 394         if (are_ex_trivially_equal(value, subsed_value))
 395                 return *this;
 396
 397         idx *i_copy = duplicate();
 398         i_copy->value = subsed_value;
 399         i_copy->clearflag(status_flags::hash_calculated);
 400         return i_copy->setflag(status_flags::dynallocated);
 401 }
 402
 403 /** Implementation of ex::diff() for an index always returns 0.
 404  *
 405  *  @see ex::diff */
 406 ex idx::derivative(const symbol & s) const
 407 {
 408         return _ex0;
 409 }
 410
 411 //////////
 412 // new virtual functions
 413 //////////
 414
 415 bool idx::is_dummy_pair_same_type(const basic & other) const
 416 {
 417         const idx &o = static_cast<const idx &>(other);
 418
 419         // Only pure symbols form dummy pairs, "2n+1" doesn't
 420         if (!is_a<symbol>(value))
 421                 return false;
 422
 423         // Value must be equal, of course
 424         if (!value.is_equal(o.value))
 425                 return false;
 426
 427         // Dimensions need not be equal but must be comparable (so we can
 428         // determine the minimum dimension of contractions)
 429         if (dim.is_equal(o.dim))
 430                 return true;
 431
 432         return is_exactly_a<numeric>(dim) || is_exactly_a<numeric>(o.dim);
 433 }
 434
 435 bool varidx::is_dummy_pair_same_type(const basic & other) const
 436 {
 437         const varidx &o = static_cast<const varidx &>(other);
 438
 439         // Variance must be opposite
 440         if (covariant == o.covariant)
 441                 return false;
 442
 443         return inherited::is_dummy_pair_same_type(other);
 444 }
 445
 446 bool spinidx::is_dummy_pair_same_type(const basic & other) const
 447 {
 448         const spinidx &o = static_cast<const spinidx &>(other);
 449
 450         // Dottedness must be the same
 451         if (dotted != o.dotted)
 452                 return false;
 453
 454         return inherited::is_dummy_pair_same_type(other);
 455 }
 456
 457
 458 //////////
 459 // non-virtual functions
 460 //////////
 461
 462 ex idx::replace_dim(const ex & new_dim) const
 463 {
 464         idx *i_copy = duplicate();
 465         i_copy->dim = new_dim;
 466         i_copy->clearflag(status_flags::hash_calculated);
 467         return i_copy->setflag(status_flags::dynallocated);
 468 }
 469
 470 ex idx::minimal_dim(const idx & other) const
 471 {
 472         return GiNaC::minimal_dim(dim, other.dim);
 473 }
 474
 475 ex varidx::toggle_variance() const
 476 {
 477         varidx *i_copy = duplicate();
 478         i_copy->covariant = !i_copy->covariant;
 479         i_copy->clearflag(status_flags::hash_calculated);
 480         return i_copy->setflag(status_flags::dynallocated);
 481 }
 482
 483 ex spinidx::toggle_dot() const
 484 {
 485         spinidx *i_copy = duplicate();
 486         i_copy->dotted = !i_copy->dotted;
 487         i_copy->clearflag(status_flags::hash_calculated);
 488         return i_copy->setflag(status_flags::dynallocated);
 489 }
 490
 491 ex spinidx::toggle_variance_dot() const
 492 {
 493         spinidx *i_copy = duplicate();
 494         i_copy->covariant = !i_copy->covariant;
 495         i_copy->dotted = !i_copy->dotted;
 496         i_copy->clearflag(status_flags::hash_calculated);
 497         return i_copy->setflag(status_flags::dynallocated);
 498 }
 499
 500 //////////
 501 // global functions
 502 //////////
 503
 504 bool is_dummy_pair(const idx & i1, const idx & i2)
 505 {
 506         // The indices must be of exactly the same type
 507         if (i1.tinfo() != i2.tinfo())
 508                 return false;
 509
 510         // Same type, let the indices decide whether they are paired
 511         return i1.is_dummy_pair_same_type(i2);
 512 }
 513
 514 bool is_dummy_pair(const ex & e1, const ex & e2)
 515 {
 516         // The expressions must be indices
 517         if (!is_a<idx>(e1) || !is_a<idx>(e2))
 518                 return false;
 519
 520         return is_dummy_pair(ex_to<idx>(e1), ex_to<idx>(e2));
 521 }
 522
 523 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
 524 {
 525         out_free.clear();
 526         out_dummy.clear();
 527
 528         // No indices? Then do nothing
 529         if (it == itend)
 530                 return;
 531
 532         // Only one index? Then it is a free one if it's not numeric
 533         if (itend - it == 1) {
 534                 if (ex_to<idx>(*it).is_symbolic())
 535                         out_free.push_back(*it);
 536                 return;
 537         }
 538
 539         // Sort index vector. This will cause dummy indices come to lie next
 540         // to each other (because the sort order is defined to guarantee this).
 541         exvector v(it, itend);
 542         shaker_sort(v.begin(), v.end(), ex_is_less(), ex_swap());
 543
 544         // Find dummy pairs and free indices
 545         it = v.begin(); itend = v.end();
 546         exvector::const_iterator last = it++;
 547         while (it != itend) {
 548                 if (is_dummy_pair(*it, *last)) {
 549                         out_dummy.push_back(*last);
 550                         it++;
 551                         if (it == itend)
 552                                 return;
 553                 } else {
 554                         if (!it->is_equal(*last) && ex_to<idx>(*last).is_symbolic())
 555                                 out_free.push_back(*last);
 556                 }
 557                 last = it++;
 558         }
 559         if (ex_to<idx>(*last).is_symbolic())
 560                 out_free.push_back(*last);
 561 }
 562
 563 ex minimal_dim(const ex & dim1, const ex & dim2)
 564 {
 565         if (dim1.is_equal(dim2) || dim1 < dim2 || (is_exactly_a<numeric>(dim1) && !is_a<numeric>(dim2)))
 566                 return dim1;
 567         else if (dim1 > dim2 || (!is_a<numeric>(dim1) && is_exactly_a<numeric>(dim2)))
 568                 return dim2;
 569         else {
 570                 std::ostringstream s;
 571                 s << "minimal_dim(): index dimensions " << dim1 << " and " << dim2 << " cannot be ordered";
 572                 throw (std::runtime_error(s.str()));
 573         }
 574 }
 575
 576 } // namespace GiNaC