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