X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=67fbf93612beaa629ed2df6a987152e56a546374;hp=610c48d6a2653bd9945ce35862b5e59e7f2f66da;hb=73081e1cccfec32a9e9dab6a5e4b063b3fa0e7fb;hpb=e7cc6a764ff67b5885d6633385fac23ccc1dc9a7

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 610c48d6..67fbf936 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -6,7 +6,7 @@
  *  computation, square-free factorization and rational function normalization. */
 
 /*
- *  GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2004 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
@@ -39,6 +39,7 @@
 #include "numeric.h"
 #include "power.h"
 #include "relational.h"
+#include "operators.h"
 #include "matrix.h"
 #include "pseries.h"
 #include "symbol.h"
@@ -88,18 +89,18 @@ static struct _stat_print {
  *  function returns for a given expression.
  *
  *  @param e  expression to search
- *  @param x  pointer to first symbol found (returned)
+ *  @param x  first symbol found (returned)
  *  @return "false" if no symbol was found, "true" otherwise */
-static bool get_first_symbol(const ex &e, const symbol *&x)
+static bool get_first_symbol(const ex &e, ex &x)
 {
-	if (is_ex_exactly_of_type(e, symbol)) {
-		x = &ex_to<symbol>(e);
+	if (is_a<symbol>(e)) {
+		x = e;
 		return true;
-	} else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) {
-		for (unsigned i=0; i<e.nops(); i++)
+	} else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
+		for (size_t i=0; i<e.nops(); i++)
 			if (get_first_symbol(e.op(i), x))
 				return true;
-	} else if (is_ex_exactly_of_type(e, power)) {
+	} else if (is_exactly_a<power>(e)) {
 		if (get_first_symbol(e.op(0), x))
 			return true;
 	}
