From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Fri, 11 Jul 2003 19:28:43 +0000 (+0000)
Subject: subs() and normal() use maps instead of lists, resulting in a huge performance
X-Git-Tag: release_1-0-15~27
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=5ef801553eb39aed7bd2df9dd1aff9d752c3ea9d

subs() and normal() use maps instead of lists, resulting in a huge performance
boost for subs()
---

diff --git a/ginac/add.h b/ginac/add.h
index f2a94ed4..20faea8b 100644
--- a/ginac/add.h
+++ b/ginac/add.h
@@ -54,7 +54,7 @@ public:
 	ex eval(int level=0) const;
 	ex evalm() const;
 	ex series(const relational & r, int order, unsigned options = 0) const;
-	ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
+	ex normal(exmap & repl, int level=0) const;
 	numeric integer_content() const;
 	ex smod(const numeric &xi) const;
 	numeric max_coefficient() const;
diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index 595d6450..a296cf9a 100644
--- a/ginac/basic.cpp
+++ b/ginac/basic.cpp
@@ -524,22 +524,19 @@ bool basic::match(const ex & pattern, lst & repl_lst) const
 }
 
 /** Helper function for subs(). Does not recurse into subexpressions. */
-ex basic::subs_one_level(const lst & ls, const lst & lr, unsigned options) const
+ex basic::subs_one_level(const exmap & m, unsigned options) const
 {
-	GINAC_ASSERT(ls.nops() == lr.nops());
-
-	lst::const_iterator its, itr;
+	exmap::const_iterator it;
 
 	if (options & subs_options::subs_no_pattern) {
-		for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
-			if (is_equal(ex_to<basic>(*its)))
-				return *itr;
-		}
+		it = m.find(*this);
+		if (it != m.end())
+			return it->second;
 	} else {
-		for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+		for (it = m.begin(); it != m.end(); ++it) {
 			lst repl_lst;
-			if (match(ex_to<basic>(*its), repl_lst))
-				return itr->subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
+			if (match(ex_to<basic>(it->first), repl_lst))
+				return it->second.subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
 		}
 	}
 
@@ -548,7 +545,7 @@ ex basic::subs_one_level(const lst & ls, const lst & lr, unsigned options) const
 
 /** Substitute a set of objects by arbitrary expressions. The ex returned
  *  will already be evaluated. */
-ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
+ex basic::subs(const exmap & m, unsigned options) const
 {
 	size_t num = nops();
 	if (num) {
@@ -556,7 +553,7 @@ ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
 		// Substitute in subexpressions
 		for (size_t i=0; i<num; i++) {
 			const ex & orig_op = op(i);
-			const ex & subsed_op = orig_op.subs(ls, lr, options);
+			const ex & subsed_op = orig_op.subs(m, options);
 			if (!are_ex_trivially_equal(orig_op, subsed_op)) {
 
 				// Something changed, clone the object
@@ -569,16 +566,16 @@ ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
 
 				// Substitute the other operands
 				for (; i<num; i++)
-					copy->let_op(i) = op(i).subs(ls, lr, options);
+					copy->let_op(i) = op(i).subs(m, options);
 
 				// Perform substitutions on the new object as a whole
-				return copy->subs_one_level(ls, lr, options);
+				return copy->subs_one_level(m, options);
 			}
 		}
 	}
 
 	// Nothing changed or no subexpressions
-	return subs_one_level(ls, lr, options);
+	return subs_one_level(m, options);
 }
 
 /** Default interface of nth derivative ex::diff(s, n).  It should be called
@@ -737,35 +734,6 @@ ex basic::expand(unsigned options) const
 
 // public
 
-/** Substitute objects in an expression (syntactic substitution) and return
- *  the result as a new expression.  There are two valid types of
- *  replacement arguments: 1) a relational like object==ex and 2) a list of
- *  relationals lst(object1==ex1,object2==ex2,...), which is converted to
- *  subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
-ex basic::subs(const ex & e, unsigned options) const
-{
-	if (e.info(info_flags::relation_equal)) {
-		return subs(lst(e.lhs()), lst(e.rhs()), options);
-	}
-	if (!e.info(info_flags::list)) {
-		throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
-	}
-
-	// Split list into two
-	lst ls;
-	lst lr;
-	GINAC_ASSERT(is_a<lst>(e));
-	for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
-		ex r = *it;
-		if (!r.info(info_flags::relation_equal)) {
-			throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
-		}
-		ls.append(r.op(0));
-		lr.append(r.op(1));
-	}
-	return subs(ls, lr, options);
-}
-
 /** Compare objects syntactically to establish canonical ordering.
  *  All compare functions return: -1 for *this less than other, 0 equal,
  *  1 greater. */
diff --git a/ginac/basic.h b/ginac/basic.h
index a68c2eee..a5af3236 100644
--- a/ginac/basic.h
+++ b/ginac/basic.h
@@ -25,6 +25,7 @@
 
 #include <cstddef> // for size_t
 #include <vector>
+#include <map>
 // CINT needs <algorithm> to work properly with <vector>
 #include <algorithm>
 
@@ -45,6 +46,7 @@ class print_context;
 template <class> class ptr;
 
 typedef std::vector<ex> exvector;
+typedef std::map<ex, ex, ex_is_less> exmap;
 
 
 /** Function object for map(). */
@@ -135,7 +137,7 @@ protected:
 public:
 
 	// substitutions
-	virtual ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+	virtual ex subs(const exmap & m, unsigned options = 0) const;
 
 	// function mapping
 	virtual ex map(map_function & f) const;
@@ -163,7 +165,7 @@ public:
 	virtual ex series(const relational & r, int order, unsigned options = 0) const;
 
 	// rational functions
-	virtual ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+	virtual ex normal(exmap & repl, int level = 0) const;
 	virtual ex to_rational(lst &repl_lst) const;
 	virtual ex to_polynomial(lst &repl_lst) const;
 
@@ -190,8 +192,7 @@ protected: // functions that should be called from class ex only
 	
 	// non-virtual functions in this class
 public:
-	ex subs(const ex & e, unsigned options = 0) const;
-	ex subs_one_level(const lst & ls, const lst & lr, unsigned options) const;
+	ex subs_one_level(const exmap & m, unsigned options) const;
 	ex diff(const symbol & s, unsigned nth = 1) const;
 	int compare(const basic & other) const;
 	bool is_equal(const basic & other) const;
