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
  25 #include "idx.h"
  26 #include "symbol.h"
  27 #include "lst.h"
  28 #include "print.h"
  29 #include "archive.h"
  30 #include "utils.h"
  31 #include "debugmsg.h"
  32
  33 namespace GiNaC {
  34
  35 GINAC_IMPLEMENT_REGISTERED_CLASS(idx, basic)
  36 GINAC_IMPLEMENT_REGISTERED_CLASS(varidx, idx)
  37
  38 //////////
  39 // default constructor, destructor, copy constructor assignment operator and helpers
  40 //////////
  41
  42 idx::idx() : inherited(TINFO_idx)
  43 {
  44         debugmsg("idx default constructor", LOGLEVEL_CONSTRUCT);
  45 }
  46
  47 varidx::varidx() : covariant(false)
  48 {
  49         debugmsg("varidx default constructor", LOGLEVEL_CONSTRUCT);
  50         tinfo_key = TINFO_varidx;
  51 }
  52
  53 void idx::copy(const idx & other)
  54 {
  55         inherited::copy(other);
  56         value = other.value;
  57         dim = other.dim;
  58 }
  59
  60 void varidx::copy(const varidx & other)
  61 {
  62         inherited::copy(other);
  63         covariant = other.covariant;
  64 }
  65
  66 DEFAULT_DESTROY(idx)
  67 DEFAULT_DESTROY(varidx)
  68
  69 //////////
  70 // other constructors
  71 //////////
  72
  73 idx::idx(const ex & v, const ex & d) : inherited(TINFO_idx), value(v), dim(d)
  74 {
  75         debugmsg("idx constructor from ex,ex", LOGLEVEL_CONSTRUCT);
  76         if (is_dim_numeric())
  77                 if (!dim.info(info_flags::posint))
  78                         throw(std::invalid_argument("dimension of space must be a positive integer"));
  79 }
  80
  81 varidx::varidx(const ex & v, const ex & d, bool cov) : inherited(v, d), covariant(cov)
  82 {
  83         debugmsg("varidx constructor from ex,ex,bool", LOGLEVEL_CONSTRUCT);
  84         tinfo_key = TINFO_varidx;
  85 }
  86
  87 //////////
  88 // archiving
  89 //////////
  90
  91 idx::idx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
  92 {
  93         debugmsg("idx constructor from archive_node", LOGLEVEL_CONSTRUCT);
  94         n.find_ex("value", value, sym_lst);
  95         n.find_ex("dim", dim, sym_lst);
  96 }
  97
  98 varidx::varidx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
  99 {
 100         debugmsg("varidx constructor from archive_node", LOGLEVEL_CONSTRUCT);
 101         n.find_bool("covariant", covariant);
 102 }
 103
 104 void idx::archive(archive_node &n) const
 105 {
 106         inherited::archive(n);
 107         n.add_ex("value", value);
 108         n.add_ex("dim", dim);
 109 }
 110
 111 void varidx::archive(archive_node &n) const
 112 {
 113         inherited::archive(n);
 114         n.add_bool("covariant", covariant);
 115 }
 116
 117 DEFAULT_UNARCHIVE(idx)
 118 DEFAULT_UNARCHIVE(varidx)
 119
 120 //////////
 121 // functions overriding virtual functions from bases classes
 122 //////////
 123
 124 void idx::print(const print_context & c, unsigned level) const
 125 {
 126         debugmsg("idx print", LOGLEVEL_PRINT);
 127
 128         if (is_of_type(c, print_tree)) {
 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_of_type(c, print_latex))
 140                         c.s << ".";
 141                 bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
 142                 if (need_parens)
 143                         c.s << "(";
 144                 value.print(c);
 145                 if (need_parens)
 146                         c.s << ")";
 147         }
 148 }
 149
 150 void varidx::print(const print_context & c, unsigned level) const
 151 {
 152         debugmsg("varidx print", LOGLEVEL_PRINT);
 153
 154         if (is_of_type(c, print_tree)) {
 155
 156                 c.s << std::string(level, ' ') << class_name()
 157                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
 158                     << (covariant ? ", covariant" : ", contravariant")
 159                     << std::endl;
 160                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
 161                 value.print(c, level + delta_indent);
 162                 dim.print(c, level + delta_indent);
 163
 164         } else {
 165
 166                 if (!is_of_type(c, print_latex)) {
 167                         if (covariant)
 168                                 c.s << ".";
 169                         else
 170                                 c.s << "~";
 171                 }
 172                 bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
 173                 if (need_parens)
 174                         c.s << "(";
 175                 value.print(c);
 176                 if (need_parens)
 177                         c.s << ")";
 178         }
 179 }
 180
 181 bool idx::info(unsigned inf) const
 182 {
 183         if (inf == info_flags::idx)
 184                 return true;
 185         return inherited::info(inf);
 186 }
 187
 188 unsigned idx::nops() const
 189 {
 190         // don't count the dimension as that is not really a sub-expression
 191         return 1;
 192 }
 193
 194 ex & idx::let_op(int i)
 195 {
 196         GINAC_ASSERT(i == 0);
 197         return value;
 198 }
 199
 200 /** Returns order relation between two indices of the same type. The order
 201  *  must be such that dummy indices lie next to each other. */
 202 int idx::compare_same_type(const basic & other) const
 203 {
 204         GINAC_ASSERT(is_of_type(other, idx));
 205         const idx &o = static_cast<const idx &>(other);
 206
 207         int cmpval = value.compare(o.value);
 208         if (cmpval)
 209                 return cmpval;
 210         return dim.compare(o.dim);
 211 }
 212
 213 int varidx::compare_same_type(const basic & other) const
 214 {
 215         GINAC_ASSERT(is_of_type(other, varidx));
 216         const varidx &o = static_cast<const varidx &>(other);
 217
 218         int cmpval = inherited::compare_same_type(other);
 219         if (cmpval)
 220                 return cmpval;
 221
 222         // Check variance last so dummy indices will end up next to each other
 223         if (covariant != o.covariant)
 224                 return covariant ? -1 : 1;
 225         return 0;
 226 }
 227
 228 ex idx::subs(const lst & ls, const lst & lr) const
 229 {
 230         GINAC_ASSERT(ls.nops() == lr.nops());
 231
 232         // First look for index substitutions
 233         for (unsigned i=0; i<ls.nops(); i++) {
 234                 if (is_equal(*(ls.op(i)).bp)) {
 235
 236                         // Substitution index->index
 237                         if (is_ex_of_type(lr.op(i), idx))
 238                                 return lr.op(i);
 239
 240                         // Otherwise substitute value
 241                         idx *i_copy = static_cast<idx *>(duplicate());
 242                         i_copy->value = lr.op(i);
 243                         i_copy->clearflag(status_flags::hash_calculated);
 244                         return i_copy->setflag(status_flags::dynallocated);
 245                 }
 246         }
 247
 248         // None, substitute objects in value (not in dimension)
 249         const ex &subsed_value = value.subs(ls, lr);
 250         if (are_ex_trivially_equal(value, subsed_value))
 251                 return *this;
 252
 253         idx *i_copy = static_cast<idx *>(duplicate());
 254         i_copy->value = subsed_value;
 255         i_copy->clearflag(status_flags::hash_calculated);
 256         return i_copy->setflag(status_flags::dynallocated);
 257 }
 258
 259 //////////
 260 // new virtual functions
 261 //////////
 262
 263 bool idx::is_dummy_pair_same_type(const basic & other) const
 264 {
 265         const idx &o = static_cast<const idx &>(other);
 266
 267         // Only pure symbols form dummy pairs, "2n+1" doesn't
 268         if (!is_ex_of_type(value, symbol))
 269                 return false;
 270
 271         // Value must be equal, of course
 272         if (!value.is_equal(o.value))
 273                 return false;
 274
 275         // Also the dimension
 276         return dim.is_equal(o.dim);
 277 }
 278
 279 bool varidx::is_dummy_pair_same_type(const basic & other) const
 280 {
 281         const varidx &o = static_cast<const varidx &>(other);
 282
 283         // Variance must be opposite
 284         if (covariant == o.covariant)
 285                 return false;
 286
 287         return inherited::is_dummy_pair_same_type(other);
 288 }
 289
 290 //////////
 291 // non-virtual functions
 292 //////////
 293
 294 ex varidx::toggle_variance(void) const
 295 {
 296         varidx *i_copy = static_cast<varidx *>(duplicate());
 297         i_copy->covariant = !i_copy->covariant;
 298         i_copy->clearflag(status_flags::hash_calculated);
 299         return i_copy->setflag(status_flags::dynallocated);
 300 }
 301
 302 //////////
 303 // global functions
 304 //////////
 305
 306 bool is_dummy_pair(const idx & i1, const idx & i2)
 307 {
 308         // The indices must be of exactly the same type
 309         if (i1.tinfo() != i2.tinfo())
 310                 return false;
 311
 312         // Same type, let the indices decide whether they are paired
 313         return i1.is_dummy_pair_same_type(i2);
 314 }
 315
 316 bool is_dummy_pair(const ex & e1, const ex & e2)
 317 {
 318         // The expressions must be indices
 319         if (!is_ex_of_type(e1, idx) || !is_ex_of_type(e2, idx))
 320                 return false;
 321
 322         return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
 323 }
 324
 325 /** Bring a vector of indices into a canonic order. Dummy indices will lie
 326  *  next to each other after the sorting. */
 327 static void sort_index_vector(exvector &v)
 328 {
 329         // Nothing to sort if less than 2 elements
 330         if (v.size() < 2)
 331                 return;
 332
 333         // Simple bubble sort algorithm should be sufficient for the small
 334         // number of indices expected
 335         exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1;
 336         while (it1 != next_to_last_idx) {
 337                 exvector::iterator it2 = it1 + 1;
 338                 while (it2 != itend) {
 339                         if (it1->compare(*it2) > 0)
 340                                 it1->swap(*it2);
 341                         it2++;
 342                 }
 343                 it1++;
 344         }
 345 }
 346
 347
 348 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
 349 {
 350         out_free.clear();
 351         out_dummy.clear();
 352
 353         // No indices? Then do nothing
 354         if (it == itend)
 355                 return;
 356
 357         // Only one index? Then it is a free one if it's not numeric
 358         if (itend - it == 1) {
 359                 if (ex_to_idx(*it).is_symbolic())
 360                         out_free.push_back(*it);
 361                 return;
 362         }
 363
 364         // Sort index vector. This will cause dummy indices come to lie next
 365         // to each other (because the sort order is defined to guarantee this).
 366         exvector v(it, itend);
 367         sort_index_vector(v);
 368
 369         // Find dummy pairs and free indices
 370         it = v.begin(); itend = v.end();
 371         exvector::const_iterator last = it++;
 372         while (it != itend) {
 373                 if (is_dummy_pair(*it, *last)) {
 374                         out_dummy.push_back(*last);
 375                         it++;
 376                         if (it == itend)
 377                                 return;
 378                 } else {
 379                         if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic())
 380                                 out_free.push_back(*last);
 381                 }
 382                 last = it++;
 383         }
 384         if (ex_to_idx(*last).is_symbolic())
 385                 out_free.push_back(*last);
 386 }
 387
 388 exvector index_set_difference(const exvector & set1, const exvector & set2)
 389 {
 390         exvector ret;
 391
 392         exvector::const_iterator ait = set1.begin(), aitend = set1.end();
 393         while (ait != aitend) {
 394                 exvector::const_iterator bit = set2.begin(), bitend = set2.end();
 395                 bool found = false;
 396                 while (bit != bitend) {
 397                         if (ait->is_equal(*bit)) {
 398                                 found = true;
 399                                 break;
 400                         }
 401                         bit++;
 402                 }
 403                 if (!found)
 404                         ret.push_back(*ait);
 405                 ait++;
 406         }
 407
 408         return ret;
 409 }
 410
 411 } // namespace GiNaC