From ea5d361d94e49ca3f3b73db8c9812ee519f0633f Mon Sep 17 00:00:00 2001 From: Chris Dams Date: Fri, 11 Nov 2005 15:05:23 +0000 Subject: [PATCH] Index renaming issues, sped up simplify_indexed, used defined NC-objects can carry indices. --- doc/tutorial/ginac.texi | 8 ++- ginac/color.h | 2 + ginac/expairseq.cpp | 130 +++++++++++++++++++++++++++++------ ginac/flags.h | 3 +- ginac/indexed.cpp | 147 ++++++++++++++++++++++++++++++---------- ginac/indexed.h | 6 +- ginac/symmetry.cpp | 4 +- ginac/tensor.h | 4 ++ 8 files changed, 243 insertions(+), 61 deletions(-) diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 06ab56fe..3cd6c617 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -4278,13 +4278,17 @@ contain the same number of elements). Using this form, you would write @end example The optional last argument to @code{subs()} is a combination of -@code{subs_options} flags. There are two options available: +@code{subs_options} flags. There are three options available: @code{subs_options::no_pattern} disables pattern matching, which makes large @code{subs()} operations significantly faster if you are not using patterns. The second option, @code{subs_options::algebraic} enables algebraic substitutions in products and powers. @ref{Pattern Matching and Advanced Substitutions}, for more information -about patterns and algebraic substitutions. +about patterns and algebraic substitutions. The third option, +@code{subs_options::no_index_renaming} disables the feature that dummy +indices are renamed if the subsitution could give a result in which a +dummy index occurs more than two times. This is sometimes necessary if +you want to use @code{subs()} to rename your dummy indices. @code{subs()} performs syntactic substitution of any complete algebraic object; it does not try to match sub-expressions as is demonstrated by the diff --git a/ginac/color.h b/ginac/color.h index a585f2e5..9d9ba26e 100644 --- a/ginac/color.h +++ b/ginac/color.h @@ -109,6 +109,7 @@ public: // non-virtual functions in this class protected: + unsigned return_type() const { return return_types::commutative; } void do_print(const print_context & c, unsigned level) const; void do_print_latex(const print_latex & c, unsigned level) const; }; @@ -125,6 +126,7 @@ public: // non-virtual functions in this class protected: + unsigned return_type() const { return return_types::commutative; } void do_print(const print_context & c, unsigned level) const; void do_print_latex(const print_latex & c, unsigned level) const; }; diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 7555bd18..9f580d9e 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -34,6 +34,7 @@ #include "archive.h" #include "operators.h" #include "utils.h" +#include "indexed.h" #if EXPAIRSEQ_USE_HASHTAB #include @@ -757,8 +758,15 @@ void expairseq::construct_from_2_ex(const ex &lh, const ex &rh) construct_from_2_ex_via_exvector(lh,rh); } else { #endif // EXPAIRSEQ_USE_HASHTAB - construct_from_2_expairseq(ex_to(lh), - ex_to(rh)); + if(is_a(lh)) + { + ex newrh=rename_dummy_indices_uniquely(lh, rh); + construct_from_2_expairseq(ex_to(lh), + ex_to(newrh)); + } + else + construct_from_2_expairseq(ex_to(lh), + ex_to(rh)); #if EXPAIRSEQ_USE_HASHTAB } #endif // EXPAIRSEQ_USE_HASHTAB @@ -1008,13 +1016,27 @@ void expairseq::make_flat(const exvector &v) seq.reserve(v.size()+noperands-nexpairseqs); // copy elements and split off numerical part + exvector dummy_indices; cit = v.begin(); while (cit!=v.end()) { if (ex_to(*cit).tinfo()==this->tinfo()) { - const expairseq &subseqref = ex_to(*cit); - combine_overall_coeff(subseqref.overall_coeff); - epvector::const_iterator cit_s = subseqref.seq.begin(); - while (cit_s!=subseqref.seq.end()) { + const expairseq *subseqref; + ex newfactor; + if(is_a(*cit)) + { + exvector dummies_of_factor = get_all_dummy_indices(*cit); + sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less()); + newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, *cit); + subseqref = &(ex_to(newfactor)); + exvector new_dummy_indices; + set_union(dummy_indices.begin(), dummy_indices.end(), dummies_of_factor.begin(), dummies_of_factor.end(), std::back_insert_iterator(new_dummy_indices), ex_is_less()); + dummy_indices.swap(new_dummy_indices); + } + else + subseqref = &ex_to(*cit); + combine_overall_coeff(subseqref->overall_coeff); + epvector::const_iterator cit_s = subseqref->seq.begin(); + while (cit_s!=subseqref->seq.end()) { seq.push_back(*cit_s); ++cit_s; } @@ -1579,6 +1601,67 @@ std::auto_ptr expairseq::evalchildren(int level) const return std::auto_ptr(0); // signalling nothing has changed } +class safe_inserter +{ + public: + safe_inserter(const ex&, const bool disable_renaming=false); + std::auto_ptr getseq(){return epv;} + void insert_old_pair(const expair &p) + { + epv->push_back(p); + } + void insert_new_pair(const expair &p, const ex &orig_ex); + private: + std::auto_ptr epv; + bool dodummies; + exvector dummy_indices; + void update_dummy_indices(const exvector&); +}; + +safe_inserter::safe_inserter(const ex&e, const bool disable_renaming) + :epv(new epvector) +{ + epv->reserve(e.nops()); + dodummies=is_a(e); + if(disable_renaming) + dodummies=false; + if(dodummies) { + dummy_indices = get_all_dummy_indices(e); + sort(dummy_indices.begin(), dummy_indices.end(), ex_is_less()); + } +} + +void safe_inserter::update_dummy_indices(const exvector &v) +{ + exvector new_dummy_indices; + set_union(dummy_indices.begin(), dummy_indices.end(), v.begin(), v.end(), + std::back_insert_iterator(new_dummy_indices), ex_is_less()); + dummy_indices.swap(new_dummy_indices); +} + +void safe_inserter::insert_new_pair(const expair &p, const ex &orig_ex) +{ + if(!dodummies) { + epv->push_back(p); + return; + } + exvector dummies_of_factor = get_all_dummy_indices(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); + 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()); + set_difference(dummy_indices.begin(), dummy_indices.end(), dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(), + std::back_insert_iterator(new_dummy_indices), ex_is_less()); + dummy_indices.swap(new_dummy_indices); + ex newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, p.rest); + update_dummy_indices(dummies_of_factor); + epv -> push_back(expair(newfactor, p.coeff)); +} /** Member-wise substitute in this sequence. * @@ -1614,22 +1697,27 @@ std::auto_ptr expairseq::subschildren(const exmap & m, unsigned option if (!are_ex_trivially_equal(orig_ex, subsed_ex)) { // Something changed, copy seq, subs and return it - std::auto_ptr s(new epvector); - s->reserve(seq.size()); + safe_inserter s(*this, options & subs_options::no_index_renaming); // Copy parts of seq which are known not to have changed - s->insert(s->begin(), seq.begin(), cit); + for(epvector::const_iterator i=seq.begin(); i!=cit; ++i) + s.insert_old_pair(*i); // Copy first changed element - s->push_back(split_ex_to_pair(subsed_ex)); + s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex); ++cit; // Copy rest while (cit != last) { - s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options))); + ex orig_ex = recombine_pair_to_ex(*cit); + ex subsed_ex = orig_ex.subs(m, options); + if(are_ex_trivially_equal(orig_ex, subsed_ex)) + s.insert_old_pair(*cit); + else + s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex); ++cit; } - return s; + return s.getseq(); } ++cit; @@ -1645,23 +1733,27 @@ std::auto_ptr expairseq::subschildren(const exmap & m, unsigned option if (!are_ex_trivially_equal(cit->rest, subsed_ex)) { // Something changed, copy seq, subs and return it - std::auto_ptr s(new epvector); - s->reserve(seq.size()); + safe_inserter s(*this, options & subs_options::no_index_renaming); // Copy parts of seq which are known not to have changed - s->insert(s->begin(), seq.begin(), cit); + for(epvector::const_iterator i=seq.begin(); i!=cit; ++i) + s.insert_old_pair(*i); // Copy first changed element - s->push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff)); + s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), cit->rest); ++cit; // Copy rest while (cit != last) { - s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options), - cit->coeff)); + const ex &orig_ex = cit->rest; + const ex &subsed_ex = cit->rest.subs(m, options); + if(are_ex_trivially_equal(orig_ex, subsed_ex)) + s.insert_old_pair(*cit); + else + s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), orig_ex); ++cit; } - return s; + return s.getseq(); } ++cit; diff --git a/ginac/flags.h b/ginac/flags.h index 90f58ad2..7c51b75d 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -43,7 +43,8 @@ public: algebraic = 0x0002, ///< enable algebraic substitutions subs_algebraic = 0x0002, // for backwards compatibility pattern_is_product = 0x0004, ///< used internally by expairseq::subschildren() - pattern_is_not_product = 0x0008 ///< used internally by expairseq::subschildren() + pattern_is_not_product = 0x0008, ///< used internally by expairseq::subschildren() + no_index_renaming = 0x0010 }; }; diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index ef904f72..64eba6d0 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -297,6 +297,9 @@ ex indexed::eval(int level) const return f * thiscontainer(v); } + if(this->tinfo()==TINFO_indexed && seq.size()==1) + return base; + // Canonicalize indices according to the symmetry properties if (seq.size() > 2) { exvector v = seq; @@ -324,6 +327,14 @@ ex indexed::thiscontainer(std::auto_ptr vp) const return indexed(ex_to(symtree), vp); } +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 { GINAC_ASSERT(seq.size() > 0); @@ -532,6 +543,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(exvector::const_iterator i=v.begin(); i!=v.end(); ++i) + if(is_exactly_a(*i)) + ++number; + return number; +} + /** Rename dummy indices in an expression. * * @param e Expression to work on @@ -540,10 +560,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) @@ -557,7 +577,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex int remaining = local_size - global_size; exvector::const_iterator 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 +595,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 @@ -704,6 +726,18 @@ 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)); + if(dummy_syms.size() < 2) + return r; + ex q=symmetrize(r, dummy_syms); + return q; +} + /** 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) @@ -864,19 +898,26 @@ contraction_done: // 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()) { - free_indices.clear(); - return _ex0; - } + 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 @@ -943,6 +984,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 +1070,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; iop(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); if (term_symm.is_zero()) continue; terms.push_back(terminfo(term, term_symm)); @@ -1318,9 +1374,9 @@ exvector get_all_dummy_indices(const ex & e) return v; } -ex rename_dummy_indices_uniquely(const ex & a, const ex & b) +ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b) { - 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; @@ -1330,15 +1386,25 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b) 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=(new symbol)->setflag(status_flags::dynallocated); + ex newidx; + if(is_exactly_a(*ip)) + newidx = (new spinidx(newsym, ex_to(*ip).get_dim(), + ex_to(*ip).is_covariant(), + ex_to(*ip).is_dotted())) + -> setflag(status_flags::dynallocated); + else if (is_exactly_a(*ip)) + newidx = (new varidx(newsym, ex_to(*ip).get_dim(), + ex_to(*ip).is_covariant())) + -> setflag(status_flags::dynallocated); + else + newidx = (new idx(newsym, ex_to(*ip).get_dim())) + -> setflag(status_flags::dynallocated); + 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; } @@ -1346,6 +1412,15 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b) } } +ex rename_dummy_indices_uniquely(const ex & a, const ex & b) +{ + exvector va = get_all_dummy_indices(a); + exvector vb = get_all_dummy_indices(b); + sort(va.begin(), va.end(), ex_is_less()); + sort(vb.begin(), vb.end(), ex_is_less()); + return rename_dummy_indices_uniquely(va, vb, b); +} + ex expand_dummy_sum(const ex & e, bool subs_idx) { ex e_expanded = e.expand(); diff --git a/ginac/indexed.h b/ginac/indexed.h index ad32c237..1ce81945 100644 --- a/ginac/indexed.h +++ b/ginac/indexed.h @@ -153,7 +153,8 @@ protected: ex derivative(const symbol & s) const; ex thiscontainer(const exvector & v) const; ex thiscontainer(std::auto_ptr vp) const; - unsigned return_type() const { return return_types::commutative; } + unsigned return_type() const; + unsigned return_type_tinfo() const { return op(0).return_type_tinfo(); } ex expand(unsigned options = 0) const; // new virtual functions which can be overridden by derived classes @@ -256,6 +257,9 @@ 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); +/** Same as above, where va and vb contain the indices of a and b and are sorted */ +ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, 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. diff --git a/ginac/symmetry.cpp b/ginac/symmetry.cpp index f29c6a38..39b194a6 100644 --- a/ginac/symmetry.cpp +++ b/ginac/symmetry.cpp @@ -427,7 +427,7 @@ static ex symm(const ex & e, exvector::const_iterator first, exvector::const_ite lst new_lst; for (unsigned i=0; i