diff --git a/ginac/container.h b/ginac/container.h
index 135adb96..9a40f2ed 100644
--- a/ginac/container.h
+++ b/ginac/container.h
@@ -275,7 +275,7 @@ public:
 	ex op(size_t i) const;
 	ex & let_op(size_t i);
 	ex eval(int level = 0) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const;
 
 protected:
 	bool is_equal_same_type(const basic & other) const;
@@ -328,7 +328,7 @@ public:
 
 protected:
 	STLT evalchildren(int level) const;
-	STLT *subschildren(const lst & ls, const lst & lr, unsigned options = 0) const;
+	STLT *subschildren(const exmap & m, unsigned options = 0) const;
 };
 
 /** Default constructor */
@@ -424,13 +424,13 @@ ex container<C>::eval(int level) const
 }
 
 template <template <class> class C>
-ex container<C>::subs(const lst & ls, const lst & lr, unsigned options) const
+ex container<C>::subs(const exmap & m, unsigned options) const
 {
-	STLT *vp = subschildren(ls, lr, options);
+	STLT *vp = subschildren(m, options);
 	if (vp)
-		return ex_to<basic>(thiscontainer(vp)).subs_one_level(ls, lr, options);
+		return ex_to<basic>(thiscontainer(vp)).subs_one_level(m, options);
 	else
-		return subs_one_level(ls, lr, options);
+		return subs_one_level(m, options);
 }
 
 /** Compare two containers of the same type. */
@@ -587,7 +587,7 @@ typename container<C>::STLT container<C>::evalchildren(int level) const
 }
 
 template <template <class> class C>
