X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=c8b710ae7f1e63690cbb0356e2e1ccc10159f1e0;hp=c899e13cd7f237580b11ae8becb5bb9474b7645c;hb=6850c16435de33396ccd505ff26b4d4574fa74cf;hpb=68fdf425abf14d016d5f95ee7b9d06a19a3c5926

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index c899e13c..c8b710ae 100644
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2008 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,12 +17,13 @@
  *
  *  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>
 #include <sstream>
 #include <stdexcept>
+#include <limits>
 
 #include "indexed.h"
 #include "idx.h"
@@ -34,102 +35,108 @@
 #include "symmetry.h"
 #include "operators.h"
 #include "lst.h"
-#include "print.h"
 #include "archive.h"
+#include "symbol.h"
 #include "utils.h"
+#include "integral.h"
+#include "matrix.h"
+#include "inifcns.h"
 
 namespace GiNaC {
 
-GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq)
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(indexed, exprseq,
+  print_func<print_context>(&indexed::do_print).
+  print_func<print_latex>(&indexed::do_print_latex).
+  print_func<print_tree>(&indexed::do_print_tree))
 
 //////////
 // default constructor
 //////////
 
-indexed::indexed() : symtree(sy_none())
+indexed::indexed() : symtree(not_symmetric())
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 }
 
 //////////
 // other constructors
 //////////
 
-indexed::indexed(const ex & b) : inherited(b), symtree(sy_none())
+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(sy_none())
+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(sy_none())
+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(sy_none())
+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(sy_none())
+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(sy_none())
+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, exvector * vp) : inherited(vp), symtree(symm)
+indexed::indexed(const symmetry & symm, std::auto_ptr<exvector> vp) : inherited(vp), symtree(symm)
 {
-	tinfo_key = TINFO_indexed;
+	tinfo_key = &indexed::tinfo_static;
 }
 
 //////////
@@ -150,7 +157,7 @@ indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
 				symtree = sy_anti();
 				break;
 			default:
-				symtree = sy_none();
+				symtree = not_symmetric();
 				break;
 		}
 		const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
@@ -169,38 +176,80 @@ DEFAULT_UNARCHIVE(indexed)
 // functions overriding virtual functions from base classes
 //////////
 
