+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)) {
+ exvector v = get_all_dummy_indices(e_expanded);
+ exvector::const_iterator it = v.begin(), itend = v.end();
+ while (it != itend) {
+ varidx nu = ex_to<varidx>(*it);
+ if (nu.is_dim_numeric()) {
+ ex en = 0;
+ for (int i=0; i < ex_to<numeric>(nu.get_dim()).to_int(); i++) {
+ if (is_a<varidx>(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<indexed>(e_expanded)) {
+ exvector v = ex_to<indexed>(e_expanded).get_dummy_indices();
+ exvector::const_iterator it = v.begin(), itend = v.end();
+ while (it != itend) {
+ varidx nu = ex_to<varidx>(*it);
+ if (nu.is_dim_numeric()) {
+ ex en = 0;
+ for (int i=0; i < ex_to<numeric>(nu.get_dim()).to_int(); i++) {
+ if (is_a<varidx>(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;
+ }