-typename container<C>::STLT *container<C>::subschildren(const lst & ls, const lst & lr, unsigned options) const
+typename container<C>::STLT *container<C>::subschildren(const exmap & m, unsigned options) const
 {
 	// returns a NULL pointer if nothing had to be substituted
 	// returns a pointer to a newly created epvector otherwise
@@ -595,7 +595,7 @@ typename container<C>::STLT *container<C>::subschildren(const lst & ls, const ls
 
 	const_iterator cit = seq.begin(), end = seq.end();
 	while (cit != end) {
-		const ex & subsed_ex = cit->subs(ls, lr, options);
+		const ex & subsed_ex = cit->subs(m, options);
 		if (!are_ex_trivially_equal(*cit, subsed_ex)) {
 
 			// copy first part of seq which hasn't changed
@@ -608,7 +608,7 @@ typename container<C>::STLT *container<C>::subschildren(const lst & ls, const ls
 
 			// copy rest
 			while (cit != end) {
-				s->push_back(cit->subs(ls, lr, options));
+				s->push_back(cit->subs(m, options));
 				++cit;
 			}
 
diff --git a/ginac/ex.cpp b/ginac/ex.cpp
index 43b93766..e53cbd33 100644
--- a/ginac/ex.cpp
+++ b/ginac/ex.cpp
@@ -117,6 +117,65 @@ bool ex::find(const ex & pattern, lst & found) const
 	return any_found;
 }
 
+/** Substitute objects in an expression (syntactic substitution) and return
+ *  the result as a new expression. */
+ex ex::subs(const lst & ls, const lst & lr, unsigned options) const
+{
+	GINAC_ASSERT(ls.nops() == lr.nops());
+
+	// Convert the lists to a map
+	exmap m;
+	for (lst::const_iterator its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+		m[*its] = *itr;
+
+		// Search for products and powers in the expressions to be substituted
+		// (for an optimization in expairseq::subs())
+		if (is_exactly_a<mul>(*its) || is_exactly_a<power>(*its))
+			options |= subs_options::pattern_is_product;
+	}
+	if (!(options & subs_options::pattern_is_product))
+		options |= subs_options::pattern_is_not_product;
+
+	return bp->subs(m, options);
+}
+
+/** Substitute objects in an expression (syntactic substitution) and return
+ *  the result as a new expression.  There are two valid types of
+ *  replacement arguments: 1) a relational like object==ex and 2) a list of
+ *  relationals lst(object1==ex1,object2==ex2,...). */
+ex ex::subs(const ex & e, unsigned options) const
+{
+	if (e.info(info_flags::relation_equal)) {
+		exmap m;
+		const ex & s = e.lhs();
+		m[s] = e.rhs();
+		if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
+			options |= subs_options::pattern_is_product;
+		return bp->subs(m, options);
+	} else if (!e.info(info_flags::list))
+		throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
+
+	// Convert the list to a map
+	exmap m;
+	GINAC_ASSERT(is_a<lst>(e));
+	for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
+		ex r = *it;
+		if (!r.info(info_flags::relation_equal))
+			throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
+		const ex & s = r.lhs();
+		m[s] = r.rhs();
+
+		// Search for products and powers in the expressions to be substituted
+		// (for an optimization in expairseq::subs())
+		if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
+			options |= subs_options::pattern_is_product;
+	}
+	if (!(options & subs_options::pattern_is_product))
+		options |= subs_options::pattern_is_not_product;
+
+	return bp->subs(m, options);
+}
+
 /** Traverse expression tree with given visitor, preorder traversal. */
 void ex::traverse_preorder(visitor & v) const
 {
diff --git a/ginac/ex.h b/ginac/ex.h
index eaa6b808..b9360416 100644
--- a/ginac/ex.h
+++ b/ginac/ex.h
@@ -131,8 +131,9 @@ public:
 	bool match(const ex & pattern, lst & repl_lst) const { return bp->match(pattern, repl_lst); }
 
 	// substitutions
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return bp->subs(ls, lr, options); }
-	ex subs(const ex & e, unsigned options = 0) const { return bp->subs(e, options); }
+	ex subs(const exmap & m, unsigned options = 0) const;
+	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+	ex subs(const ex & e, unsigned options = 0) const;
 
 	// function mapping
 	ex map(map_function & f) const { return bp->map(f); }
@@ -429,12 +430,6 @@ inline ex series(const ex & thisex, const ex & r, int order, unsigned options =
 inline bool match(const ex & thisex, const ex & pattern, lst & repl_lst)
 { return thisex.match(pattern, repl_lst); }
 
-inline ex subs(const ex & thisex, const ex & e, unsigned options = 0)
-{ return thisex.subs(e, options); }
-
-inline ex subs(const ex & thisex, const lst & ls, const lst & lr, unsigned options = 0)
-{ return thisex.subs(ls, lr, options); }
-
 inline ex simplify_indexed(const ex & thisex)
 { return thisex.simplify_indexed(); }
 
@@ -497,13 +492,27 @@ struct ex_swap : public std::binary_function<ex, ex, void> {
 	void operator() (ex &lh, ex &rh) const { lh.swap(rh); }
 };
 
+inline ex ex::subs(const exmap & m, unsigned options) const
+{
+	return bp->subs(m, options);
+}
+
+inline ex subs(const ex & thisex, const exmap & m, unsigned options = 0)
+{ return thisex.subs(m, options); }
+
+inline ex subs(const ex & thisex, const lst & ls, const lst & lr, unsigned options = 0)
+{ return thisex.subs(ls, lr, options); }
+
+inline ex subs(const ex & thisex, const ex & e, unsigned options = 0)
+{ return thisex.subs(e, options); }
+
 
 /* Convert function pointer to function object suitable for map(). */
 class pointer_to_map_function : public map_function {
 protected:
 	ex (*ptr)(const ex &);
 public:
-	explicit pointer_to_map_function(ex (*x)(const ex &)) : ptr(x) {}
+	explicit pointer_to_map_function(ex x(const ex &)) : ptr(x) {}
 	ex operator()(const ex & e) { return ptr(e); }
 };
 
@@ -513,7 +522,7 @@ protected:
 	ex (*ptr)(const ex &, T1);
 	T1 arg1;
 public:
-	explicit pointer_to_map_function_1arg(ex (*x)(const ex &, T1), T1 a1) : ptr(x), arg1(a1) {}
+	explicit pointer_to_map_function_1arg(ex x(const ex &, T1), T1 a1) : ptr(x), arg1(a1) {}
 	ex operator()(const ex & e) { return ptr(e, arg1); }
 };
 
@@ -524,7 +533,7 @@ protected:
 	T1 arg1;
 	T2 arg2;
 public:
-	explicit pointer_to_map_function_2args(ex (*x)(const ex &, T1, T2), T1 a1, T2 a2) : ptr(x), arg1(a1), arg2(a2) {}
+	explicit pointer_to_map_function_2args(ex x(const ex &, T1, T2), T1 a1, T2 a2) : ptr(x), arg1(a1), arg2(a2) {}
 	ex operator()(const ex & e) { return ptr(e, arg1, arg2); }
 };
 
@@ -536,11 +545,11 @@ protected:
 	T2 arg2;
 	T3 arg3;
 public:
-	explicit pointer_to_map_function_3args(ex (*x)(const ex &, T1, T2, T3), T1 a1, T2 a2, T3 a3) : ptr(x), arg1(a1), arg2(a2), arg3(a3) {}
+	explicit pointer_to_map_function_3args(ex x(const ex &, T1, T2, T3), T1 a1, T2 a2, T3 a3) : ptr(x), arg1(a1), arg2(a2), arg3(a3) {}
 	ex operator()(const ex & e) { return ptr(e, arg1, arg2, arg3); }
 };
 
-inline ex ex::map(ex (*f)(const ex & e)) const
+inline ex ex::map(ex f(const ex &)) const
 {
 	pointer_to_map_function fcn(f);
 	return bp->map(fcn);
diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index f8675b76..5abdbc04 100644
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -386,15 +386,15 @@ found:		;
 	return inherited::match(pattern, repl_lst);
 }
 
-ex expairseq::subs(const lst &ls, const lst &lr, unsigned options) const
+ex expairseq::subs(const exmap & m, unsigned options) const
 {
-	epvector *vp = subschildren(ls, lr, options);
+	epvector *vp = subschildren(m, options);
 	if (vp)
 		return ex_to<basic>(thisexpairseq(vp, overall_coeff));
 	else if ((options & subs_options::subs_algebraic) && is_exactly_a<mul>(*this))
-		return static_cast<const mul *>(this)->algebraic_subs_mul(ls, lr, options);
+		return static_cast<const mul *>(this)->algebraic_subs_mul(m, options);
 	else
-		return subs_one_level(ls, lr, options);
+		return subs_one_level(m, options);
 }
 
 // protected
@@ -1546,29 +1546,32 @@ epvector * expairseq::evalchildren(int level) const
  *  @see expairseq::subs()
  *  @return pointer to epvector containing pairs after application of subs,
  *    or NULL pointer if no members were changed. */
-epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned options) const
+epvector * expairseq::subschildren(const exmap & m, unsigned options) const
 {
-	GINAC_ASSERT(ls.nops()==lr.nops());
-
-	// The substitution is "complex" when any of the objects to be substituted
-	// is a product or power. In this case we have to recombine the pairs
-	// because the numeric coefficients may be part of the search pattern.
-	bool complex_subs = false;
-	for (lst::const_iterator it = ls.begin(); it != ls.end(); ++it) {
-		if (is_exactly_a<mul>(*it) || is_exactly_a<power>(*it)) {
-			complex_subs = true;
-			break;
+	// When any of the objects to be substituted is a product or power
+	// we have to recombine the pairs because the numeric coefficients may
+	// be part of the search pattern.
+	if (!(options & (subs_options::pattern_is_product | subs_options::pattern_is_not_product))) {
+
+		// Search the list of substitutions and cache our findings
+		for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
+			if (is_exactly_a<mul>(it->first) || is_exactly_a<power>(it->first)) {
+				options |= subs_options::pattern_is_product;
+				break;
+			}
 		}
+		if (!(options & subs_options::pattern_is_product))
+			options |= subs_options::pattern_is_not_product;
 	}
 
-	if (complex_subs) {
+	if (options & subs_options::pattern_is_product) {
 
 		// Substitute in the recombined pairs
 		epvector::const_iterator cit = seq.begin(), last = seq.end();
 		while (cit != last) {
 
 			const ex &orig_ex = recombine_pair_to_ex(*cit);
-			const ex &subsed_ex = orig_ex.subs(ls, lr, options);
+			const ex &subsed_ex = orig_ex.subs(m, options);
 			if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
 
 				// Something changed, copy seq, subs and return it
@@ -1584,7 +1587,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned option
 
 				// Copy rest
 				while (cit != last) {
-					s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, options)));
+					s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options)));
 					++cit;
 				}
 				return s;
