ginac/idx.cpp

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