ginac/idx.cpp

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