@@ -1599,7 +1602,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned option
 		epvector::const_iterator cit = seq.begin(), last = seq.end();
 		while (cit != last) {
 
-			const ex &subsed_ex = cit->rest.subs(ls, lr, options);
+			const ex &subsed_ex = cit->rest.subs(m, options);
 			if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
 			
 				// Something changed, copy seq, subs and return it
@@ -1615,7 +1618,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned option
 
 				// Copy rest
 				while (cit != last) {
-					s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, options),
+					s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options),
 					                                           cit->coeff));
 					++cit;
 				}
diff --git a/ginac/expairseq.h b/ginac/expairseq.h
index f1086353..b6f9bc7e 100644
--- a/ginac/expairseq.h
+++ b/ginac/expairseq.h
@@ -80,7 +80,7 @@ public:
 	ex to_rational(lst &repl_lst) const;
 	ex to_polynomial(lst &repl_lst) const;
 	bool match(const ex & pattern, lst & repl_lst) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const;
 protected:
 	bool is_equal_same_type(const basic & other) const;
 	unsigned return_type() const;
@@ -146,7 +146,7 @@ protected:
 	bool is_canonical() const;
 	epvector * expandchildren(unsigned options) const;
 	epvector * evalchildren(int level) const;
-	epvector * subschildren(const lst & ls, const lst & lr, unsigned options = 0) const;
+	epvector * subschildren(const exmap & m, unsigned options = 0) const;
 	
 // member variables
 	
diff --git a/ginac/flags.h b/ginac/flags.h
index 04c430c4..c14e5840 100644
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -29,9 +29,8 @@ namespace GiNaC {
 class expand_options {
 public:
 	enum {
-		expand_trigonometric = 0x0001,
-		expand_indexed = 0x0002,
-		expand_function_args = 0x0004
+		expand_indexed = 0x0001,      ///< expands (a+b).i to a.i+b.i
+		expand_function_args = 0x0002 ///< expands the arguments of functions
 	};
 };
 
@@ -39,8 +38,12 @@ public:
 class subs_options {
 public:
 	enum {
-		subs_no_pattern = 0x0001,
-		subs_algebraic = 0x0002
+		// this should be called "no_pattern"...
+		subs_no_pattern = 0x0001,        ///< disable pattern matching
+		// this should be called "algebraic"...
+		subs_algebraic = 0x0002,         ///< enable algebraic substitutions
+        pattern_is_product = 0x0004,     ///< used internally by expairseq::subschildren()
+        pattern_is_not_product = 0x0008  ///< used internally by expairseq::subschildren()
 	};
 };
 
@@ -163,10 +166,10 @@ public:
 class status_flags {
 public:
 	enum {
-		dynallocated    = 0x0001,       ///< Heap-allocated (i.e. created by new if we want to be clever and bypass the stack, @see ex::construct_from_basic() )
-		evaluated       = 0x0002,       ///< .eval() has already done its job
-		expanded        = 0x0004,       ///< .expand(0) has already done its job (other expand() options ignore this flag)
-		hash_calculated = 0x0008        ///< .calchash() has already done its job
+		dynallocated    = 0x0001, ///< heap-allocated (i.e. created by new if we want to be clever and bypass the stack, @see ex::construct_from_basic() )
+		evaluated       = 0x0002, ///< .eval() has already done its job
+		expanded        = 0x0004, ///< .expand(0) has already done its job (other expand() options ignore this flag)
+		hash_calculated = 0x0008  ///< .calchash() has already done its job
 	};
 };
 
diff --git a/ginac/idx.cpp b/ginac/idx.cpp
index 265c1075..cc5f52c3 100644
--- a/ginac/idx.cpp
+++ b/ginac/idx.cpp
@@ -348,29 +348,25 @@ ex idx::evalf(int level) const
 	return *this;
 }
 
-ex idx::subs(const lst & ls, const lst & lr, unsigned options) const
+ex idx::subs(const exmap & m, unsigned options) const
 {
-	GINAC_ASSERT(ls.nops() == lr.nops());
-
 	// First look for index substitutions
-	lst::const_iterator its, itr;
-	for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
-		if (is_equal(ex_to<basic>(*its))) {
-
-			// Substitution index->index
-			if (is_a<idx>(*itr))
-				return *itr;
-
-			// Otherwise substitute value
-			idx *i_copy = duplicate();
-			i_copy->value = *itr;
-			i_copy->clearflag(status_flags::hash_calculated);
-			return i_copy->setflag(status_flags::dynallocated);
-		}
+	exmap::const_iterator it = m.find(*this);
+	if (it != m.end()) {
+
+		// Substitution index->index
+		if (is_a<idx>(it->second))
+			return it->second;
+
+		// Otherwise substitute value
+		idx *i_copy = duplicate();
+		i_copy->value = it->second;
+		i_copy->clearflag(status_flags::hash_calculated);
+		return i_copy->setflag(status_flags::dynallocated);
 	}
 
 	// None, substitute objects in value (not in dimension)
-	const ex &subsed_value = value.subs(ls, lr, options);
+	const ex &subsed_value = value.subs(m, options);
 	if (are_ex_trivially_equal(value, subsed_value))
 		return *this;
 
diff --git a/ginac/idx.h b/ginac/idx.h
index bb74ddac..b6e63e07 100644
--- a/ginac/idx.h
+++ b/ginac/idx.h
@@ -54,7 +54,7 @@ public:
 	ex op(size_t i) const;
 	ex map(map_function & f) const;
 	ex evalf(int level = 0) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const;
 
 protected:
 	ex derivative(const symbol & s) const;
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 990da839..cdb875d3 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -229,14 +229,14 @@ ex matrix::eval(int level) const
 											   status_flags::evaluated);
 }
 
