X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=2835a904ec767d1208596266bc2f0b15ba97dd4c;hp=ef904f7222bfc798420dee5cdb6bed6373426a26;hb=25b4d894b5b2f7f56fb49e62a7767ab969d5249f;hpb=0633ed8082961673eedc092689e06fa39d6bc322 diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index ef904f72..2835a904 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's indexed expressions. */ /* - * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2016 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,10 +20,6 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#include -#include -#include - #include "indexed.h" #include "idx.h" #include "add.h" @@ -39,6 +35,12 @@ #include "utils.h" #include "integral.h" #include "matrix.h" +#include "inifcns.h" + +#include +#include +#include +#include namespace GiNaC { @@ -53,96 +55,83 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(indexed, exprseq, indexed::indexed() : symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; } ////////// // other constructors ////////// -indexed::indexed(const ex & b) : inherited(b), symtree(not_symmetric()) +indexed::indexed(const ex & b) : inherited{b}, symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(not_symmetric()) +indexed::indexed(const ex & b, const ex & i1) : inherited{b, i1}, symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(not_symmetric()) +indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited{b, i1, i2}, symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(not_symmetric()) +indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited{b, i1, i2, i3}, symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(not_symmetric()) +indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited{b, i1, i2, i3, i4}, symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(symm) +indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited{b, i1, i2}, symtree(symm) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(symm) +indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited{b, i1, i2, i3}, symtree(symm) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(symm) +indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited{b, i1, i2, i3, i4}, symtree(symm) { - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(not_symmetric()) +indexed::indexed(const ex & b, const exvector & v) : inherited{b}, symtree(not_symmetric()) { seq.insert(seq.end(), v.begin(), v.end()); - tinfo_key = TINFO_indexed; validate(); } -indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm) +indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited{b}, symtree(symm) { seq.insert(seq.end(), v.begin(), v.end()); - tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm) { - tinfo_key = TINFO_indexed; } -indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm) +indexed::indexed(const symmetry & symm, const exvector & v) : inherited(v), symtree(symm) { - tinfo_key = TINFO_indexed; } -indexed::indexed(const symmetry & symm, std::auto_ptr vp) : inherited(vp), symtree(symm) +indexed::indexed(const symmetry & symm, exvector && v) : inherited(std::move(v)), symtree(symm) { - tinfo_key = TINFO_indexed; } ////////// // archiving ////////// -indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) +void indexed::read_archive(const archive_node &n, lst &sym_lst) { + inherited::read_archive(n, sym_lst); if (!n.find_ex("symmetry", symtree, sym_lst)) { // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property unsigned symm = 0; @@ -161,6 +150,7 @@ indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) const_cast(ex_to(symtree)).validate(seq.size() - 1); } } +GINAC_BIND_UNARCHIVER(indexed); void indexed::archive(archive_node &n) const { @@ -168,8 +158,6 @@ void indexed::archive(archive_node &n) const n.add_ex("symmetry", symtree); } -DEFAULT_UNARCHIVE(indexed) - ////////// // functions overriding virtual functions from base classes ////////// @@ -178,7 +166,7 @@ void indexed::printindices(const print_context & c, unsigned level) const { if (seq.size() > 1) { - exvector::const_iterator it=seq.begin() + 1, itend = seq.end(); + auto it = seq.begin() + 1, itend = seq.end(); if (is_a(c)) { @@ -255,12 +243,6 @@ bool indexed::info(unsigned inf) const return inherited::info(inf); } -struct idx_is_not : public std::binary_function { - bool operator() (const ex & e, unsigned inf) const { - return !(ex_to(e).get_value().info(inf)); - } -}; - bool indexed::all_index_values_are(unsigned inf) const { // No indices? Then no property can be fulfilled @@ -268,7 +250,8 @@ bool indexed::all_index_values_are(unsigned inf) const return false; // Check all indices - return find_if(seq.begin() + 1, seq.end(), bind2nd(idx_is_not(), inf)) == seq.end(); + return find_if(seq.begin() + 1, seq.end(), + [inf](const ex & e) { return !(ex_to(e).get_value().info(inf)); }) == seq.end(); } int indexed::compare_same_type(const basic & other) const @@ -277,12 +260,8 @@ int indexed::compare_same_type(const basic & other) const return inherited::compare_same_type(other); } -ex indexed::eval(int level) const +ex indexed::eval() const { - // First evaluate children, then we will end up here again - if (level > 1) - return indexed(ex_to(symtree), evalchildren(level)); - const ex &base = seq[0]; // If the base object is 0, the whole object is 0 @@ -297,12 +276,15 @@ ex indexed::eval(int level) const return f * thiscontainer(v); } + if((typeid(*this) == typeid(indexed)) && seq.size()==1) + return base; + // Canonicalize indices according to the symmetry properties if (seq.size() > 2) { exvector v = seq; GINAC_ASSERT(is_exactly_a(symtree)); int sig = canonicalize(v.begin() + 1, ex_to(symtree)); - if (sig != INT_MAX) { + if (sig != std::numeric_limits::max()) { // Something has changed while sorting indices, more evaluations later if (sig == 0) return _ex0; @@ -314,14 +296,36 @@ ex indexed::eval(int level) const return ex_to(base).eval_indexed(*this); } +ex indexed::real_part() const +{ + if(op(0).info(info_flags::real)) + return *this; + return real_part_function(*this).hold(); +} + +ex indexed::imag_part() const +{ + if(op(0).info(info_flags::real)) + return 0; + return imag_part_function(*this).hold(); +} + ex indexed::thiscontainer(const exvector & v) const { return indexed(ex_to(symtree), v); } -ex indexed::thiscontainer(std::auto_ptr vp) const +ex indexed::thiscontainer(exvector && v) const { - return indexed(ex_to(symtree), vp); + return indexed(ex_to(symtree), std::move(v)); +} + +unsigned indexed::return_type() const +{ + if(is_a(op(0))) + return return_types::commutative; + else + return op(0).return_type(); } ex indexed::expand(unsigned options) const @@ -364,7 +368,7 @@ ex indexed::expand(unsigned options) const void indexed::validate() const { GINAC_ASSERT(seq.size() > 0); - exvector::const_iterator it = seq.begin() + 1, itend = seq.end(); + auto it = seq.begin() + 1, itend = seq.end(); while (it != itend) { if (!is_a(*it)) throw(std::invalid_argument("indices of indexed object must be of type idx")); @@ -441,7 +445,7 @@ exvector indexed::get_dummy_indices(const indexed & other) const bool indexed::has_dummy_index_for(const ex & i) const { - exvector::const_iterator it = seq.begin() + 1, itend = seq.end(); + auto it = seq.begin() + 1, itend = seq.end(); while (it != itend) { if (is_dummy_pair(*it, i)) return true; @@ -509,22 +513,6 @@ struct is_summation_idx : public std::unary_function { } }; -exvector power::get_free_indices() const -{ - // Get free indices of basis - exvector basis_indices = basis.get_free_indices(); - - if (exponent.info(info_flags::even)) { - // If the exponent is an even number, then any "free" index that - // forms a dummy pair with itself is actually a summation index - exvector really_free; - std::remove_copy_if(basis_indices.begin(), basis_indices.end(), - std::back_inserter(really_free), is_summation_idx()); - return really_free; - } else - return basis_indices; -} - exvector integral::get_free_indices() const { if (a.get_free_indices().size() || b.get_free_indices().size()) @@ -532,6 +520,15 @@ exvector integral::get_free_indices() const return f.get_free_indices(); } +template size_t number_of_type(const exvector&v) +{ + size_t number = 0; + for (auto & it : v) + if (is_exactly_a(it)) + ++number; + return number; +} + /** Rename dummy indices in an expression. * * @param e Expression to work on @@ -540,10 +537,10 @@ exvector integral::get_free_indices() const * @param global_dummy_indices The set of dummy indices that have appeared * before and which we would like to use in "e", too. This gets updated * by the function */ -static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices) +template static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices) { - size_t global_size = global_dummy_indices.size(), - local_size = local_dummy_indices.size(); + size_t global_size = number_of_type(global_dummy_indices), + local_size = number_of_type(local_dummy_indices); // Any local dummy indices at all? if (local_size == 0) @@ -555,9 +552,9 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex // to the global set size_t old_global_size = global_size; int remaining = local_size - global_size; - exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end(); + auto it = local_dummy_indices.begin(), itend = local_dummy_indices.end(); while (it != itend && remaining > 0) { - if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(op0_is_equal(), *it)) == global_dummy_indices.end()) { + if (is_exactly_a(*it) && find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(idx_is_equal_ignore_dim(), *it)) == global_dummy_indices.end()) { global_dummy_indices.push_back(*it); global_size++; remaining--; @@ -575,11 +572,13 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex exvector local_syms, global_syms; local_syms.reserve(local_size); global_syms.reserve(local_size); - for (size_t i=0; i(local_dummy_indices[i])) + local_syms.push_back(local_dummy_indices[i].op(0)); shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap()); - for (size_t i=0; i(global_dummy_indices[i])) + global_syms.push_back(global_dummy_indices[i].op(0)); shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap()); // Remove common indices @@ -618,10 +617,58 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector { bool something_changed = false; + // Find dummy symbols that occur twice in the same indexed object. + exvector local_var_dummies; + local_var_dummies.reserve(e.nops()/2); + for (size_t i=1; i(e.op(i))) + continue; + for (size_t j=i+1; jop(0)) { + variant_dummy_indices.erase(k); + break; + } + } + break; + } + } + } + + // In the case where a dummy symbol occurs twice in the same indexed object + // we try all possibilities of raising/lowering and keep the least one in + // the sense of ex_is_less. + ex optimal_e = e; + size_t numpossibs = 1 << local_var_dummies.size(); + for (size_t i=0; i(curr_idx).toggle_variance(); + m[curr_idx] = curr_toggle; + m[curr_toggle] = curr_idx; + } + try_e = e.subs(m, subs_options::no_pattern); + } + if(ex_is_less()(try_e, optimal_e)) + { optimal_e = try_e; + something_changed = true; + } + } + e = optimal_e; + + if (!is_a(e)) + return true; + + exvector seq = ex_to(e).seq; + // If a dummy index is encountered for the first time in the // product, pull it up, otherwise, pull it down - exvector::const_iterator it2, it2start, it2end; - for (it2start = ex_to(e).seq.begin(), it2end = ex_to(e).seq.end(), it2 = it2start + 1; it2 != it2end; ++it2) { + for (auto it2 = seq.begin()+1, it2end = seq.end(); it2 != it2end; ++it2) { if (!is_exactly_a(*it2)) continue; @@ -629,14 +676,20 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector for (vit = variant_dummy_indices.begin(), vitend = variant_dummy_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to(*it2).is_covariant()) { - e = e.subs(lst( - *it2 == ex_to(*it2).toggle_variance(), - ex_to(*it2).toggle_variance() == *it2 - ), subs_options::no_pattern); + /* + * N.B. we don't want to use + * + * e = e.subs(lst{ + * *it2 == ex_to(*it2).toggle_variance(), + * ex_to(*it2).toggle_variance() == *it2 + * }, subs_options::no_pattern); + * + * since this can trigger non-trivial repositioning of indices, + * e.g. due to non-trivial symmetry properties of e, thus + * invalidating iterators + */ + *it2 = ex_to(*it2).toggle_variance(); something_changed = true; - it2 = ex_to(e).seq.begin() + (it2 - it2start); - it2start = ex_to(e).seq.begin(); - it2end = ex_to(e).seq.end(); } moved_indices.push_back(*vit); variant_dummy_indices.erase(vit); @@ -647,11 +700,8 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector for (vit = moved_indices.begin(), vitend = moved_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to(*it2).is_contravariant()) { - e = e.subs(*it2 == ex_to(*it2).toggle_variance(), subs_options::no_pattern); + *it2 = ex_to(*it2).toggle_variance(); something_changed = true; - it2 = ex_to(e).seq.begin() + (it2 - it2start); - it2start = ex_to(e).seq.begin(); - it2end = ex_to(e).seq.end(); } goto next_index; } @@ -660,6 +710,9 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector next_index: ; } + if (something_changed) + e = ex_to(e).thiscontainer(seq); + return something_changed; } @@ -704,6 +757,21 @@ static void product_to_exvector(const ex & e, exvector & v, bool & non_commutati } } +template ex idx_symmetrization(const ex& r,const exvector& local_dummy_indices) +{ exvector dummy_syms; + dummy_syms.reserve(r.nops()); + for (auto & it : local_dummy_indices) + if(is_exactly_a(it)) + dummy_syms.push_back(it.op(0)); + if(dummy_syms.size() < 2) + return r; + ex q=symmetrize(r, dummy_syms); + return q; +} + +// Forward declaration needed in absence of friend injection, C.f. [namespace.memdef]: +ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp); + /** Simplify product of indexed expressions (commutative, noncommutative and * simple squares), return list of free indices. */ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp) @@ -718,6 +786,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du // Perform contractions bool something_changed = false; + bool has_nonsymmetric = false; GINAC_ASSERT(v.size() > 1); exvector::iterator it1, itend = v.end(), next_to_last = itend - 1; for (it1 = v.begin(); it1 != next_to_last; it1++) { @@ -727,6 +796,7 @@ try_again: continue; bool first_noncommutative = (it1->return_type() != return_types::commutative); + bool first_nonsymmetric = ex_to(ex_to(*it1).get_symmetry()).has_nonsymmetric(); // Indexed factor found, get free indices and look for contraction // candidates @@ -756,20 +826,12 @@ try_again: // At least one dummy index, is it a defined scalar product? bool contracted = false; - if (free.empty()) { - - // Find minimal dimension of all indices of both factors - exvector::const_iterator dit = ex_to(*it1).seq.begin() + 1, ditend = ex_to(*it1).seq.end(); - ex dim = ex_to(*dit).get_dim(); - ++dit; - for (; dit != ditend; ++dit) { - dim = minimal_dim(dim, ex_to(*dit).get_dim()); - } - dit = ex_to(*it2).seq.begin() + 1; - ditend = ex_to(*it2).seq.end(); - for (; dit != ditend; ++dit) { - dim = minimal_dim(dim, ex_to(*dit).get_dim()); - } + if (free.empty() && it1->nops()==2 && it2->nops()==2) { + + ex dim = minimal_dim( + ex_to(it1->op(1)).get_dim(), + ex_to(it2->op(1)).get_dim() + ); // User-defined scalar product? if (sp.is_defined(*it1, *it2, dim)) { @@ -801,8 +863,16 @@ contraction_done: // Non-commutative products are always re-expanded to give // eval_ncmul() the chance to re-order and canonicalize // the product - ex r = (non_commutative ? ex(ncmul(v, true)) : ex(mul(v))); - return simplify_indexed(r, free_indices, dummy_indices, sp); + bool is_a_product = (is_exactly_a(*it1) || is_exactly_a(*it1)) && + (is_exactly_a(*it2) || is_exactly_a(*it2)); + ex r = (non_commutative ? ex(ncmul(std::move(v))) : ex(mul(std::move(v)))); + + // If new expression is a product we can call this function again, + // otherwise we need to pass argument to simplify_indexed() to be expanded + if (is_a_product) + return simplify_indexed_product(r, free_indices, dummy_indices, sp); + else + return simplify_indexed(r, free_indices, dummy_indices, sp); } // Both objects may have new indices now or they might @@ -811,6 +881,11 @@ contraction_done: something_changed = true; goto try_again; } + else if (!has_nonsymmetric && + (first_nonsymmetric || + ex_to(ex_to(*it2).get_symmetry()).has_nonsymmetric())) { + has_nonsymmetric = true; + } } } @@ -857,26 +932,35 @@ contraction_done: ex r; if (something_changed) - r = non_commutative ? ex(ncmul(v, true)) : ex(mul(v)); + r = non_commutative ? ex(ncmul(std::move(v))) : ex(mul(std::move(v))); else r = e; // The result should be symmetric with respect to exchange of dummy // indices, so if the symmetrization vanishes, the whole expression is // zero. This detects things like eps.i.j.k * p.j * p.k = 0. - if (local_dummy_indices.size() >= 2) { - exvector dummy_syms; - dummy_syms.reserve(local_dummy_indices.size()); - for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it) - dummy_syms.push_back(it->op(0)); - if (symmetrize(r, dummy_syms).is_zero()) { + if (has_nonsymmetric) { + ex q = idx_symmetrization(r, local_dummy_indices); + if (q.is_zero()) { + free_indices.clear(); + return _ex0; + } + q = idx_symmetrization(q, local_dummy_indices); + if (q.is_zero()) { + free_indices.clear(); + return _ex0; + } + q = idx_symmetrization(q, local_dummy_indices); + if (q.is_zero()) { free_indices.clear(); return _ex0; } } // Dummy index renaming - r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); + r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); + r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); + r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); // Product of indexed object with a scalar? if (is_exactly_a(r) && r.nops() == 2 @@ -893,7 +977,7 @@ public: terminfo(const ex & orig_, const ex & symm_) : orig(orig_), symm(symm_) {} ex orig; /**< original term */ - ex symm; /**< symmtrized term */ + ex symm; /**< symmetrized term */ }; class terminfo_is_less { @@ -943,6 +1027,17 @@ public: } }; +bool hasindex(const ex &x, const ex &sym) +{ + if(is_a(x) && x.op(0)==sym) + return true; + else + for(size_t i=0; i(e_expanded, dummy_indices, local_dummy_indices); + e_expanded = rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices); + e_expanded = rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices); + return e_expanded; } // Simplification of sum = sum of simplifications, check consistency of @@ -1015,18 +1113,19 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi if (num_terms_orig < 2 || dummy_indices.size() < 2) return sum; - // Yes, construct vector of all dummy index symbols - exvector dummy_syms; - dummy_syms.reserve(dummy_indices.size()); - for (exvector::const_iterator it = dummy_indices.begin(); it != dummy_indices.end(); ++it) - dummy_syms.push_back(it->op(0)); - // Chop the sum into terms and symmetrize each one over the dummy // indices std::vector terms; for (size_t i=0; i(term, dummy_indices_of_term); + term_symm = idx_symmetrization(term_symm, dummy_indices_of_term); + term_symm = idx_symmetrization(term_symm, dummy_indices_of_term); if (term_symm.is_zero()) continue; terms.push_back(terminfo(term, term_symm)); @@ -1039,7 +1138,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi std::vector terms_pass2; for (std::vector::const_iterator i=terms.begin(); i!=terms.end(); ) { size_t num = 1; - std::vector::const_iterator j = i + 1; + auto j = i + 1; while (j != terms.end() && j->symm == i->symm) { num++; j++; @@ -1054,13 +1153,13 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi // Chop the symmetrized terms into subterms std::vector sy; - for (std::vector::const_iterator i=terms_pass2.begin(); i!=terms_pass2.end(); ++i) { - if (is_exactly_a(i->symm)) { - size_t num = i->symm.nops(); + for (auto & i : terms_pass2) { + if (is_exactly_a(i.symm)) { + size_t num = i.symm.nops(); for (size_t j=0; jsymm.op(j), i->orig, num)); + sy.push_back(symminfo(i.symm.op(j), i.orig, num)); } else - sy.push_back(symminfo(i->symm, i->orig, 1)); + sy.push_back(symminfo(i.symm, i.orig, 1)); } // Sort by symmetrized subterms @@ -1069,10 +1168,10 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi // Combine equal symmetrized subterms std::vector sy_pass2; exvector result; - for (std::vector::const_iterator i=sy.begin(); i!=sy.end(); ) { + for (auto i=sy.begin(); i!=sy.end(); ) { // Combine equal terms - std::vector::const_iterator j = i + 1; + auto j = i + 1; if (j != sy.end() && j->symmterm == i->symmterm) { // More than one term, collect the coefficients @@ -1105,7 +1204,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi // How many symmetrized terms of this original term are left? size_t num = 1; - std::vector::const_iterator j = i + 1; + auto j = i + 1; while (j != sy_pass2.end() && j->orig == i->orig) { num++; j++; @@ -1129,7 +1228,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi } // Add all resulting terms - ex sum_symm = (new add(result))->setflag(status_flags::dynallocated); + ex sum_symm = dynallocate(result); if (sum_symm.is_zero()) free_indices.clear(); return sum_symm; @@ -1257,9 +1356,9 @@ void scalar_products::add(const ex & v1, const ex & v2, const ex & dim, const ex void scalar_products::add_vectors(const lst & l, const ex & dim) { // Add all possible pairs of products - for (lst::const_iterator it1 = l.begin(); it1 != l.end(); ++it1) - for (lst::const_iterator it2 = l.begin(); it2 != l.end(); ++it2) - add(*it1, *it2, *it1 * *it2); + for (auto & it1 : l) + for (auto & it2 : l) + add(it1, it2, it1 * it2); } void scalar_products::clear() @@ -1282,14 +1381,52 @@ ex scalar_products::evaluate(const ex & v1, const ex & v2, const ex & dim) const void scalar_products::debugprint() const { std::cerr << "map size=" << spm.size() << std::endl; - spmap::const_iterator i = spm.begin(), end = spm.end(); - while (i != end) { - const spmapkey & k = i->first; + for (auto & it : spm) { + const spmapkey & k = it.first; std::cerr << "item key="; k.debugprint(); - std::cerr << ", value=" << i->second << std::endl; - ++i; + std::cerr << ", value=" << it.second << std::endl; + } +} + +exvector get_all_dummy_indices_safely(const ex & e) +{ + if (is_a(e)) + return ex_to(e).get_dummy_indices(); + else if (is_a(e) && e.op(1)==2) { + return e.op(0).get_free_indices(); + } + else if (is_a(e) || is_a(e)) { + exvector dummies; + exvector free_indices; + for (std::size_t i = 0; i < e.nops(); ++i) { + exvector dummies_of_factor = get_all_dummy_indices_safely(e.op(i)); + dummies.insert(dummies.end(), dummies_of_factor.begin(), + dummies_of_factor.end()); + exvector free_of_factor = e.op(i).get_free_indices(); + free_indices.insert(free_indices.begin(), free_of_factor.begin(), + free_of_factor.end()); + } + exvector free_out, dummy_out; + find_free_and_dummy(free_indices.begin(), free_indices.end(), free_out, + dummy_out); + dummies.insert(dummies.end(), dummy_out.begin(), dummy_out.end()); + return dummies; + } + else if(is_a(e)) { + exvector result; + for(std::size_t i = 0; i < e.nops(); ++i) { + exvector dummies_of_term = get_all_dummy_indices_safely(e.op(i)); + sort(dummies_of_term.begin(), dummies_of_term.end()); + exvector new_vec; + set_union(result.begin(), result.end(), dummies_of_term.begin(), + dummies_of_term.end(), std::back_inserter(new_vec), + ex_is_less()); + result.swap(new_vec); + } + return result; } + return exvector(); } /** Returns all dummy indices from the exvector */ @@ -1298,13 +1435,13 @@ exvector get_all_dummy_indices(const ex & e) exvector p; bool nc; product_to_exvector(e, p, nc); - exvector::const_iterator ip = p.begin(), ipend = p.end(); + auto ip = p.begin(), ipend = p.end(); exvector v, v1; while (ip != ipend) { if (is_a(*ip)) { v1 = ex_to(*ip).get_dummy_indices(); v.insert(v.end(), v1.begin(), v1.end()); - exvector::const_iterator ip1 = ip+1; + auto ip1 = ip + 1; while (ip1 != ipend) { if (is_a(*ip1)) { v1 = ex_to(*ip).get_dummy_indices(ex_to(*ip1)); @@ -1318,78 +1455,119 @@ exvector get_all_dummy_indices(const ex & e) return v; } -ex rename_dummy_indices_uniquely(const ex & a, const ex & b) +lst rename_dummy_indices_uniquely(const exvector & va, const exvector & vb) { - exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), common_indices; + exvector common_indices; set_intersection(va.begin(), va.end(), vb.begin(), vb.end(), std::back_insert_iterator(common_indices), ex_is_less()); if (common_indices.empty()) { - return b; + return lst{lst{}, lst{}}; } else { exvector new_indices, old_indices; old_indices.reserve(2*common_indices.size()); new_indices.reserve(2*common_indices.size()); exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end(); while (ip != ipend) { - if (is_a(*ip)) { - varidx mu((new symbol)->setflag(status_flags::dynallocated), ex_to(*ip).get_dim(), ex_to(*ip).is_covariant()); - old_indices.push_back(*ip); - new_indices.push_back(mu); + ex newsym = dynallocate(); + ex newidx; + if(is_exactly_a(*ip)) + newidx = dynallocate(newsym, ex_to(*ip).get_dim(), + ex_to(*ip).is_covariant(), + ex_to(*ip).is_dotted()); + else if (is_exactly_a(*ip)) + newidx = dynallocate(newsym, ex_to(*ip).get_dim(), + ex_to(*ip).is_covariant()); + else + newidx = dynallocate(newsym, ex_to(*ip).get_dim()); + old_indices.push_back(*ip); + new_indices.push_back(newidx); + if(is_a(*ip)) { old_indices.push_back(ex_to(*ip).toggle_variance()); - new_indices.push_back(mu.toggle_variance()); - } else { - old_indices.push_back(*ip); - new_indices.push_back(idx((new symbol)->setflag(status_flags::dynallocated), ex_to(*ip).get_dim())); + new_indices.push_back(ex_to(newidx).toggle_variance()); } ++ip; } - return b.subs(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()), subs_options::no_pattern); + return lst{lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end())}; } } +ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b) +{ + lst indices_subs = rename_dummy_indices_uniquely(va, vb); + return (indices_subs.op(0).nops()>0 ? b.subs(ex_to(indices_subs.op(0)), ex_to(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming) : b); +} + +ex rename_dummy_indices_uniquely(const ex & a, const ex & b) +{ + exvector va = get_all_dummy_indices_safely(a); + if (va.size() > 0) { + exvector vb = get_all_dummy_indices_safely(b); + if (vb.size() > 0) { + sort(va.begin(), va.end(), ex_is_less()); + sort(vb.begin(), vb.end(), ex_is_less()); + lst indices_subs = rename_dummy_indices_uniquely(va, vb); + if (indices_subs.op(0).nops() > 0) + return b.subs(ex_to(indices_subs.op(0)), ex_to(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming); + } + } + return b; +} + +ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va) +{ + if (va.size() > 0) { + exvector vb = get_all_dummy_indices_safely(b); + if (vb.size() > 0) { + sort(vb.begin(), vb.end(), ex_is_less()); + lst indices_subs = rename_dummy_indices_uniquely(va, vb); + if (indices_subs.op(0).nops() > 0) { + if (modify_va) { + for (auto & i : ex_to(indices_subs.op(1))) + va.push_back(i); + exvector uncommon_indices; + set_difference(vb.begin(), vb.end(), indices_subs.op(0).begin(), indices_subs.op(0).end(), std::back_insert_iterator(uncommon_indices), ex_is_less()); + for (auto & ip : uncommon_indices) + va.push_back(ip); + sort(va.begin(), va.end(), ex_is_less()); + } + return b.subs(ex_to(indices_subs.op(0)), ex_to(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming); + } + } + } + return b; +} + ex expand_dummy_sum(const ex & e, bool subs_idx) { ex e_expanded = e.expand(); pointer_to_map_function_1arg fcn(expand_dummy_sum, subs_idx); if (is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded)) { return e_expanded.map(fcn); - } else if (is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded)) { - exvector v = get_all_dummy_indices(e_expanded); - exvector::const_iterator it = v.begin(), itend = v.end(); - while (it != itend) { - varidx nu = ex_to(*it); - if (nu.is_dim_numeric()) { - ex en = 0; - for (int i=0; i < ex_to(nu.get_dim()).to_int(); i++) { - if (is_a(nu) && !subs_idx) { - en += e_expanded.subs(lst(nu == varidx(i, nu.get_dim(), true), nu.toggle_variance() == varidx(i, nu.get_dim()))); - } else { - en += e_expanded.subs(lst(nu == idx(i, nu.get_dim()), nu.toggle_variance() == idx(i, nu.get_dim()))); - } - } - return expand_dummy_sum(en, subs_idx); - } - ++it; - } - return e; - } else if (is_a(e_expanded)) { - exvector v = ex_to(e_expanded).get_dummy_indices(); - exvector::const_iterator it = v.begin(), itend = v.end(); - while (it != itend) { - varidx nu = ex_to(*it); - if (nu.is_dim_numeric()) { + } else if (is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded)) { + exvector v; + if (is_a(e_expanded)) + v = ex_to(e_expanded).get_dummy_indices(); + else + v = get_all_dummy_indices(e_expanded); + ex result = e_expanded; + for (const auto & nu : v) { + if (ex_to(nu).get_dim().info(info_flags::nonnegint)) { + int idim = ex_to(ex_to(nu).get_dim()).to_int(); ex en = 0; - for (int i=0; i < ex_to(nu.get_dim()).to_int(); i++) { - if (is_a(nu) && !subs_idx) { - en += e_expanded.subs(lst(nu == varidx(i, nu.get_dim(), true), nu.toggle_variance() == varidx(i, nu.get_dim()))); + for (int i=0; i < idim; i++) { + if (subs_idx && is_a(nu)) { + ex other = ex_to(nu).toggle_variance(); + en += result.subs(lst{ + nu == idx(i, idim), + other == idx(i, idim) + }); } else { - en += e_expanded.subs(lst(nu == idx(i, nu.get_dim()), nu.toggle_variance() == idx(i, nu.get_dim()))); + en += result.subs( nu.op(0) == i ); } } - return expand_dummy_sum(en, subs_idx); - } - ++it; + result = en; + } } - return e; + return result; } else { return e; }