From 9ee1dc214b3d03747db93463dac9b3f0c63e1501 Mon Sep 17 00:00:00 2001
From: Chris Dams <Chris.Dams@mi.infn.it>
Date: Mon, 7 Aug 2006 16:13:43 +0000
Subject: [PATCH] Bug in expand_dummy_sum is fixed.

---
 doc/tutorial/ginac.texi |  6 +----
 ginac/container.h       |  1 +
 ginac/indexed.cpp       | 55 +++++++++++++++++------------------------
 ginac/indexed.h         |  2 +-
 ginac/inifcns.cpp       |  4 +--
 5 files changed, 27 insertions(+), 41 deletions(-)

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 7235248a..815064a9 100644
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -2674,11 +2674,7 @@ arithmetic class, you just pass it to @code{simplify_indexed()}):
 The @code{scalar_products} object @code{sp} acts as a storage for the
 scalar products added to it with the @code{.add()} method. This method
 takes three arguments: the two expressions of which the scalar product is
-taken, and the expression to replace it with. After @code{sp.add(A, B, 0)},
-@code{simplify_indexed()} will replace all scalar products of indexed
-objects that have the symbols @code{A} and @code{B} as base expressions
-with the single value 0. The number, type and dimension of the indices
-don't matter; @samp{A~mu~nu*B.mu.nu} would also be replaced by 0.
+taken, and the expression to replace it with.
 
 @cindex @code{expand()}
 The example above also illustrates a feature of the @code{expand()} method:
diff --git a/ginac/container.h b/ginac/container.h
index 21e7d4ce..57749cd5 100644
--- a/ginac/container.h
+++ b/ginac/container.h
@@ -129,6 +129,7 @@ template <template <class T, class = std::allocator<T> > class C>
 class container : public basic, public container_storage<C> {
 	GINAC_DECLARE_REGISTERED_CLASS(container, basic)
 
+protected:
 	typedef typename container_storage<C>::STLT STLT;
 
 public:
diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 15455980..05c3b90a 100644
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -1395,44 +1395,33 @@ ex expand_dummy_sum(const ex & e, bool subs_idx)
 	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()) {
+	} 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 < 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())));
+				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 += e_expanded.subs(lst(nu == idx(i, nu.get_dim()), nu.toggle_variance() == idx(i, nu.get_dim())));
+						en += result.subs( nu.op(0) == i );
 					}
 				}
-				return expand_dummy_sum(en, subs_idx);
-			} 
-			++it;
+				result = en;
+			}
 		}
-		return e;
+		return result;
 	} else {
 		return e;
 	}
diff --git a/ginac/indexed.h b/ginac/indexed.h
index ad32c237..8288441f 100644
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -258,7 +258,7 @@ ex rename_dummy_indices_uniquely(const ex & a, const ex & b);
 
 /** This function returns the given expression with expanded sums
  *  for all dummy index summations, where the dimensionality of 
- *  the dummy index is numeric.
+ *  the dummy index is a nonnegative integer.
  *  Optionally all indices with a variance will be substituted by 
  *  indices with the corresponding numeric values without variance.
  *
diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index e7ab635a..256d1315 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -589,7 +589,7 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
 	
 	// syntax checks
 	if (!eqns.info(info_flags::list)) {
-		throw(std::invalid_argument("lsolve(): 1st argument must be a list"));
+		throw(std::invalid_argument("lsolve(): 1st argument must be a list or an equation"));
 	}
 	for (size_t i=0; i<eqns.nops(); i++) {
 		if (!eqns.op(i).info(info_flags::relation_equal)) {
@@ -597,7 +597,7 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
 		}
 	}
 	if (!symbols.info(info_flags::list)) {
-		throw(std::invalid_argument("lsolve(): 2nd argument must be a list"));
+		throw(std::invalid_argument("lsolve(): 2nd argument must be a list or a symbol"));
 	}
 	for (size_t i=0; i<symbols.nops(); i++) {
 		if (!symbols.op(i).info(info_flags::symbol)) {
-- 
2.50.0