X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=64eba6d08cf1cccf6dbf455cbcce4c01356dd3e7;hp=ef904f7222bfc798420dee5cdb6bed6373426a26;hb=ea5d361d94e49ca3f3b73db8c9812ee519f0633f;hpb=6715b097ae926fecfc62e53d4af38b7217634908 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();