X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=bd41f44b0b38a271d243564c945d674bea9acaa3;hp=34fbc65f93122eeee823b6cb3e0da8138719b1f6;hb=ef2ed9a711c416722c34548015b4ae3621948c1c;hpb=695f6ae955ec530cded8f21efd5569df39447f76

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 34fbc65f..bd41f44b 100644
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
 #include <iostream>
@@ -35,8 +35,11 @@
 #include "operators.h"
 #include "lst.h"
 #include "archive.h"
+#include "symbol.h"
 #include "utils.h"
 #include "integral.h"
+#include "matrix.h"
+#include "inifcns.h"
 
 namespace GiNaC {
 
@@ -51,7 +54,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(indexed, exprseq,
 
 indexed::indexed() : symtree(not_symmetric())
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 }
 
 //////////
@@ -60,79 +63,79 @@ indexed::indexed() : symtree(not_symmetric())
 
 indexed::indexed(const ex & b) : inherited(b), symtree(not_symmetric())
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(not_symmetric())
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(not_symmetric())
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(not_symmetric())
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(not_symmetric())
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(symm)
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(symm)
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(symm)
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(not_symmetric())
 {
 	seq.insert(seq.end(), v.begin(), v.end());
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm)
 {
 	seq.insert(seq.end(), v.begin(), v.end());
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 	validate();
 }
 
 indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm)
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 }
 
 indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm)
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 }
 
 indexed::indexed(const symmetry & symm, std::auto_ptr<exvector> vp) : inherited(vp), symtree(symm)
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 }
 
 //////////
@@ -295,6 +298,9 @@ ex indexed::eval(int level) const
 		return f * thiscontainer(v);
 	}
 
+	if(this->tinfo()==&indexed::tinfo_static && seq.size()==1)
+		return base;
+
 	// Canonicalize indices according to the symmetry properties
 	if (seq.size() > 2) {
 		exvector v = seq;
@@ -312,6 +318,20 @@ ex indexed::eval(int level) const
 	return ex_to<basic>(base).eval_indexed(*this);
 }
 
+ex indexed::real_part() const
+{
+	if(op(0).info(info_flags::real))
+		return *this;
+	return real_part_function(*this).hold();
+}
+
+ex indexed::imag_part() const
+{
+	if(op(0).info(info_flags::real))
+		return 0;
+	return imag_part_function(*this).hold();
+}
+
 ex indexed::thiscontainer(const exvector & v) const
 {
 	return indexed(ex_to<symmetry>(symtree), v);
@@ -322,6 +342,14 @@ ex indexed::thiscontainer(std::auto_ptr<exvector> vp) const
 	return indexed(ex_to<symmetry>(symtree), vp);
 }
 
+unsigned indexed::return_type() const
+{
+	if(is_a<matrix>(op(0)))
+		return return_types::commutative;
+	else
+		return op(0).return_type();
+}
+
 ex indexed::expand(unsigned options) const
 {
 	GINAC_ASSERT(seq.size() > 0);
@@ -507,22 +535,6 @@ struct is_summation_idx : public std::unary_function<ex, bool> {
 	}
 };
 
-exvector power::get_free_indices() const
-{
-	// Get free indices of basis
-	exvector basis_indices = basis.get_free_indices();
-
-	if (exponent.info(info_flags::even)) {
-		// If the exponent is an even number, then any "free" index that
-		// forms a dummy pair with itself is actually a summation index
-		exvector really_free;
-		std::remove_copy_if(basis_indices.begin(), basis_indices.end(),
-		                    std::back_inserter(really_free), is_summation_idx());
-		return really_free;
-	} else
-		return basis_indices;
-}
-
 exvector integral::get_free_indices() const
 {
 	if (a.get_free_indices().size() || b.get_free_indices().size())
@@ -530,6 +542,15 @@ exvector integral::get_free_indices() const
 	return f.get_free_indices();
 }
 
