X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=e32f4c7a103f047cc09a5332665b7ef0e6c91ef8;hp=56fd0b6b4080ebaafa72845978bfea39d5198ae9;hb=e93f20fd0deb9b45d2163cbec247e5181f021a28;hpb=c4c7b40ea71dd20ecc8942f9cab50c4bb9766624

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 56fd0b6b..e32f4c7a 100644
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's indexed expressions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2002 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
@@ -20,8 +20,8 @@
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
+#include <iostream>
 #include <stdexcept>
-#include <algorithm>
 
 #include "indexed.h"
 #include "idx.h"
@@ -29,24 +29,23 @@
 #include "mul.h"
 #include "ncmul.h"
 #include "power.h"
+#include "relational.h"
 #include "symmetry.h"
 #include "lst.h"
 #include "print.h"
 #include "archive.h"
 #include "utils.h"
-#include "debugmsg.h"
 
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq)
 
 //////////
-// default constructor, destructor, copy constructor assignment operator and helpers
+// default ctor, dtor, copy ctor, assignment operator and helpers
 //////////
 
 indexed::indexed() : symtree(sy_none())
 {
-	debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 }
 
@@ -64,63 +63,54 @@ DEFAULT_DESTROY(indexed)
 
 indexed::indexed(const ex & b) : inherited(b), symtree(sy_none())
 {
-	debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(sy_none())
 {
-	debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(sy_none())
 {
-	debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	validate();
 }
 
 indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(sy_none())
 {
-	debugmsg("indexed constructor from ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	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())
 {
-	debugmsg("indexed constructor from ex,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	validate();
 }
 
 indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(symm)
 {
-	debugmsg("indexed constructor from ex,symmetry,ex,ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	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)
 {
-	debugmsg("indexed constructor from ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	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)
 {
-	debugmsg("indexed constructor from ex,symmetry,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 	validate();
 }
 
 indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_none())
 {
-	debugmsg("indexed constructor from ex,exvector", LOGLEVEL_CONSTRUCT);
 	seq.insert(seq.end(), v.begin(), v.end());
 	tinfo_key = TINFO_indexed;
 	validate();
@@ -128,7 +118,6 @@ indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_no
 
 indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm)
 {
-	debugmsg("indexed constructor from ex,symmetry,exvector", LOGLEVEL_CONSTRUCT);
 	seq.insert(seq.end(), v.begin(), v.end());
 	tinfo_key = TINFO_indexed;
 	validate();
@@ -136,19 +125,16 @@ indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inhe
 
 indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm)
 {
-	debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 }
 
 indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm)
 {
-	debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 }
 
 indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(symm)
 {
-	debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT);
 	tinfo_key = TINFO_indexed;
 }
 
@@ -158,7 +144,6 @@ indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(
 
 indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
-	debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT);
 	if (!n.find_ex("symmetry", symtree, sym_lst)) {
 		// GiNaC versions <= 0.9.0 had an unsigned "symmetry" property
 		unsigned symm = 0;
@@ -174,7 +159,7 @@ indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_l
 				symtree = sy_none();
 				break;
 		}
-		ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+		const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
 	}
 }
 
@@ -192,10 +177,9 @@ DEFAULT_UNARCHIVE(indexed)
 
 void indexed::print(const print_context & c, unsigned level) const
 {
-	debugmsg("indexed print", LOGLEVEL_PRINT);
 	GINAC_ASSERT(seq.size() > 0);
 
-	if (is_of_type(c, print_tree)) {
+	if (is_a<print_tree>(c)) {
 
 		c.s << std::string(level, ' ') << class_name()
 		    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
@@ -207,21 +191,19 @@ void indexed::print(const print_context & c, unsigned level) const
 
 	} else {
 
-		bool is_tex = is_of_type(c, print_latex);
+		bool is_tex = is_a<print_latex>(c);
 		const ex & base = seq[0];
-		bool need_parens = is_ex_exactly_of_type(base, add) || is_ex_exactly_of_type(base, mul)
-		                || is_ex_exactly_of_type(base, ncmul) || is_ex_exactly_of_type(base, power)
-		                || is_ex_of_type(base, indexed);
+
+		if (precedence() <= level)
+			c.s << (is_tex ? "{(" : "(");
 		if (is_tex)
 			c.s << "{";
-		if (need_parens)
-			c.s << "(";
-		base.print(c);
-		if (need_parens)
-			c.s << ")";
+		base.print(c, precedence());
 		if (is_tex)
 			c.s << "}";
 		printindices(c, level);
+		if (precedence() <= level)
+			c.s << (is_tex ? ")}" : ")");
 	}
 }
 
@@ -250,7 +232,7 @@ bool indexed::all_index_values_are(unsigned inf) const
 
 int indexed::compare_same_type(const basic & other) const
 {
-	GINAC_ASSERT(is_of_type(other, indexed));
+	GINAC_ASSERT(is_a<indexed>(other));
 	return inherited::compare_same_type(other);
 }
 
@@ -264,7 +246,7 @@ ex indexed::eval(int level) const
 
 	// If the base object is 0, the whole object is 0
 	if (base.is_zero())
-		return _ex0();
+		return _ex0;
 
 	// If the base object is a product, pull out the numeric factor
 	if (is_ex_exactly_of_type(base, mul) && is_ex_exactly_of_type(base.op(base.nops() - 1), numeric)) {
@@ -277,36 +259,18 @@ ex indexed::eval(int level) const
 	// Canonicalize indices according to the symmetry properties
 	if (seq.size() > 2) {
 		exvector v = seq;
-		GINAC_ASSERT(is_ex_exactly_of_type(symtree, symmetry));
+		GINAC_ASSERT(is_exactly_a<symmetry>(symtree));
 		int sig = canonicalize(v.begin() + 1, ex_to<symmetry>(symtree));
 		if (sig != INT_MAX) {
 			// Something has changed while sorting indices, more evaluations later
 			if (sig == 0)
-				return _ex0();
+				return _ex0;
 			return ex(sig) * thisexprseq(v);
 		}
 	}
 
 	// Let the class of the base object perform additional evaluations
-	return base.bp->eval_indexed(*this);
-}
-
-int indexed::degree(const ex & s) const
-{
-	return is_equal(*s.bp) ? 1 : 0;
-}
-
-int indexed::ldegree(const ex & s) const
-{
-	return is_equal(*s.bp) ? 1 : 0;
-}
-
-ex indexed::coeff(const ex & s, int n) const
-{
-	if (is_equal(*s.bp))
-		return n==1 ? _ex1() : _ex0();
-	else
-		return n==0 ? ex(*this) : _ex0();
+	return ex_to<basic>(base).eval_indexed(*this);
 }
 
 ex indexed::thisexprseq(const exvector & v) const
@@ -327,7 +291,7 @@ ex indexed::expand(unsigned options) const
 
 		// expand_indexed expands (a+b).i -> a.i + b.i
 		const ex & base = seq[0];
-		ex sum = _ex0();
+		ex sum = _ex0;
 		for (unsigned i=0; i<base.nops(); i++) {
 			exvector s = seq;
 			s[0] = base.op(i);
@@ -355,7 +319,7 @@ void indexed::printindices(const print_context & c, unsigned level) const
 
 		exvector::const_iterator it=seq.begin() + 1, itend = seq.end();
 
-		if (is_of_type(c, print_latex)) {
+		if (is_a<print_latex>(c)) {
 
 			// TeX output: group by variance
 			bool first = true;
@@ -363,9 +327,10 @@ void indexed::printindices(const print_context & c, unsigned level) const
 
 			while (it != itend) {
 				bool cur_covariant = (is_ex_of_type(*it, varidx) ? ex_to<varidx>(*it).is_covariant() : true);
-				if (first || cur_covariant != covariant) {
+				if (first || cur_covariant != covariant) { // Variance changed
+					// The empty {} prevents indices from ending up on top of each other
 					if (!first)
-						c.s << "}";
+						c.s << "}{}";
 					covariant = cur_covariant;
 					if (covariant)
 						c.s << "_{";
@@ -406,7 +371,7 @@ void indexed::validate(void) const
 	if (!symtree.is_zero()) {
 		if (!is_ex_exactly_of_type(symtree, symmetry))
 			throw(std::invalid_argument("symmetry of indexed object must be of type symmetry"));
-		ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+		const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
 	}
 }
 
@@ -415,7 +380,7 @@ void indexed::validate(void) const
  *  @see ex::diff */
 ex indexed::derivative(const symbol & s) const
 {
-	return _ex0();
+	return _ex0;
 }
 
 //////////
@@ -526,7 +491,7 @@ exvector power::get_free_indices(void) const
 
 /** Rename dummy indices in an expression.
  *
- *  @param e Expression to be worked on
+ *  @param e Expression to work on
  *  @param local_dummy_indices The set of dummy indices that appear in the
  *    expression "e"
  *  @param global_dummy_indices The set of dummy indices that have appeared
@@ -556,7 +521,6 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
 			}
 			it++;
 		}
-		shaker_sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less(), ex_swap());
 
 		// If this is the first set of local indices, do nothing
 		if (old_global_size == 0)
@@ -571,6 +535,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
 	shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
 	for (unsigned i=0; i<global_size; 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
 	exlist local_uniq, global_uniq;
@@ -587,6 +552,80 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex
 	}
 }
 
+/** Given a set of indices, extract those of class varidx. */
+static void find_variant_indices(const exvector & v, exvector & variant_indices)
+{
+	exvector::const_iterator it1, itend;
+	for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
+		if (is_exactly_a<varidx>(*it1))
+			variant_indices.push_back(*it1);
+	}
+}
+
+/** Raise/lower dummy indices in a single indexed objects to canonicalize their
+ *  variance.
+ *
+ *  @param e Object to work on
+ *  @param variant_dummy_indices The set of indices that might need repositioning (will be changed by this function)
+ *  @param moved_indices The set of indices that have been repositioned (will be changed by this function)
+ *  @return true if 'e' was changed */
+bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector & moved_indices)
+{
+	bool something_changed = false;
+
+	// 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) {
+		if (!is_exactly_a<varidx>(*it2))
+			continue;
+
+		exvector::iterator vit, vitend;
+		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
+					));
+					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);
+				goto next_index;
+			}
+		}
+
+		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());
+					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;
+			}
+		}
+
+next_index: ;
+	}
+
+	return something_changed;
+}
+
+/* Ordering that only compares the base expressions of indexed objects. */
+struct ex_base_is_less : public std::binary_function<ex, ex, bool> {
+	bool operator() (const ex &lh, const ex &rh) const
+	{
+		return (is_a<indexed>(lh) ? lh.op(0) : lh).compare(is_a<indexed>(rh) ? rh.op(0) : rh) < 0;
+	}
+};
+
 /** 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)
@@ -601,13 +640,13 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du
 
 	if (is_ex_exactly_of_type(e, power)) {
 		// We only get called for simple squares, split a^2 -> a*a
-		GINAC_ASSERT(e.op(1).is_equal(_ex2()));
+		GINAC_ASSERT(e.op(1).is_equal(_ex2));
 		v.push_back(e.op(0));
 		v.push_back(e.op(0));
 	} else {
 		for (unsigned i=0; i<e.nops(); i++) {
 			ex f = e.op(i);
-			if (is_ex_exactly_of_type(f, power) && f.op(1).is_equal(_ex2())) {
+			if (is_ex_exactly_of_type(f, power) && f.op(1).is_equal(_ex2)) {
 				v.push_back(f.op(0));
 	            v.push_back(f.op(0));
 			} else if (is_ex_exactly_of_type(f, ncmul)) {
@@ -663,41 +702,18 @@ try_again:
 			if (free.empty()) {
 				if (sp.is_defined(*it1, *it2)) {
 					*it1 = sp.evaluate(*it1, *it2);
-					*it2 = _ex1();
+					*it2 = _ex1;
 					goto contraction_done;
 				}
 			}
 
-			// Contraction of symmetric with antisymmetric object is zero
-			if (num_dummies > 1
-			 && ex_to<symmetry>(ex_to<indexed>(*it1).symtree).has_symmetry()
-			 && ex_to<symmetry>(ex_to<indexed>(*it2).symtree).has_symmetry()) {
-
-				// Check all pairs of dummy indices
-				for (unsigned idx1=0; idx1<num_dummies-1; idx1++) {
-					for (unsigned idx2=idx1+1; idx2<num_dummies; idx2++) {
-
-						// Try and swap the index pair and check whether the
-						// relative sign changed
-						lst subs_lst(dummy[idx1].op(0), dummy[idx2].op(0)), repl_lst(dummy[idx2].op(0), dummy[idx1].op(0));
-						ex swapped1 = it1->subs(subs_lst, repl_lst);
-						ex swapped2 = it2->subs(subs_lst, repl_lst);
-						if (it1->is_equal(swapped1) && it2->is_equal(-swapped2)
-						 || it1->is_equal(-swapped1) && it2->is_equal(swapped2)) {
-							free_indices.clear();
-							return _ex0();
-						}
-					}
-				}
-			}
-
 			// Try to contract the first one with the second one
-			contracted = it1->op(0).bp->contract_with(it1, it2, v);
+			contracted = ex_to<basic>(it1->op(0)).contract_with(it1, it2, v);
 			if (!contracted) {
 
 				// That didn't work; maybe the second object knows how to
 				// contract itself with the first one
-				contracted = it2->op(0).bp->contract_with(it2, it1, v);
+				contracted = ex_to<basic>(it2->op(0)).contract_with(it2, it1, v);
 			}
 			if (contracted) {
 contraction_done:
@@ -727,8 +743,7 @@ contraction_done:
 	// Find free indices (concatenate them all and call find_free_and_dummy())
 	// and all dummy indices that appear
 	exvector un, individual_dummy_indices;
-	it1 = v.begin(); itend = v.end();
-	while (it1 != itend) {
+	for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
 		exvector free_indices_of_factor;
 		if (is_ex_of_type(*it1, indexed)) {
 			exvector dummy_indices_of_factor;
@@ -737,29 +752,127 @@ contraction_done:
 		} else
 			free_indices_of_factor = it1->get_free_indices();
 		un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
-		it1++;
 	}
 	exvector local_dummy_indices;
 	find_free_and_dummy(un, free_indices, local_dummy_indices);
 	local_dummy_indices.insert(local_dummy_indices.end(), individual_dummy_indices.begin(), individual_dummy_indices.end());
 
+	// Filter out the dummy indices with variance
+	exvector variant_dummy_indices;
+	find_variant_indices(local_dummy_indices, variant_dummy_indices);
+
+	// Any indices with variance present at all?
+	if (!variant_dummy_indices.empty()) {
+
+		// Yes, bring the product into a canonical order that only depends on
+		// the base expressions of indexed objects
+		if (!non_commutative)
+			std::sort(v.begin(), v.end(), ex_base_is_less());
+
+		exvector moved_indices;
+
+		// Iterate over all indexed objects in the product
+		for (it1 = v.begin(), itend = v.end(); it1 != itend; ++it1) {
+			if (!is_ex_of_type(*it1, indexed))
+				continue;
+
+			if (reposition_dummy_indices(*it1, variant_dummy_indices, moved_indices))
+				something_changed = true;
+		}
+	}
+
 	ex r;
 	if (something_changed)
 		r = non_commutative ? ex(ncmul(v, true)) : ex(mul(v));
 	else
 		r = e;
 
+	// 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) {
+		lst dummy_syms;
+		for (int i=0; i<local_dummy_indices.size(); i++)
+			dummy_syms.append(local_dummy_indices[i].op(0));
+		if (r.symmetrize(dummy_syms).is_zero()) {
+			free_indices.clear();
+			return _ex0;
+		}
+	}
+
 	// Dummy index renaming
 	r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
 
 	// Product of indexed object with a scalar?
 	if (is_ex_exactly_of_type(r, mul) && r.nops() == 2
 	 && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed))
-		return r.op(0).op(0).bp->scalar_mul_indexed(r.op(0), ex_to<numeric>(r.op(1)));
+		return ex_to<basic>(r.op(0).op(0)).scalar_mul_indexed(r.op(0), ex_to<numeric>(r.op(1)));
 	else
 		return r;
 }
 
+/** This structure stores the original and symmetrized versions of terms
+ *  obtained during the simplification of sums. */
+class symminfo {
+public:
+	symminfo() {}
+	~symminfo() {}
+
+	symminfo(const ex & symmterm_, const ex & orig_)
+	{
+		if (is_a<mul>(orig_) && is_a<numeric>(orig_.op(orig_.nops()-1))) {
+			ex tmp = orig_.op(orig_.nops()-1);
+			orig = orig_ / tmp;
+		} else 
+			orig = orig_;
+
+		if (is_a<mul>(symmterm_) && is_a<numeric>(symmterm_.op(symmterm_.nops()-1))) {
+			coeff = symmterm_.op(symmterm_.nops()-1);
+			symmterm = symmterm_ / coeff;
+		} else {
+			coeff = 1;
+			symmterm = symmterm_;
+		}
+	}
+
+	symminfo(const symminfo & other)
+	{
+		symmterm = other.symmterm;
+		coeff = other.coeff;
+		orig = other.orig;
+	}
+
+	const symminfo & operator=(const symminfo & other)
+	{
+		if (this != &other) {
+			symmterm = other.symmterm;
+			coeff = other.coeff;
+			orig = other.orig;
+		}
+		return *this;
+	}
+
+	ex symmterm;
+	ex coeff;
+	ex orig;
+};
+
+class symminfo_is_less {
+public:
+	bool operator() (const symminfo & si1, const symminfo & si2)
+	{
+		int comp = si1.symmterm.compare(si2.symmterm);
+		if (comp < 0) return true;
+		if (comp > 0) return false;
+		comp = si1.orig.compare(si2.orig);
+		if (comp < 0) return true;
+		if (comp > 0) return false;
+		comp = si1.coeff.compare(si2.coeff);
+		if (comp < 0) 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)
 {
@@ -767,11 +880,27 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
 	ex e_expanded = e.expand();
 
 	// Simplification of single indexed object: just find the free indices
-	// and perform dummy index renaming
+	// and perform dummy index renaming/repositioning
 	if (is_ex_of_type(e_expanded, indexed)) {
+
+		// Find the dummy indices
 		const indexed &i = ex_to<indexed>(e_expanded);
 		exvector local_dummy_indices;
 		find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, local_dummy_indices);
+
+		// Filter out the dummy indices with variance
+		exvector variant_dummy_indices;
+		find_variant_indices(local_dummy_indices, variant_dummy_indices);
+
+		// Any indices with variance present at all?
+		if (!variant_dummy_indices.empty()) {
+
+			// Yes, reposition them
+			exvector moved_indices;
+			reposition_dummy_indices(e_expanded, variant_dummy_indices, moved_indices);
+		}
+
+		// Rename the dummy indices
 		return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);
 	}
 