@@ -118,8 +119,8 @@ static bool get_first_symbol(const ex &e, const symbol *&x)
  *
  *  @see get_symbol_stats */
 struct sym_desc {
-	/** Pointer to symbol */
-	const symbol *sym;
+	/** Reference to symbol */
+	ex sym;
 
 	/** Highest degree of symbol in polynomial "a" */
 	int deg_a;
@@ -137,7 +138,7 @@ struct sym_desc {
 	int max_deg;
 
 	/** Maximum number of terms of leading coefficient of symbol in both polynomials */
-	int max_lcnops;
+	size_t max_lcnops;
 
 	/** Commparison operator for sorting */
 	bool operator<(const sym_desc &x) const
@@ -153,11 +154,11 @@ struct sym_desc {
 typedef std::vector<sym_desc> sym_desc_vec;
 
 // Add symbol the sym_desc_vec (used internally by get_symbol_stats())
-static void add_symbol(const symbol *s, sym_desc_vec &v)
+static void add_symbol(const ex &s, sym_desc_vec &v)
 {
 	sym_desc_vec::const_iterator it = v.begin(), itend = v.end();
 	while (it != itend) {
-		if (it->sym->compare(*s) == 0)  // If it's already in there, don't add it a second time
+		if (it->sym.is_equal(s))  // If it's already in there, don't add it a second time
 			return;
 		++it;
 	}
@@ -169,12 +170,12 @@ static void add_symbol(const symbol *s, sym_desc_vec &v)
 // Collect all symbols of an expression (used internally by get_symbol_stats())
 static void collect_symbols(const ex &e, sym_desc_vec &v)
 {
-	if (is_ex_exactly_of_type(e, symbol)) {
-		add_symbol(&ex_to<symbol>(e), v);
-	} else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) {
-		for (unsigned i=0; i<e.nops(); i++)
+	if (is_a<symbol>(e)) {
+		add_symbol(e, v);
+	} else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
+		for (size_t i=0; i<e.nops(); i++)
 			collect_symbols(e.op(i), v);
-	} else if (is_ex_exactly_of_type(e, power)) {
+	} else if (is_exactly_a<power>(e)) {
 		collect_symbols(e.op(0), v);
 	}
 }
@@ -197,23 +198,24 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 	collect_symbols(b.eval(), v);
 	sym_desc_vec::iterator it = v.begin(), itend = v.end();
 	while (it != itend) {
-		int deg_a = a.degree(*(it->sym));
-		int deg_b = b.degree(*(it->sym));
+		int deg_a = a.degree(it->sym);
+		int deg_b = b.degree(it->sym);
 		it->deg_a = deg_a;
 		it->deg_b = deg_b;
 		it->max_deg = std::max(deg_a, deg_b);
-		it->max_lcnops = std::max(a.lcoeff(*(it->sym)).nops(), b.lcoeff(*(it->sym)).nops());
-		it->ldeg_a = a.ldegree(*(it->sym));
-		it->ldeg_b = b.ldegree(*(it->sym));
+		it->max_lcnops = std::max(a.lcoeff(it->sym).nops(), b.lcoeff(it->sym).nops());
+		it->ldeg_a = a.ldegree(it->sym);
+		it->ldeg_b = b.ldegree(it->sym);
 		++it;
 	}
 	std::sort(v.begin(), v.end());
+
 #if 0
 	std::clog << "Symbols:\n";
 	it = v.begin(); itend = v.end();
 	while (it != itend) {
-		std::clog << " " << *it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
-		std::clog << "  lcoeff_a=" << a.lcoeff(*(it->sym)) << ", lcoeff_b=" << b.lcoeff(*(it->sym)) << endl;
+		std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
+		std::clog << "  lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << endl;
 		++it;
 	}
 #endif
@@ -230,18 +232,18 @@ static numeric lcmcoeff(const ex &e, const numeric &l)
 {
 	if (e.info(info_flags::rational))
 		return lcm(ex_to<numeric>(e).denom(), l);
-	else if (is_ex_exactly_of_type(e, add)) {
+	else if (is_exactly_a<add>(e)) {
 		numeric c = _num1;
-		for (unsigned i=0; i<e.nops(); i++)
+		for (size_t i=0; i<e.nops(); i++)
 			c = lcmcoeff(e.op(i), c);
 		return lcm(c, l);
-	} else if (is_ex_exactly_of_type(e, mul)) {
+	} else if (is_exactly_a<mul>(e)) {
 		numeric c = _num1;
-		for (unsigned i=0; i<e.nops(); i++)
+		for (size_t i=0; i<e.nops(); i++)
 			c *= lcmcoeff(e.op(i), _num1);
 		return lcm(c, l);
-	} else if (is_ex_exactly_of_type(e, power)) {
-		if (is_ex_exactly_of_type(e.op(0), symbol))
+	} else if (is_exactly_a<power>(e)) {
+		if (is_a<symbol>(e.op(0)))
 			return l;
 		else
 			return pow(lcmcoeff(e.op(0), l), ex_to<numeric>(e.op(1)));
@@ -268,25 +270,25 @@ static numeric lcm_of_coefficients_denominators(const ex &e)
  *  @param lcm  LCM to multiply in */
 static ex multiply_lcm(const ex &e, const numeric &lcm)
 {
-	if (is_ex_exactly_of_type(e, mul)) {
-		unsigned num = e.nops();
+	if (is_exactly_a<mul>(e)) {
+		size_t num = e.nops();
 		exvector v; v.reserve(num + 1);
 		numeric lcm_accum = _num1;
-		for (unsigned i=0; i<e.nops(); i++) {
+		for (size_t i=0; i<num; i++) {
 			numeric op_lcm = lcmcoeff(e.op(i), _num1);
 			v.push_back(multiply_lcm(e.op(i), op_lcm));
 			lcm_accum *= op_lcm;
 		}
 		v.push_back(lcm / lcm_accum);
 		return (new mul(v))->setflag(status_flags::dynallocated);
-	} else if (is_ex_exactly_of_type(e, add)) {
-		unsigned num = e.nops();
+	} else if (is_exactly_a<add>(e)) {
+		size_t num = e.nops();
 		exvector v; v.reserve(num);
-		for (unsigned i=0; i<num; i++)
+		for (size_t i=0; i<num; i++)
 			v.push_back(multiply_lcm(e.op(i), lcm));
 		return (new add(v))->setflag(status_flags::dynallocated);
-	} else if (is_ex_exactly_of_type(e, power)) {
-		if (is_ex_exactly_of_type(e.op(0), symbol))
+	} else if (is_exactly_a<power>(e)) {
+		if (is_a<symbol>(e.op(0)))
 			return e * lcm;
 		else
 			return pow(multiply_lcm(e.op(0), lcm.power(ex_to<numeric>(e.op(1)).inverse())), e.op(1));
@@ -296,43 +298,45 @@ static ex multiply_lcm(const ex &e, const numeric &lcm)
 
 
 /** Compute the integer content (= GCD of all numeric coefficients) of an
- *  expanded polynomial.
+ *  expanded polynomial. For a polynomial with rational coefficients, this
+ *  returns g/l where g is the GCD of the coefficients' numerators and l
+ *  is the LCM of the coefficients' denominators.
  *
- *  @param e  expanded polynomial
  *  @return integer content */
-numeric ex::integer_content(void) const
+numeric ex::integer_content() const
 {
-	GINAC_ASSERT(bp!=0);
 	return bp->integer_content();
 }
 
-numeric basic::integer_content(void) const
+numeric basic::integer_content() const
 {
 	return _num1;
 }
 
-numeric numeric::integer_content(void) const
+numeric numeric::integer_content() const
 {
 	return abs(*this);
 }
 
-numeric add::integer_content(void) const
+numeric add::integer_content() const
 {
 	epvector::const_iterator it = seq.begin();
 	epvector::const_iterator itend = seq.end();
-	numeric c = _num0;
+	numeric c = _num0, l = _num1;
 	while (it != itend) {
 		GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
 		GINAC_ASSERT(is_exactly_a<numeric>(it->coeff));
-		c = gcd(ex_to<numeric>(it->coeff), c);
+		c = gcd(ex_to<numeric>(it->coeff).numer(), c);
+		l = lcm(ex_to<numeric>(it->coeff).denom(), l);
 		it++;
 	}
 	GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
-	c = gcd(ex_to<numeric>(overall_coeff),c);
-	return c;
+	c = gcd(ex_to<numeric>(overall_coeff).numer(), c);
+	l = lcm(ex_to<numeric>(overall_coeff).denom(), l);
+	return c/l;
 }
 
-numeric mul::integer_content(void) const
+numeric mul::integer_content() const
 {
 #ifdef DO_GINAC_ASSERT
 	epvector::const_iterator it = seq.begin();
@@ -360,11 +364,11 @@ numeric mul::integer_content(void) const
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return quotient of a and b in Q[x] */
-ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
+ex quo(const ex &a, const ex &b, const ex &x, bool check_args)
 {
 	if (b.is_zero())
 		throw(std::overflow_error("quo: division by zero"));
-	if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric))
+	if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b))
 		return a / b;
 #if FAST_COMPARE
 	if (a.is_equal(b))
@@ -380,7 +384,7 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
 	int bdeg = b.degree(x);
 	int rdeg = r.degree(x);
 	ex blcoeff = b.expand().coeff(x, bdeg);
-	bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
+	bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
 	exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
 	while (rdeg >= bdeg) {
 		ex term, rcoeff = r.coeff(x, rdeg);
@@ -410,12 +414,12 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return remainder of a(x) and b(x) in Q[x] */
-ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
+ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
 {
 	if (b.is_zero())
 		throw(std::overflow_error("rem: division by zero"));
-	if (is_ex_exactly_of_type(a, numeric)) {
-		if  (is_ex_exactly_of_type(b, numeric))
+	if (is_exactly_a<numeric>(a)) {
+		if  (is_exactly_a<numeric>(b))
 			return _ex0;
 		else
 			return a;
@@ -434,7 +438,7 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
 	int bdeg = b.degree(x);
 	int rdeg = r.degree(x);
 	ex blcoeff = b.expand().coeff(x, bdeg);
-	bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
+	bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
 	while (rdeg >= bdeg) {
 		ex term, rcoeff = r.coeff(x, rdeg);
 		if (blcoeff_is_numeric)
@@ -459,32 +463,32 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
  *  @param a rational function in x
  *  @param x a is a function of x
  *  @return decomposed function. */
-ex decomp_rational(const ex &a, const symbol &x)
+ex decomp_rational(const ex &a, const ex &x)
 {
 	ex nd = numer_denom(a);
 	ex numer = nd.op(0), denom = nd.op(1);
 	ex q = quo(numer, denom, x);
-	if (is_ex_exactly_of_type(q, fail))
+	if (is_exactly_a<fail>(q))
 		return a;
 	else
 		return q + rem(numer, denom, x) / denom;
 }
 
 
-/** Pseudo-remainder of polynomials a(x) and b(x) in Z[x].
+/** Pseudo-remainder of polynomials a(x) and b(x) in Q[x].
  *
  *  @param a  first polynomial in x (dividend)
  *  @param b  second polynomial in x (divisor)
  *  @param x  a and b are polynomials in x
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
- *  @return pseudo-remainder of a(x) and b(x) in Z[x] */
-ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
+ *  @return pseudo-remainder of a(x) and b(x) in Q[x] */
+ex prem(const ex &a, const ex &b, const ex &x, bool check_args)
 {
 	if (b.is_zero())
 		throw(std::overflow_error("prem: division by zero"));
-	if (is_ex_exactly_of_type(a, numeric)) {
-		if (is_ex_exactly_of_type(b, numeric))
+	if (is_exactly_a<numeric>(a)) {
+		if (is_exactly_a<numeric>(b))
 			return _ex0;
 		else
 			return b;
@@ -523,20 +527,20 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
 }
 
 
-/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Z[x].
+/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Q[x].
  *
  *  @param a  first polynomial in x (dividend)
  *  @param b  second polynomial in x (divisor)
  *  @param x  a and b are polynomials in x
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
- *  @return sparse pseudo-remainder of a(x) and b(x) in Z[x] */
-ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args)
+ *  @return sparse pseudo-remainder of a(x) and b(x) in Q[x] */
+ex sprem(const ex &a, const ex &b, const ex &x, bool check_args)
 {
 	if (b.is_zero())
 		throw(std::overflow_error("prem: division by zero"));
-	if (is_ex_exactly_of_type(a, numeric)) {
-		if (is_ex_exactly_of_type(b, numeric))
+	if (is_exactly_a<numeric>(a)) {
+		if (is_exactly_a<numeric>(b))
 			return _ex0;
 		else
 			return b;
@@ -590,10 +594,10 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 		q = _ex0;
 		return true;
 	}
-	if (is_ex_exactly_of_type(b, numeric)) {
+	if (is_exactly_a<numeric>(b)) {
 		q = a / b;
 		return true;
-	} else if (is_ex_exactly_of_type(a, numeric))
+	} else if (is_exactly_a<numeric>(a))
 		return false;
 #if FAST_COMPARE
 	if (a.is_equal(b)) {
@@ -606,7 +610,7 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 		throw(std::invalid_argument("divide: arguments must be polynomials over the rationals"));
 
 	// Find first symbol
-	const symbol *x;
+	ex x;
 	if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
 		throw(std::invalid_argument("invalid expression in divide()"));
 
@@ -616,26 +620,26 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 		q = _ex0;
 		return true;
 	}
-	int bdeg = b.degree(*x);
-	int rdeg = r.degree(*x);
-	ex blcoeff = b.expand().coeff(*x, bdeg);
-	bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
+	int bdeg = b.degree(x);
+	int rdeg = r.degree(x);
+	ex blcoeff = b.expand().coeff(x, bdeg);
+	bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
 	exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
 	while (rdeg >= bdeg) {
-		ex term, rcoeff = r.coeff(*x, rdeg);
+		ex term, rcoeff = r.coeff(x, rdeg);
 		if (blcoeff_is_numeric)
 			term = rcoeff / blcoeff;
 		else
 			if (!divide(rcoeff, blcoeff, term, false))
 				return false;
-		term *= power(*x, rdeg - bdeg);
+		term *= power(x, rdeg - bdeg);
 		v.push_back(term);
 		r -= (term * b).expand();
 		if (r.is_zero()) {
 			q = (new add(v))->setflag(status_flags::dynallocated);
 			return true;
 		}
-		rdeg = r.degree(*x);
+		rdeg = r.degree(x);
 	}
 	return false;
 }
@@ -664,7 +668,7 @@ typedef std::map<ex2, exbool, ex2_less> ex2_exbool_remember;
 /** Exact polynomial division of a(X) by b(X) in Z[X].
  *  This functions works like divide() but the input and output polynomials are
  *  in Z[X] instead of Q[X] (i.e. they have integer coefficients). Unlike
- *  divide(), it doesn´t check whether the input polynomials really are integer
+ *  divide(), it doesn't check whether the input polynomials really are integer
  *  polynomials, so be careful of what you pass in. Also, you have to run
  *  get_symbol_stats() over the input polynomials before calling this function
  *  and pass an iterator to the first element of the sym_desc vector. This
@@ -686,8 +690,8 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 		q = a;
 		return true;
 	}
-	if (is_ex_exactly_of_type(a, numeric)) {
-		if (is_ex_exactly_of_type(b, numeric)) {
+	if (is_exactly_a<numeric>(a)) {
+		if (is_exactly_a<numeric>(b)) {
 			q = a / b;
 			return q.info(info_flags::integer);
 		} else
@@ -711,10 +715,10 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 #endif
 
 	// Main symbol
-	const symbol *x = var->sym;
+	const ex &x = var->sym;
 
 	// Compare degrees
-	int adeg = a.degree(*x), bdeg = b.degree(*x);
+	int adeg = a.degree(x), bdeg = b.degree(x);
 	if (bdeg > adeg)
 		return false;
 
@@ -729,12 +733,12 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 	numeric point = _num0;
 	ex c;
 	for (i=0; i<=adeg; i++) {
-		ex bs = b.subs(*x == point);
+		ex bs = b.subs(x == point, subs_options::no_pattern);
 		while (bs.is_zero()) {
 			point += _num1;
-			bs = b.subs(*x == point);
+			bs = b.subs(x == point, subs_options::no_pattern);
 		}
-		if (!divide_in_z(a.subs(*x == point), bs, c, var+1))
+		if (!divide_in_z(a.subs(x == point, subs_options::no_pattern), bs, c, var+1))
 			return false;
 		alpha.push_back(point);
 		u.push_back(c);
@@ -764,9 +768,9 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 	// Convert from Newton form to standard form
 	c = v[adeg];
 	for (k=adeg-1; k>=0; k--)
-		c = c * (*x - alpha[k]) + v[k];
+		c = c * (x - alpha[k]) + v[k];
 
-	if (c.degree(*x) == (adeg - bdeg)) {
+	if (c.degree(x) == (adeg - bdeg)) {
 		q = c.expand();
 		return true;
 	} else
@@ -780,13 +784,13 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 		return true;
 	int rdeg = adeg;
 	ex eb = b.expand();
-	ex blcoeff = eb.coeff(*x, bdeg);
+	ex blcoeff = eb.coeff(x, bdeg);
 	exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
 	while (rdeg >= bdeg) {
-		ex term, rcoeff = r.coeff(*x, rdeg);
+		ex term, rcoeff = r.coeff(x, rdeg);
 		if (!divide_in_z(rcoeff, blcoeff, term, var+1))
 			break;
-		term = (term * power(*x, rdeg - bdeg)).expand();
+		term = (term * power(x, rdeg - bdeg)).expand();
 		v.push_back(term);
 		r -= (term * eb).expand();
 		if (r.is_zero()) {
@@ -796,7 +800,7 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 #endif
 			return true;
 		}
-		rdeg = r.degree(*x);
+		rdeg = r.degree(x);
 	}
 #if USE_REMEMBER
 	dr_remember[ex2(a, b)] = exbool(q, false);
@@ -812,21 +816,21 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
  */
 
 /** Compute unit part (= sign of leading coefficient) of a multivariate
- *  polynomial in Z[x]. The product of unit part, content part, and primitive
+ *  polynomial in Q[x]. The product of unit part, content part, and primitive
  *  part is the polynomial itself.
  *
  *  @param x  variable in which to compute the unit part
  *  @return unit part
  *  @see ex::content, ex::primpart */
-ex ex::unit(const symbol &x) const
+ex ex::unit(const ex &x) const
 {
 	ex c = expand().lcoeff(x);
-	if (is_ex_exactly_of_type(c, numeric))
+	if (is_exactly_a<numeric>(c))
 		return c < _ex0 ? _ex_1 : _ex1;
 	else {
-		const symbol *y;
+		ex y;
 		if (get_first_symbol(c, y))
-			return c.unit(*y);
+			return c.unit(y);
 		else
 			throw(std::invalid_argument("invalid expression in unit()"));
 	}
@@ -834,23 +838,24 @@ ex ex::unit(const symbol &x) const
 
 
 /** Compute content part (= unit normal GCD of all coefficients) of a
- *  multivariate polynomial in Z[x].  The product of unit part, content part,
+ *  multivariate polynomial in Q[x]. The product of unit part, content part,
  *  and primitive part is the polynomial itself.
  *
  *  @param x  variable in which to compute the content part
  *  @return content part
  *  @see ex::unit, ex::primpart */
-ex ex::content(const symbol &x) const
+ex ex::content(const ex &x) const
 {
 	if (is_zero())
 		return _ex0;
-	if (is_ex_exactly_of_type(*this, numeric))
+	if (is_exactly_a<numeric>(*this))
 		return info(info_flags::negative) ? -*this : *this;
 	ex e = expand();
 	if (e.is_zero())
 		return _ex0;
 
-	// First, try the integer content
+	// First, divide out the integer content (which we can calculate very efficiently).
+	// If the leading coefficient of the quotient is an integer, we are done.
 	ex c = e.integer_content();
 	ex r = e / c;
 	ex lcoeff = r.lcoeff(x);
@@ -858,60 +863,60 @@ ex ex::content(const symbol &x) const
 		return c;
 
 	// GCD of all coefficients
-	int deg = e.degree(x);
-	int ldeg = e.ldegree(x);
+	int deg = r.degree(x);
+	int ldeg = r.ldegree(x);
 	if (deg == ldeg)
-		return e.lcoeff(x) / e.unit(x);
-	c = _ex0;
+		return lcoeff * c;
+	ex cont = _ex0;
 	for (int i=ldeg; i<=deg; i++)
-		c = gcd(e.coeff(x, i), c, NULL, NULL, false);
-	return c;
+		cont = gcd(r.coeff(x, i), cont, NULL, NULL, false);
+	return cont * c;
 }
 
 
-/** Compute primitive part of a multivariate polynomial in Z[x].
+/** Compute primitive part of a multivariate polynomial in Q[x].
  *  The product of unit part, content part, and primitive part is the
  *  polynomial itself.
  *
  *  @param x  variable in which to compute the primitive part
  *  @return primitive part
  *  @see ex::unit, ex::content */
-ex ex::primpart(const symbol &x) const
+ex ex::primpart(const ex &x) const
 {
 	if (is_zero())
 		return _ex0;
-	if (is_ex_exactly_of_type(*this, numeric))
+	if (is_exactly_a<numeric>(*this))
 		return _ex1;
 
 	ex c = content(x);
 	if (c.is_zero())
 		return _ex0;
 	ex u = unit(x);
-	if (is_ex_exactly_of_type(c, numeric))
+	if (is_exactly_a<numeric>(c))
 		return *this / (c * u);
 	else
 		return quo(*this, c * u, x, false);
 }
 
 
-/** Compute primitive part of a multivariate polynomial in Z[x] when the
+/** Compute primitive part of a multivariate polynomial in Q[x] when the
  *  content part is already known. This function is faster in computing the
  *  primitive part than the previous function.
  *
  *  @param x  variable in which to compute the primitive part
  *  @param c  previously computed content part
  *  @return primitive part */
-ex ex::primpart(const symbol &x, const ex &c) const
+ex ex::primpart(const ex &x, const ex &c) const
 {
 	if (is_zero())
 		return _ex0;
 	if (c.is_zero())
 		return _ex0;
-	if (is_ex_exactly_of_type(*this, numeric))
+	if (is_exactly_a<numeric>(*this))
 		return _ex1;
 
 	ex u = unit(x);
-	if (is_ex_exactly_of_type(c, numeric))
+	if (is_exactly_a<numeric>(c))
 		return *this / (c * u);
 	else
 		return quo(*this, c * u, x, false);
@@ -922,206 +927,6 @@ ex ex::primpart(const symbol &x, const ex &c) const
  *  GCD of multivariate polynomials
  */
 
-/** Compute GCD of polynomials in Q[X] using the Euclidean algorithm (not
- *  really suited for multivariate GCDs). This function is only provided for
- *  testing purposes.
- *
- *  @param a  first multivariate polynomial
- *  @param b  second multivariate polynomial
- *  @param x  pointer to symbol (main variable) in which to compute the GCD in
- *  @return the GCD as a new expression
- *  @see gcd */
-
-static ex eu_gcd(const ex &a, const ex &b, const symbol *x)
-{
-//std::clog << "eu_gcd(" << a << "," << b << ")\n";
-
-	// Sort c and d so that c has higher degree
-	ex c, d;
-	int adeg = a.degree(*x), bdeg = b.degree(*x);
-	if (adeg >= bdeg) {
-		c = a;
-		d = b;
-	} else {
-		c = b;
-		d = a;
-	}
-
-	// Normalize in Q[x]
-	c = c / c.lcoeff(*x);
-	d = d / d.lcoeff(*x);
-
-	// Euclidean algorithm
-	ex r;
-	for (;;) {
-//std::clog << " d = " << d << endl;
-		r = rem(c, d, *x, false);
-		if (r.is_zero())
-			return d / d.lcoeff(*x);
-		c = d;
-		d = r;
-	}
-}
-
-
-/** Compute GCD of multivariate polynomials using the Euclidean PRS algorithm
- *  with pseudo-remainders ("World's Worst GCD Algorithm", staying in Z[X]).
- *  This function is only provided for testing purposes.
- *
- *  @param a  first multivariate polynomial
- *  @param b  second multivariate polynomial
- *  @param x  pointer to symbol (main variable) in which to compute the GCD in
- *  @return the GCD as a new expression
- *  @see gcd */
-
-static ex euprem_gcd(const ex &a, const ex &b, const symbol *x)
-{
-//std::clog << "euprem_gcd(" << a << "," << b << ")\n";
-
-	// Sort c and d so that c has higher degree
-	ex c, d;
-	int adeg = a.degree(*x), bdeg = b.degree(*x);
-	if (adeg >= bdeg) {
-		c = a;
-		d = b;
-	} else {
-		c = b;
-		d = a;
-	}
-
-	// Calculate GCD of contents
-	ex gamma = gcd(c.content(*x), d.content(*x), NULL, NULL, false);
-
-	// Euclidean algorithm with pseudo-remainders
-	ex r;
-	for (;;) {
-//std::clog << " d = " << d << endl;
-		r = prem(c, d, *x, false);
-		if (r.is_zero())
-			return d.primpart(*x) * gamma;
-		c = d;
-		d = r;
-	}
-}
-
-
-/** Compute GCD of multivariate polynomials using the primitive Euclidean
- *  PRS algorithm (complete content removal at each step). This function is
- *  only provided for testing purposes.
- *
- *  @param a  first multivariate polynomial
- *  @param b  second multivariate polynomial
- *  @param x  pointer to symbol (main variable) in which to compute the GCD in
- *  @return the GCD as a new expression
- *  @see gcd */
-
-static ex peu_gcd(const ex &a, const ex &b, const symbol *x)
-{
-//std::clog << "peu_gcd(" << a << "," << b << ")\n";
-
-	// Sort c and d so that c has higher degree
-	ex c, d;
-	int adeg = a.degree(*x), bdeg = b.degree(*x);
-	int ddeg;
-	if (adeg >= bdeg) {
-		c = a;
-		d = b;
-		ddeg = bdeg;
-	} else {
-		c = b;
-		d = a;
-		ddeg = adeg;
-	}
-
-	// Remove content from c and d, to be attached to GCD later
-	ex cont_c = c.content(*x);
-	ex cont_d = d.content(*x);
-	ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
-	if (ddeg == 0)
-		return gamma;
-	c = c.primpart(*x, cont_c);
-	d = d.primpart(*x, cont_d);
-
-	// Euclidean algorithm with content removal
-	ex r;
-	for (;;) {
-//std::clog << " d = " << d << endl;
-		r = prem(c, d, *x, false);
-		if (r.is_zero())
-			return gamma * d;
-		c = d;
-		d = r.primpart(*x);
-	}
-}
-
-
-/** Compute GCD of multivariate polynomials using the reduced PRS algorithm.
- *  This function is only provided for testing purposes.
- *
- *  @param a  first multivariate polynomial
- *  @param b  second multivariate polynomial
- *  @param x  pointer to symbol (main variable) in which to compute the GCD in
- *  @return the GCD as a new expression
- *  @see gcd */
-
-static ex red_gcd(const ex &a, const ex &b, const symbol *x)
-{
-//std::clog << "red_gcd(" << a << "," << b << ")\n";
-
-	// Sort c and d so that c has higher degree
-	ex c, d;
-	int adeg = a.degree(*x), bdeg = b.degree(*x);
-	int cdeg, ddeg;
-	if (adeg >= bdeg) {
-		c = a;
-		d = b;
-		cdeg = adeg;
-		ddeg = bdeg;
-	} else {
-		c = b;
-		d = a;
-		cdeg = bdeg;
-		ddeg = adeg;
-	}
-
-	// Remove content from c and d, to be attached to GCD later
-	ex cont_c = c.content(*x);
-	ex cont_d = d.content(*x);
-	ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
-	if (ddeg == 0)
-		return gamma;
-	c = c.primpart(*x, cont_c);
-	d = d.primpart(*x, cont_d);
-
-	// First element of divisor sequence
-	ex r, ri = _ex1;
-	int delta = cdeg - ddeg;
-
-	for (;;) {
-		// Calculate polynomial pseudo-remainder
-//std::clog << " d = " << d << endl;
-		r = prem(c, d, *x, false);
-		if (r.is_zero())
-			return gamma * d.primpart(*x);
-		c = d;
-		cdeg = ddeg;
-
-		if (!divide(r, pow(ri, delta), d, false))
-			throw(std::runtime_error("invalid expression in red_gcd(), division failed"));
-		ddeg = d.degree(*x);
-		if (ddeg == 0) {
-			if (is_ex_exactly_of_type(r, numeric))
-				return gamma;
-			else
-				return gamma * r.primpart(*x);
-		}
-
-		ri = c.expand().lcoeff(*x);
-		delta = cdeg - ddeg;
-	}
-}
-
-
 /** Compute GCD of multivariate polynomials using the subresultant PRS
  *  algorithm. This function is used internally by gcd().
  *
@@ -1133,13 +938,12 @@ static ex red_gcd(const ex &a, const ex &b, const symbol *x)
 
 static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
 {
-//std::clog << "sr_gcd(" << a << "," << b << ")\n";
 #if STATISTICS
 	sr_gcd_called++;
 #endif
 
 	// The first symbol is our main variable
-	const symbol &x = *(var->sym);
+	const ex &x = var->sym;
 
 	// Sort c and d so that c has higher degree
 	ex c, d;
@@ -1165,34 +969,31 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
 		return gamma;
 	c = c.primpart(x, cont_c);
 	d = d.primpart(x, cont_d);
-//std::clog << " content " << gamma << " removed, continuing with sr_gcd(" << c << "," << d << ")\n";
 
 	// First element of subresultant sequence
 	ex r = _ex0, ri = _ex1, psi = _ex1;
 	int delta = cdeg - ddeg;
 
 	for (;;) {
+
 		// Calculate polynomial pseudo-remainder
-//std::clog << " start of loop, psi = " << psi << ", calculating pseudo-remainder...\n";
-//std::clog << " d = " << d << endl;
 		r = prem(c, d, x, false);
 		if (r.is_zero())
 			return gamma * d.primpart(x);
+
 		c = d;
 		cdeg = ddeg;
-//std::clog << " dividing...\n";
 		if (!divide_in_z(r, ri * pow(psi, delta), d, var))
 			throw(std::runtime_error("invalid expression in sr_gcd(), division failed"));
 		ddeg = d.degree(x);
 		if (ddeg == 0) {
-			if (is_ex_exactly_of_type(r, numeric))
+			if (is_exactly_a<numeric>(r))
 				return gamma;
 			else
 				return gamma * r.primpart(x);
 		}
 
 		// Next element of subresultant sequence
-//std::clog << " calculating next subresultant...\n";
 		ri = c.expand().lcoeff(x);
 		if (delta == 1)
 			psi = ri;
@@ -1206,28 +1007,26 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
 /** Return maximum (absolute value) coefficient of a polynomial.
  *  This function is used internally by heur_gcd().
  *
- *  @param e  expanded multivariate polynomial
  *  @return maximum coefficient
  *  @see heur_gcd */
-numeric ex::max_coefficient(void) const
+numeric ex::max_coefficient() const
 {
-	GINAC_ASSERT(bp!=0);
 	return bp->max_coefficient();
 }
 
 /** Implementation ex::max_coefficient().
  *  @see heur_gcd */
-numeric basic::max_coefficient(void) const
+numeric basic::max_coefficient() const
 {
 	return _num1;
 }
 
-numeric numeric::max_coefficient(void) const
+numeric numeric::max_coefficient() const
 {
 	return abs(*this);
 }
 
-numeric add::max_coefficient(void) const
+numeric add::max_coefficient() const
 {
 	epvector::const_iterator it = seq.begin();
 	epvector::const_iterator itend = seq.end();
@@ -1244,7 +1043,7 @@ numeric add::max_coefficient(void) const
 	return cur_max;
 }
 
-numeric mul::max_coefficient(void) const
+numeric mul::max_coefficient() const
 {
 #ifdef DO_GINAC_ASSERT
 	epvector::const_iterator it = seq.begin();
@@ -1313,7 +1112,7 @@ ex mul::smod(const numeric &xi) const
 
 
 /** xi-adic polynomial interpolation */
-static ex interpolate(const ex &gamma, const numeric &xi, const symbol &x, int degree_hint = 1)
+static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint = 1)
 {
 	exvector g; g.reserve(degree_hint);
 	ex e = gamma;
@@ -1346,7 +1145,6 @@ class gcdheu_failed {};
  *  @exception gcdheu_failed() */
 static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
 {
-//std::clog << "heur_gcd(" << a << "," << b << ")\n";
 #if STATISTICS
 	heur_gcd_called++;
 #endif
@@ -1356,7 +1154,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 		return (new fail())->setflag(status_flags::dynallocated);
 
 	// GCD of two numeric values -> CLN
-	if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) {
+	if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
 		numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
 		if (ca)
 			*ca = ex_to<numeric>(a) / g;
@@ -1366,7 +1164,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 	}
 
 	// The first symbol is our main variable
-	const symbol &x = *(var->sym);
+	const ex &x = var->sym;
 
 	// Remove integer content
 	numeric gc = gcd(a.integer_content(), b.integer_content());
@@ -1387,14 +1185,13 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 	// 6 tries maximum
 	for (int t=0; t<6; t++) {
 		if (xi.int_length() * maxdeg > 100000) {
-//std::clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << std::endl;
 			throw gcdheu_failed();
 		}
 
 		// Apply evaluation homomorphism and calculate GCD
 		ex cp, cq;
-		ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), &cp, &cq, var+1).expand();
-		if (!is_ex_exactly_of_type(gamma, fail)) {
+		ex gamma = heur_gcd(p.subs(x == xi, subs_options::no_pattern), q.subs(x == xi, subs_options::no_pattern), &cp, &cq, var+1).expand();
+		if (!is_exactly_a<fail>(gamma)) {
 
 			// Reconstruct polynomial from GCD of mapped polynomials
 			ex g = interpolate(gamma, xi, x, maxdeg);
@@ -1407,39 +1204,11 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 			if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
 				g *= gc;
 				ex lc = g.lcoeff(x);
-				if (is_ex_exactly_of_type(lc, numeric) && ex_to<numeric>(lc).is_negative())
+				if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
 					return -g;
 				else
 					return g;
 			}
-#if 0
-			cp = interpolate(cp, xi, x);
-			if (divide_in_z(cp, p, g, var)) {
-				if (divide_in_z(g, q, cb ? *cb : dummy, var)) {
-					g *= gc;
-					if (ca)
-						*ca = cp;
-					ex lc = g.lcoeff(x);
-					if (is_ex_exactly_of_type(lc, numeric) && ex_to<numeric>(lc).is_negative())
-						return -g;
-					else
-						return g;
-				}
-			}
-			cq = interpolate(cq, xi, x);
-			if (divide_in_z(cq, q, g, var)) {
-				if (divide_in_z(g, p, ca ? *ca : dummy, var)) {
-					g *= gc;
-					if (cb)
-						*cb = cq;
-					ex lc = g.lcoeff(x);
-					if (is_ex_exactly_of_type(lc, numeric) && ex_to<numeric>(lc).is_negative())
-						return -g;
-					else
-						return g;
-				}
-			}
-#endif
 		}
 
 		// Next evaluation point
@@ -1454,18 +1223,19 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
  *
  *  @param a  first multivariate polynomial
  *  @param b  second multivariate polynomial
+ *  @param ca pointer to expression that will receive the cofactor of a, or NULL
+ *  @param cb pointer to expression that will receive the cofactor of b, or NULL
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return the GCD as a new expression */
 ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
 {
-//std::clog << "gcd(" << a << "," << b << ")\n";
 #if STATISTICS
 	gcd_called++;
 #endif
 
 	// GCD of numerics -> CLN
-	if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) {
+	if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
 		numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
 		if (ca || cb) {
 			if (g.is_zero()) {
@@ -1489,15 +1259,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
 	}
 
 	// Partially factored cases (to avoid expanding large expressions)
-	if (is_ex_exactly_of_type(a, mul)) {
-		if (is_ex_exactly_of_type(b, mul) && b.nops() > a.nops())
+	if (is_exactly_a<mul>(a)) {
+		if (is_exactly_a<mul>(b) && b.nops() > a.nops())
 			goto factored_b;
 factored_a:
-		unsigned num = a.nops();
+		size_t num = a.nops();
 		exvector g; g.reserve(num);
 		exvector acc_ca; acc_ca.reserve(num);
 		ex part_b = b;
-		for (unsigned i=0; i<num; i++) {
+		for (size_t i=0; i<num; i++) {
 			ex part_ca, part_cb;
 			g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args));
 			acc_ca.push_back(part_ca);
@@ -1508,15 +1278,15 @@ factored_a:
 		if (cb)
 			*cb = part_b;
 		return (new mul(g))->setflag(status_flags::dynallocated);
-	} else if (is_ex_exactly_of_type(b, mul)) {
-		if (is_ex_exactly_of_type(a, mul) && a.nops() > b.nops())
+	} else if (is_exactly_a<mul>(b)) {
+		if (is_exactly_a<mul>(a) && a.nops() > b.nops())
 			goto factored_a;
 factored_b:
-		unsigned num = b.nops();
+		size_t num = b.nops();
 		exvector g; g.reserve(num);
 		exvector acc_cb; acc_cb.reserve(num);
 		ex part_a = a;
-		for (unsigned i=0; i<num; i++) {
+		for (size_t i=0; i<num; i++) {
 			ex part_ca, part_cb;
 			g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args));
 			acc_cb.push_back(part_cb);
@@ -1531,9 +1301,9 @@ factored_b:
 
 #if FAST_COMPARE
 	// Input polynomials of the form poly^n are sometimes also trivial
-	if (is_ex_exactly_of_type(a, power)) {
+	if (is_exactly_a<power>(a)) {
 		ex p = a.op(0);
-		if (is_ex_exactly_of_type(b, power)) {
+		if (is_exactly_a<power>(b)) {
 			if (p.is_equal(b.op(0))) {
 				// a = p^n, b = p^m, gcd = p^min(n, m)
 				ex exp_a = a.op(1), exp_b = b.op(1);
@@ -1561,7 +1331,7 @@ factored_b:
 				return p;
 			}
 		}
-	} else if (is_ex_exactly_of_type(b, power)) {
+	} else if (is_exactly_a<power>(b)) {
 		ex p = b.op(0);
 		if (p.is_equal(a)) {
 			// a = p, b = p^n, gcd = p
@@ -1613,7 +1383,7 @@ factored_b:
 
 	// The symbol with least degree is our main variable
 	sym_desc_vec::const_iterator var = sym_stats.begin();
-	const symbol &x = *(var->sym);
+	const ex &x = var->sym;
 
 	// Cancel trivial common factor
 	int ldeg_a = var->ldeg_a;
@@ -1621,20 +1391,17 @@ factored_b:
 	int min_ldeg = std::min(ldeg_a,ldeg_b);
 	if (min_ldeg > 0) {
 		ex common = power(x, min_ldeg);
-//std::clog << "trivial common factor " << common << std::endl;
 		return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
 	}
 
 	// Try to eliminate variables
 	if (var->deg_a == 0) {
-//std::clog << "eliminating variable " << x << " from b" << std::endl;
 		ex c = bex.content(x);
 		ex g = gcd(aex, c, ca, cb, false);
 		if (cb)
 			*cb *= bex.unit(x) * bex.primpart(x, c);
 		return g;
 	} else if (var->deg_b == 0) {
-//std::clog << "eliminating variable " << x << " from a" << std::endl;
 		ex c = aex.content(x);
 		ex g = gcd(c, bex, ca, cb, false);
 		if (ca)
@@ -1642,25 +1409,17 @@ factored_b:
 		return g;
 	}
 
-	ex g;
-#if 1
 	// Try heuristic algorithm first, fall back to PRS if that failed
+	ex g;
 	try {
 		g = heur_gcd(aex, bex, ca, cb, var);
 	} catch (gcdheu_failed) {
 		g = fail();
 	}
-	if (is_ex_exactly_of_type(g, fail)) {
-//std::clog << "heuristics failed" << std::endl;
+	if (is_exactly_a<fail>(g)) {
 #if STATISTICS
 		heur_gcd_failed++;
 #endif
-#endif
-//		g = heur_gcd(aex, bex, ca, cb, var);
-//		g = eu_gcd(aex, bex, &x);
-//		g = euprem_gcd(aex, bex, &x);
-//		g = peu_gcd(aex, bex, &x);
-//		g = red_gcd(aex, bex, &x);
 		g = sr_gcd(aex, bex, var);
 		if (g.is_equal(_ex1)) {
 			// Keep cofactors factored if possible
@@ -1674,7 +1433,6 @@ factored_b:
 			if (cb)
 				divide(bex, g, *cb, false);
 		}
-#if 1
 	} else {
 		if (g.is_equal(_ex1)) {
 			// Keep cofactors factored if possible
@@ -1684,7 +1442,7 @@ factored_b:
 				*cb = b;
 		}
 	}
-#endif
+
 	return g;
 }
 
@@ -1698,7 +1456,7 @@ factored_b:
  *  @return the LCM as a new expression */
 ex lcm(const ex &a, const ex &b, bool check_args)
 {
-	if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric))
+	if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b))
 		return lcm(ex_to<numeric>(a), ex_to<numeric>(b));
 	if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
 		throw(std::invalid_argument("lcm: arguments must be polynomials over the rationals"));
@@ -1741,10 +1499,11 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
 	return res;
 }
 
+
 /** Compute a square-free factorization of a multivariate polynomial in Q[X].
  *
  *  @param a  multivariate polynomial over Q[X]
- *  @param x  lst of variables to factor in, may be left empty for autodetection
+ *  @param l  lst of variables to factor in, may be left empty for autodetection
  *  @return   a square-free factorization of \p a.
  *
  * \note
@@ -1778,7 +1537,7 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
  */
 ex sqrfree(const ex &a, const lst &l)
 {
-	if (is_a<numeric>(a) ||     // algorithm does not trap a==0
+	if (is_exactly_a<numeric>(a) ||     // algorithm does not trap a==0
 	    is_a<symbol>(a))        // shortcut
 		return a;
 
@@ -1791,7 +1550,7 @@ ex sqrfree(const ex &a, const lst &l)
 		get_symbol_stats(a, _ex0, sdv);
 		sym_desc_vec::const_iterator it = sdv.begin(), itend = sdv.end();
 		while (it != itend) {
-			args.append(*it->sym);
+			args.append(it->sym);
 			++it;
 		}
 	} else {
@@ -1799,7 +1558,7 @@ ex sqrfree(const ex &a, const lst &l)
 	}
 
 	// Find the symbol to factor in at this stage
-	if (!is_ex_of_type(args.op(0), symbol))
+	if (!is_a<symbol>(args.op(0)))
 		throw (std::runtime_error("sqrfree(): invalid factorization variable"));
 	const symbol &x = ex_to<symbol>(args.op(0));
 
@@ -1808,7 +1567,7 @@ ex sqrfree(const ex &a, const lst &l)
 	const ex tmp = multiply_lcm(a,lcm);
 
 	// find the factors
-	exvector factors = sqrfree_yun(tmp,x);
+	exvector factors = sqrfree_yun(tmp, x);
 
 	// construct the next list of symbols with the first element popped
 	lst newargs = args;
@@ -1842,6 +1601,7 @@ ex sqrfree(const ex &a, const lst &l)
 	return result *	lcm.inverse();
 }
 
+
 /** Compute square-free partial fraction decomposition of rational function
  *  a(x).
  *
@@ -1863,15 +1623,15 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 	// Factorize denominator and compute cofactors
 	exvector yun = sqrfree_yun(denom, x);
 //clog << "yun factors: " << exprseq(yun) << endl;
-	unsigned num_yun = yun.size();
+	size_t num_yun = yun.size();
 	exvector factor; factor.reserve(num_yun);
 	exvector cofac; cofac.reserve(num_yun);
-	for (unsigned i=0; i<num_yun; i++) {
+	for (size_t i=0; i<num_yun; i++) {
 		if (!yun[i].is_equal(_ex1)) {
-			for (unsigned j=0; j<=i; j++) {
+			for (size_t j=0; j<=i; j++) {
 				factor.push_back(pow(yun[i], j+1));
 				ex prod = _ex1;
-				for (unsigned k=0; k<num_yun; k++) {
+				for (size_t k=0; k<num_yun; k++) {
 					if (k == i)
 						prod *= pow(yun[k], i-j);
 					else
@@ -1881,7 +1641,7 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 			}
 		}
 	}
-	unsigned num_factors = factor.size();
+	size_t num_factors = factor.size();
 //clog << "factors  : " << exprseq(factor) << endl;
 //clog << "cofactors: " << exprseq(cofac) << endl;
 
@@ -1890,7 +1650,7 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 	matrix sys(max_denom_deg + 1, num_factors);
 	matrix rhs(max_denom_deg + 1, 1);
 	for (int i=0; i<=max_denom_deg; i++) {
-		for (unsigned j=0; j<num_factors; j++)
+		for (size_t j=0; j<num_factors; j++)
 			sys(i, j) = cofac[j].coeff(x, i);
 		rhs(i, 0) = red_numer.coeff(x, i);
 	}
@@ -1899,13 +1659,13 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 
 	// Solve resulting linear system
 	matrix vars(num_factors, 1);
-	for (unsigned i=0; i<num_factors; i++)
+	for (size_t i=0; i<num_factors; i++)
 		vars(i, 0) = symbol();
 	matrix sol = sys.solve(vars, rhs);
 
 	// Sum up decomposed fractions
 	ex sum = 0;
-	for (unsigned i=0; i<num_factors; i++)
+	for (size_t i=0; i<num_factors; i++)
 		sum += sol(i, 0) / factor[i];
 
 	return red_poly + sum;
@@ -1926,45 +1686,44 @@ 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, exmap & rev_lookup)
 {
-	// Expression already in repl_lst? Then return the assigned symbol
-	for (unsigned i=0; i<repl_lst.nops(); i++)
-		if (repl_lst.op(i).is_equal(e))
-			return sym_lst.op(i);
+	// Expression already replaced? Then return the assigned symbol
+	exmap::const_iterator it = rev_lookup.find(e);
+	if (it != rev_lookup.end())
+		return it->second;
 	
 	// 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, subs_options::no_pattern);
+	repl.insert(std::make_pair(es, e_replaced));
+	rev_lookup.insert(std::make_pair(e_replaced, es));
 	return es;
 }
 
 /** Create a symbol for replacing the expression "e" (or return a previously
- *  assigned symbol). An expression of the form "symbol == expression" is added
- *  to repl_lst and the symbol is returned.
- *  @see basic::to_rational */
-static ex replace_with_symbol(const ex &e, lst &repl_lst)
-{
-	// Expression already in repl_lst? Then return the assigned symbol
-	for (unsigned i=0; i<repl_lst.nops(); i++)
-		if (repl_lst.op(i).op(1).is_equal(e))
-			return repl_lst.op(i).op(0);
+ *  assigned symbol). The symbol and expression are appended to repl, and the
+ *  symbol is returned.
+ *  @see basic::to_rational
+ *  @see basic::to_polynomial */
+static ex replace_with_symbol(const ex & e, exmap & repl)
+{
+	// Expression already replaced? 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(repl_lst);
-	repl_lst.append(es == e_replaced);
+	ex es = (new symbol)->setflag(status_flags::dynallocated);
+	ex e_replaced = e.subs(repl, subs_options::no_pattern);
+	repl.insert(std::make_pair(es, e_replaced));
 	return es;
 }
 
@@ -1979,18 +1738,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, exmap & rev_lookup, 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, rev_lookup), _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, rev_lookup), _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, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
 		}
 	}
 }
@@ -1998,7 +1757,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, exmap & rev_lookup, int level) const
 {
 	return (new lst(*this, _ex1))->setflag(status_flags::dynallocated);
 }
@@ -2008,19 +1767,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, exmap & rev_lookup, 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, rev_lookup);
 	} 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, rev_lookup);
+		ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, rev_lookup);
+		numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup);
 	}
 
 	// Denominator is always a real integer (see numeric::denom())
@@ -2067,13 +1826,20 @@ static ex frac_cancel(const ex &n, const ex &d)
 
 	// Make denominator unit normal (i.e. coefficient of first symbol
 	// as defined by get_first_symbol() is made positive)
-	const symbol *x;
-	if (get_first_symbol(den, x)) {
-		GINAC_ASSERT(is_exactly_a<numeric>(den.unit(*x)));
-		if (ex_to<numeric>(den.unit(*x)).is_negative()) {
+	if (is_exactly_a<numeric>(den)) {
+		if (ex_to<numeric>(den).is_negative()) {
 			num *= _ex_1;
 			den *= _ex_1;
 		}
+	} else {
+		ex x;
+		if (get_first_symbol(den, x)) {
+			GINAC_ASSERT(is_exactly_a<numeric>(den.unit(x)));
+			if (ex_to<numeric>(den.unit(x)).is_negative()) {
+				num *= _ex_1;
+				den *= _ex_1;
+			}
+		}
 	}
 
 	// Return result as list
@@ -2085,10 +1851,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, exmap & rev_lookup, 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, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
 	else if (level == -max_recursion_level)
 		throw(std::runtime_error("max recursion level reached"));
 
@@ -2098,12 +1864,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, rev_lookup, 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, rev_lookup, level-1);
 	nums.push_back(n.op(0));
 	dens.push_back(n.op(1));
 	GINAC_ASSERT(nums.size() == dens.size());
@@ -2144,10 +1910,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, exmap & rev_lookup, 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, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
 	else if (level == -max_recursion_level)
 		throw(std::runtime_error("max recursion level reached"));
 
@@ -2157,12 +1923,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, rev_lookup, 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, rev_lookup, level-1);
 	num.push_back(n.op(0));
 	den.push_back(n.op(1));
 
@@ -2172,20 +1938,20 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
 }
 
 
-/** Implementation of ex::normal() for powers. It normalizes the basis,
+/** Implementation of ex::normal([B) for powers. It normalizes the basis,
  *  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, exmap & rev_lookup, 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, rev_lookup), _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, rev_lookup, level-1);
+	ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, level-1);
 	n_exponent = n_exponent.op(0) / n_exponent.op(1);
 
 	if (n_exponent.info(info_flags::integer)) {
@@ -2206,34 +1972,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, rev_lookup), _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, rev_lookup)))->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, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
 			}
-
-		} else {	// n_exponent not numeric
-
-			// (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);
 		}
 	}
+
+	// (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), repl, rev_lookup), _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, exmap & rev_lookup, int level) const
 {
 	epvector newseq;
 	epvector::const_iterator i = seq.begin(), end = seq.end();
@@ -2244,7 +2008,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, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
 }
 
 
@@ -2262,14 +2026,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, rev_lookup;
 
-	ex e = bp->normal(sym_lst, repl_lst, level);
+	ex e = bp->normal(repl, rev_lookup, 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, subs_options::no_pattern);
 
 	// Convert {numerator, denominator} form back to fraction
 	return e.op(0) / e.op(1);
@@ -2281,18 +2045,18 @@ ex ex::normal(int level) const
  *
  *  @see ex::normal
  *  @return numerator */
-ex ex::numer(void) const
+ex ex::numer() const
 {
-	lst sym_lst, repl_lst;
+	exmap repl, rev_lookup;
 
-	ex e = bp->normal(sym_lst, repl_lst, 0);
+	ex e = bp->normal(repl, rev_lookup, 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, subs_options::no_pattern);
 }
 
 /** Get denominator of an expression. If the expression is not of the normal
@@ -2301,18 +2065,18 @@ ex ex::numer(void) const
  *
  *  @see ex::normal
  *  @return denominator */
-ex ex::denom(void) const
+ex ex::denom() const
 {
-	lst sym_lst, repl_lst;
+	exmap repl, rev_lookup;
 
-	ex e = bp->normal(sym_lst, repl_lst, 0);
+	ex e = bp->normal(repl, rev_lookup, 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, subs_options::no_pattern);
 }
 
 /** Get numerator and denominator of an expression. If the expresison is not
@@ -2321,44 +2085,103 @@ ex ex::denom(void) const
  *
  *  @see ex::normal
  *  @return a list [numerator, denominator] */
-ex ex::numer_denom(void) const
+ex ex::numer_denom() const
 {
-	lst sym_lst, repl_lst;
+	exmap repl, rev_lookup;
 
-	ex e = bp->normal(sym_lst, repl_lst, 0);
+	ex e = bp->normal(repl, rev_lookup, 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, subs_options::no_pattern);
 }
 
 
 /** Rationalization of non-rational functions.
- *  This function converts a general expression to a rational polynomial
+ *  This function converts a general expression to a rational function
  *  by replacing all non-rational subexpressions (like non-rational numbers,
  *  non-integer powers or functions like sin(), cos() etc.) to temporary
  *  symbols. This makes it possible to use functions like gcd() and divide()
  *  on non-rational functions by applying to_rational() on the arguments,
  *  calling the desired function and re-substituting the temporary symbols
  *  in the result. To make the last step possible, all temporary symbols and
- *  their associated expressions are collected in the list specified by the
- *  repl_lst parameter in the form {symbol == expression}, ready to be passed
- *  as an argument to ex::subs().
+ *  their associated expressions are collected in the map specified by the
+ *  repl parameter, ready to be passed as an argument to ex::subs().
  *
- *  @param repl_lst collects a list of all temporary symbols and their replacements
+ *  @param repl collects all temporary symbols and their replacements
  *  @return rationalized expression */
-ex basic::to_rational(lst &repl_lst) const
+ex ex::to_rational(exmap & repl) const
+{
+	return bp->to_rational(repl);
+}
+
+// GiNaC 1.1 compatibility function
+ex ex::to_rational(lst & repl_lst) const
+{
+	// Convert lst to exmap
+	exmap m;
+	for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
+		m.insert(std::make_pair(it->op(0), it->op(1)));
+
+	ex ret = bp->to_rational(m);
+
+	// Convert exmap back to lst
+	repl_lst.remove_all();
+	for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
+		repl_lst.append(it->first == it->second);
+
+	return ret;
+}
+
+ex ex::to_polynomial(exmap & repl) const
+{
+	return bp->to_polynomial(repl);
+}
+
+// GiNaC 1.1 compatibility function
+ex ex::to_polynomial(lst & repl_lst) const
+{
+	// Convert lst to exmap
+	exmap m;
+	for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
+		m.insert(std::make_pair(it->op(0), it->op(1)));
+
+	ex ret = bp->to_polynomial(m);
+
+	// Convert exmap back to lst
+	repl_lst.remove_all();
+	for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
+		repl_lst.append(it->first == it->second);
+
+	return ret;
+}
+
+/** Default implementation of ex::to_rational(). This replaces the object with
+ *  a temporary symbol. */
+ex basic::to_rational(exmap & repl) const
+{
+	return replace_with_symbol(*this, repl);
+}
+
+ex basic::to_polynomial(exmap & repl) const
 {
-	return replace_with_symbol(*this, repl_lst);
+	return replace_with_symbol(*this, repl);
 }
 
 
 /** Implementation of ex::to_rational() for symbols. This returns the
  *  unmodified symbol. */
-ex symbol::to_rational(lst &repl_lst) const
+ex symbol::to_rational(exmap & repl) const
+{
+	return *this;
+}
+
+/** Implementation of ex::to_polynomial() for symbols. This returns the
+ *  unmodified symbol. */
+ex symbol::to_polynomial(exmap & repl) const
 {
 	return *this;
 }
@@ -2367,17 +2190,35 @@ ex symbol::to_rational(lst &repl_lst) const
 /** Implementation of ex::to_rational() for a numeric. It splits complex
  *  numbers into re+I*im and replaces I and non-rational real numbers with a
  *  temporary symbol. */
-ex numeric::to_rational(lst &repl_lst) const
+ex numeric::to_rational(exmap & repl) const
 {
 	if (is_real()) {
 		if (!is_rational())
-			return replace_with_symbol(*this, repl_lst);
+			return replace_with_symbol(*this, repl);
+	} else { // complex
+		numeric re = real();
+		numeric im = imag();
+		ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
+		ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
+		return re_ex + im_ex * replace_with_symbol(I, repl);
+	}
+	return *this;
+}
+
+/** Implementation of ex::to_polynomial() for a numeric. It splits complex
+ *  numbers into re+I*im and replaces I and non-integer real numbers with a
+ *  temporary symbol. */
+ex numeric::to_polynomial(exmap & repl) const
+{
+	if (is_real()) {
+		if (!is_integer())
+			return replace_with_symbol(*this, repl);
 	} else { // complex
 		numeric re = real();
 		numeric im = imag();
-		ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst);
-		ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst);
-		return re_ex + im_ex * replace_with_symbol(I, repl_lst);
+		ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl);
+		ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl);
+		return re_ex + im_ex * replace_with_symbol(I, repl);
 	}
 	return *this;
 }
@@ -2385,26 +2226,54 @@ ex numeric::to_rational(lst &repl_lst) const
 
 /** Implementation of ex::to_rational() for powers. It replaces non-integer
  *  powers by temporary symbols. */
-ex power::to_rational(lst &repl_lst) const
+ex power::to_rational(exmap & repl) const
 {
 	if (exponent.info(info_flags::integer))
-		return power(basis.to_rational(repl_lst), exponent);
+		return power(basis.to_rational(repl), exponent);
 	else
-		return replace_with_symbol(*this, repl_lst);
+		return replace_with_symbol(*this, repl);
+}
+
+/** Implementation of ex::to_polynomial() for powers. It replaces non-posint
+ *  powers by temporary symbols. */
+ex power::to_polynomial(exmap & repl) const
+{
+	if (exponent.info(info_flags::posint))
+		return power(basis.to_rational(repl), exponent);
+	else
+		return replace_with_symbol(*this, repl);
 }
 
 
 /** Implementation of ex::to_rational() for expairseqs. */
-ex expairseq::to_rational(lst &repl_lst) const
+ex expairseq::to_rational(exmap & repl) const
+{
+	epvector s;
+	s.reserve(seq.size());
+	epvector::const_iterator i = seq.begin(), end = seq.end();
+	while (i != end) {
+		s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_rational(repl)));
+		++i;
+	}
+	ex oc = overall_coeff.to_rational(repl);
+	if (oc.info(info_flags::numeric))
+		return thisexpairseq(s, overall_coeff);
+	else
+		s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1));
+	return thisexpairseq(s, default_overall_coeff());
+}
+
+/** Implementation of ex::to_polynomial() for expairseqs. */
+ex expairseq::to_polynomial(exmap & repl) const
 {
 	epvector s;
 	s.reserve(seq.size());
 	epvector::const_iterator i = seq.begin(), end = seq.end();
 	while (i != end) {
-		s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_rational(repl_lst)));
+		s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_polynomial(repl)));
 		++i;
 	}
