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