-void indexed::print(const print_context & c, unsigned level) const
+void indexed::printindices(const print_context & c, unsigned level) const
 {
-	GINAC_ASSERT(seq.size() > 0);
-
-	if (is_a<print_tree>(c)) {
+	if (seq.size() > 1) {
 
-		c.s << std::string(level, ' ') << class_name()
-		    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
-		    << ", " << seq.size()-1 << " indices"
-		    << ", symmetry=" << symtree << std::endl;
-		unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
-		seq[0].print(c, level + delta_indent);
-		printindices(c, level + delta_indent);
+		exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
 
-	} else {
+		if (is_a<print_latex>(c)) {
 
-		bool is_tex = is_a<print_latex>(c);
-		const ex & base = seq[0];
+			// TeX output: group by variance
+			bool first = true;
+			bool covariant = true;
 
-		if (precedence() <= level)
-			c.s << (is_tex ? "{(" : "(");
-		if (is_tex)
-			c.s << "{";
-		base.print(c, precedence());
-		if (is_tex)
+			while (it != itend) {
+				bool cur_covariant = (is_a<varidx>(*it) ? ex_to<varidx>(*it).is_covariant() : true);
+				if (first || cur_covariant != covariant) { // Variance changed
+					// The empty {} prevents indices from ending up on top of each other
+					if (!first)
+						c.s << "}{}";
+					covariant = cur_covariant;
+					if (covariant)
+						c.s << "_{";
+					else
+						c.s << "^{";
+				}
+				it->print(c, level);
+				c.s << " ";
+				first = false;
+				it++;
+			}
 			c.s << "}";
-		printindices(c, level);
-		if (precedence() <= level)
-			c.s << (is_tex ? ")}" : ")");
+
+		} else {
+
+			// Ordinary output
+			while (it != itend) {
+				it->print(c, level);
+				it++;
+			}
+		}
 	}
 }
 
+void indexed::print_indexed(const print_context & c, const char *openbrace, const char *closebrace, unsigned level) const
+{
+	if (precedence() <= level)
+		c.s << openbrace << '(';
+	c.s << openbrace;
+	seq[0].print(c, precedence());
+	c.s << closebrace;
+	printindices(c, level);
+	if (precedence() <= level)
+		c.s << ')' << closebrace;
+}
+
+void indexed::do_print(const print_context & c, unsigned level) const
+{
+	print_indexed(c, "", "", level);
+}
+
+void indexed::do_print_latex(const print_latex & c, unsigned level) const
+{
+	print_indexed(c, "{", "}", level);
+}
+
+void indexed::do_print_tree(const print_tree & c, unsigned level) const
+{
+	c.s << std::string(level, ' ') << class_name() << " @" << this
+	    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
+	    << ", " << seq.size()-1 << " indices"
+	    << ", symmetry=" << symtree << std::endl;
+	seq[0].print(c, level + c.delta_indent);
+	printindices(c, level + c.delta_indent);
+}
+
 bool indexed::info(unsigned inf) const
 {
 	if (inf == info_flags::indexed) return true;
@@ -250,12 +299,15 @@ 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;
 		GINAC_ASSERT(is_exactly_a<symmetry>(symtree));
 		int sig = canonicalize(v.begin() + 1, ex_to<symmetry>(symtree));
-		if (sig != INT_MAX) {
+		if (sig != std::numeric_limits<int>::max()) {
 			// Something has changed while sorting indices, more evaluations later
 			if (sig == 0)
 				return _ex0;
@@ -267,34 +319,60 @@ 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);
 }
 
-ex indexed::thiscontainer(exvector * vp) const
+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);
 
-	if ((options & expand_options::expand_indexed) && is_exactly_a<add>(seq[0])) {
-
-		// expand_indexed expands (a+b).i -> a.i + b.i
-		const ex & base = seq[0];
-		ex sum = _ex0;
-		for (size_t i=0; i<base.nops(); i++) {
+	if (options & expand_options::expand_indexed) {
+		ex newbase = seq[0].expand(options);
+		if (is_exactly_a<add>(newbase)) {
+			ex sum = _ex0;
+			for (size_t i=0; i<newbase.nops(); i++) {
+				exvector s = seq;
+				s[0] = newbase.op(i);
+				sum += thiscontainer(s).expand(options);
+			}
+			return sum;
+		}
+		if (!are_ex_trivially_equal(newbase, seq[0])) {
 			exvector s = seq;
-			s[0] = base.op(i);
-			sum += thiscontainer(s).expand();
+			s[0] = newbase;
+			return ex_to<indexed>(thiscontainer(s)).inherited::expand(options);
 		}
-		return sum;
-
-	} else
-		return inherited::expand(options);
+	}
+	return inherited::expand(options);
 }
 
 //////////
@@ -307,48 +385,6 @@ ex indexed::expand(unsigned options) const
 // non-virtual functions in this class
 //////////
 
-void indexed::printindices(const print_context & c, unsigned level) const
-{
-	if (seq.size() > 1) {
-
-		exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
-
-		if (is_a<print_latex>(c)) {
-
-			// TeX output: group by variance
-			bool first = true;
-			bool covariant = true;
-
-			while (it != itend) {
-				bool cur_covariant = (is_a<varidx>(*it) ? ex_to<varidx>(*it).is_covariant() : true);
-				if (first || cur_covariant != covariant) { // Variance changed
-					// The empty {} prevents indices from ending up on top of each other
-					if (!first)
-						c.s << "}{}";
-					covariant = cur_covariant;
-					if (covariant)
-						c.s << "_{";
-					else
-						c.s << "^{";
-				}
-				it->print(c, level);
-				c.s << " ";
-				first = false;
-				it++;
-			}
-			c.s << "}";
-
-		} else {
-
-			// Ordinary output
-			while (it != itend) {
-				it->print(c, level);
-				it++;
-			}
-		}
-	}
-}
-
 /** Check whether all indices are of class idx and validate the symmetry
  *  tree. This function is used internally to make sure that all constructed
  *  indexed objects really carry indices and not some other classes. */
@@ -493,10 +529,27 @@ exvector ncmul::get_free_indices() const
 	return free_indices;
 }
 