-	ex oc = overall_coeff.to_rational(repl_lst);
+	ex oc = overall_coeff.to_polynomial(repl);
 	if (oc.info(info_flags::numeric))
 		return thisexpairseq(s, overall_coeff);
 	else
@@ -2413,4 +2282,135 @@ ex expairseq::to_rational(lst &repl_lst) const
 }
 
 
+/** Remove the common factor in the terms of a sum 'e' by calculating the GCD,
+ *  and multiply it into the expression 'factor' (which needs to be initialized
+ *  to 1, unless you're accumulating factors). */
+static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
+{
+	if (is_exactly_a<add>(e)) {
+
+		size_t num = e.nops();
+		exvector terms; terms.reserve(num);
+		ex gc;
+
+		// Find the common GCD
+		for (size_t i=0; i<num; i++) {
+			ex x = e.op(i).to_polynomial(repl);
+
+			if (is_exactly_a<add>(x) || is_exactly_a<mul>(x)) {
+				ex f = 1;
+				x = find_common_factor(x, f, repl);
+				x *= f;
+			}
+
+			if (i == 0)
+				gc = x;
+			else
+				gc = gcd(gc, x);
+
+			terms.push_back(x);
+		}
+
+		if (gc.is_equal(_ex1))
+			return e;
+
+		// The GCD is the factor we pull out
+		factor *= gc;
+
+		// Now divide all terms by the GCD
+		for (size_t i=0; i<num; i++) {
+			ex x;
+
+			// Try to avoid divide() because it expands the polynomial
+			ex &t = terms[i];
+			if (is_exactly_a<mul>(t)) {
+				for (size_t j=0; j<t.nops(); j++) {
+					if (t.op(j).is_equal(gc)) {
+						exvector v; v.reserve(t.nops());
+						for (size_t k=0; k<t.nops(); k++) {
+							if (k == j)
+								v.push_back(_ex1);
+							else
+								v.push_back(t.op(k));
+						}
+						t = (new mul(v))->setflag(status_flags::dynallocated);
+						goto term_done;
+					}
+				}
+			}
+
+			divide(t, gc, x);
+			t = x;
+term_done:	;
+		}
+		return (new add(terms))->setflag(status_flags::dynallocated);
+
+	} else if (is_exactly_a<mul>(e)) {
+
+		size_t num = e.nops();
+		exvector v; v.reserve(num);
+
+		for (size_t i=0; i<num; i++)
+			v.push_back(find_common_factor(e.op(i), factor, repl));
+
+		return (new mul(v))->setflag(status_flags::dynallocated);
+
+	} else if (is_exactly_a<power>(e)) {
+
+		return e.to_polynomial(repl);
+
+	} else
+		return e;
+}
+
+
+/** Collect common factors in sums. This converts expressions like
+ *  'a*(b*x+b*y)' to 'a*b*(x+y)'. */
+ex collect_common_factors(const ex & e)
+{
+	if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
+
+		exmap repl;
+		ex factor = 1;
+		ex r = find_common_factor(e, factor, repl);
+		return factor.subs(repl, subs_options::no_pattern) * r.subs(repl, subs_options::no_pattern);
+
+	} else
+		return e;
+}
+
+
+/** Resultant of two expressions e1,e2 with respect to symbol s.
+ *  Method: Compute determinant of Sylvester matrix of e1,e2,s.  */
+ex resultant(const ex & e1, const ex & e2, const ex & s)
+{
+	const ex ee1 = e1.expand();
+	const ex ee2 = e2.expand();
+	if (!ee1.info(info_flags::polynomial) ||
+	    !ee2.info(info_flags::polynomial))
+		throw(std::runtime_error("resultant(): arguments must be polynomials"));
+
+	const int h1 = ee1.degree(s);
+	const int l1 = ee1.ldegree(s);
+	const int h2 = ee2.degree(s);
+	const int l2 = ee2.ldegree(s);
+
+	const int msize = h1 + h2;
+	matrix m(msize, msize);
+
+	for (int l = h1; l >= l1; --l) {
+		const ex e = ee1.coeff(s, l);
+		for (int k = 0; k < h2; ++k)
+			m(k, k+h1-l) = e;
+	}
+	for (int l = h2; l >= l2; --l) {
+		const ex e = ee2.coeff(s, l);
+		for (int k = 0; k < h1; ++k)
+			m(k+h2, k+h2-l) = e;
+	}
+
+	return m.determinant();
+}
+
+
 } // namespace GiNaC