+template<class T> 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<T>(*i))
+			++number;
+	return number;
+}
+
 /** Rename dummy indices in an expression.
  *
  *  @param e Expression to work on
@@ -538,10 +559,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<class T> 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<T>(global_dummy_indices),
+	       local_size = number_of_type<T>(local_dummy_indices);
 
 	// Any local dummy indices at all?
 	if (local_size == 0)
@@ -555,7 +576,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<T>(*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--;
@@ -573,11 +594,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_size; i++)
-		local_syms.push_back(local_dummy_indices[i].op(0));
+	for (size_t i=0; local_syms.size()!=local_size; i++)
+		if(is_exactly_a<T>(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<local_size; i++) // don't use more global symbols than necessary
-		global_syms.push_back(global_dummy_indices[i].op(0));
+	for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary
+		if(is_exactly_a<T>(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
@@ -669,16 +692,15 @@ struct ex_base_is_less : public std::binary_function<ex, ex, bool> {
 	}
 };
 
-/** 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<ncmul>(e);
+	non_commutative = is_exactly_a<ncmul>(e);
 
 	// Collect factors in an exvector, store squares twice
-	exvector v;
 	v.reserve(e.nops() * 2);
 
 	if (is_exactly_a<power>(e)) {
@@ -691,7 +713,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du
 			ex f = e.op(i);
 			if (is_exactly_a<power>(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<ncmul>(f)) {
 				// Noncommutative factor found, split it as well
 				non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later
@@ -701,6 +723,34 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du
 				v.push_back(f);
 		}
 	}
+}
+
+template<class T> 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<T>(*it))
+				dummy_syms.push_back(it->op(0));
+	if(dummy_syms.size() < 2)
+		return r;
+	ex q=symmetrize(r, dummy_syms);
+	return q;
+}
+
+// Forward declaration needed in absence of friend injection, C.f. [namespace.memdef]:
+ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp);
+
+/** 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;
@@ -742,20 +792,12 @@ try_again:
 
 			// At least one dummy index, is it a defined scalar product?
 			bool contracted = false;
-			if (free.empty()) {
-
-				// Find minimal dimension of all indices of both factors
-				exvector::const_iterator dit = ex_to<indexed>(*it1).seq.begin() + 1, ditend = ex_to<indexed>(*it1).seq.end();
-				ex dim = ex_to<idx>(*dit).get_dim();
-				++dit;
-				for (; dit != ditend; ++dit) {
-					dim = minimal_dim(dim, ex_to<idx>(*dit).get_dim());
-				}
-				dit = ex_to<indexed>(*it2).seq.begin() + 1;
-				ditend = ex_to<indexed>(*it2).seq.end();
-				for (; dit != ditend; ++dit) {
-					dim = minimal_dim(dim, ex_to<idx>(*dit).get_dim());
-				}
+			if (free.empty() && it1->nops()==2 && it2->nops()==2) {
+
+				ex dim = minimal_dim(
+					ex_to<idx>(it1->op(1)).get_dim(),
+					ex_to<idx>(it2->op(1)).get_dim()
+				);
 
 				// User-defined scalar product?
 				if (sp.is_defined(*it1, *it2, dim)) {
@@ -850,19 +892,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<idx>(r, local_dummy_indices);
+	if (q.is_zero()) {
+		free_indices.clear();
+		return _ex0;
+	}
+	q = idx_symmetrization<varidx>(q, local_dummy_indices);
+	if (q.is_zero()) {
+		free_indices.clear();
+		return _ex0;
+	}
+	q = idx_symmetrization<spinidx>(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<idx>(r, dummy_indices, local_dummy_indices);
+	r = rename_dummy_indices<varidx>(r, dummy_indices, local_dummy_indices);
+	r = rename_dummy_indices<spinidx>(r, dummy_indices, local_dummy_indices);
 
 	// Product of indexed object with a scalar?
 	if (is_exactly_a<mul>(r) && r.nops() == 2
@@ -929,6 +978,17 @@ public:
 	}
 };
 
+bool hasindex(const ex &x, const ex &sym)
+{	
+	if(is_a<idx>(x) && x.op(0)==sym)
+		return true;
+	else
+		for(size_t i=0; i<x.nops(); ++i)
+			if(hasindex(x.op(i), sym))
+				return true;
+	return false;
+}
+
 /** Simplify indexed expression, return list of free indices. */
 ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
 {
@@ -957,7 +1017,10 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
 		}
 
 		// Rename the dummy indices
-		return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);
+		e_expanded = rename_dummy_indices<idx>(e_expanded, dummy_indices, local_dummy_indices);
+		e_expanded = rename_dummy_indices<varidx>(e_expanded, dummy_indices, local_dummy_indices);
+		e_expanded = rename_dummy_indices<spinidx>(e_expanded, dummy_indices, local_dummy_indices);
+		return e_expanded;
 	}
 
 	// Simplification of sum = sum of simplifications, check consistency of
@@ -1001,18 +1064,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<terminfo> terms;
 		for (size_t i=0; i<sum.nops(); i++) {
 			const ex & term = sum.op(i);
-			ex term_symm = symmetrize(term, dummy_syms);
+			exvector dummy_indices_of_term;
+			dummy_indices_of_term.reserve(dummy_indices.size());
+			for(exvector::iterator i=dummy_indices.begin(); i!=dummy_indices.end(); ++i)
+				if(hasindex(term,i->op(0)))
+					dummy_indices_of_term.push_back(*i);
+			ex term_symm = idx_symmetrization<idx>(term, dummy_indices_of_term);
+			term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices_of_term);
+			term_symm = idx_symmetrization<spinidx>(term_symm, dummy_indices_of_term);
 			if (term_symm.is_zero())
 				continue;
 			terms.push_back(terminfo(term, term_symm));
@@ -1278,4 +1342,195 @@ void scalar_products::debugprint() const
 	}
 }
 