-exvector power::get_free_indices() const
+struct is_summation_idx : public std::unary_function<ex, bool> {
+	bool operator()(const ex & e)
+	{
+		return is_dummy_pair(e, e);
+	}
+};
+
+exvector integral::get_free_indices() const
+{
+	if (a.get_free_indices().size() || b.get_free_indices().size())
+		throw (std::runtime_error("integral::get_free_indices: boundary values should not have free indices"));
+	return f.get_free_indices();
+}
+
+template<class T> size_t number_of_type(const exvector&v)
 {
-	// Return free indices of basis
-	return basis.get_free_indices();
+	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.
@@ -507,10 +560,10 @@ exvector power::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)
@@ -524,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<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--;
@@ -542,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_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
@@ -560,7 +615,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
 	else {
 		while (global_uniq.size() > local_uniq.size())
 			global_uniq.pop_back();
-		return e.subs(lst(local_uniq.begin(), local_uniq.end()), lst(global_uniq.begin(), global_uniq.end()));
+		return e.subs(lst(local_uniq.begin(), local_uniq.end()), lst(global_uniq.begin(), global_uniq.end()), subs_options::no_pattern);
 	}
 }
 
@@ -585,10 +640,60 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 {
 	bool something_changed = false;
 
+	// Find dummy symbols that occur twice in the same indexed object.
+	exvector local_var_dummies;
+	local_var_dummies.reserve(e.nops()/2);
+	for (size_t i=1; i<e.nops(); ++i) {
+		if (!is_a<varidx>(e.op(i)))
+			continue;
+		for (size_t j=i+1; j<e.nops(); ++j) {
+			if (is_dummy_pair(e.op(i), e.op(j))) {
+				local_var_dummies.push_back(e.op(i));
+				for (exvector::iterator k = variant_dummy_indices.begin();
+						k!=variant_dummy_indices.end(); ++k) {
+					if (e.op(i).op(0) == k->op(0)) {
+						variant_dummy_indices.erase(k);
+						break;
+					}
+				}
+				break;
+			}
+		}
+	}
+
+	// In the case where a dummy symbol occurs twice in the same indexed object
+	// we try all posibilities of raising/lowering and keep the least one in
+	// the sense of ex_is_less.
+	ex optimal_e = e;
+	size_t numpossibs = 1 << local_var_dummies.size();
+	for (size_t i=0; i<numpossibs; ++i) {
+		ex try_e = e;
+		for (size_t j=0; j<local_var_dummies.size(); ++j) {
+			exmap m;
+			if (1<<j & i) {
+				ex curr_idx = local_var_dummies[j];
+				ex curr_toggle = ex_to<varidx>(curr_idx).toggle_variance();
+				m[curr_idx] = curr_toggle;
+				m[curr_toggle] = curr_idx;
+			}
+			try_e = e.subs(m, subs_options::no_pattern);
+		}
+		if(ex_is_less()(try_e, optimal_e))
+		{	optimal_e = try_e;
+			something_changed = true;
+		}
+	}
+	e = optimal_e;
+
+	if (!is_a<indexed>(e))
+		return true;
+
+	exvector seq = ex_to<indexed>(e).seq;
+
 	// If a dummy index is encountered for the first time in the
 	// product, pull it up, otherwise, pull it down
-	exvector::const_iterator it2, it2start, it2end;
-	for (it2start = ex_to<indexed>(e).seq.begin(), it2end = ex_to<indexed>(e).seq.end(), it2 = it2start + 1; it2 != it2end; ++it2) {
+	for (exvector::iterator it2 = seq.begin()+1, it2end = seq.end();
+			it2 != it2end; ++it2) {
 		if (!is_exactly_a<varidx>(*it2))
 			continue;
 
@@ -596,14 +701,20 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 		for (vit = variant_dummy_indices.begin(), vitend = variant_dummy_indices.end(); vit != vitend; ++vit) {
 			if (it2->op(0).is_equal(vit->op(0))) {
 				if (ex_to<varidx>(*it2).is_covariant()) {
-					e = e.subs(lst(
-						*it2 == ex_to<varidx>(*it2).toggle_variance(),
-						ex_to<varidx>(*it2).toggle_variance() == *it2
-					));
+					/*
+					 * N.B. we don't want to use
+					 *
+					 *  e = e.subs(lst(
+					 *  *it2 == ex_to<varidx>(*it2).toggle_variance(),
+					 *  ex_to<varidx>(*it2).toggle_variance() == *it2
+					 *  ), subs_options::no_pattern);
+					 *
+					 * since this can trigger non-trivial repositioning of indices,
+					 * e.g. due to non-trivial symmetry properties of e, thus
+					 * invalidating iterators
+					 */
+					*it2 = ex_to<varidx>(*it2).toggle_variance();
 					something_changed = true;
-					it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start);
-					it2start = ex_to<indexed>(e).seq.begin();
-					it2end = ex_to<indexed>(e).seq.end();
 				}
 				moved_indices.push_back(*vit);
 				variant_dummy_indices.erase(vit);
@@ -614,11 +725,8 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 		for (vit = moved_indices.begin(), vitend = moved_indices.end(); vit != vitend; ++vit) {
 			if (it2->op(0).is_equal(vit->op(0))) {
 				if (ex_to<varidx>(*it2).is_contravariant()) {
-					e = e.subs(*it2 == ex_to<varidx>(*it2).toggle_variance());
+					*it2 = ex_to<varidx>(*it2).toggle_variance();
 					something_changed = true;
-					it2 = ex_to<indexed>(e).seq.begin() + (it2 - it2start);
-					it2start = ex_to<indexed>(e).seq.begin();
-					it2end = ex_to<indexed>(e).seq.end();
 				}
 				goto next_index;
 			}
@@ -627,6 +735,9 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector
 next_index: ;
 	}
 
+	if (something_changed)
+		e = ex_to<indexed>(e).thiscontainer(seq);
+
 	return something_changed;
 }
 
@@ -638,16 +749,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)) {
@@ -660,7 +770,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
@@ -670,6 +780,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;
@@ -711,20 +849,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)) {
@@ -819,19 +949,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
@@ -898,6 +1035,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)
 {
@@ -926,7 +1074,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
@@ -970,18 +1121,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));
@@ -1105,8 +1257,9 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
  *  performs contraction of dummy indices where possible and checks whether
  *  the free indices in sums are consistent.
  *
+ *  @param options Simplification options (currently unused)
  *  @return simplified expression */
-ex ex::simplify_indexed() const
+ex ex::simplify_indexed(unsigned options) const
 {
 	exvector free_indices, dummy_indices;
 	scalar_products sp;
@@ -1119,8 +1272,9 @@ ex ex::simplify_indexed() const
  *  scalar products by known values if desired.
  *
  *  @param sp Scalar products to be replaced automatically
+ *  @param options Simplification options (currently unused)
  *  @return simplified expression */
-ex ex::simplify_indexed(const scalar_products & sp) const
+ex ex::simplify_indexed(const scalar_products & sp, unsigned options) const
 {
 	exvector free_indices, dummy_indices;
 	return GiNaC::simplify_indexed(*this, free_indices, dummy_indices, sp);
@@ -1245,4 +1399,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