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