ginac/idx.cpp

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