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