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 constructor
  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 //////////
  59 // other constructors
  60 //////////
  61
  62 idx::idx(const ex & v, const ex & d) : inherited(TINFO_idx), value(v), dim(d)
  63 {
  64         if (is_dim_numeric())
  65                 if (!dim.info(info_flags::posint))
  66                         throw(std::invalid_argument("dimension of space must be a positive integer"));
  67 }
  68
  69 varidx::varidx(const ex & v, const ex & d, bool cov) : inherited(v, d), covariant(cov)
  70 {
  71         tinfo_key = TINFO_varidx;
  72 }
  73
  74 spinidx::spinidx(const ex & v, const ex & d, bool cov, bool dot) : inherited(v, d, cov), dotted(dot)
  75 {
  76         tinfo_key = TINFO_spinidx;
  77 }
  78
  79 //////////
  80 // archiving
  81 //////////
  82
  83 idx::idx(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
  84 {
  85         n.find_ex("value", value, sym_lst);
  86         n.find_ex("dim", dim, sym_lst);
  87 }
  88
  89 varidx::varidx(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
  90 {
  91         n.find_bool("covariant", covariant);
  92 }
  93
  94 spinidx::spinidx(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
  95 {
  96         n.find_bool("dotted", dotted);
  97 }
  98
  99 void idx::archive(archive_node &n) const
 100 {
 101         inherited::archive(n);
 102         n.add_ex("value", value);
 103         n.add_ex("dim", dim);
 104 }
 105
 106 void varidx::archive(archive_node &n) const
 107 {
 108         inherited::archive(n);
 109         n.add_bool("covariant", covariant);
 110 }
 111
 112 void spinidx::archive(archive_node &n) const
 113 {
 114         inherited::archive(n);
 115         n.add_bool("dotted", dotted);
 116 }
 117
 118 DEFAULT_UNARCHIVE(idx)
 119 DEFAULT_UNARCHIVE(varidx)
 120 DEFAULT_UNARCHIVE(spinidx)
 121
 122 //////////
 123 // functions overriding virtual functions from base classes
 124 //////////
 125
 126 void idx::print(const print_context & c, unsigned level) const
 127 {
 128         if (is_a<print_tree>(c)) {
 129
 130                 c.s << std::string(level, ' ') << class_name()
 131                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 132                     << std::endl;
 133                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
 134                 value.print(c, level + delta_indent);
 135                 dim.print(c, level + delta_indent);
 136
 137         } else {
 138
 139                 if (is_a<print_latex>(c))
 140                         c.s << "{";
 141                 else
 142                         c.s << ".";
 143                 bool need_parens = !(is_exactly_a<numeric>(value) || is_a<symbol>(value));
 144                 if (need_parens)
 145                         c.s << "(";
 146                 value.print(c);
 147                 if (need_parens)
 148                         c.s << ")";
 149                 if (c.options & print_options::print_index_dimensions) {
 150                         c.s << "[";
 151                         dim.print(c);
 152                         c.s << "]";
 153                 }
 154                 if (is_a<print_latex>(c))
 155                         c.s << "}";
 156         }
 157 }
 158
 159 void varidx::print(const print_context & c, unsigned level) const
 160 {
 161         if (is_a<print_tree>(c)) {
 162
 163                 c.s << std::string(level, ' ') << class_name()
 164                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 165                     << (covariant ? ", covariant" : ", contravariant")
 166                     << std::endl;
 167                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
 168                 value.print(c, level + delta_indent);
 169                 dim.print(c, level + delta_indent);
 170
 171         } else {
 172                 if (is_a<print_latex>(c))
 173                         c.s << "{";
 174                 else {
 175                         if (covariant)
 176                                 c.s << ".";
 177                         else
 178                                 c.s << "~";
 179                 }
 180                 bool need_parens = !(is_exactly_a<numeric>(value) || is_a<symbol>(value));
 181                 if (need_parens)
 182                         c.s << "(";
 183                 value.print(c);
 184                 if (need_parens)
 185                         c.s << ")";
 186                 if (c.options & print_options::print_index_dimensions) {
 187                         c.s << "[";
 188                         dim.print(c);
 189                         c.s << "]";
 190                 }
 191                 if (is_a<print_latex>(c))
 192                         c.s << "}";
 193         }
 194 }
 195
 196 void spinidx::print(const print_context & c, unsigned level) const
 197 {
 198         if (is_a<print_tree>(c)) {
 199
 200                 c.s << std::string(level, ' ') << class_name()
 201                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 202                     << (covariant ? ", covariant" : ", contravariant")
 203                     << (dotted ? ", dotted" : ", undotted")
 204                     << std::endl;
 205                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
 206                 value.print(c, level + delta_indent);
 207                 dim.print(c, level + delta_indent);
 208
 209         } else {
 210
 211                 bool is_tex = is_a<print_latex>(c);
 212                 if (is_tex) {
 213                         if (covariant)
 214                                 c.s << "_{";
 215                         else
 216                                 c.s << "^{";
 217                 } else {
 218                         if (covariant)
 219                                 c.s << ".";
 220                         else
 221                                 c.s << "~";
 222                 }
 223                 if (dotted) {
 224                         if (is_tex)
 225                                 c.s << "\\dot{";
 226                         else
 227                                 c.s << "*";
 228                 }
 229                 bool need_parens = !(is_exactly_a<numeric>(value) || is_a<symbol>(value));
 230                 if (need_parens)
 231                         c.s << "(";
 232                 value.print(c);
 233                 if (need_parens)
 234                         c.s << ")";
 235                 if (is_tex && dotted)
 236                         c.s << "}";
 237                 if (is_tex)
 238                         c.s << "}";
 239         }
 240 }
 241
 242 bool idx::info(unsigned inf) const
 243 {
 244         if (inf == info_flags::idx)
 245                 return true;
 246         return inherited::info(inf);
 247 }
 248
 249 size_t idx::nops() const
 250 {
 251         // don't count the dimension as that is not really a sub-expression
 252         return 1;
 253 }
 254
 255 ex idx::op(size_t i) const
 256 {
 257         GINAC_ASSERT(i == 0);
 258         return value;
 259 }
 260
 261 ex idx::map(map_function & f) const
 262 {
 263         idx *copy = duplicate();
 264         copy->setflag(status_flags::dynallocated);
 265         copy->clearflag(status_flags::hash_calculated);
 266         copy->value = f(value);
 267         return *copy;
 268 }
 269
 270 /** Returns order relation between two indices of the same type. The order
 271  *  must be such that dummy indices lie next to each other. */
 272 int idx::compare_same_type(const basic & other) const
 273 {
 274         GINAC_ASSERT(is_a<idx>(other));
 275         const idx &o = static_cast<const idx &>(other);
 276
 277         int cmpval = value.compare(o.value);
 278         if (cmpval)
 279                 return cmpval;
 280         return dim.compare(o.dim);
 281 }
 282
 283 bool idx::match_same_type(const basic & other) const
 284 {
 285         GINAC_ASSERT(is_a<idx>(other));
 286         const idx &o = static_cast<const idx &>(other);
 287
 288         return dim.is_equal(o.dim);
 289 }
 290
 291 int varidx::compare_same_type(const basic & other) const
 292 {
 293         GINAC_ASSERT(is_a<varidx>(other));
 294         const varidx &o = static_cast<const varidx &>(other);
 295
 296         int cmpval = inherited::compare_same_type(other);
 297         if (cmpval)
 298                 return cmpval;
 299
 300         // Check variance last so dummy indices will end up next to each other
 301         if (covariant != o.covariant)
 302                 return covariant ? -1 : 1;
 303
 304         return 0;
 305 }
 306
 307 bool varidx::match_same_type(const basic & other) const
 308 {
 309         GINAC_ASSERT(is_a<varidx>(other));
 310         const varidx &o = static_cast<const varidx &>(other);
 311
 312         if (covariant != o.covariant)
 313                 return false;
 314
 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 unsigned idx::calchash() const
 345 {
 346         unsigned v = golden_ratio_hash(tinfo());
 347         v = rotate_left(v);
 348         v ^= value.gethash();
 349         v = rotate_left(v);
 350         v ^= dim.gethash();
 351
 352         // Store calculated hash value only if object is already evaluated
 353         if (flags & status_flags::evaluated) {
 354                 setflag(status_flags::hash_calculated);
 355                 hashvalue = v;
 356         }
 357
 358         return v;
 359 }
 360
 361 /** By default, basic::evalf would evaluate the index value but we don't want
 362  *  a.1 to become a.(1.0). */
 363 ex idx::evalf(int level) const
 364 {
 365         return *this;
 366 }
 367
 368 ex idx::subs(const exmap & m, unsigned options) const
 369 {
 370         // First look for index substitutions
 371         exmap::const_iterator it = m.find(*this);
 372         if (it != m.end()) {
 373
 374                 // Substitution index->index
 375                 if (is_a<idx>(it->second))
 376                         return it->second;
 377
 378                 // Otherwise substitute value
 379                 idx *i_copy = duplicate();
 380                 i_copy->value = it->second;
 381                 i_copy->clearflag(status_flags::hash_calculated);
 382                 return i_copy->setflag(status_flags::dynallocated);
 383         }
 384
 385         // None, substitute objects in value (not in dimension)
 386         const ex &subsed_value = value.subs(m, options);
 387         if (are_ex_trivially_equal(value, subsed_value))
 388                 return *this;
 389
 390         idx *i_copy = 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 = 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() const
 469 {
 470         varidx *i_copy = 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() const
 477 {
 478         spinidx *i_copy = 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() const
 485 {
 486         spinidx *i_copy = 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