-ex matrix::subs(const lst & ls, const lst & lr, unsigned options) const
-{	
+ex matrix::subs(const exmap & mp, unsigned options) const
+{
 	exvector m2(row * col);
 	for (unsigned r=0; r<row; ++r)
 		for (unsigned c=0; c<col; ++c)
-			m2[r*col+c] = m[r*col+c].subs(ls, lr, options);
+			m2[r*col+c] = m[r*col+c].subs(mp, options);
 
-	return matrix(row, col, m2).subs_one_level(ls, lr, options);
+	return matrix(row, col, m2).subs_one_level(mp, options);
 }
 
 // protected
diff --git a/ginac/matrix.h b/ginac/matrix.h
index d0d9d7d1..f951850c 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -49,7 +49,7 @@ public:
 	ex & let_op(size_t i);
 	ex eval(int level=0) const;
 	ex evalm() const {return *this;}
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const;
 	ex eval_indexed(const basic & i) const;
 	ex add_indexed(const ex & self, const ex & other) const;
 	ex scalar_mul_indexed(const ex & self, const numeric & other) const;
diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index 8fe73998..2c845651 100644
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -558,27 +558,26 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
 	return true;
 }
 
-ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const
+ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
 {	
 	std::vector<bool> subsed(seq.size(), false);
 	exvector subsresult(seq.size());
 
-	lst::const_iterator its, itr;
-	for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+	for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
 
-		if (is_exactly_a<mul>(*its)) {
+		if (is_exactly_a<mul>(it->first)) {
 
 			int nummatches = std::numeric_limits<int>::max();
 			std::vector<bool> currsubsed(seq.size(), false);
 			bool succeed = true;
 			lst repls;
 
-			for (size_t j=0; j<its->nops(); j++) {
+			for (size_t j=0; j<it->first.nops(); j++) {
 				bool found=false;
 				for (size_t k=0; k<nops(); k++) {
 					if (currsubsed[k] || subsed[k])
 						continue;
-					if (tryfactsubs(op(k), its->op(j), nummatches, repls)) {
+					if (tryfactsubs(op(k), it->first.op(j), nummatches, repls)) {
 						currsubsed[k] = true;
 						found = true;
 						break;
@@ -599,7 +598,7 @@ ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) con
 						subsresult[j] = op(j);
 					else {
 						foundfirstsubsedfactor = true;
-						subsresult[j] = op(j) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches);
+						subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches);
 					}
 					subsed[j] = true;
 				}
@@ -611,9 +610,9 @@ ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) con
 			lst repls;
 
 			for (size_t j=0; j<this->nops(); j++) {
-				if (!subsed[j] && tryfactsubs(op(j), *its, nummatches, repls)) {
+				if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)) {
 					subsed[j] = true;
-					subsresult[j] = op(j) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches);
+					subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches);
 				}
 			}
 		}
@@ -627,7 +626,7 @@ ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) con
 		}
 	}
 	if (!subsfound)
-		return subs_one_level(ls, lr, options | subs_options::subs_algebraic);
+		return subs_one_level(m, options | subs_options::subs_algebraic);
 
 	exvector ev; ev.reserve(nops());
 	for (size_t i=0; i<nops(); i++) {
diff --git a/ginac/mul.h b/ginac/mul.h
index 01fb2056..0ec57202 100644
--- a/ginac/mul.h
+++ b/ginac/mul.h
@@ -57,7 +57,7 @@ public:
 	ex evalf(int level=0) const;
 	ex evalm() const;
 	ex series(const relational & s, int order, unsigned options = 0) const;
-	ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+	ex normal(exmap & repl, int level = 0) const;
 	numeric integer_content() const;
 	ex smod(const numeric &xi) const;
 	numeric max_coefficient() const;
@@ -85,7 +85,7 @@ protected:
 	
 	// non-virtual functions in this class
 public:
-	ex algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const;
+	ex algebraic_subs_mul(const exmap & m, unsigned options) const;
 protected:
 	epvector * expandchildren(unsigned options) const;
 };
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index cd2ad8ec..3a8a82fc 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1681,25 +1681,22 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 
 
 /** Create a symbol for replacing the expression "e" (or return a previously
- *  assigned symbol). The symbol is appended to sym_lst and returned, the
- *  expression is appended to repl_lst.
+ *  assigned symbol). The symbol and expression are appended to repl, for
+ *  a later application of subs().
  *  @see ex::normal */
-static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst)
+static ex replace_with_symbol(const ex & e, exmap & repl)
 {
-	// Expression already in repl_lst? Then return the assigned symbol
-	lst::const_iterator its, itr;
-	for (its = sym_lst.begin(), itr = repl_lst.begin(); itr != repl_lst.end(); ++its, ++itr)
-		if (itr->is_equal(e))
-			return *its;
+	// Expression already in repl? Then return the assigned symbol
+	for (exmap::const_iterator it = repl.begin(); it != repl.end(); ++it)
+		if (it->second.is_equal(e))
+			return it->first;
 	
 	// Otherwise create new symbol and add to list, taking care that the
-	// replacement expression doesn't contain symbols from the sym_lst
+	// replacement expression doesn't itself contain symbols from repl,
 	// because subs() is not recursive
-	symbol s;
-	ex es(s);
-	ex e_replaced = e.subs(sym_lst, repl_lst);
-	sym_lst.append(es);
-	repl_lst.append(e_replaced);
+	ex es = (new symbol)->setflag(status_flags::dynallocated);
+	ex e_replaced = e.subs(repl);
+	repl[es] = e_replaced;
 	return es;
 }
 
@@ -1708,7 +1705,7 @@ static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst)
  *  to repl_lst and the symbol is returned.
  *  @see basic::to_rational
  *  @see basic::to_polynomial */