@@ -779,7 +908,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
 	// free indices in each term
 	if (is_ex_exactly_of_type(e_expanded, add)) {
 		bool first = true;
-		ex sum = _ex0();
+		ex sum = _ex0;
 		free_indices.clear();
 
 		for (unsigned i=0; i<e_expanded.nops(); i++) {
@@ -794,20 +923,64 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
 					if (!indices_consistent(free_indices, free_indices_of_term))
 						throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
 					if (is_ex_of_type(sum, indexed) && is_ex_of_type(term, indexed))
-						sum = sum.op(0).bp->add_indexed(sum, term);
+						sum = ex_to<basic>(sum.op(0)).add_indexed(sum, term);
 					else
 						sum += term;
 				}
 			}
 		}
 
+		// If the sum turns out to be zero, we are finished
+		if (sum.is_zero()) {
+			free_indices.clear();
+			return sum;
+		}
+
+		// Symmetrizing over the dummy indices may cancel terms
+		int num_terms_orig = (is_a<add>(sum) ? sum.nops() : 1);
+		if (num_terms_orig > 1 && dummy_indices.size() >= 2) {
+
+			// Construct list of all dummy index symbols
+			lst dummy_syms;
+			for (int i=0; i<dummy_indices.size(); i++)
+				dummy_syms.append(dummy_indices[i].op(0));
+
+			// Symmetrize each term separately and store the resulting
+			// terms in a list of symminfo structures
+			std::vector<symminfo> v;
+			for (int i=0; i<sum.nops(); i++) {
+				ex sum_symm = sum.op(i).symmetrize(dummy_syms);
+				if (is_a<add>(sum_symm))
+					for (int j=0; j<sum_symm.nops(); j++)
+						v.push_back(symminfo(sum_symm.op(j), sum.op(i)));
+				else
+					v.push_back(symminfo(sum_symm, sum.op(i)));
+			}
+
+			// Now add up all the unsymmetrized versions of the terms that
+			// did not cancel out in the symmetrization
+			exvector result;
+			std::sort(v.begin(), v.end(), symminfo_is_less());
+			for (std::vector<symminfo>::iterator i=v.begin(); i!=v.end(); ) {
+				std::vector<symminfo>::iterator j = i;
+				for (j++; j!=v.end() && i->symmterm == j->symmterm; j++) ;
+				for (std::vector<symminfo>::iterator k=i; k!=j; k++)
+					result.push_back((k->coeff)*(i->orig));
+				i = j;
+			}
+			ex sum_symm = (new add(result))->setflag(status_flags::dynallocated);
+			if (sum_symm.is_zero())
+				free_indices.clear();
+			return sum_symm;
+		}
+
 		return sum;
 	}
 
 	// Simplification of products
 	if (is_ex_exactly_of_type(e_expanded, mul)
 	 || is_ex_exactly_of_type(e_expanded, ncmul)
-	 || (is_ex_exactly_of_type(e_expanded, power) && is_ex_of_type(e_expanded.op(0), indexed) && e_expanded.op(1).is_equal(_ex2())))
+	 || (is_ex_exactly_of_type(e_expanded, power) && is_ex_of_type(e_expanded.op(0), indexed) && e_expanded.op(1).is_equal(_ex2)))
 		return simplify_indexed_product(e_expanded, free_indices, dummy_indices, sp);
 
 	// Cannot do anything