X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=d5f7bf85eccaeaa6fc92736a3530afc4939dfe46;hp=0360d3416fc9f5c372b933efd4ffba427807f006;hb=9593ce33c14b7ff535d113f8a825f4c42ca81912;hpb=1602530f716ba1d425a0667b897182b99c374823 diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 0360d341..d5f7bf85 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's indexed expressions. */ /* - * GiNaC Copyright (C) 1999-2009 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2015 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 @@ -117,11 +117,11 @@ indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), sym { } -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) { } -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) { } @@ -166,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)) { @@ -243,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 @@ -256,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 @@ -324,9 +319,9 @@ 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 @@ -377,7 +372,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")); @@ -454,7 +449,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; @@ -532,8 +527,8 @@ exvector integral::get_free_indices() const template size_t number_of_type(const exvector&v) { size_t number = 0; - for(exvector::const_iterator i=v.begin(); i!=v.end(); ++i) - if(is_exactly_a(*i)) + for (auto & it : v) + if (is_exactly_a(it)) ++number; return number; } @@ -561,7 +556,7 @@ template static ex rename_dummy_indices(const ex & e, exvector & global // 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 (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); @@ -635,8 +630,7 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector for (size_t j=i+1; jop(0)) { variant_dummy_indices.erase(k); break; @@ -648,7 +642,7 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector } // In the case where a dummy symbol occurs twice in the same indexed object - // we try all posibilities of raising/lowering and keep the least one in + // 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(); @@ -678,8 +672,7 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector // If a dummy index is encountered for the first time in the // product, pull it up, otherwise, pull it down - for (exvector::iterator it2 = seq.begin()+1, it2end = seq.end(); - it2 != it2end; ++it2) { + for (auto it2 = seq.begin()+1, it2end = seq.end(); it2 != it2end; ++it2) { if (!is_exactly_a(*it2)) continue; @@ -771,9 +764,9 @@ 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 (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it) - if(is_exactly_a(*it)) - dummy_syms.push_back(it->op(0)); + 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); @@ -797,6 +790,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++) { @@ -806,6 +800,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 @@ -872,8 +867,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 @@ -882,6 +885,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; + } } } @@ -928,27 +936,29 @@ 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. - 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; + 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 @@ -971,7 +981,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 { @@ -1114,9 +1124,9 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi const ex & term = sum.op(i); exvector dummy_indices_of_term; dummy_indices_of_term.reserve(dummy_indices.size()); - for(exvector::iterator i=dummy_indices.begin(); i!=dummy_indices.end(); ++i) - if(hasindex(term,i->op(0))) - dummy_indices_of_term.push_back(*i); + for (auto & i : dummy_indices) + if (hasindex(term,i.op(0))) + dummy_indices_of_term.push_back(i); ex term_symm = idx_symmetrization(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); @@ -1132,7 +1142,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++; @@ -1147,13 +1157,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 @@ -1162,10 +1172,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 @@ -1198,7 +1208,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++; @@ -1350,9 +1360,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() @@ -1375,13 +1385,11 @@ 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; } } @@ -1431,13 +1439,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)); @@ -1520,15 +1528,12 @@ ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va) lst indices_subs = rename_dummy_indices_uniquely(va, vb); if (indices_subs.op(0).nops() > 0) { if (modify_va) { - for (lst::const_iterator i = ex_to(indices_subs.op(1)).begin(); i != ex_to(indices_subs.op(1)).end(); ++i) - va.push_back(*i); + 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()); - exvector::const_iterator ip = uncommon_indices.begin(), ipend = uncommon_indices.end(); - while (ip != ipend) { - va.push_back(*ip); - ++ip; - } + 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); @@ -1551,8 +1556,7 @@ ex expand_dummy_sum(const ex & e, bool subs_idx) else v = get_all_dummy_indices(e_expanded); ex result = e_expanded; - for(exvector::const_iterator it=v.begin(); it!=v.end(); ++it) { - ex nu = *it; + 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;