ginac/idx.cpp

   1 /** @file idx.cpp
   2  *
   3  *  Implementation of GiNaC's indices. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2004 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  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_latex inherited from idx
  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(TINFO_idx) {}
  57
  58 varidx::varidx() : covariant(false)
  59 {
  60         tinfo_key = TINFO_varidx;
  61 }
  62
  63 spinidx::spinidx() : dotted(false)
  64 {
  65         tinfo_key = TINFO_spinidx;
  66 }
  67
  68 //////////
  69 // other constructors
  70 //////////
  71
  72 idx::idx(const ex & v, const ex & d) : inherited(TINFO_idx), 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 = TINFO_varidx;
  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 = TINFO_spinidx;
  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         idx *copy = duplicate();
 246         copy->setflag(status_flags::dynallocated);
 247         copy->clearflag(status_flags::hash_calculated);
 248         copy->value = f(value);
 249         return *copy;
 250 }
 251
 252 /** Returns order relation between two indices of the same type. The order
 253  *  must be such that dummy indices lie next to each other. */
 254 int idx::compare_same_type(const basic & other) const
 255 {
 256         GINAC_ASSERT(is_a<idx>(other));
 257         const idx &o = static_cast<const idx &>(other);
 258
 259         int cmpval = value.compare(o.value);
 260         if (cmpval)
 261                 return cmpval;
 262         return dim.compare(o.dim);
 263 }
 264
 265 bool idx::match_same_type(const basic & other) const
 266 {
 267         GINAC_ASSERT(is_a<idx>(other));
 268         const idx &o = static_cast<const idx &>(other);
 269
 270         return dim.is_equal(o.dim);
 271 }
 272
 273 int varidx::compare_same_type(const basic & other) const
 274 {
 275         GINAC_ASSERT(is_a<varidx>(other));
 276         const varidx &o = static_cast<const varidx &>(other);
 277
 278         int cmpval = inherited::compare_same_type(other);
 279         if (cmpval)
 280                 return cmpval;
 281
 282         // Check variance last so dummy indices will end up next to each other
 283         if (covariant != o.covariant)
 284                 return covariant ? -1 : 1;
 285
 286         return 0;
 287 }
 288
 289 bool varidx::match_same_type(const basic & other) const
 290 {
 291         GINAC_ASSERT(is_a<varidx>(other));
 292         const varidx &o = static_cast<const varidx &>(other);
 293
 294         if (covariant != o.covariant)
 295                 return false;
 296
 297         return inherited::match_same_type(other);
 298 }
 299
 300 int spinidx::compare_same_type(const basic & other) const
 301 {
 302         GINAC_ASSERT(is_a<spinidx>(other));
 303         const spinidx &o = static_cast<const spinidx &>(other);
 304
 305         // Check dottedness first so dummy indices will end up next to each other
 306         if (dotted != o.dotted)
 307                 return dotted ? -1 : 1;
 308
 309         int cmpval = inherited::compare_same_type(other);
 310         if (cmpval)
 311                 return cmpval;
 312
 313         return 0;
 314 }
 315
 316 bool spinidx::match_same_type(const basic & other) const
 317 {
 318         GINAC_ASSERT(is_a<spinidx>(other));
 319         const spinidx &o = static_cast<const spinidx &>(other);
 320
 321         if (dotted != o.dotted)
 322                 return false;
 323         return inherited::match_same_type(other);
 324 }
 325
 326 unsigned idx::calchash() const
 327 {
 328         // NOTE: The code in simplify_indexed() assumes that canonically
 329         // ordered sequences of indices have the two members of dummy index
 330         // pairs lying next to each other. The hash values for indices must
 331         // be devised accordingly. The easiest (only?) way to guarantee the
 332         // desired ordering is to make indices with the same value have equal
 333         // hash keys. That is, the hash values must not depend on the index
 334         // dimensions or other attributes (variance etc.).
 335         // The compare_same_type() methods will take care of the rest.
 336         unsigned v = golden_ratio_hash(tinfo());
 337         v = rotate_left(v);
 338         v ^= value.gethash();
 339
 340         // Store calculated hash value only if object is already evaluated
 341         if (flags & status_flags::evaluated) {
 342                 setflag(status_flags::hash_calculated);
 343                 hashvalue = v;
 344         }
 345
 346         return v;
 347 }
 348
 349 /** By default, basic::evalf would evaluate the index value but we don't want
 350  *  a.1 to become a.(1.0). */
 351 ex idx::evalf(int level) const
 352 {
 353         return *this;
 354 }
 355
 356 ex idx::subs(const exmap & m, unsigned options) const
 357 {
 358         // First look for index substitutions
 359         exmap::const_iterator it = m.find(*this);
 360         if (it != m.end()) {
 361
 362                 // Substitution index->index
 363                 if (is_a<idx>(it->second))
 364                         return it->second;
 365
 366                 // Otherwise substitute value
 367                 idx *i_copy = duplicate();
 368                 i_copy->value = it->second;
 369                 i_copy->clearflag(status_flags::hash_calculated);
 370                 return i_copy->setflag(status_flags::dynallocated);
 371         }
 372
 373         // None, substitute objects in value (not in dimension)
 374         const ex &subsed_value = value.subs(m, options);
 375         if (are_ex_trivially_equal(value, subsed_value))
 376                 return *this;
 377
 378         idx *i_copy = duplicate();
 379         i_copy->value = subsed_value;
 380         i_copy->clearflag(status_flags::hash_calculated);
 381         return i_copy->setflag(status_flags::dynallocated);
 382 }
 383
 384 /** Implementation of ex::diff() for an index always returns 0.
 385  *
 386  *  @see ex::diff */
 387 ex idx::derivative(const symbol & s) const
 388 {
 389         return _ex0;
 390 }
 391
 392 //////////
 393 // new virtual functions
 394 //////////
 395
 396 bool idx::is_dummy_pair_same_type(const basic & other) const
 397 {
 398         const idx &o = static_cast<const idx &>(other);
 399
 400         // Only pure symbols form dummy pairs, "2n+1" doesn't
 401         if (!is_a<symbol>(value))
 402                 return false;
 403
 404         // Value must be equal, of course
 405         if (!value.is_equal(o.value))
 406                 return false;
 407
 408         // Dimensions need not be equal but must be comparable (so we can
 409         // determine the minimum dimension of contractions)
 410         if (dim.is_equal(o.dim))
 411                 return true;
 412
 413         return (dim < o.dim || dim > o.dim || (is_exactly_a<numeric>(dim) && is_a<symbol>(o.dim)) || (is_a<symbol>(dim) && is_exactly_a<numeric>(o.dim)));
 414 }
 415
 416 bool varidx::is_dummy_pair_same_type(const basic & other) const
 417 {
 418         const varidx &o = static_cast<const varidx &>(other);
 419
 420         // Variance must be opposite
 421         if (covariant == o.covariant)
 422                 return false;
 423
 424         return inherited::is_dummy_pair_same_type(other);
 425 }
 426
 427 bool spinidx::is_dummy_pair_same_type(const basic & other) const
 428 {
 429         const spinidx &o = static_cast<const spinidx &>(other);
 430
 431         // Dottedness must be the same
 432         if (dotted != o.dotted)
 433                 return false;
 434
 435         return inherited::is_dummy_pair_same_type(other);
 436 }
 437
 438
 439 //////////
 440 // non-virtual functions
 441 //////////
 442
 443 ex idx::replace_dim(const ex & new_dim) const
 444 {
 445         idx *i_copy = duplicate();
 446         i_copy->dim = new_dim;
 447         i_copy->clearflag(status_flags::hash_calculated);
 448         return i_copy->setflag(status_flags::dynallocated);
 449 }
 450
 451 ex idx::minimal_dim(const idx & other) const
 452 {
 453         return GiNaC::minimal_dim(dim, other.dim);
 454 }
 455
 456 ex varidx::toggle_variance() const
 457 {
 458         varidx *i_copy = duplicate();
 459         i_copy->covariant = !i_copy->covariant;
 460         i_copy->clearflag(status_flags::hash_calculated);
 461         return i_copy->setflag(status_flags::dynallocated);
 462 }
 463
 464 ex spinidx::toggle_dot() const
 465 {
 466         spinidx *i_copy = duplicate();
 467         i_copy->dotted = !i_copy->dotted;
 468         i_copy->clearflag(status_flags::hash_calculated);
 469         return i_copy->setflag(status_flags::dynallocated);
 470 }
 471
 472 ex spinidx::toggle_variance_dot() const
 473 {
 474         spinidx *i_copy = duplicate();
 475         i_copy->covariant = !i_copy->covariant;
 476         i_copy->dotted = !i_copy->dotted;
 477         i_copy->clearflag(status_flags::hash_calculated);
 478         return i_copy->setflag(status_flags::dynallocated);
 479 }
 480
 481 //////////
 482 // global functions
 483 //////////
 484
 485 bool is_dummy_pair(const idx & i1, const idx & i2)
 486 {
 487         // The indices must be of exactly the same type
 488         if (i1.tinfo() != i2.tinfo())
 489                 return false;
 490
 491         // Same type, let the indices decide whether they are paired
 492         return i1.is_dummy_pair_same_type(i2);
 493 }
 494
 495 bool is_dummy_pair(const ex & e1, const ex & e2)
 496 {
 497         // The expressions must be indices
 498         if (!is_a<idx>(e1) || !is_a<idx>(e2))
 499                 return false;
 500
 501         return is_dummy_pair(ex_to<idx>(e1), ex_to<idx>(e2));
 502 }
 503
 504 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
 505 {
 506         out_free.clear();
 507         out_dummy.clear();
 508
 509         // No indices? Then do nothing
 510         if (it == itend)
 511                 return;
 512
 513         // Only one index? Then it is a free one if it's not numeric
 514         if (itend - it == 1) {
 515                 if (ex_to<idx>(*it).is_symbolic())
 516                         out_free.push_back(*it);
 517                 return;
 518         }
 519
 520         // Sort index vector. This will cause dummy indices come to lie next
 521         // to each other (because the sort order is defined to guarantee this).
 522         exvector v(it, itend);
 523         shaker_sort(v.begin(), v.end(), ex_is_less(), ex_swap());
 524
 525         // Find dummy pairs and free indices
 526         it = v.begin(); itend = v.end();
 527         exvector::const_iterator last = it++;
 528         while (it != itend) {
 529                 if (is_dummy_pair(*it, *last)) {
 530                         out_dummy.push_back(*last);
 531                         it++;
 532                         if (it == itend)
 533                                 return;
 534                 } else {
 535                         if (!it->is_equal(*last) && ex_to<idx>(*last).is_symbolic())
 536                                 out_free.push_back(*last);
 537                 }
 538                 last = it++;
 539         }
 540         if (ex_to<idx>(*last).is_symbolic())
 541                 out_free.push_back(*last);
 542 }
 543
 544 ex minimal_dim(const ex & dim1, const ex & dim2)
 545 {
 546         if (dim1.is_equal(dim2) || dim1 < dim2 || (is_exactly_a<numeric>(dim1) && is_a<symbol>(dim2)))
 547                 return dim1;
 548         else if (dim1 > dim2 || (is_a<symbol>(dim1) && is_exactly_a<numeric>(dim2)))
 549                 return dim2;
 550         else {
 551                 std::ostringstream s;
 552                 s << "minimal_dim(): index dimensions " << dim1 << " and " << dim2 << " cannot be ordered";
 553                 throw (std::runtime_error(s.str()));
 554         }
 555 }
 556
 557 } // namespace GiNaC