X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=fb34df9502a861c35f32cd739f01c4410e605ed6;hp=4eeb2db13600980050b2d7f7c91159e9e0044511;hb=745d46b8492878030cc0decd1cb9bf2de2de03ff;hpb=0a1b35cf1e59c9e3aae33de8febaa1c8f4bbe630 diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 4eeb2db1..fb34df95 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -1,9 +1,9 @@ /** @file indexed.cpp * - * Implementation of GiNaC's index carrying objects. */ + * Implementation of GiNaC's indexed expressions. */ /* - * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -20,16 +20,19 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -#include +#include #include "indexed.h" -#include "ex.h" #include "idx.h" +#include "add.h" +#include "mul.h" +#include "ncmul.h" +#include "power.h" +#include "archive.h" +#include "utils.h" #include "debugmsg.h" -#ifndef NO_GINAC_NAMESPACE namespace GiNaC { -#endif // ndef NO_GINAC_NAMESPACE GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq) @@ -37,98 +40,119 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq) // default constructor, destructor, copy constructor assignment operator and helpers ////////// -// public +indexed::indexed() : symmetry(unknown) +{ + debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; +} -indexed::indexed() +void indexed::copy(const indexed & other) { - debugmsg("indexed default constructor",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; + inherited::copy(other); + symmetry = other.symmetry; } -indexed::~indexed() +void indexed::destroy(bool call_parent) { - debugmsg("indexed destructor",LOGLEVEL_DESTRUCT); - destroy(0); + if (call_parent) + inherited::destroy(call_parent); } -indexed::indexed(indexed const & other) +////////// +// other constructors +////////// + +indexed::indexed(const ex & b) : inherited(b), symmetry(unknown) { - debugmsg("indexed copy constructor",LOGLEVEL_CONSTRUCT); - copy (other); + debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -indexed const & indexed::operator=(indexed const & other) +indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symmetry(unknown) { - debugmsg("indexed operator=",LOGLEVEL_ASSIGNMENT); - if (this != &other) { - destroy(1); - copy(other); - } - return *this; + debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -// protected +indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(unknown) +{ + debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); +} -void indexed::copy(indexed const & other) +indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(unknown) { - inherited::copy(other); + debugmsg("indexed constructor from ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -void indexed::destroy(bool call_parent) +indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(unknown) { - if (call_parent) { - inherited::destroy(call_parent); - } + debugmsg("indexed constructor from ex,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -////////// -// other constructors -////////// +indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(symm) +{ + debugmsg("indexed constructor from ex,symmetry,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); +} -// public +indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(symm) +{ + debugmsg("indexed constructor from ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); +} -indexed::indexed(ex const & i1) : inherited(i1) +indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(symm) { - debugmsg("indexed constructor from ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,symmetry,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -indexed::indexed(ex const & i1, ex const & i2) : inherited(i1,i2) +indexed::indexed(const ex & b, const exvector & v) : inherited(b), symmetry(unknown) { - debugmsg("indexed constructor from ex,ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,exvector", LOGLEVEL_CONSTRUCT); + seq.insert(seq.end(), v.begin(), v.end()); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -indexed::indexed(ex const & i1, ex const & i2, ex const & i3) - : inherited(i1,i2,i3) +indexed::indexed(const ex & b, symmetry_type symm, const exvector & v) : inherited(b), symmetry(symm) { - debugmsg("indexed constructor from ex,ex,ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from ex,symmetry,exvector", LOGLEVEL_CONSTRUCT); + seq.insert(seq.end(), v.begin(), v.end()); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -indexed::indexed(ex const & i1, ex const & i2, ex const & i3, ex const & i4) - : inherited(i1,i2,i3,i4) +indexed::indexed(symmetry_type symm, const exprseq & es) : inherited(es), symmetry(symm) { - debugmsg("indexed constructor from ex,ex,ex,ex",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -indexed::indexed(exvector const & iv) : inherited(iv) +indexed::indexed(symmetry_type symm, const exvector & v, bool discardable) : inherited(v, discardable), symmetry(symm) { - debugmsg("indexed constructor from exvector",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } -indexed::indexed(exvector * ivp) : inherited(ivp) +indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(symm) { - debugmsg("indexed constructor from exvector *",LOGLEVEL_CONSTRUCT); - tinfo_key=TINFO_indexed; - GINAC_ASSERT(all_of_type_idx()); + debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_indexed; + assert_all_indices_of_type_idx(); } ////////// @@ -138,118 +162,218 @@ indexed::indexed(exvector * ivp) : inherited(ivp) /** Construct object from archive_node. */ indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) { - debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT); + debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT); + unsigned int symm; + if (!(n.find_unsigned("symmetry", symm))) + throw (std::runtime_error("unknown indexed symmetry type in archive")); } /** Unarchive the object. */ ex indexed::unarchive(const archive_node &n, const lst &sym_lst) { - return (new indexed(n, sym_lst))->setflag(status_flags::dynallocated); + return (new indexed(n, sym_lst))->setflag(status_flags::dynallocated); } /** Archive the object. */ void indexed::archive(archive_node &n) const { - inherited::archive(n); + inherited::archive(n); + n.add_unsigned("symmetry", symmetry); } ////////// // functions overriding virtual functions from bases classes ////////// -// public - -basic * indexed::duplicate() const +void indexed::printraw(std::ostream & os) const { - debugmsg("indexed duplicate",LOGLEVEL_DUPLICATE); - return new indexed(*this); -} + debugmsg("indexed printraw", LOGLEVEL_PRINT); + GINAC_ASSERT(seq.size() > 0); -void indexed::printraw(ostream & os) const -{ - debugmsg("indexed printraw",LOGLEVEL_PRINT); - os << "indexed(indices="; - printrawindices(os); - os << ",hash=" << hashvalue << ",flags=" << flags << ")"; + os << class_name() << "("; + seq[0].printraw(os); + os << ",indices="; + printrawindices(os); + os << ",hash=" << hashvalue << ",flags=" << flags << ")"; } -void indexed::printtree(ostream & os, unsigned indent) const +void indexed::printtree(std::ostream & os, unsigned indent) const { - debugmsg("indexed printtree",LOGLEVEL_PRINT); - os << string(indent,' ') << "indexed: " << seq.size() << " indices"; - os << ",hash=" << hashvalue << ",flags=" << flags << endl; - printtreeindices(os,indent); -} + debugmsg("indexed printtree", LOGLEVEL_PRINT); + GINAC_ASSERT(seq.size() > 0); -void indexed::print(ostream & os, unsigned upper_precedence) const -{ - debugmsg("indexed print",LOGLEVEL_PRINT); - os << "UNNAMEDINDEX"; - printindices(os); + os << std::string(indent, ' ') << class_name() << ", " << seq.size()-1 << " indices"; + os << ",hash=" << hashvalue << ",flags=" << flags << std::endl; + printtreeindices(os, indent); } -void indexed::printcsrc(ostream & os, unsigned type, - unsigned upper_precedence) const +void indexed::print(std::ostream & os, unsigned upper_precedence) const { - debugmsg("indexed print csrc",LOGLEVEL_PRINT); - print(os,upper_precedence); + debugmsg("indexed print", LOGLEVEL_PRINT); + GINAC_ASSERT(seq.size() > 0); + + const ex & base = seq[0]; + bool need_parens = is_ex_exactly_of_type(base, add) || is_ex_exactly_of_type(base, mul) + || is_ex_exactly_of_type(base, ncmul) || is_ex_exactly_of_type(base, power); + if (need_parens) + os << "("; + os << base; + if (need_parens) + os << ")"; + printindices(os); } bool indexed::info(unsigned inf) const { - if (inf==info_flags::indexed) return true; - if (inf==info_flags::has_indices) return seq.size()!=0; - return inherited::info(inf); + if (inf == info_flags::indexed) return true; + if (inf == info_flags::has_indices) return seq.size() > 1; + return inherited::info(inf); } -exvector indexed::get_indices(void) const +bool indexed::all_index_values_are(unsigned inf) const { - return seq; + // No indices? Then no property can be fulfilled + if (seq.size() < 2) + return false; - /* - idxvector filtered_indices; - filtered_indices.reserve(indices.size()); - for (idxvector::const_iterator cit=indices.begin(); cit!=indices.end(); ++cit) { - if ((*cit).get_type()==t) { - filtered_indices.push_back(*cit); - } - } - return filtered_indices; - */ + // Check all indices + exvector::const_iterator it = seq.begin() + 1, itend = seq.end(); + while (it != itend) { + GINAC_ASSERT(is_ex_of_type(*it, idx)); + if (!ex_to_idx(*it).get_value().info(inf)) + return false; + it++; + } + return true; } -// protected - -int indexed::compare_same_type(basic const & other) const +int indexed::compare_same_type(const basic & other) const { - GINAC_ASSERT(is_of_type(other,indexed)); - return inherited::compare_same_type(other); + GINAC_ASSERT(is_of_type(other, indexed)); + return inherited::compare_same_type(other); } -bool indexed::is_equal_same_type(basic const & other) const -{ - GINAC_ASSERT(is_of_type(other,indexed)); - return inherited::is_equal_same_type(other); -} +// The main difference between sort_index_vector() and canonicalize_indices() +// is that the latter takes the symmetry of the object into account. Once we +// implement mixed symmetries, canonicalize_indices() will only be able to +// reorder index pairs with known symmetry properties, while sort_index_vector() +// always sorts the whole vector. -unsigned indexed::return_type(void) const -{ - return return_types::noncommutative; -} - -unsigned indexed::return_type_tinfo(void) const -{ - return tinfo_key; +/** Bring a vector of indices into a canonic order (don't care about the + * symmetry of the objects carrying the indices). Dummy indices will lie + * next to each other after the sorting. + * + * @param v Index vector to be sorted */ +static void sort_index_vector(exvector &v) +{ + // Nothing to sort if less than 2 elements + if (v.size() < 2) + return; + + // Simple bubble sort algorithm should be sufficient for the small + // number of indices expected + exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1; + while (it1 != next_to_last_idx) { + exvector::iterator it2 = it1 + 1; + while (it2 != itend) { + if (it1->compare(*it2) > 0) + it1->swap(*it2); + it2++; + } + it1++; + } +} + +/** Bring a vector of indices into a canonic order. This operation only makes + * sense if the object carrying these indices is either symmetric or totally + * antisymmetric with respect to the indices. + * + * @param itbegin Start of index vector + * @param itend End of index vector + * @param antisymm Whether the object is antisymmetric + * @return the sign introduced by the reordering of the indices if the object + * is antisymmetric (or 0 if two equal indices are encountered). For + * symmetric objects, this is always +1. If the index vector was + * already in a canonic order this function returns INT_MAX. */ +static int canonicalize_indices(exvector::iterator itbegin, exvector::iterator itend, bool antisymm) +{ + bool something_changed = false; + int sig = 1; + + // Simple bubble sort algorithm should be sufficient for the small + // number of indices expected + exvector::iterator it1 = itbegin, next_to_last_idx = itend - 1; + while (it1 != next_to_last_idx) { + exvector::iterator it2 = it1 + 1; + while (it2 != itend) { + int cmpval = it1->compare(*it2); + if (cmpval == 1) { + it1->swap(*it2); + something_changed = true; + if (antisymm) + sig = -sig; + } else if (cmpval == 0 && antisymm) { + something_changed = true; + sig = 0; + } + it2++; + } + it1++; + } + + return something_changed ? sig : INT_MAX; +} + +ex indexed::eval(int level) const +{ + // First evaluate children, then we will end up here again + if (level > 1) + return indexed(symmetry, evalchildren(level)); + + // Canonicalize indices according to the symmetry properties + if (seq.size() > 2 && (symmetry != unknown && symmetry != mixed)) { + exvector v = seq; + int sig = canonicalize_indices(v.begin() + 1, v.end(), symmetry == antisymmetric); + if (sig != INT_MAX) { + // Something has changed while sorting indices, more evaluations later + if (sig == 0) + return _ex0(); + return ex(sig) * thisexprseq(v); + } + } + + // Let the class of the base object perform additional evaluations + return op(0).bp->eval_indexed(*this); +} + +ex indexed::thisexprseq(const exvector & v) const +{ + return indexed(symmetry, v); } -ex indexed::thisexprseq(exvector const & v) const +ex indexed::thisexprseq(exvector * vp) const { - return indexed(v); + return indexed(symmetry, vp); } -ex indexed::thisexprseq(exvector * vp) const +ex indexed::expand(unsigned options) const { - return indexed(vp); + GINAC_ASSERT(seq.size() > 0); + + if ((options & expand_options::expand_indexed) && is_ex_exactly_of_type(seq[0], add)) { + + // expand_indexed expands (a+b).i -> a.i + b.i + const ex & base = seq[0]; + ex sum = _ex0(); + for (unsigned i=0; i 1) { + exvector::const_iterator it=seq.begin() + 1, itend = seq.end(); + while (it != itend) { + it->printraw(os); + it++; + if (it != itend) + os << ","; + } + } } -void indexed::printtreeindices(ostream & os, unsigned indent) const +void indexed::printtreeindices(std::ostream & os, unsigned indent) const { - if (seq.size()!=0) { - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - os << string(indent+delta_indent,' '); - (*cit).printraw(os); - os << endl; - } - } + if (seq.size() > 1) { + exvector::const_iterator it=seq.begin() + 1, itend = seq.end(); + while (it != itend) { + os << std::string(indent + delta_indent, ' '); + it->printraw(os); + os << std::endl; + it++; + } + } } -void indexed::printindices(ostream & os) const +void indexed::printindices(std::ostream & os) const { - if (seq.size()!=0) { - if (seq.size()>1) { - os << "{"; - } - exvector::const_iterator last=seq.end()-1; - exvector::const_iterator cit=seq.begin(); - for (; cit!=last; ++cit) { - (*cit).print(os); - os << ","; - } - (*cit).print(os); - if (seq.size()>1) { - os << "}"; - } - } + if (seq.size() > 1) { + exvector::const_iterator it=seq.begin() + 1, itend = seq.end(); + while (it != itend) { + it->print(os); + it++; + } + } } -bool indexed::all_of_type_idx(void) const +/** Check whether all indices are of class idx. This function is used + * internally to make sure that all constructed indexed objects really + * carry indices and not some other classes. */ +void indexed::assert_all_indices_of_type_idx(void) const { - // used only inside of ASSERTs - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - if (!is_ex_of_type(*cit,idx)) return false; - } - return true; + GINAC_ASSERT(seq.size() > 0); + exvector::const_iterator it = seq.begin() + 1, itend = seq.end(); + while (it != itend) { + if (!is_ex_of_type(*it, idx)) + throw(std::invalid_argument("indices of indexed object must be of type idx")); + it++; + } } ////////// -// static member variables +// global functions ////////// -// none +/** Given a vector of indices, split them into two vectors, one containing + * the free indices, the other containing the dummy indices. */ +static void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy) +{ + out_free.clear(); + out_dummy.clear(); + + // No indices? Then do nothing + if (it == itend) + return; + + // Only one index? Then it is a free one if it's not numeric + if (itend - it == 1) { + if (ex_to_idx(*it).is_symbolic()) + out_free.push_back(*it); + return; + } + + // Sort index vector. This will cause dummy indices come to lie next + // to each other (because the sort order is defined to guarantee this). + exvector v(it, itend); + sort_index_vector(v); + + // Find dummy pairs and free indices + it = v.begin(); itend = v.end(); + exvector::const_iterator last = it++; + while (it != itend) { + if (is_dummy_pair(*it, *last)) { + out_dummy.push_back(*last); + it++; + if (it == itend) + return; + } else { + if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic()) + out_free.push_back(*last); + } + last = it++; + } + if (ex_to_idx(*last).is_symbolic()) + out_free.push_back(*last); +} + +/** Check whether two sorted index vectors are consistent (i.e. equal). */ +static bool indices_consistent(const exvector & v1, const exvector & v2) +{ + // Number of indices must be the same + if (v1.size() != v2.size()) + return false; + + // And also the indices themselves + exvector::const_iterator ait = v1.begin(), aitend = v1.end(), + bit = v2.begin(), bitend = v2.end(); + while (ait != aitend) { + if (!ait->is_equal(*bit)) + return false; + ait++; bit++; + } + return true; +} + +exvector indexed::get_dummy_indices(void) const +{ + exvector free_indices, dummy_indices; + find_free_and_dummy(seq.begin() + 1, seq.end(), free_indices, dummy_indices); + return dummy_indices; +} + +exvector indexed::get_free_indices(void) const +{ + exvector free_indices, dummy_indices; + find_free_and_dummy(seq.begin() + 1, seq.end(), free_indices, dummy_indices); + return free_indices; +} + +exvector add::get_free_indices(void) const +{ + exvector free_indices; + for (unsigned i=0; i a*a + GINAC_ASSERT(e.op(1).is_equal(_ex2())); + v.push_back(e.op(0)); + v.push_back(e.op(0)); + } else { + for (int i=0; i 1); + exvector::iterator it1, itend = v.end(), next_to_last = itend - 1; + for (it1 = v.begin(); it1 != next_to_last; it1++) { + +try_again: + if (!is_ex_of_type(*it1, indexed)) + continue; + + // Indexed factor found, look for contraction candidates + exvector::iterator it2; + for (it2 = it1 + 1; it2 != itend; it2++) { + + if (!is_ex_of_type(*it2, indexed)) + continue; + + // Check whether the two factors share dummy indices + exvector un(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end()); + un.insert(un.end(), ex_to_indexed(*it2).seq.begin() + 1, ex_to_indexed(*it2).seq.end()); + exvector free, dummy; + find_free_and_dummy(un.begin(), un.end(), free, dummy); + if (dummy.size() == 0) + continue; + + // At least one dummy index, is it a defined scalar product? + if (free.size() == 0) { + if (sp.is_defined(*it1, *it2)) { + *it1 = sp.evaluate(*it1, *it2); + *it2 = _ex1(); + something_changed = true; + goto try_again; + } + } + + // Try to contract the first one with the second one + bool contracted = it1->op(0).bp->contract_with(it1, it2, v); + if (!contracted) { + + // That didn't work; maybe the second object knows how to + // contract itself with the first one + contracted = it2->op(0).bp->contract_with(it2, it1, v); + } + if (contracted) { + something_changed = true; + + // Both objects may have new indices now or they might + // even not be indexed objects any more, so we have to + // start over + goto try_again; + } + } + } + + // Find free indices (concatenate them all and call find_free_and_dummy()) + exvector un, dummy_indices; + it1 = v.begin(); itend = v.end(); + while (it1 != itend) { + if (is_ex_of_type(*it1, indexed)) { + const indexed & o = ex_to_indexed(*it1); + un.insert(un.end(), o.seq.begin() + 1, o.seq.end()); + } + it1++; + } + find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices); + + if (something_changed) { + if (non_commutative) + return ncmul(v); + else + return mul(v); + } else + return e; +} + +/** Simplify indexed expression, return list of free indices. */ +ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products & sp) +{ + // Expand the expression + ex e_expanded = e.expand(); + + // Simplification of single indexed object: just find the free indices + if (is_ex_of_type(e_expanded, indexed)) { + const indexed &i = ex_to_indexed(e_expanded); + exvector dummy_indices; + find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, dummy_indices); + return e_expanded; + } + + // Simplification of sum = sum of simplifications, check consistency of + // free indices in each term + if (is_ex_exactly_of_type(e_expanded, add)) { + ex sum = _ex0(); + + for (unsigned i=0; isecond; +} + +void scalar_products::debugprint(void) const +{ + std::cerr << "map size=" << spm.size() << std::endl; + for (spmap::const_iterator cit=spm.begin(); cit!=spm.end(); ++cit) { + const spmapkey & k = cit->first; + std::cerr << "item key=(" << k.first << "," << k.second; + std::cerr << "), value=" << cit->second << std::endl; + } +} + +/** Make key from object pair. */ +spmapkey scalar_products::make_key(const ex & v1, const ex & v2) +{ + // If indexed, extract base objects + ex s1 = is_ex_of_type(v1, indexed) ? v1.op(0) : v1; + ex s2 = is_ex_of_type(v2, indexed) ? v2.op(0) : v2; + + // Enforce canonical order in pair + if (s1.compare(s2) > 0) + return spmapkey(s2, s1); + else + return spmapkey(s1, s2); +} -#ifndef NO_GINAC_NAMESPACE } // namespace GiNaC -#endif // ndef NO_GINAC_NAMESPACE