From 0633ed8082961673eedc092689e06fa39d6bc322 Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Thu, 19 May 2005 18:10:40 +0000 Subject: [PATCH] Fixed bug in expanding expressions containing dummy indices. [V.Kisil] --- ginac/indexed.cpp | 129 +++++++++++++++++++++++++++++++++++++++++++--- ginac/indexed.h | 17 ++++++ ginac/mul.cpp | 14 ++--- ginac/ncmul.cpp | 22 ++++++-- ginac/power.cpp | 11 +++- 5 files changed, 177 insertions(+), 16 deletions(-) diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 6d74b09a..ef904f72 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -35,8 +35,10 @@ #include "operators.h" #include "lst.h" #include "archive.h" +#include "symbol.h" #include "utils.h" #include "integral.h" +#include "matrix.h" namespace GiNaC { @@ -669,16 +671,15 @@ struct ex_base_is_less : public std::binary_function { } }; -/** 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) +/* An auxiliary function used by simplify_indexed() and expand_dummy_sum() + * It returns an exvector of factors from the supplied product */ +static void product_to_exvector(const ex & e, exvector & v, bool & non_commutative) { // Remember whether the product was commutative or noncommutative // (because we chop it into factors and need to reassemble later) - bool non_commutative = is_exactly_a(e); + non_commutative = is_exactly_a(e); // Collect factors in an exvector, store squares twice - exvector v; v.reserve(e.nops() * 2); if (is_exactly_a(e)) { @@ -691,7 +692,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du ex f = e.op(i); if (is_exactly_a(f) && f.op(1).is_equal(_ex2)) { v.push_back(f.op(0)); - v.push_back(f.op(0)); + v.push_back(f.op(0)); } else if (is_exactly_a(f)) { // Noncommutative factor found, split it as well non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later @@ -701,6 +702,19 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du v.push_back(f); } } +} + +/** 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) +{ + // Collect factors in an exvector + exvector v; + + // Remember whether the product was commutative or noncommutative + // (because we chop it into factors and need to reassemble later) + bool non_commutative; + product_to_exvector(e, v, non_commutative); // Perform contractions bool something_changed = false; @@ -1278,4 +1292,107 @@ void scalar_products::debugprint() const } } +/** Returns all dummy indices from the exvector */ +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(); + 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; + while (ip1 != ipend) { + if (is_a(*ip1)) { + v1 = ex_to(*ip).get_dummy_indices(ex_to(*ip1)); + v.insert(v.end(), v1.begin(), v1.end()); + } + ++ip1; + } + } + ++ip; + } + return v; +} + +ex rename_dummy_indices_uniquely(const ex & a, const ex & b) +{ + exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), 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; + } 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); + 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())); + } + ++ip; + } + return b.subs(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()), subs_options::no_pattern); + } +} + +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()) { + 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 { + return e; + } +} + } // namespace GiNaC diff --git a/ginac/indexed.h b/ginac/indexed.h index 7329e87d..ad32c237 100644 --- a/ginac/indexed.h +++ b/ginac/indexed.h @@ -250,6 +250,23 @@ template<> inline bool is_exactly_a(const basic & obj) return obj.tinfo()==TINFO_indexed; } +/** Returns all dummy indices from the expression */ +exvector get_all_dummy_indices(const ex & e); + +/** Returns b with all dummy indices, which are common with a, renamed */ +ex rename_dummy_indices_uniquely(const ex & a, const ex & b); + +/** This function returns the given expression with expanded sums + * for all dummy index summations, where the dimensionality of + * the dummy index is numeric. + * Optionally all indices with a variance will be substituted by + * indices with the corresponding numeric values without variance. + * + * @param e the given expression + * @param subs_idx indicates if variance of dummy indixes should be neglected + */ +ex expand_dummy_sum(const ex & e, bool subs_idx = false); + } // namespace GiNaC #endif // ndef __GINAC_INDEXED_H__ diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 924493d2..fea5b346 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -30,6 +30,7 @@ #include "power.h" #include "operators.h" #include "matrix.h" +#include "indexed.h" #include "lst.h" #include "archive.h" #include "utils.h" @@ -904,11 +905,12 @@ ex mul::expand(unsigned options) const for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) { // Don't push_back expairs which might have a rest that evaluates to a numeric, // since that would violate an invariant of expairseq: - const ex rest = (new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated); - if (is_exactly_a(rest)) + const ex rest = (new mul(i1->rest, rename_dummy_indices_uniquely(i1->rest, i2->rest)))->setflag(status_flags::dynallocated); + if (is_exactly_a(rest)) { oc += ex_to(rest).mul(ex_to(i1->coeff).mul(ex_to(i2->coeff))); - else + } else { distrseq.push_back(expair(rest, ex_to(i1->coeff).mul_dyn(ex_to(i2->coeff)))); + } } tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated); } @@ -934,11 +936,11 @@ ex mul::expand(unsigned options) const for (size_t i=0; isetflag(status_flags::dynallocated); - if (can_be_further_expanded(term)) + if (can_be_further_expanded(term)) { distrseq.push_back(term.expand()); - else { + } else { if (options == 0) ex_to(term).setflag(status_flags::expanded); distrseq.push_back(term); diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp index 4424253c..eeadf4e7 100644 --- a/ginac/ncmul.cpp +++ b/ginac/ncmul.cpp @@ -30,6 +30,7 @@ #include "mul.h" #include "matrix.h" #include "archive.h" +#include "indexed.h" #include "utils.h" namespace GiNaC { @@ -167,10 +168,25 @@ ex ncmul::expand(unsigned options) const intvector k(number_of_adds); + /* Rename indices in the static members of the product */ + exvector expanded_seq_mod; + size_t j = 0; + size_t i; + for (i=0; i setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0))); diff --git a/ginac/power.cpp b/ginac/power.cpp index bd276c56..79297685 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -853,8 +853,17 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr { GINAC_ASSERT(n.is_integer()); - if (n.is_zero()) + if (n.is_zero()) { return _ex1; + } + + // Leave it to multiplication since dummy indices have to be renamed + if (get_all_dummy_indices(m).size() > 0) { + ex result = m; + for (int i=1; i < n.to_int(); i++) + result *= rename_dummy_indices_uniquely(m,m); + return result; + } epvector distrseq; distrseq.reserve(m.seq.size()); -- 2.45.0