+exvector get_all_dummy_indices_safely(const ex & e)
+{
+	if (is_a<indexed>(e))
+		return ex_to<indexed>(e).get_dummy_indices();
+	else if (is_a<power>(e) && e.op(1)==2) {
+		return e.op(0).get_free_indices();
+	}	
+	else if (is_a<mul>(e) || is_a<ncmul>(e)) {
+		exvector dummies;
+		exvector free_indices;
+		for (int i=0; i<e.nops(); ++i) {
+			exvector dummies_of_factor = get_all_dummy_indices_safely(e.op(i));
+			dummies.insert(dummies.end(), dummies_of_factor.begin(),
+				dummies_of_factor.end());
+			exvector free_of_factor = e.op(i).get_free_indices();
+			free_indices.insert(free_indices.begin(), free_of_factor.begin(),
+				free_of_factor.end());
+		}
+		exvector free_out, dummy_out;
+		find_free_and_dummy(free_indices.begin(), free_indices.end(), free_out,
+			dummy_out);
+		dummies.insert(dummies.end(), dummy_out.begin(), dummy_out.end());
+		return dummies;
+	}
+	else if(is_a<add>(e)) {
+		exvector result;
+		for(int i=0; i<e.nops(); ++i) {
+			exvector dummies_of_term = get_all_dummy_indices_safely(e.op(i));
+			sort(dummies_of_term.begin(), dummies_of_term.end());
+			exvector new_vec;
+			set_union(result.begin(), result.end(), dummies_of_term.begin(),
+				dummies_of_term.end(), std::back_inserter<exvector>(new_vec),
+				ex_is_less());
+			result.swap(new_vec);
+		}
+		return result;
+	}
+	return exvector();
+}
+
+/** 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;
+}
+
+lst rename_dummy_indices_uniquely(const exvector & va, const exvector & vb)
+{
+	exvector 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 lst(lst(), lst());
+	} 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) {
+			ex newsym=(new symbol)->setflag(status_flags::dynallocated);
+			ex newidx;
+			if(is_exactly_a<spinidx>(*ip))
+				newidx = (new spinidx(newsym, ex_to<spinidx>(*ip).get_dim(),
+						ex_to<spinidx>(*ip).is_covariant(),
+						ex_to<spinidx>(*ip).is_dotted()))
+					-> setflag(status_flags::dynallocated);
+			else if (is_exactly_a<varidx>(*ip))
+				newidx = (new varidx(newsym, ex_to<varidx>(*ip).get_dim(),
+						ex_to<varidx>(*ip).is_covariant()))
+					-> setflag(status_flags::dynallocated);
+			else
+				newidx = (new idx(newsym, ex_to<idx>(*ip).get_dim()))
+					-> setflag(status_flags::dynallocated);
+			old_indices.push_back(*ip);
+			new_indices.push_back(newidx);
+			if(is_a<varidx>(*ip)) {
+				old_indices.push_back(ex_to<varidx>(*ip).toggle_variance());
+				new_indices.push_back(ex_to<varidx>(newidx).toggle_variance());
+			}
+			++ip;
+		}
+		return lst(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()));
+	}
+}
+
+ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)
+{
+	lst indices_subs = rename_dummy_indices_uniquely(va, vb);
+	return (indices_subs.op(0).nops()>0 ? b.subs(ex_to<lst>(indices_subs.op(0)), ex_to<lst>(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming) : b);
+}
+
+ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+{
+	exvector va = get_all_dummy_indices_safely(a);
+	if (va.size() > 0) {
+		exvector vb = get_all_dummy_indices_safely(b);
+		if (vb.size() > 0) {
+			sort(va.begin(), va.end(), ex_is_less());
+			sort(vb.begin(), vb.end(), ex_is_less());
+			lst indices_subs = rename_dummy_indices_uniquely(va, vb);
+			if (indices_subs.op(0).nops() > 0)
+				return b.subs(ex_to<lst>(indices_subs.op(0)), ex_to<lst>(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming);
+		}
+	}
+	return b;
+}
+
+ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va)
+{
+	if (va.size() > 0) {
+		exvector vb = get_all_dummy_indices_safely(b);
+		if (vb.size() > 0) {
+			sort(vb.begin(), vb.end(), ex_is_less());
+			lst indices_subs = rename_dummy_indices_uniquely(va, vb);
+			if (indices_subs.op(0).nops() > 0) {
+				if (modify_va) {
+					for (lst::const_iterator i = ex_to<lst>(indices_subs.op(1)).begin(); i != ex_to<lst>(indices_subs.op(1)).end(); ++i)
+						va.push_back(*i);
+					exvector uncommon_indices;
+					set_difference(vb.begin(), vb.end(), indices_subs.op(0).begin(), indices_subs.op(0).end(), std::back_insert_iterator<exvector>(uncommon_indices), ex_is_less());
+					exvector::const_iterator ip = uncommon_indices.begin(), ipend = uncommon_indices.end();
+					while (ip != ipend) {
+						va.push_back(*ip);
+						++ip;
+					}
+					sort(va.begin(), va.end(), ex_is_less());
+				}
+				return b.subs(ex_to<lst>(indices_subs.op(0)), ex_to<lst>(indices_subs.op(1)), subs_options::no_pattern|subs_options::no_index_renaming);
+			}
+		}
+	}
+	return b;
+}
+
+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;
+	}
+}
+
 } // namespace GiNaC