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