+/** 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<indexed>(*ip)) {
+ v1 = ex_to<indexed>(*ip).get_dummy_indices();
+ v.insert(v.end(), v1.begin(), v1.end());
+ exvector::const_iterator ip1 = ip+1;
+ while (ip1 != ipend) {
+ if (is_a<indexed>(*ip1)) {
+ v1 = ex_to<indexed>(*ip).get_dummy_indices(ex_to<indexed>(*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<exvector>(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<varidx>(*ip)) {
+ varidx mu((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim(), ex_to<varidx>(*ip).is_covariant());
+ old_indices.push_back(*ip);
+ new_indices.push_back(mu);
+ old_indices.push_back(ex_to<varidx>(*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<varidx>(*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<bool> fcn(expand_dummy_sum, subs_idx);
+ if (is_a<add>(e_expanded) || is_a<lst>(e_expanded) || is_a<matrix>(e_expanded)) {
+ return e_expanded.map(fcn);
+ } else if (is_a<ncmul>(e_expanded) || is_a<mul>(e_expanded) || is_a<power>(e_expanded) || is_a<indexed>(e_expanded)) {
+ exvector v;
+ if (is_a<indexed>(e_expanded))
+ v = ex_to<indexed>(e_expanded).get_dummy_indices();
+ 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;
+ if (ex_to<idx>(nu).get_dim().info(info_flags::nonnegint)) {
+ int idim = ex_to<numeric>(ex_to<idx>(nu).get_dim()).to_int();
+ ex en = 0;
+ for (int i=0; i < idim; i++) {
+ if (subs_idx && is_a<varidx>(nu)) {
+ ex other = ex_to<varidx>(nu).toggle_variance();
+ en += result.subs(lst(
+ nu == idx(i, idim),
+ other == idx(i, idim)
+ ));
+ } else {
+ en += result.subs( nu.op(0) == i );
+ }
+ }
+ result = en;
+ }
+ }
+ return result;
+ } else {
+ return e;
+ }
+}
+