ginac/idx.cpp

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