-static ex replace_with_symbol(const ex &e, lst &repl_lst)
+static ex replace_with_symbol(const ex & e, lst & repl_lst)
 {
 	// Expression already in repl_lst? Then return the assigned symbol
 	for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
@@ -1716,10 +1713,9 @@ static ex replace_with_symbol(const ex &e, lst &repl_lst)
 			return it->op(0);
 	
 	// Otherwise create new symbol and add to list, taking care that the
-	// replacement expression doesn't contain symbols from the sym_lst
+	// replacement expression doesn't itself contain symbols from the repl_lst,
 	// because subs() is not recursive
-	symbol s;
-	ex es(s);
+	ex es = (new symbol)->setflag(status_flags::dynallocated);
 	ex e_replaced = e.subs(repl_lst);
 	repl_lst.append(es == e_replaced);
 	return es;
@@ -1736,18 +1732,18 @@ struct normal_map_function : public map_function {
 /** Default implementation of ex::normal(). It normalizes the children and
  *  replaces the object with a temporary symbol.
  *  @see ex::normal */
-ex basic::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex basic::normal(exmap & repl, int level) const
 {
 	if (nops() == 0)
-		return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+		return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
 	else {
 		if (level == 1)
-			return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+			return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
 		else if (level == -max_recursion_level)
 			throw(std::runtime_error("max recursion level reached"));
 		else {
 			normal_map_function map_normal(level - 1);
-			return (new lst(replace_with_symbol(map(map_normal), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+			return (new lst(replace_with_symbol(map(map_normal), repl), _ex1))->setflag(status_flags::dynallocated);
 		}
 	}
 }
@@ -1755,7 +1751,7 @@ ex basic::normal(lst &sym_lst, lst &repl_lst, int level) const
 
 /** Implementation of ex::normal() for symbols. This returns the unmodified symbol.
  *  @see ex::normal */
-ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex symbol::normal(exmap & repl, int level) const
 {
 	return (new lst(*this, _ex1))->setflag(status_flags::dynallocated);
 }
@@ -1765,19 +1761,19 @@ ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  into re+I*im and replaces I and non-rational real numbers with a temporary
  *  symbol.
  *  @see ex::normal */
-ex numeric::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex numeric::normal(exmap & repl, int level) const
 {
 	numeric num = numer();
 	ex numex = num;
 
 	if (num.is_real()) {
 		if (!num.is_integer())
-			numex = replace_with_symbol(numex, sym_lst, repl_lst);
+			numex = replace_with_symbol(numex, repl);
 	} else { // complex
 		numeric re = num.real(), im = num.imag();
-		ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst);
-		ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst);
-		numex = re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst);
+		ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
+		ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
+		numex = re_ex + im_ex * replace_with_symbol(I, repl);
 	}
 
 	// Denominator is always a real integer (see numeric::denom())
@@ -1849,10 +1845,10 @@ static ex frac_cancel(const ex &n, const ex &d)
 /** Implementation of ex::normal() for a sum. It expands terms and performs
  *  fractional addition.
  *  @see ex::normal */
-ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex add::normal(exmap & repl, int level) const
 {
 	if (level == 1)
-		return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+		return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
 	else if (level == -max_recursion_level)
 		throw(std::runtime_error("max recursion level reached"));
 
@@ -1862,12 +1858,12 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
 	dens.reserve(seq.size()+1);
 	epvector::const_iterator it = seq.begin(), itend = seq.end();
 	while (it != itend) {
-		ex n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(sym_lst, repl_lst, level-1);
+		ex n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, level-1);
 		nums.push_back(n.op(0));
 		dens.push_back(n.op(1));
 		it++;
 	}
-	ex n = ex_to<numeric>(overall_coeff).normal(sym_lst, repl_lst, level-1);
+	ex n = ex_to<numeric>(overall_coeff).normal(repl, level-1);
 	nums.push_back(n.op(0));
 	dens.push_back(n.op(1));
 	GINAC_ASSERT(nums.size() == dens.size());
@@ -1908,10 +1904,10 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
 /** Implementation of ex::normal() for a product. It cancels common factors
  *  from fractions.
  *  @see ex::normal() */
-ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex mul::normal(exmap & repl, int level) const
 {
 	if (level == 1)
-		return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+		return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
 	else if (level == -max_recursion_level)
 		throw(std::runtime_error("max recursion level reached"));
 
@@ -1921,12 +1917,12 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
 	ex n;
 	epvector::const_iterator it = seq.begin(), itend = seq.end();
 	while (it != itend) {
-		n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(sym_lst, repl_lst, level-1);
+		n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, level-1);
 		num.push_back(n.op(0));
 		den.push_back(n.op(1));
 		it++;
 	}
-	n = ex_to<numeric>(overall_coeff).normal(sym_lst, repl_lst, level-1);
+	n = ex_to<numeric>(overall_coeff).normal(repl, level-1);
 	num.push_back(n.op(0));
 	den.push_back(n.op(1));
 
@@ -1940,16 +1936,16 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  distributes integer exponents to numerator and denominator, and replaces
  *  non-integer powers by temporary symbols.
  *  @see ex::normal */
-ex power::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex power::normal(exmap & repl, int level) const
 {
 	if (level == 1)
-		return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+		return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
 	else if (level == -max_recursion_level)
 		throw(std::runtime_error("max recursion level reached"));
 
 	// Normalize basis and exponent (exponent gets reassembled)
-	ex n_basis = ex_to<basic>(basis).normal(sym_lst, repl_lst, level-1);
-	ex n_exponent = ex_to<basic>(exponent).normal(sym_lst, repl_lst, level-1);
+	ex n_basis = ex_to<basic>(basis).normal(repl, level-1);
+	ex n_exponent = ex_to<basic>(exponent).normal(repl, level-1);
 	n_exponent = n_exponent.op(0) / n_exponent.op(1);
 
 	if (n_exponent.info(info_flags::integer)) {
@@ -1970,32 +1966,32 @@ ex power::normal(lst &sym_lst, lst &repl_lst, int level) const
 		if (n_exponent.info(info_flags::positive)) {
 
 			// (a/b)^x -> {sym((a/b)^x), 1}
-			return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+			return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
 
 		} else if (n_exponent.info(info_flags::negative)) {
 
 			if (n_basis.op(1).is_equal(_ex1)) {
 
 				// a^-x -> {1, sym(a^x)}
-				return (new lst(_ex1, replace_with_symbol(power(n_basis.op(0), -n_exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated);
+				return (new lst(_ex1, replace_with_symbol(power(n_basis.op(0), -n_exponent), repl)))->setflag(status_flags::dynallocated);
 
 			} else {
 
 				// (a/b)^-x -> {sym((b/a)^x), 1}
-				return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+				return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
 			}
 		}
 	}
 
 	// (a/b)^x -> {sym((a/b)^x, 1}
-	return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+	return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
 }
 
 
 /** Implementation of ex::normal() for pseries. It normalizes each coefficient
  *  and replaces the series by a temporary symbol.
  *  @see ex::normal */
-ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex pseries::normal(exmap & repl, int level) const
 {
 	epvector newseq;
 	epvector::const_iterator i = seq.begin(), end = seq.end();
@@ -2006,7 +2002,7 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
 		++i;
 	}
 	ex n = pseries(relational(var,point), newseq);
-	return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+	return (new lst(replace_with_symbol(n, repl), _ex1))->setflag(status_flags::dynallocated);
 }
 
 
@@ -2024,14 +2020,14 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  @return normalized expression */
 ex ex::normal(int level) const
 {
-	lst sym_lst, repl_lst;
+	exmap repl;
 
-	ex e = bp->normal(sym_lst, repl_lst, level);
+	ex e = bp->normal(repl, level);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
-	if (sym_lst.nops() > 0)
-		e = e.subs(sym_lst, repl_lst);
+	if (!repl.empty())
+		e = e.subs(repl);
 
 	// Convert {numerator, denominator} form back to fraction
 	return e.op(0) / e.op(1);
@@ -2045,16 +2041,16 @@ ex ex::normal(int level) const
  *  @return numerator */
 ex ex::numer() const
 {
-	lst sym_lst, repl_lst;
+	exmap repl;
 
-	ex e = bp->normal(sym_lst, repl_lst, 0);
+	ex e = bp->normal(repl, 0);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
-	if (sym_lst.nops() > 0)
-		return e.op(0).subs(sym_lst, repl_lst);
-	else
+	if (repl.empty())
 		return e.op(0);
+	else
+		return e.op(0).subs(repl);
 }
 
 /** Get denominator of an expression. If the expression is not of the normal
@@ -2065,16 +2061,16 @@ ex ex::numer() const
  *  @return denominator */
 ex ex::denom() const
 {
-	lst sym_lst, repl_lst;
+	exmap repl;
 
-	ex e = bp->normal(sym_lst, repl_lst, 0);
+	ex e = bp->normal(repl, 0);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
-	if (sym_lst.nops() > 0)
-		return e.op(1).subs(sym_lst, repl_lst);
-	else
+	if (repl.empty())
 		return e.op(1);
+	else
+		return e.op(1).subs(repl);
 }
 
 /** Get numerator and denominator of an expression. If the expresison is not
@@ -2085,16 +2081,16 @@ ex ex::denom() const
  *  @return a list [numerator, denominator] */
 ex ex::numer_denom() const
 {
-	lst sym_lst, repl_lst;
+	exmap repl;
 
-	ex e = bp->normal(sym_lst, repl_lst, 0);
+	ex e = bp->normal(repl, 0);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
-	if (sym_lst.nops() > 0)
-		return e.subs(sym_lst, repl_lst);
-	else
+	if (repl.empty())
 		return e;
+	else
+		return e.subs(repl);
 }
 
 
diff --git a/ginac/numeric.h b/ginac/numeric.h
index 33498f3a..a4f4447b 100644
--- a/ginac/numeric.h
+++ b/ginac/numeric.h
@@ -104,8 +104,8 @@ public:
 	bool has(const ex &other) const;
 	ex eval(int level = 0) const;
 	ex evalf(int level = 0) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return subs_one_level(ls, lr, options); } // overwrites basic::subs() for performance reasons
-	ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const { return subs_one_level(m, options); } // overwrites basic::subs() for performance reasons
+	ex normal(exmap & repl, int level = 0) const;
 	ex to_rational(lst &repl_lst) const;
 	ex to_polynomial(lst &repl_lst) const;
 	numeric integer_content() const;
diff --git a/ginac/power.cpp b/ginac/power.cpp
index 0d24a04f..a30182ef 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -518,27 +518,26 @@ ex power::evalm() const
 // from mul.cpp
 extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
 
-ex power::subs(const lst & ls, const lst & lr, unsigned options) const
+ex power::subs(const exmap & m, unsigned options) const
 {	
-	const ex &subsed_basis = basis.subs(ls, lr, options);
-	const ex &subsed_exponent = exponent.subs(ls, lr, options);
+	const ex &subsed_basis = basis.subs(m, options);
+	const ex &subsed_exponent = exponent.subs(m, options);
 
 	if (!are_ex_trivially_equal(basis, subsed_basis)
 	 || !are_ex_trivially_equal(exponent, subsed_exponent)) 
-		return power(subsed_basis, subsed_exponent).subs_one_level(ls, lr, options);
+		return power(subsed_basis, subsed_exponent).subs_one_level(m, options);
 
 	if (!(options & subs_options::subs_algebraic))
-		return subs_one_level(ls, lr, options);
+		return subs_one_level(m, options);
 
-	lst::const_iterator its, itr;
-	for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+	for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
 		int nummatches = std::numeric_limits<int>::max();
 		lst repls;
-		if (tryfactsubs(*this, *its, nummatches, repls))
-			return (ex_to<basic>((*this) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(ls, lr, options);
+		if (tryfactsubs(*this, it->first, nummatches, repls))
+			return (ex_to<basic>((*this) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(m, options);
 	}
 
-	return subs_one_level(ls, lr, options);
+	return subs_one_level(m, options);
 }
 
 ex power::eval_ncmul(const exvector & v) const
diff --git a/ginac/power.h b/ginac/power.h
index 4e1df9fb..a1f5a23a 100644
--- a/ginac/power.h
+++ b/ginac/power.h
@@ -61,8 +61,8 @@ public:
 	ex evalf(int level=0) const;
 	ex evalm() const;
 	ex series(const relational & s, int order, unsigned options = 0) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
-	ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const;
+	ex normal(exmap & repl, int level = 0) const;
 	ex to_rational(lst &repl_lst) const;
 	ex to_polynomial(lst &repl_lst) const;
 	exvector get_free_indices() const;
diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index e9887055..22d7b777 100644
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -396,13 +396,13 @@ ex pseries::evalf(int level) const
 	return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
 }
 
-ex pseries::subs(const lst & ls, const lst & lr, unsigned options) const
+ex pseries::subs(const exmap & m, unsigned options) const
 {
 	// If expansion variable is being substituted, convert the series to a
 	// polynomial and do the substitution there because the result might
 	// no longer be a power series
-	if (ls.has(var))
-		return convert_to_poly(true).subs(ls, lr, options);
+	if (m.find(var) != m.end())
+		return convert_to_poly(true).subs(m, options);
 	
 	// Otherwise construct a new series with substituted coefficients and
 	// expansion point
@@ -410,10 +410,10 @@ ex pseries::subs(const lst & ls, const lst & lr, unsigned options) const
 	newseq.reserve(seq.size());
 	epvector::const_iterator it = seq.begin(), itend = seq.end();
 	while (it != itend) {
-		newseq.push_back(expair(it->rest.subs(ls, lr, options), it->coeff));
+		newseq.push_back(expair(it->rest.subs(m, options), it->coeff));
 		++it;
 	}
-	return (new pseries(relational(var,point.subs(ls, lr, options)), newseq))->setflag(status_flags::dynallocated);
+	return (new pseries(relational(var,point.subs(m, options)), newseq))->setflag(status_flags::dynallocated);
 }
 
 /** Implementation of ex::expand() for a power series.  It expands all the
diff --git a/ginac/pseries.h b/ginac/pseries.h
index 17ed024d..5dbdc66e 100644
--- a/ginac/pseries.h
+++ b/ginac/pseries.h
@@ -53,8 +53,8 @@ public:
 	ex eval(int level=0) const;
 	ex evalf(int level=0) const;
 	ex series(const relational & r, int order, unsigned options = 0) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
-	ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const;
+	ex normal(exmap & repl, int level = 0) const;
 	ex expand(unsigned options = 0) const;
 protected:
 	ex derivative(const symbol & s) const;
diff --git a/ginac/relational.cpp b/ginac/relational.cpp
index c45dd1ef..453b4c00 100644
--- a/ginac/relational.cpp
+++ b/ginac/relational.cpp
@@ -178,15 +178,15 @@ ex relational::eval(int level) const
 	return (new relational(lh.eval(level-1),rh.eval(level-1),o))->setflag(status_flags::dynallocated | status_flags::evaluated);
 }
 
-ex relational::subs(const lst & ls, const lst & lr, unsigned options) const
+ex relational::subs(const exmap & m, unsigned options) const
 {
-	const ex & subsed_lh = lh.subs(ls, lr, options);
-	const ex & subsed_rh = rh.subs(ls, lr, options);
+	const ex & subsed_lh = lh.subs(m, options);
+	const ex & subsed_rh = rh.subs(m, options);
 
 	if (!are_ex_trivially_equal(lh, subsed_lh) || !are_ex_trivially_equal(rh, subsed_rh))
-		return relational(subsed_lh, subsed_rh, o).subs_one_level(ls, lr, options);
+		return relational(subsed_lh, subsed_rh, o).subs_one_level(m, options);
 	else
-		return subs_one_level(ls, lr, options);
+		return subs_one_level(m, options);
 }
 
 ex relational::eval_ncmul(const exvector & v) const
diff --git a/ginac/relational.h b/ginac/relational.h
index 591cc3b7..cbe51468 100644
--- a/ginac/relational.h
+++ b/ginac/relational.h
@@ -57,7 +57,7 @@ public:
 	size_t nops() const;
 	ex op(size_t i) const;
 	ex map(map_function & f) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const;
 	ex eval(int level=0) const;
 
 protected:
diff --git a/ginac/structure.h b/ginac/structure.h
index 23ae93c0..41b243ac 100644
--- a/ginac/structure.h
+++ b/ginac/structure.h
@@ -160,7 +160,7 @@ protected:
 public:
 
 	// substitutions
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return inherited::subs(ls, lr, options); }
+	ex subs(const exmap & m, unsigned options = 0) const { return inherited::subs(m, options); }
 
 	// function mapping
 	ex map(map_function & f) const { return inherited::map(f); }
@@ -181,7 +181,7 @@ public:
 	ex series(const relational & r, int order, unsigned options = 0) const { return inherited::series(r, order, options); }
 
 	// rational functions
-	ex normal(lst & sym_lst, lst & repl_lst, int level = 0) const { return inherited::normal(sym_lst, repl_lst, level); }
+	ex normal(exmap & repl, int level = 0) const { return inherited::normal(repl, level); }
 	ex to_rational(lst & repl_lst) const { return inherited::to_rational(repl_lst); }
 	ex to_polynomial(lst & repl_lst) const { return inherited::to_polynomial(repl_lst); }
 
diff --git a/ginac/symbol.h b/ginac/symbol.h
index 8ff63130..8cd738b5 100644
--- a/ginac/symbol.h
+++ b/ginac/symbol.h
@@ -67,8 +67,8 @@ public:
 	ex eval(int level = 0) const;
 	ex evalf(int level = 0) const { return *this; } // overwrites basic::evalf() for performance reasons
 	ex series(const relational & s, int order, unsigned options = 0) const;
-	ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return subs_one_level(ls, lr, options); } // overwrites basic::subs() for performance reasons
-	ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+	ex subs(const exmap & m, unsigned options = 0) const { return subs_one_level(m, options); } // overwrites basic::subs() for performance reasons
+	ex normal(exmap & repl, int level = 0) const;
 	ex to_rational(lst &repl_lst) const;
 	ex to_polynomial(lst &repl_lst) const;
 protected:
diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
index d78eff66..f8e5635a 100644
--- a/ginac/tensor.cpp
+++ b/ginac/tensor.cpp
@@ -171,7 +171,10 @@ ex tensdelta::eval_indexed(const basic & i) const
 	// dimension
 	if (!i1.get_dim().is_equal(i2.get_dim())) {
 		ex min_dim = i1.minimal_dim(i2);
-		return i.subs(lst(i1 == i1.replace_dim(min_dim), i2 == i2.replace_dim(min_dim)));
+		exmap m;
+		m[i1] = i1.replace_dim(min_dim);
+		m[i2] = i2.replace_dim(min_dim);
+		return i.subs(m);
 	}
 
 	// Trace of delta tensor is the (effective) dimension of the space
@@ -212,7 +215,10 @@ ex tensmetric::eval_indexed(const basic & i) const
 	// dimension
 	if (!i1.get_dim().is_equal(i2.get_dim())) {
 		ex min_dim = i1.minimal_dim(i2);
-		return i.subs(lst(i1 == i1.replace_dim(min_dim), i2 == i2.replace_dim(min_dim)));
+		exmap m;
+		m[i1] = i1.replace_dim(min_dim);
+		m[i2] = i2.replace_dim(min_dim);
+		return i.subs(m);
 	}
 
 	// A metric tensor with one covariant and one contravariant index gets