ginac/idx.cpp

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