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