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, const 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, const lst &sym_lst) : inherited(n, sym_lst)
 113 {
 114         n.find_bool("covariant", covariant);
 115 }
 116
 117 spinidx::spinidx(const archive_node &n, const 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 unsigned 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::let_op(int i)
 279 {
 280         GINAC_ASSERT(i == 0);
 281         return value;
 282 }
 283
 284 /** Returns order relation between two indices of the same type. The order
 285  *  must be such that dummy indices lie next to each other. */
 286 int idx::compare_same_type(const basic & other) const
 287 {
 288         GINAC_ASSERT(is_a<idx>(other));
 289         const idx &o = static_cast<const idx &>(other);
 290
 291         int cmpval = value.compare(o.value);
 292         if (cmpval)
 293                 return cmpval;
 294         return dim.compare(o.dim);
 295 }
 296
 297 bool idx::match_same_type(const basic & other) const
 298 {
 299         GINAC_ASSERT(is_a<idx>(other));
 300         const idx &o = static_cast<const idx &>(other);
 301
 302         return dim.is_equal(o.dim);
 303 }
 304
 305 int varidx::compare_same_type(const basic & other) const
 306 {
 307         GINAC_ASSERT(is_a<varidx>(other));
 308         const varidx &o = static_cast<const varidx &>(other);
 309
 310         int cmpval = inherited::compare_same_type(other);
 311         if (cmpval)
 312                 return cmpval;
 313
 314         // Check variance last so dummy indices will end up next to each other
 315         if (covariant != o.covariant)
 316                 return covariant ? -1 : 1;
 317
 318         return 0;
 319 }
 320
 321 bool varidx::match_same_type(const basic & other) const
 322 {
 323         GINAC_ASSERT(is_a<varidx>(other));
 324         const varidx &o = static_cast<const varidx &>(other);
 325
 326         if (covariant != o.covariant)
 327                 return false;
 328
 329         return inherited::match_same_type(other);
 330 }
 331
 332 int spinidx::compare_same_type(const basic & other) const
 333 {
 334         GINAC_ASSERT(is_a<spinidx>(other));
 335         const spinidx &o = static_cast<const spinidx &>(other);
 336
 337         // Check dottedness first so dummy indices will end up next to each other
 338         if (dotted != o.dotted)
 339                 return dotted ? -1 : 1;
 340
 341         int cmpval = inherited::compare_same_type(other);
 342         if (cmpval)
 343                 return cmpval;
 344
 345         return 0;
 346 }
 347
 348 bool spinidx::match_same_type(const basic & other) const
 349 {
 350         GINAC_ASSERT(is_a<spinidx>(other));
 351         const spinidx &o = static_cast<const spinidx &>(other);
 352
 353         if (dotted != o.dotted)
 354                 return false;
 355         return inherited::match_same_type(other);
 356 }
 357
 358 /** By default, basic::evalf would evaluate the index value but we don't want
 359  *  a.1 to become a.(1.0). */
 360 ex idx::evalf(int level) const
 361 {
 362         return *this;
 363 }
 364
 365 ex idx::subs(const lst & ls, const lst & lr, unsigned options) const
 366 {
 367         GINAC_ASSERT(ls.nops() == lr.nops());
 368
 369         // First look for index substitutions
 370         for (unsigned i=0; i<ls.nops(); i++) {
 371                 if (is_equal(ex_to<basic>(ls.op(i)))) {
 372
 373                         // Substitution index->index
 374                         if (is_a<idx>(lr.op(i)))
 375                                 return lr.op(i);
 376
 377                         // Otherwise substitute value
 378                         idx *i_copy = static_cast<idx *>(duplicate());
 379                         i_copy->value = lr.op(i);
 380                         i_copy->clearflag(status_flags::hash_calculated);
 381                         return i_copy->setflag(status_flags::dynallocated);
 382                 }
 383         }
 384
 385         // None, substitute objects in value (not in dimension)
 386         const ex &subsed_value = value.subs(ls, lr, options);
 387         if (are_ex_trivially_equal(value, subsed_value))
 388                 return *this;
 389
 390         idx *i_copy = static_cast<idx *>(duplicate());
 391         i_copy->value = subsed_value;
 392         i_copy->clearflag(status_flags::hash_calculated);
 393         return i_copy->setflag(status_flags::dynallocated);
 394 }
 395
 396 /** Implementation of ex::diff() for an index always returns 0.
 397  *
 398  *  @see ex::diff */
 399 ex idx::derivative(const symbol & s) const
 400 {
 401         return _ex0;
 402 }
 403
 404 //////////
 405 // new virtual functions
 406 //////////
 407
 408 bool idx::is_dummy_pair_same_type(const basic & other) const
 409 {
 410         const idx &o = static_cast<const idx &>(other);
 411
 412         // Only pure symbols form dummy pairs, "2n+1" doesn't
 413         if (!is_a<symbol>(value))
 414                 return false;
 415
 416         // Value must be equal, of course
 417         if (!value.is_equal(o.value))
 418                 return false;
 419
 420         // Dimensions need not be equal but must be comparable (so we can
 421         // determine the minimum dimension of contractions)
 422         if (dim.is_equal(o.dim))
 423                 return true;
 424
 425         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)));
 426 }
 427
 428 bool varidx::is_dummy_pair_same_type(const basic & other) const
 429 {
 430         const varidx &o = static_cast<const varidx &>(other);
 431
 432         // Variance must be opposite
 433         if (covariant == o.covariant)
 434                 return false;
 435
 436         return inherited::is_dummy_pair_same_type(other);
 437 }
 438
 439 bool spinidx::is_dummy_pair_same_type(const basic & other) const
 440 {
 441         const spinidx &o = static_cast<const spinidx &>(other);
 442
 443         // Dottedness must be the same
 444         if (dotted != o.dotted)
 445                 return false;
 446
 447         return inherited::is_dummy_pair_same_type(other);
 448 }
 449
 450
 451 //////////
 452 // non-virtual functions
 453 //////////
 454
 455 ex idx::replace_dim(const ex & new_dim) const
 456 {
 457         idx *i_copy = static_cast<idx *>(duplicate());
 458         i_copy->dim = new_dim;
 459         i_copy->clearflag(status_flags::hash_calculated);
 460         return i_copy->setflag(status_flags::dynallocated);
 461 }
 462
 463 ex idx::minimal_dim(const idx & other) const
 464 {
 465         return GiNaC::minimal_dim(dim, other.dim);
 466 }
 467
 468 ex varidx::toggle_variance(void) const
 469 {
 470         varidx *i_copy = static_cast<varidx *>(duplicate());
 471         i_copy->covariant = !i_copy->covariant;
 472         i_copy->clearflag(status_flags::hash_calculated);
 473         return i_copy->setflag(status_flags::dynallocated);
 474 }
 475
 476 ex spinidx::toggle_dot(void) const
 477 {
 478         spinidx *i_copy = static_cast<spinidx *>(duplicate());
 479         i_copy->dotted = !i_copy->dotted;
 480         i_copy->clearflag(status_flags::hash_calculated);
 481         return i_copy->setflag(status_flags::dynallocated);
 482 }
 483
 484 ex spinidx::toggle_variance_dot(void) const
 485 {
 486         spinidx *i_copy = static_cast<spinidx *>(duplicate());
 487         i_copy->covariant = !i_copy->covariant;
 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 //////////
 494 // global functions
 495 //////////
 496
 497 bool is_dummy_pair(const idx & i1, const idx & i2)
 498 {
 499         // The indices must be of exactly the same type
 500         if (i1.tinfo() != i2.tinfo())
 501                 return false;
 502
 503         // Same type, let the indices decide whether they are paired
 504         return i1.is_dummy_pair_same_type(i2);
 505 }
 506
 507 bool is_dummy_pair(const ex & e1, const ex & e2)
 508 {
 509         // The expressions must be indices
 510         if (!is_a<idx>(e1) || !is_a<idx>(e2))
 511                 return false;
 512
 513         return is_dummy_pair(ex_to<idx>(e1), ex_to<idx>(e2));
 514 }
 515
 516 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
 517 {
 518         out_free.clear();
 519         out_dummy.clear();
 520
 521         // No indices? Then do nothing
 522         if (it == itend)
 523                 return;
 524
 525         // Only one index? Then it is a free one if it's not numeric
 526         if (itend - it == 1) {
 527                 if (ex_to<idx>(*it).is_symbolic())
 528                         out_free.push_back(*it);
 529                 return;
 530         }
 531
 532         // Sort index vector. This will cause dummy indices come to lie next
 533         // to each other (because the sort order is defined to guarantee this).
 534         exvector v(it, itend);
 535         shaker_sort(v.begin(), v.end(), ex_is_less(), ex_swap());
 536
 537         // Find dummy pairs and free indices
 538         it = v.begin(); itend = v.end();
 539         exvector::const_iterator last = it++;
 540         while (it != itend) {
 541                 if (is_dummy_pair(*it, *last)) {
 542                         out_dummy.push_back(*last);
 543                         it++;
 544                         if (it == itend)
 545                                 return;
 546                 } else {
 547                         if (!it->is_equal(*last) && ex_to<idx>(*last).is_symbolic())
 548                                 out_free.push_back(*last);
 549                 }
 550                 last = it++;
 551         }
 552         if (ex_to<idx>(*last).is_symbolic())
 553                 out_free.push_back(*last);
 554 }
 555
 556 ex minimal_dim(const ex & dim1, const ex & dim2)
 557 {
 558         if (dim1.is_equal(dim2) || dim1 < dim2 || (is_exactly_a<numeric>(dim1) && is_a<symbol>(dim2)))
 559                 return dim1;
 560         else if (dim1 > dim2 || (is_a<symbol>(dim1) && is_exactly_a<numeric>(dim2)))
 561                 return dim2;
 562         else {
 563                 std::ostringstream s;
 564                 s << "minimal_dim(): index dimensions " << dim1 << " and " << dim2 << " cannot be ordered";
 565                 throw (std::runtime_error(s.str()));
 566         }
 567 }
 568
 569 } // namespace GiNaC