From 22d004fe5759658a8232736e40b49202987e9407 Mon Sep 17 00:00:00 2001 From: Chris Dams Date: Fri, 16 Dec 2005 17:19:12 +0000 Subject: [PATCH] Improvements w.r.t. automatic renaming of dummy indices. --- ginac/expairseq.cpp | 6 +-- ginac/indexed.cpp | 104 +++++++++++++++++++++++++------------------- ginac/indexed.h | 4 ++ ginac/power.h | 1 - 4 files changed, 66 insertions(+), 49 deletions(-) diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 9f580d9e..8cf8d741 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -1626,7 +1626,7 @@ safe_inserter::safe_inserter(const ex&e, const bool disable_renaming) if(disable_renaming) dodummies=false; if(dodummies) { - dummy_indices = get_all_dummy_indices(e); + dummy_indices = get_all_dummy_indices_safely(e); sort(dummy_indices.begin(), dummy_indices.end(), ex_is_less()); } } @@ -1645,13 +1645,13 @@ void safe_inserter::insert_new_pair(const expair &p, const ex &orig_ex) epv->push_back(p); return; } - exvector dummies_of_factor = get_all_dummy_indices(p.rest); + exvector dummies_of_factor = get_all_dummy_indices_safely(p.rest); if(dummies_of_factor.size() == 0) { epv->push_back(p); return; } sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less()); - exvector dummies_of_orig_ex = get_all_dummy_indices(orig_ex); + exvector dummies_of_orig_ex = get_all_dummy_indices_safely(orig_ex); sort(dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(), ex_is_less()); exvector new_dummy_indices; new_dummy_indices.reserve(dummy_indices.size()); diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 0c1c904f..d39692c4 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -520,22 +520,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()) @@ -793,31 +777,21 @@ try_again: // At least one dummy index, is it a defined scalar product? bool contracted = false; - if (free.empty()) { - - try { - // 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) { - // User-defined scalar product? - if (sp.is_defined(*it1, *it2, dim)) { + ex dim = minimal_dim( + ex_to(it1->op(1)).get_dim(), + ex_to(it2->op(1)).get_dim() + ); - // Yes, substitute it - *it1 = sp.evaluate(*it1, *it2, dim); - *it2 = _ex1; - goto contraction_done; - } - } catch (const std::runtime_error&) {} + // User-defined scalar product? + if (sp.is_defined(*it1, *it2, dim)) { + + // Yes, substitute it + *it1 = sp.evaluate(*it1, *it2, dim); + *it2 = _ex1; + goto contraction_done; + } } // Try to contract the first one with the second one @@ -1353,6 +1327,46 @@ void scalar_products::debugprint() const } } +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 (int i=0; i(e)) { + exvector result; + for(int i=0; i(new_vec), + ex_is_less()); + result.swap(new_vec); + } + return result; + } + return exvector(); +} + /** Returns all dummy indices from the exvector */ exvector get_all_dummy_indices(const ex & e) { @@ -1420,20 +1434,20 @@ lst rename_dummy_indices_uniquely(const exvector & va, const exvector & vb) 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((lst)indices_subs.op(0), (lst)indices_subs.op(1), subs_options::no_pattern) : b); + 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(a); + exvector va = get_all_dummy_indices_safely(a); if (va.size() > 0) { - exvector vb = get_all_dummy_indices(b); + 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((lst)indices_subs.op(0), (lst)indices_subs.op(1), subs_options::no_pattern); + 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; @@ -1442,7 +1456,7 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b) ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va) { if (va.size() > 0) { - exvector vb = get_all_dummy_indices(b); + 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); @@ -1459,7 +1473,7 @@ ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va) } sort(va.begin(), va.end(), ex_is_less()); } - return b.subs((lst)indices_subs.op(0), (lst)indices_subs.op(1), subs_options::no_pattern); + return b.subs(ex_to(indices_subs.op(0)), ex_to(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming); } } } diff --git a/ginac/indexed.h b/ginac/indexed.h index 8f2779ff..0da2c421 100644 --- a/ginac/indexed.h +++ b/ginac/indexed.h @@ -254,6 +254,10 @@ template<> inline bool is_exactly_a(const basic & obj) /** Returns all dummy indices from the expression */ exvector get_all_dummy_indices(const ex & e); +/** More reliable version of the form. The former assumes that e is an + * expanded epxression. */ +exvector get_all_dummy_indices_safely(const ex & e); + /** Returns b with all dummy indices, which are listed in va, renamed * if modify_va is set to TRUE all dummy indices of b will be appended to va */ ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va = false); diff --git a/ginac/power.h b/ginac/power.h index 2a2c5004..d93bf981 100644 --- a/ginac/power.h +++ b/ginac/power.h @@ -65,7 +65,6 @@ public: ex normal(exmap & repl, exmap & rev_lookup, int level = 0) const; ex to_rational(exmap & repl) const; ex to_polynomial(exmap & repl) const; - exvector get_free_indices() const; ex conjugate() const; protected: ex derivative(const symbol & s) const; -- 2.44.0