X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=a80e15084c3affee06491acdba6013aaae0ee2ee;hp=64c34ac290037348d835420d49ff1b2bf019000a;hb=c3a7dda76e171bbd1e1418b3dab56df49920f721;hpb=22abfbe8c78e339188096a5bf749a7c2d4f0a368

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 64c34ac2..a80e1508 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-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2020 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
@@ -23,9 +23,6 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <algorithm>
-#include <map>
-
 #include "normal.h"
 #include "basic.h"
 #include "ex.h"
@@ -44,6 +41,10 @@
 #include "pseries.h"
 #include "symbol.h"
 #include "utils.h"
+#include "polynomial/chinrem_gcd.h"
+
+#include <algorithm>
+#include <map>
 
 namespace GiNaC {
 
@@ -119,6 +120,11 @@ static bool get_first_symbol(const ex &e, ex &x)
  *
  *  @see get_symbol_stats */
 struct sym_desc {
+	/** Initialize symbol, leave other variables uninitialized */
+	sym_desc(const ex& s)
+	  : sym(s), deg_a(0), deg_b(0), ldeg_a(0), ldeg_b(0), max_deg(0), max_lcnops(0)
+	{ }
+
 	/** Reference to symbol */
 	ex sym;
 
@@ -140,7 +146,7 @@ struct sym_desc {
 	/** Maximum number of terms of leading coefficient of symbol in both polynomials */
 	size_t max_lcnops;
 
-	/** Commparison operator for sorting */
+	/** Comparison operator for sorting */
 	bool operator<(const sym_desc &x) const
 	{
 		if (max_deg == x.max_deg)
@@ -156,15 +162,11 @@ 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 ex &s, sym_desc_vec &v)
 {
-	sym_desc_vec::const_iterator it = v.begin(), itend = v.end();
-	while (it != itend) {
-		if (it->sym.is_equal(s))  // If it's already in there, don't add it a second time
+	for (auto & it : v)
+		if (it.sym.is_equal(s))  // If it's already in there, don't add it a second time
 			return;
-		++it;
-	}
-	sym_desc d;
-	d.sym = s;
-	v.push_back(d);
+
+	v.push_back(sym_desc(s));
 }
 
 // Collect all symbols of an expression (used internally by get_symbol_stats())
@@ -194,28 +196,26 @@ static void collect_symbols(const ex &e, sym_desc_vec &v)
  *  @param v  vector of sym_desc structs (filled in) */
 static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 {
-	collect_symbols(a.eval(), v);   // eval() to expand assigned symbols
-	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);
-		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;
+	collect_symbols(a, v);
+	collect_symbols(b, v);
+	for (auto & it : v) {
+		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);
 	}
 	std::sort(v.begin(), v.end());
 
 #if 0
 	std::clog << "Symbols:\n";
-	it = v.begin(); itend = v.end();
+	auto 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 << std::endl;
+		std::clog << "  lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << std::endl;
 		++it;
 	}
 #endif
@@ -270,9 +270,15 @@ 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 (lcm.is_equal(*_num1_p))
+		// e * 1 -> e;
+		return e;
+
 	if (is_exactly_a<mul>(e)) {
+		// (a*b*...)*lcm -> (a*lcma)*(b*lcmb)*...*(lcm/(lcma*lcmb*...))
 		size_t num = e.nops();
-		exvector v; v.reserve(num + 1);
+		exvector v;
+		v.reserve(num + 1);
 		numeric lcm_accum = *_num1_p;
 		for (size_t i=0; i<num; i++) {
 			numeric op_lcm = lcmcoeff(e.op(i), *_num1_p);
@@ -280,20 +286,26 @@ static ex multiply_lcm(const ex &e, const numeric &lcm)
 			lcm_accum *= op_lcm;
 		}
 		v.push_back(lcm / lcm_accum);
-		return (new mul(v))->setflag(status_flags::dynallocated);
+		return dynallocate<mul>(v);
 	} else if (is_exactly_a<add>(e)) {
+		// (a+b+...)*lcm -> a*lcm+b*lcm+...
 		size_t num = e.nops();
-		exvector v; v.reserve(num);
+		exvector v;
+		v.reserve(num);
 		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);
+		return dynallocate<add>(v);
 	} 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));
-	} else
-		return e * lcm;
+		if (!is_a<symbol>(e.op(0))) {
+			// (b^e)*lcm -> (b*lcm^(1/e))^e if lcm^(1/e) ∈ ℚ (i.e. not a float)
+			// but not for symbolic b, as evaluation would undo this again
+			numeric root_of_lcm = lcm.power(ex_to<numeric>(e.op(1)).inverse());
+			if (root_of_lcm.is_rational())
+				return pow(multiply_lcm(e.op(0), root_of_lcm), e.op(1));
+		}
+	}
+	// can't recurse down into e
+	return dynallocate<mul>(e, lcm);
 }
 
 
@@ -320,15 +332,12 @@ numeric numeric::integer_content() const
 
 numeric add::integer_content() const
 {
-	epvector::const_iterator it = seq.begin();
-	epvector::const_iterator itend = seq.end();
 	numeric c = *_num0_p, l = *_num1_p;
-	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).numer(), c);
-		l = lcm(ex_to<numeric>(it->coeff).denom(), l);
-		it++;
+	for (auto & it : seq) {
+		GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
+		GINAC_ASSERT(is_exactly_a<numeric>(it.coeff));
+		c = gcd(ex_to<numeric>(it.coeff).numer(), c);
+		l = lcm(ex_to<numeric>(it.coeff).denom(), l);
 	}
 	GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
 	c = gcd(ex_to<numeric>(overall_coeff).numer(), c);
@@ -339,11 +348,8 @@ numeric add::integer_content() const
 numeric mul::integer_content() const
 {
 #ifdef DO_GINAC_ASSERT
-	epvector::const_iterator it = seq.begin();
-	epvector::const_iterator itend = seq.end();
-	while (it != itend) {
-		GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*it)));
-		++it;
+	for (auto & it : seq) {
+		GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
 	}
 #endif // def DO_GINAC_ASSERT
 	GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
@@ -392,16 +398,16 @@ ex quo(const ex &a, const ex &b, const ex &x, bool check_args)
 			term = rcoeff / blcoeff;
 		else {
 			if (!divide(rcoeff, blcoeff, term, false))
-				return (new fail())->setflag(status_flags::dynallocated);
+				return dynallocate<fail>();
 		}
-		term *= power(x, rdeg - bdeg);
+		term *= pow(x, rdeg - bdeg);
 		v.push_back(term);
 		r -= (term * b).expand();
 		if (r.is_zero())
 			break;
 		rdeg = r.degree(x);
 	}
-	return (new add(v))->setflag(status_flags::dynallocated);
+	return dynallocate<add>(v);
 }
 
 
@@ -445,9 +451,9 @@ ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
 			term = rcoeff / blcoeff;
 		else {
 			if (!divide(rcoeff, blcoeff, term, false))
-				return (new fail())->setflag(status_flags::dynallocated);
+				return dynallocate<fail>();
 		}
-		term *= power(x, rdeg - bdeg);
+		term *= pow(x, rdeg - bdeg);
 		r -= (term * b).expand();
 		if (r.is_zero())
 			break;
@@ -507,23 +513,23 @@ ex prem(const ex &a, const ex &b, const ex &x, bool check_args)
 		if (bdeg == 0)
 			eb = _ex0;
 		else
-			eb -= blcoeff * power(x, bdeg);
+			eb -= blcoeff * pow(x, bdeg);
 	} else
 		blcoeff = _ex1;
 
 	int delta = rdeg - bdeg + 1, i = 0;
 	while (rdeg >= bdeg && !r.is_zero()) {
 		ex rlcoeff = r.coeff(x, rdeg);
-		ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand();
+		ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand();
 		if (rdeg == 0)
 			r = _ex0;
 		else
-			r -= rlcoeff * power(x, rdeg);
+			r -= rlcoeff * pow(x, rdeg);
 		r = (blcoeff * r).expand() - term;
 		rdeg = r.degree(x);
 		i++;
 	}
-	return power(blcoeff, delta - i) * r;
+	return pow(blcoeff, delta - i) * r;
 }
 
 
@@ -559,17 +565,17 @@ ex sprem(const ex &a, const ex &b, const ex &x, bool check_args)
 		if (bdeg == 0)
 			eb = _ex0;
 		else
-			eb -= blcoeff * power(x, bdeg);
+			eb -= blcoeff * pow(x, bdeg);
 	} else
 		blcoeff = _ex1;
 
 	while (rdeg >= bdeg && !r.is_zero()) {
 		ex rlcoeff = r.coeff(x, rdeg);
-		ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand();
+		ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand();
 		if (rdeg == 0)
 			r = _ex0;
 		else
-			r -= rlcoeff * power(x, rdeg);
+			r -= rlcoeff * pow(x, rdeg);
 		r = (blcoeff * r).expand() - term;
 		rdeg = r.degree(x);
 	}
@@ -614,6 +620,73 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 	if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
 		throw(std::invalid_argument("invalid expression in divide()"));
 
+	// Try to avoid expanding partially factored expressions.
+	if (is_exactly_a<mul>(b)) {
+	// Divide sequentially by each term
+		ex rem_new, rem_old = a;
+		for (size_t i=0; i < b.nops(); i++) {
+			if (! divide(rem_old, b.op(i), rem_new, false))
+				return false;
+			rem_old = rem_new;
+		}
+		q = rem_new;
+		return true;
+	} else if (is_exactly_a<power>(b)) {
+		const ex& bb(b.op(0));
+		int exp_b = ex_to<numeric>(b.op(1)).to_int();
+		ex rem_new, rem_old = a;
+		for (int i=exp_b; i>0; i--) {
+			if (! divide(rem_old, bb, rem_new, false))
+				return false;
+			rem_old = rem_new;
+		}
+		q = rem_new;
+		return true;
+	} 
+	
+	if (is_exactly_a<mul>(a)) {
+		// Divide sequentially each term. If some term in a is divisible 
+		// by b we are done... and if not, we can't really say anything.
+		size_t i;
+		ex rem_i;
+		bool divisible_p = false;
+		for (i=0; i < a.nops(); ++i) {
+			if (divide(a.op(i), b, rem_i, false)) {
+				divisible_p = true;
+				break;
+			}
+		}
+		if (divisible_p) {
+			exvector resv;
+			resv.reserve(a.nops());
+			for (size_t j=0; j < a.nops(); j++) {
+				if (j==i)
+					resv.push_back(rem_i);
+				else
+					resv.push_back(a.op(j));
+			}
+			q = dynallocate<mul>(resv);
+			return true;
+		}
+	} else if (is_exactly_a<power>(a)) {
+		// The base itself might be divisible by b, in that case we don't
+		// need to expand a
+		const ex& ab(a.op(0));
+		int a_exp = ex_to<numeric>(a.op(1)).to_int();
+		ex rem_i;
+		if (divide(ab, b, rem_i, false)) {
+			q = rem_i * pow(ab, a_exp - 1);
+			return true;
+		}
+// code below is commented-out because it leads to a significant slowdown
+//		for (int i=2; i < a_exp; i++) {
+//			if (divide(power(ab, i), b, rem_i, false)) {
+//				q = rem_i*power(ab, a_exp - i);
+//				return true;
+//			}
+//		} // ... so we *really* need to expand expression.
+	}
+	
 	// Polynomial long division (recursive)
 	ex r = a.expand();
 	if (r.is_zero()) {
@@ -632,11 +705,11 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 		else
 			if (!divide(rcoeff, blcoeff, term, false))
 				return false;
-		term *= power(x, rdeg - bdeg);
+		term *= pow(x, rdeg - bdeg);
 		v.push_back(term);
 		r -= (term * b).expand();
 		if (r.is_zero()) {
-			q = (new add(v))->setflag(status_flags::dynallocated);
+			q = dynallocate<add>(v);
 			return true;
 		}
 		rdeg = r.degree(x);
@@ -714,6 +787,31 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 	}
 #endif
 
+	if (is_exactly_a<power>(b)) {
+		const ex& bb(b.op(0));
+		ex qbar = a;
+		int exp_b = ex_to<numeric>(b.op(1)).to_int();
+		for (int i=exp_b; i>0; i--) {
+			if (!divide_in_z(qbar, bb, q, var))
+				return false;
+			qbar = q;
+		}
+		return true;
+	}
+
+	if (is_exactly_a<mul>(b)) {
+		ex qbar = a;
+		for (const auto & it : b) {
+			sym_desc_vec sym_stats;
+			get_symbol_stats(a, it, sym_stats);
+			if (!divide_in_z(qbar, it, q, sym_stats.begin()))
+				return false;
+
+			qbar = q;
+		}
+		return true;
+	}
+
 	// Main symbol
 	const ex &x = var->sym;
 
@@ -790,11 +888,11 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 		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 * pow(x, rdeg - bdeg)).expand();
 		v.push_back(term);
 		r -= (term * eb).expand();
 		if (r.is_zero()) {
-			q = (new add(v))->setflag(status_flags::dynallocated);
+			q = dynallocate<add>(v);
 #if USE_REMEMBER
 			dr_remember[ex2(a, b)] = exbool(q, true);
 #endif
@@ -868,7 +966,7 @@ ex ex::content(const ex &x) const
 		return lcoeff * c / lcoeff.unit(x);
 	ex cont = _ex0;
 	for (int i=ldeg; i<=deg; i++)
-		cont = gcd(r.coeff(x, i), cont, NULL, NULL, false);
+		cont = gcd(r.coeff(x, i), cont, nullptr, nullptr, false);
 	return cont * c;
 }
 
@@ -1008,7 +1106,7 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
 	// 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);
+	ex gamma = gcd(cont_c, cont_d, nullptr, nullptr, false);
 	if (ddeg == 0)
 		return gamma;
 	c = c.primpart(x, cont_c);
@@ -1072,17 +1170,14 @@ numeric numeric::max_coefficient() const
 
 numeric add::max_coefficient() const
 {
-	epvector::const_iterator it = seq.begin();
-	epvector::const_iterator itend = seq.end();
 	GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
 	numeric cur_max = abs(ex_to<numeric>(overall_coeff));
-	while (it != itend) {
+	for (auto & it : seq) {
 		numeric a;
-		GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
-		a = abs(ex_to<numeric>(it->coeff));
+		GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
+		a = abs(ex_to<numeric>(it.coeff));
 		if (a > cur_max)
 			cur_max = a;
-		it++;
 	}
 	return cur_max;
 }
@@ -1090,11 +1185,8 @@ numeric add::max_coefficient() const
 numeric mul::max_coefficient() const
 {
 #ifdef DO_GINAC_ASSERT
-	epvector::const_iterator it = seq.begin();
-	epvector::const_iterator itend = seq.end();
-	while (it != itend) {
-		GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*it)));
-		it++;
+	for (auto & it : seq) {
+		GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
 	}
 #endif // def DO_GINAC_ASSERT
 	GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
@@ -1122,36 +1214,30 @@ ex add::smod(const numeric &xi) const
 {
 	epvector newseq;
 	newseq.reserve(seq.size()+1);
-	epvector::const_iterator it = seq.begin();
-	epvector::const_iterator itend = seq.end();
-	while (it != itend) {
-		GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
-		numeric coeff = GiNaC::smod(ex_to<numeric>(it->coeff), xi);
+	for (auto & it : seq) {
+		GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
+		numeric coeff = GiNaC::smod(ex_to<numeric>(it.coeff), xi);
 		if (!coeff.is_zero())
-			newseq.push_back(expair(it->rest, coeff));
-		it++;
+			newseq.push_back(expair(it.rest, coeff));
 	}
 	GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
 	numeric coeff = GiNaC::smod(ex_to<numeric>(overall_coeff), xi);
-	return (new add(newseq,coeff))->setflag(status_flags::dynallocated);
+	return dynallocate<add>(std::move(newseq), coeff);
 }
 
 ex mul::smod(const numeric &xi) const
 {
 #ifdef DO_GINAC_ASSERT
-	epvector::const_iterator it = seq.begin();
-	epvector::const_iterator itend = seq.end();
-	while (it != itend) {
-		GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*it)));
-		it++;
+	for (auto & it : seq) {
+		GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
 	}
 #endif // def DO_GINAC_ASSERT
-	mul * mulcopyp = new mul(*this);
+	mul & mulcopy = dynallocate<mul>(*this);
 	GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
-	mulcopyp->overall_coeff = GiNaC::smod(ex_to<numeric>(overall_coeff),xi);
-	mulcopyp->clearflag(status_flags::evaluated);
-	mulcopyp->clearflag(status_flags::hash_calculated);
-	return mulcopyp->setflag(status_flags::dynallocated);
+	mulcopy.overall_coeff = GiNaC::smod(ex_to<numeric>(overall_coeff),xi);
+	mulcopy.clearflag(status_flags::evaluated);
+	mulcopy.clearflag(status_flags::hash_calculated);
+	return mulcopy;
 }
 
 
@@ -1163,10 +1249,10 @@ static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degre
 	numeric rxi = xi.inverse();
 	for (int i=0; !e.is_zero(); i++) {
 		ex gi = e.smod(xi);
-		g.push_back(gi * power(x, i));
+		g.push_back(gi * pow(x, i));
 		e = (e - gi) * rxi;
 	}
-	return (new add(g))->setflag(status_flags::dynallocated);
+	return dynallocate<add>(g);
 }
 
 /** Exception thrown by heur_gcd() to signal failure. */
@@ -1177,17 +1263,19 @@ class gcdheu_failed {};
  *  polynomials and an iterator to the first element of the sym_desc vector
  *  passed in. This function is used internally by gcd().
  *
- *  @param a  first multivariate polynomial (expanded)
- *  @param b  second multivariate polynomial (expanded)
- *  @param ca  cofactor of polynomial a (returned), NULL to suppress
+ *  @param a  first integer multivariate polynomial (expanded)
+ *  @param b  second integer multivariate polynomial (expanded)
+ *  @param ca  cofactor of polynomial a (returned), nullptr to suppress
  *             calculation of cofactor
- *  @param cb  cofactor of polynomial b (returned), NULL to suppress
+ *  @param cb  cofactor of polynomial b (returned), nullptr to suppress
  *             calculation of cofactor
  *  @param var iterator to first element of vector of sym_desc structs
- *  @return the GCD as a new expression
+ *  @param res the GCD (returned)
+ *  @return true if GCD was computed, false otherwise.
  *  @see gcd
  *  @exception gcdheu_failed() */
-static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
+static bool heur_gcd_z(ex& res, const ex &a, const ex &b, ex *ca, ex *cb,
+	               sym_desc_vec::const_iterator var)
 {
 #if STATISTICS
 	heur_gcd_called++;
@@ -1195,7 +1283,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 
 	// Algorithm only works for non-vanishing input polynomials
 	if (a.is_zero() || b.is_zero())
-		return (new fail())->setflag(status_flags::dynallocated);
+		return false;
 
 	// GCD of two numeric values -> CLN
 	if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
@@ -1204,7 +1292,8 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 			*ca = ex_to<numeric>(a) / g;
 		if (cb)
 			*cb = ex_to<numeric>(b) / g;
-		return g;
+		res = g;
+		return true;
 	}
 
 	// The first symbol is our main variable
@@ -1234,9 +1323,13 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 
 		// Apply evaluation homomorphism and calculate GCD
 		ex cp, cq;
-		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)) {
-
+		ex gamma;
+		bool found = heur_gcd_z(gamma,
+			                p.subs(x == xi, subs_options::no_pattern),
+			                q.subs(x == xi, subs_options::no_pattern),
+				        &cp, &cq, var+1);
+		if (found) {
+			gamma = gamma.expand();
 			// Reconstruct polynomial from GCD of mapped polynomials
 			ex g = interpolate(gamma, xi, x, maxdeg);
 
@@ -1247,29 +1340,96 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 			ex dummy;
 			if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
 				g *= gc;
-				return g;
+				res = g;
+				return true;
 			}
 		}
 
 		// Next evaluation point
 		xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011));
 	}
-	return (new fail())->setflag(status_flags::dynallocated);
+	return false;
+}
+
+/** Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
+ *  get_symbol_stats() must have been called previously with the input
+ *  polynomials and an iterator to the first element of the sym_desc vector
+ *  passed in. This function is used internally by gcd().
+ *
+ *  @param a  first rational multivariate polynomial (expanded)
+ *  @param b  second rational multivariate polynomial (expanded)
+ *  @param ca  cofactor of polynomial a (returned), nullptr to suppress
+ *             calculation of cofactor
+ *  @param cb  cofactor of polynomial b (returned), nullptr to suppress
+ *             calculation of cofactor
+ *  @param var iterator to first element of vector of sym_desc structs
+ *  @param res the GCD (returned)
+ *  @return true if GCD was computed, false otherwise.
+ *  @see heur_gcd_z
+ *  @see gcd
+ */
+static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb,
+	             sym_desc_vec::const_iterator var)
+{
+	if (a.info(info_flags::integer_polynomial) && 
+	    b.info(info_flags::integer_polynomial)) {
+		try {
+			return heur_gcd_z(res, a, b, ca, cb, var);
+		} catch (gcdheu_failed) {
+			return false;
+		}
+	}
+
+	// convert polynomials to Z[X]
+	const numeric a_lcm = lcm_of_coefficients_denominators(a);
+	const numeric ab_lcm = lcmcoeff(b, a_lcm);
+
+	const ex ai = a*ab_lcm;
+	const ex bi = b*ab_lcm;
+	if (!ai.info(info_flags::integer_polynomial))
+		throw std::logic_error("heur_gcd: not an integer polynomial [1]");
+
+	if (!bi.info(info_flags::integer_polynomial))
+		throw std::logic_error("heur_gcd: not an integer polynomial [2]");
+
+	bool found = false;
+	try {
+		found = heur_gcd_z(res, ai, bi, ca, cb, var);
+	} catch (gcdheu_failed) {
+		return false;
+	}
+	
+	// GCD is not unique, it's defined up to a unit (i.e. invertible
+	// element). If the coefficient ring is a field, every its element is
+	// invertible, so one can multiply the polynomial GCD with any element
+	// of the coefficient field. We use this ambiguity to make cofactors
+	// integer polynomials.
+	if (found)
+		res /= ab_lcm;
+	return found;
 }
 
 
+// gcd helper to handle partially factored polynomials (to avoid expanding
+// large expressions). At least one of the arguments should be a power.
+static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb);
+
+// gcd helper to handle partially factored polynomials (to avoid expanding
+// large expressions). At least one of the arguments should be a product.
+static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb);
+
 /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X)
  *  and b(X) in Z[X]. Optionally also compute the cofactors of a and b,
  *  defined by a = ca * gcd(a, b) and b = cb * gcd(a, b).
  *
  *  @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 ca pointer to expression that will receive the cofactor of a, or nullptr
+ *  @param cb pointer to expression that will receive the cofactor of b, or nullptr
  *  @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)
+ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
 {
 #if STATISTICS
 	gcd_called++;
@@ -1300,153 +1460,17 @@ 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_exactly_a<mul>(a)) {
-		if (is_exactly_a<mul>(b) && b.nops() > a.nops())
-			goto factored_b;
-factored_a:
-		size_t num = a.nops();
-		exvector g; g.reserve(num);
-		exvector acc_ca; acc_ca.reserve(num);
-		ex part_b = b;
-		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);
-			part_b = part_cb;
-		}
-		if (ca)
-			*ca = (new mul(acc_ca))->setflag(status_flags::dynallocated);
-		if (cb)
-			*cb = part_b;
-		return (new mul(g))->setflag(status_flags::dynallocated);
-	} else if (is_exactly_a<mul>(b)) {
-		if (is_exactly_a<mul>(a) && a.nops() > b.nops())
-			goto factored_a;
-factored_b:
-		size_t num = b.nops();
-		exvector g; g.reserve(num);
-		exvector acc_cb; acc_cb.reserve(num);
-		ex part_a = a;
-		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);
-			part_a = part_ca;
-		}
-		if (ca)
-			*ca = part_a;
-		if (cb)
-			*cb = (new mul(acc_cb))->setflag(status_flags::dynallocated);
-		return (new mul(g))->setflag(status_flags::dynallocated);
-	}
-
+	if (!(options & gcd_options::no_part_factored)) {
+		if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b))
+			return gcd_pf_mul(a, b, ca, cb);
 #if FAST_COMPARE
-	// Input polynomials of the form poly^n are sometimes also trivial
-	if (is_exactly_a<power>(a)) {
-		ex p = a.op(0);
-		const ex& exp_a = a.op(1);
-		if (is_exactly_a<power>(b)) {
-			ex pb = b.op(0);
-			const ex& exp_b = b.op(1);
-			if (p.is_equal(pb)) {
-				// a = p^n, b = p^m, gcd = p^min(n, m)
-				if (exp_a < exp_b) {
-					if (ca)
-						*ca = _ex1;
-					if (cb)
-						*cb = power(p, exp_b - exp_a);
-					return power(p, exp_a);
-				} else {
-					if (ca)
-						*ca = power(p, exp_a - exp_b);
-					if (cb)
-						*cb = _ex1;
-					return power(p, exp_b);
-				}
-			} else {
-				ex p_co, pb_co;
-				ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args);
-				if (p_gcd.is_equal(_ex1)) {
-					// a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==>
-					// gcd(a,b) = 1
-					if (ca)
-						*ca = a;
-					if (cb)
-						*cb = b;
-					return _ex1;
-					// XXX: do I need to check for p_gcd = -1?
-				} else {
-					// there are common factors:
-					// a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
-					// gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
-					if (exp_a < exp_b) {
-						return power(p_gcd, exp_a)*
-							gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
-					} else {
-						return power(p_gcd, exp_b)*
-							gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
-					}
-				} // p_gcd.is_equal(_ex1)
-			} // p.is_equal(pb)
-
-		} else {
-			if (p.is_equal(b)) {
-				// a = p^n, b = p, gcd = p
-				if (ca)
-					*ca = power(p, a.op(1) - 1);
-				if (cb)
-					*cb = _ex1;
-				return p;
-			} 
-
-			ex p_co, bpart_co;
-			ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
-
-			if (p_gcd.is_equal(_ex1)) {
-				// a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
-				if (ca)
-					*ca = a;
-				if (cb)
-					*cb = b;
-				return _ex1;
-			} else {
-				// a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
-				return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
-			}
-		} // is_exactly_a<power>(b)
-
-	} else if (is_exactly_a<power>(b)) {
-		ex p = b.op(0);
-		if (p.is_equal(a)) {
-			// a = p, b = p^n, gcd = p
-			if (ca)
-				*ca = _ex1;
-			if (cb)
-				*cb = power(p, b.op(1) - 1);
-			return p;
-		}
-
-		ex p_co, apart_co;
-		const ex& exp_b(b.op(1));
-		ex p_gcd = gcd(a, p, &apart_co, &p_co, false);
-		if (p_gcd.is_equal(_ex1)) {
-			// b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1
-			if (ca)
-				*ca = a;
-			if (cb)
-				*cb = b;
-			return _ex1;
-		} else {
-			// there are common factors:
-			// a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
-
-			return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false);
-		} // p_gcd.is_equal(_ex1)
-	}
+		if (is_exactly_a<power>(a) || is_exactly_a<power>(b))
+			return gcd_pf_pow(a, b, ca, cb);
 #endif
+	}
 
 	// Some trivial cases
-	ex aex = a.expand(), bex = b.expand();
+	ex aex = a.expand();
 	if (aex.is_zero()) {
 		if (ca)
 			*ca = _ex0;
@@ -1454,6 +1478,7 @@ factored_b:
 			*cb = _ex1;
 		return b;
 	}
+	ex bex = b.expand();
 	if (bex.is_zero()) {
 		if (ca)
 			*ca = _ex1;
@@ -1478,11 +1503,70 @@ factored_b:
 	}
 #endif
 
+	if (is_a<symbol>(aex)) {
+		if (! bex.subs(aex==_ex0, subs_options::no_pattern).is_zero()) {
+			if (ca)
+				*ca = a;
+			if (cb)
+				*cb = b;
+			return _ex1;
+		}
+	}
+
+	if (is_a<symbol>(bex)) {
+		if (! aex.subs(bex==_ex0, subs_options::no_pattern).is_zero()) {
+			if (ca)
+				*ca = a;
+			if (cb)
+				*cb = b;
+			return _ex1;
+		}
+	}
+
+	if (is_exactly_a<numeric>(aex)) {
+		numeric bcont = bex.integer_content();
+		numeric g = gcd(ex_to<numeric>(aex), bcont);
+		if (ca)
+			*ca = ex_to<numeric>(aex)/g;
+		if (cb)
+			*cb = bex/g;
+		return g;
+	}
+
+	if (is_exactly_a<numeric>(bex)) {
+		numeric acont = aex.integer_content();
+		numeric g = gcd(ex_to<numeric>(bex), acont);
+		if (ca)
+			*ca = aex/g;
+		if (cb)
+			*cb = ex_to<numeric>(bex)/g;
+		return g;
+	}
+
 	// Gather symbol statistics
 	sym_desc_vec sym_stats;
 	get_symbol_stats(a, b, sym_stats);
 
-	// The symbol with least degree is our main variable
+	// The symbol with least degree which is contained in both polynomials
+	// is our main variable
+	auto vari = sym_stats.begin();
+	while ((vari != sym_stats.end()) && 
+	       (((vari->ldeg_b == 0) && (vari->deg_b == 0)) ||
+	        ((vari->ldeg_a == 0) && (vari->deg_a == 0))))
+		vari++;
+
+	// No common symbols at all, just return 1:
+	if (vari == sym_stats.end()) {
+		// N.B: keep cofactors factored
+		if (ca)
+			*ca = a;
+		if (cb)
+			*cb = b;
+		return _ex1;
+	}
+	// move symbol contained only in one of the polynomials to the end:
+	rotate(sym_stats.begin(), vari, sym_stats.end());
+
 	sym_desc_vec::const_iterator var = sym_stats.begin();
 	const ex &x = var->sym;
 
@@ -1491,19 +1575,19 @@ factored_b:
 	int ldeg_b = var->ldeg_b;
 	int min_ldeg = std::min(ldeg_a,ldeg_b);
 	if (min_ldeg > 0) {
-		ex common = power(x, min_ldeg);
+		ex common = pow(x, min_ldeg);
 		return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
 	}
 
 	// Try to eliminate variables
-	if (var->deg_a == 0) {
+	if (var->deg_a == 0 && var->deg_b != 0 ) {
 		ex bex_u, bex_c, bex_p;
 		bex.unitcontprim(x, bex_u, bex_c, bex_p);
 		ex g = gcd(aex, bex_c, ca, cb, false);
 		if (cb)
 			*cb *= bex_u * bex_p;
 		return g;
-	} else if (var->deg_b == 0) {
+	} else if (var->deg_b == 0 && var->deg_a != 0) {
 		ex aex_u, aex_c, aex_p;
 		aex.unitcontprim(x, aex_u, aex_c, aex_p);
 		ex g = gcd(aex_c, bex, ca, cb, false);
@@ -1514,41 +1598,171 @@ factored_b:
 
 	// 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_exactly_a<fail>(g)) {
+	if (!(options & gcd_options::no_heur_gcd)) {
+		bool found = heur_gcd(g, aex, bex, ca, cb, var);
+		if (found) {
+			// heur_gcd have already computed cofactors...
+			if (g.is_equal(_ex1)) {
+				// ... but we want to keep them factored if possible.
+				if (ca)
+					*ca = a;
+				if (cb)
+					*cb = b;
+			}
+			return g;
+		}
 #if STATISTICS
-		heur_gcd_failed++;
+		else {
+			heur_gcd_failed++;
+		}
 #endif
+	}
+	if (options & gcd_options::use_sr_gcd) {
 		g = sr_gcd(aex, bex, var);
-		if (g.is_equal(_ex1)) {
-			// Keep cofactors factored if possible
+	} else {
+		exvector vars;
+		for (std::size_t n = sym_stats.size(); n-- != 0; )
+			vars.push_back(sym_stats[n].sym);
+		g = chinrem_gcd(aex, bex, vars);
+	}
+
+	if (g.is_equal(_ex1)) {
+		// Keep cofactors factored if possible
+		if (ca)
+			*ca = a;
+		if (cb)
+			*cb = b;
+	} else {
+		if (ca)
+			divide(aex, g, *ca, false);
+		if (cb)
+			divide(bex, g, *cb, false);
+	}
+	return g;
+}
+
+// gcd helper to handle partially factored polynomials (to avoid expanding
+// large expressions). Both arguments should be powers.
+static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb)
+{
+	ex p = a.op(0);
+	const ex& exp_a = a.op(1);
+	ex pb = b.op(0);
+	const ex& exp_b = b.op(1);
+
+	// a = p^n, b = p^m, gcd = p^min(n, m)
+	if (p.is_equal(pb)) {
+		if (exp_a < exp_b) {
 			if (ca)
-				*ca = a;
+				*ca = _ex1;
 			if (cb)
-				*cb = b;
+				*cb = pow(p, exp_b - exp_a);
+			return pow(p, exp_a);
 		} else {
 			if (ca)
-				divide(aex, g, *ca, false);
+				*ca = pow(p, exp_a - exp_b);
 			if (cb)
-				divide(bex, g, *cb, false);
+				*cb = _ex1;
+			return pow(p, exp_b);
 		}
-	} else {
-		if (g.is_equal(_ex1)) {
-			// Keep cofactors factored if possible
+	}
+
+	ex p_co, pb_co;
+	ex p_gcd = gcd(p, pb, &p_co, &pb_co, false);
+	// a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> gcd(a,b) = 1
+	if (p_gcd.is_equal(_ex1)) {
 			if (ca)
 				*ca = a;
 			if (cb)
 				*cb = b;
+			return _ex1;
+	}
+
+	// there are common factors:
+	// a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
+	// gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
+	if (exp_a < exp_b) {
+		ex pg =  gcd(pow(p_co, exp_a), pow(p_gcd, exp_b-exp_a)*pow(pb_co, exp_b), ca, cb, false);
+		return pow(p_gcd, exp_a)*pg;
+	} else {
+		ex pg = gcd(pow(p_gcd, exp_a - exp_b)*pow(p_co, exp_a), pow(pb_co, exp_b), ca, cb, false);
+		return pow(p_gcd, exp_b)*pg;
+	}
+}
+
+static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
+{
+	if (is_exactly_a<power>(a) && is_exactly_a<power>(b))
+		return gcd_pf_pow_pow(a, b, ca, cb);
+
+	if (is_exactly_a<power>(b) && (! is_exactly_a<power>(a)))
+		return gcd_pf_pow(b, a, cb, ca);
+
+	GINAC_ASSERT(is_exactly_a<power>(a));
+
+	ex p = a.op(0);
+	const ex& exp_a = a.op(1);
+	if (p.is_equal(b)) {
+		// a = p^n, b = p, gcd = p
+		if (ca)
+			*ca = pow(p, exp_a - 1);
+		if (cb)
+			*cb = _ex1;
+		return p;
+	}
+	if (is_a<symbol>(p)) {
+		// Cancel trivial common factor
+		int ldeg_a = ex_to<numeric>(exp_a).to_int();
+		int ldeg_b = b.ldegree(p);
+		int min_ldeg = std::min(ldeg_a, ldeg_b);
+		if (min_ldeg > 0) {
+			ex common = pow(p, min_ldeg);
+			return gcd(pow(p, ldeg_a - min_ldeg), (b / common).expand(), ca, cb, false) * common;
 		}
 	}
 
-	return g;
+	ex p_co, bpart_co;
+	ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
+
+	if (p_gcd.is_equal(_ex1)) {
+		// a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
+		if (ca)
+			*ca = a;
+		if (cb)
+			*cb = b;
+		return _ex1;
+	}
+	// a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+	ex rg = gcd(pow(p_gcd, exp_a-1)*pow(p_co, exp_a), bpart_co, ca, cb, false);
+	return p_gcd*rg;
 }
 
+static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb)
+{
+	if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b)
+		                 && (b.nops() >  a.nops()))
+		return gcd_pf_mul(b, a, cb, ca);
+
+	if (is_exactly_a<mul>(b) && (!is_exactly_a<mul>(a)))
+		return gcd_pf_mul(b, a, cb, ca);
+
+	GINAC_ASSERT(is_exactly_a<mul>(a));
+	size_t num = a.nops();
+	exvector g; g.reserve(num);
+	exvector acc_ca; acc_ca.reserve(num);
+	ex part_b = b;
+	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, false));
+		acc_ca.push_back(part_ca);
+		part_b = part_cb;
+	}
+	if (ca)
+		*ca = dynallocate<mul>(acc_ca);
+	if (cb)
+		*cb = part_b;
+	return dynallocate<mul>(g);
+}
 
 /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
  *
@@ -1578,34 +1792,47 @@ ex lcm(const ex &a, const ex &b, bool check_args)
  *  Yun's algorithm.  Used internally by sqrfree().
  *
  *  @param a  multivariate polynomial over Z[X], treated here as univariate
- *            polynomial in x.
+ *            polynomial in x (needs not be expanded).
  *  @param x  variable to factor in
- *  @return   vector of factors sorted in ascending degree */
-static exvector sqrfree_yun(const ex &a, const symbol &x)
+ *  @return   vector of expairs (factor, exponent), sorted by exponents */
+static epvector sqrfree_yun(const ex &a, const symbol &x)
 {
-	exvector res;
 	ex w = a;
 	ex z = w.diff(x);
 	ex g = gcd(w, z);
+	if (g.is_zero()) {
+		return epvector{};
+	}
 	if (g.is_equal(_ex1)) {
-		res.push_back(a);
-		return res;
+		return epvector{expair(a, _ex1)};
 	}
-	ex y;
+	epvector results;
+	ex exponent = _ex0;
 	do {
 		w = quo(w, g, x);
-		y = quo(z, g, x);
-		z = y - w.diff(x);
+		if (w.is_zero()) {
+			return results;
+		}
+		z = quo(z, g, x) - w.diff(x);
+		exponent = exponent + 1;
+		if (w.is_equal(x)) {
+			// shortcut for x^n with n ∈ ℕ
+			exponent += quo(z, w.diff(x), x);
+			results.push_back(expair(w, exponent));
+			break;
+		}
 		g = gcd(w, z);
-		res.push_back(g);
+		if (!g.is_equal(_ex1)) {
+			results.push_back(expair(g, exponent));
+		}
 	} while (!z.is_zero());
-	return res;
+	return results;
 }
 
 
 /** Compute a square-free factorization of a multivariate polynomial in Q[X].
  *
- *  @param a  multivariate polynomial over Q[X]
+ *  @param a  multivariate polynomial over Q[X] (needs not be expanded)
  *  @param l  lst of variables to factor in, may be left empty for autodetection
  *  @return   a square-free factorization of \p a.
  *
@@ -1640,8 +1867,8 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
  */
 ex sqrfree(const ex &a, const lst &l)
 {
-	if (is_exactly_a<numeric>(a) ||     // algorithm does not trap a==0
-	    is_a<symbol>(a))        // shortcut
+	if (is_exactly_a<numeric>(a) ||
+	    is_a<symbol>(a))        // shortcuts
 		return a;
 
 	// If no lst of variables to factorize in was specified we have to
@@ -1651,11 +1878,8 @@ ex sqrfree(const ex &a, const lst &l)
 	if (l.nops()==0) {
 		sym_desc_vec sdv;
 		get_symbol_stats(a, _ex0, sdv);
-		sym_desc_vec::const_iterator it = sdv.begin(), itend = sdv.end();
-		while (it != itend) {
-			args.append(it->sym);
-			++it;
-		}
+		for (auto & it : sdv)
+			args.append(it.sym);
 	} else {
 		args = l;
 	}
@@ -1670,37 +1894,32 @@ ex sqrfree(const ex &a, const lst &l)
 	const ex tmp = multiply_lcm(a,lcm);
 
 	// find the factors
-	exvector factors = sqrfree_yun(tmp, x);
+	epvector factors = sqrfree_yun(tmp, x);
 
-	// construct the next list of symbols with the first element popped
-	lst newargs = args;
-	newargs.remove_first();
+	// remove symbol x and proceed recursively with the remaining symbols
+	args.remove_first();
 
 	// recurse down the factors in remaining variables
-	if (newargs.nops()>0) {
-		exvector::iterator i = factors.begin();
-		while (i != factors.end()) {
-			*i = sqrfree(*i, newargs);
-			++i;
-		}
+	if (args.nops()>0) {
+		for (auto & it : factors)
+			it.rest = sqrfree(it.rest, args);
 	}
 
 	// Done with recursion, now construct the final result
 	ex result = _ex1;
-	exvector::const_iterator it = factors.begin(), itend = factors.end();
-	for (int p = 1; it!=itend; ++it, ++p)
-		result *= power(*it, p);
+	for (auto & it : factors)
+		result *= pow(it.rest, it.coeff);
 
 	// Yun's algorithm does not account for constant factors.  (For univariate
 	// polynomials it works only in the monic case.)  We can correct this by
 	// inserting what has been lost back into the result.  For completeness
 	// we'll also have to recurse down that factor in the remaining variables.
-	if (newargs.nops()>0)
-		result *= sqrfree(quo(tmp, result, x), newargs);
+	if (args.nops()>0)
+		result *= sqrfree(quo(tmp, result, x), args);
 	else
 		result *= quo(tmp, result, x);
 
-	// Put in the reational overall factor again and return
+	// Put in the rational overall factor again and return
 	return result *	lcm.inverse();
 }
 
@@ -1724,24 +1943,21 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 //clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl;
 
 	// Factorize denominator and compute cofactors
-	exvector yun = sqrfree_yun(denom, x);
-//clog << "yun factors: " << exprseq(yun) << endl;
-	size_t num_yun = yun.size();
-	exvector factor; factor.reserve(num_yun);
-	exvector cofac; cofac.reserve(num_yun);
-	for (size_t i=0; i<num_yun; i++) {
-		if (!yun[i].is_equal(_ex1)) {
-			for (size_t j=0; j<=i; j++) {
-				factor.push_back(pow(yun[i], j+1));
-				ex prod = _ex1;
-				for (size_t k=0; k<num_yun; k++) {
-					if (k == i)
-						prod *= pow(yun[k], i-j);
-					else
-						prod *= pow(yun[k], k+1);
-				}
-				cofac.push_back(prod.expand());
+	epvector yun = sqrfree_yun(denom, x);
+	size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
+	exvector factor, cofac;
+	for (size_t i=0; i<yun.size(); i++) {
+		numeric i_exponent = ex_to<numeric>(yun[i].coeff);
+		for (size_t j=0; j<i_exponent; j++) {
+			factor.push_back(pow(yun[i].rest, j+1));
+			ex prod = _ex1;
+			for (size_t k=0; k<yun.size(); k++) {
+				if (yun[k].coeff == i_exponent)
+					prod *= pow(yun[k].rest, i_exponent-1-j);
+				else
+					prod *= pow(yun[k].rest, yun[k].coeff);
 			}
+			cofac.push_back(prod.expand());
 		}
 	}
 	size_t num_factors = factor.size();
@@ -1794,16 +2010,18 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
  *  @see ex::normal */
 static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
 {
+	// Since the repl contains replaced expressions we should search for them
+	ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
 	// Expression already replaced? Then return the assigned symbol
-	exmap::const_iterator it = rev_lookup.find(e);
+	auto it = rev_lookup.find(e_replaced);
 	if (it != rev_lookup.end())
 		return it->second;
-	
+
 	// Otherwise create new symbol and add to list, taking care that the
 	// replacement expression doesn't itself contain symbols from repl,
 	// because subs() is not recursive
-	ex es = (new symbol)->setflag(status_flags::dynallocated);
-	ex e_replaced = e.subs(repl, subs_options::no_pattern);
+	ex es = dynallocate<symbol>();
 	repl.insert(std::make_pair(es, e_replaced));
 	rev_lookup.insert(std::make_pair(e_replaced, es));
 	return es;
@@ -1816,16 +2034,18 @@ static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
  *  @see basic::to_polynomial */
 static ex replace_with_symbol(const ex & e, exmap & repl)
 {
+	// Since the repl contains replaced expressions we should search for them
+	ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
 	// 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;
-	
+	for (auto & it : repl)
+		if (it.second.is_equal(e_replaced))
+			return it.first;
+
 	// Otherwise create new symbol and add to list, taking care that the
 	// replacement expression doesn't itself contain symbols from repl,
 	// because subs() is not recursive
-	ex es = (new symbol)->setflag(status_flags::dynallocated);
-	ex e_replaced = e.subs(repl, subs_options::no_pattern);
+	ex es = dynallocate<symbol>();
 	repl.insert(std::make_pair(es, e_replaced));
 	return es;
 }
@@ -1833,36 +2053,27 @@ static ex replace_with_symbol(const ex & e, exmap & repl)
 
 /** Function object to be applied by basic::normal(). */
 struct normal_map_function : public map_function {
-	int level;
-	normal_map_function(int l) : level(l) {}
-	ex operator()(const ex & e) { return normal(e, level); }
+	ex operator()(const ex & e) override { return normal(e); }
 };
 
 /** Default implementation of ex::normal(). It normalizes the children and
  *  replaces the object with a temporary symbol.
  *  @see ex::normal */
-ex basic::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex basic::normal(exmap & repl, exmap & rev_lookup) const
 {
 	if (nops() == 0)
-		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, 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), repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
-		}
-	}
+		return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup), _ex1});
+
+	normal_map_function map_normal;
+	return dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_lookup), _ex1});
 }
 
 
 /** Implementation of ex::normal() for symbols. This returns the unmodified symbol.
  *  @see ex::normal */
-ex symbol::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex symbol::normal(exmap & repl, exmap & rev_lookup) const
 {
-	return (new lst(*this, _ex1))->setflag(status_flags::dynallocated);
+	return dynallocate<lst>({*this, _ex1});
 }
 
 
@@ -1870,7 +2081,7 @@ ex symbol::normal(exmap & repl, exmap & rev_lookup, 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(exmap & repl, exmap & rev_lookup, int level) const
+ex numeric::normal(exmap & repl, exmap & rev_lookup) const
 {
 	numeric num = numer();
 	ex numex = num;
@@ -1886,7 +2097,7 @@ ex numeric::normal(exmap & repl, exmap & rev_lookup, int level) const
 	}
 
 	// Denominator is always a real integer (see numeric::denom())
-	return (new lst(numex, denom()))->setflag(status_flags::dynallocated);
+	return dynallocate<lst>({numex, denom()});
 }
 
 
@@ -1904,11 +2115,11 @@ static ex frac_cancel(const ex &n, const ex &d)
 
 	// Handle trivial case where denominator is 1
 	if (den.is_equal(_ex1))
-		return (new lst(num, den))->setflag(status_flags::dynallocated);
+		return dynallocate<lst>({num, den});
 
 	// Handle special cases where numerator or denominator is 0
 	if (num.is_zero())
-		return (new lst(num, _ex1))->setflag(status_flags::dynallocated);
+		return dynallocate<lst>({num, _ex1});
 	if (den.expand().is_zero())
 		throw(std::overflow_error("frac_cancel: division by zero in frac_cancel"));
 
@@ -1947,32 +2158,25 @@ static ex frac_cancel(const ex &n, const ex &d)
 
 	// Return result as list
 //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl;
-	return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated);
+	return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom()});
 }
 
 
 /** Implementation of ex::normal() for a sum. It expands terms and performs
  *  fractional addition.
  *  @see ex::normal */
-ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex add::normal(exmap & repl, exmap & rev_lookup) const
 {
-	if (level == 1)
-		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 children and split each one into numerator and denominator
 	exvector nums, dens;
 	nums.reserve(seq.size()+1);
 	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(repl, rev_lookup, level-1);
+	for (auto & it : seq) {
+		ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup);
 		nums.push_back(n.op(0));
 		dens.push_back(n.op(1));
-		it++;
 	}
-	ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, level-1);
+	ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup);
 	nums.push_back(n.op(0));
 	dens.push_back(n.op(1));
 	GINAC_ASSERT(nums.size() == dens.size());
@@ -1982,8 +2186,8 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
 //std::clog << "add::normal uses " << nums.size() << " summands:\n";
 
 	// Add fractions sequentially
-	exvector::const_iterator num_it = nums.begin(), num_itend = nums.end();
-	exvector::const_iterator den_it = dens.begin(), den_itend = dens.end();
+	auto num_it = nums.begin(), num_itend = nums.end();
+	auto den_it = dens.begin(), den_itend = dens.end();
 //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
 	ex num = *num_it++, den = *den_it++;
 	while (num_it != num_itend) {
@@ -1996,7 +2200,7 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
 			num_it++; den_it++;
 		}
 
-		// Additiion of two fractions, taking advantage of the fact that
+		// Addition of two fractions, taking advantage of the fact that
 		// the heuristic GCD algorithm computes the cofactors at no extra cost
 		ex co_den1, co_den2;
 		ex g = gcd(den, next_den, &co_den1, &co_den2, false);
@@ -2013,31 +2217,23 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
 /** Implementation of ex::normal() for a product. It cancels common factors
  *  from fractions.
  *  @see ex::normal() */
-ex mul::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex mul::normal(exmap & repl, exmap & rev_lookup) const
 {
-	if (level == 1)
-		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 children, separate into numerator and denominator
 	exvector num; num.reserve(seq.size());
 	exvector den; den.reserve(seq.size());
 	ex n;
-	epvector::const_iterator it = seq.begin(), itend = seq.end();
-	while (it != itend) {
-		n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, rev_lookup, level-1);
+	for (auto & it : seq) {
+		n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup);
 		num.push_back(n.op(0));
 		den.push_back(n.op(1));
-		it++;
 	}
-	n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, level-1);
+	n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup);
 	num.push_back(n.op(0));
 	den.push_back(n.op(1));
 
 	// Perform fraction cancellation
-	return frac_cancel((new mul(num))->setflag(status_flags::dynallocated),
-	                   (new mul(den))->setflag(status_flags::dynallocated));
+	return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
 }
 
 
@@ -2045,16 +2241,11 @@ ex mul::normal(exmap & repl, exmap & rev_lookup, int level) const
  *  distributes integer exponents to numerator and denominator, and replaces
  *  non-integer powers by temporary symbols.
  *  @see ex::normal */
-ex power::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex power::normal(exmap & repl, exmap & rev_lookup) const
 {
-	if (level == 1)
-		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(repl, rev_lookup, level-1);
-	ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, level-1);
+	ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup);
+	ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup);
 	n_exponent = n_exponent.op(0) / n_exponent.op(1);
 
 	if (n_exponent.info(info_flags::integer)) {
@@ -2062,12 +2253,12 @@ ex power::normal(exmap & repl, exmap & rev_lookup, int level) const
 		if (n_exponent.info(info_flags::positive)) {
 
 			// (a/b)^n -> {a^n, b^n}
-			return (new lst(power(n_basis.op(0), n_exponent), power(n_basis.op(1), n_exponent)))->setflag(status_flags::dynallocated);
+			return dynallocate<lst>({pow(n_basis.op(0), n_exponent), pow(n_basis.op(1), n_exponent)});
 
 		} else if (n_exponent.info(info_flags::negative)) {
 
 			// (a/b)^-n -> {b^n, a^n}
-			return (new lst(power(n_basis.op(1), -n_exponent), power(n_basis.op(0), -n_exponent)))->setflag(status_flags::dynallocated);
+			return dynallocate<lst>({pow(n_basis.op(1), -n_exponent), pow(n_basis.op(0), -n_exponent)});
 		}
 
 	} else {
@@ -2075,43 +2266,41 @@ ex power::normal(exmap & repl, exmap & rev_lookup, 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), repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
+			return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1});
 
 		} 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), repl, rev_lookup)))->setflag(status_flags::dynallocated);
+				return dynallocate<lst>({_ex1, replace_with_symbol(pow(n_basis.op(0), -n_exponent), repl, rev_lookup)});
 
 			} 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), repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
+				return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup), _ex1});
 			}
 		}
 	}
 
 	// (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);
+	return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1});
 }
 
 
 /** Implementation of ex::normal() for pseries. It normalizes each coefficient
  *  and replaces the series by a temporary symbol.
  *  @see ex::normal */
-ex pseries::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex pseries::normal(exmap & repl, exmap & rev_lookup) const
 {
 	epvector newseq;
-	epvector::const_iterator i = seq.begin(), end = seq.end();
-	while (i != end) {
-		ex restexp = i->rest.normal();
+	for (auto & it : seq) {
+		ex restexp = it.rest.normal();
 		if (!restexp.is_zero())
-			newseq.push_back(expair(restexp, i->coeff));
-		++i;
+			newseq.push_back(expair(restexp, it.coeff));
 	}
-	ex n = pseries(relational(var,point), newseq);
-	return (new lst(replace_with_symbol(n, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated);
+	ex n = pseries(relational(var,point), std::move(newseq));
+	return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup), _ex1});
 }
 
 
@@ -2125,13 +2314,12 @@ ex pseries::normal(exmap & repl, exmap & rev_lookup, int level) const
  *  expression can be treated as a rational function). normal() is applied
  *  recursively to arguments of functions etc.
  *
- *  @param level maximum depth of recursion
  *  @return normalized expression */
-ex ex::normal(int level) const
+ex ex::normal() const
 {
 	exmap repl, rev_lookup;
 
-	ex e = bp->normal(repl, rev_lookup, level);
+	ex e = bp->normal(repl, rev_lookup);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
@@ -2152,7 +2340,7 @@ ex ex::numer() const
 {
 	exmap repl, rev_lookup;
 
-	ex e = bp->normal(repl, rev_lookup, 0);
+	ex e = bp->normal(repl, rev_lookup);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
@@ -2172,7 +2360,7 @@ ex ex::denom() const
 {
 	exmap repl, rev_lookup;
 
-	ex e = bp->normal(repl, rev_lookup, 0);
+	ex e = bp->normal(repl, rev_lookup);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
@@ -2182,7 +2370,7 @@ ex ex::denom() const
 		return e.op(1).subs(repl, subs_options::no_pattern);
 }
 
-/** Get numerator and denominator of an expression. If the expresison is not
+/** Get numerator and denominator of an expression. If the expression is not
  *  of the normal form "numerator/denominator", it is first converted to this
  *  form and then a list [numerator, denominator] is returned.
  *
@@ -2192,7 +2380,7 @@ ex ex::numer_denom() const
 {
 	exmap repl, rev_lookup;
 
-	ex e = bp->normal(repl, rev_lookup, 0);
+	ex e = bp->normal(repl, rev_lookup);
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
@@ -2221,47 +2409,11 @@ 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
@@ -2332,7 +2484,7 @@ ex numeric::to_polynomial(exmap & repl) const
 ex power::to_rational(exmap & repl) const
 {
 	if (exponent.info(info_flags::integer))
-		return power(basis.to_rational(repl), exponent);
+		return pow(basis.to_rational(repl), exponent);
 	else
 		return replace_with_symbol(*this, repl);
 }
@@ -2342,7 +2494,18 @@ ex power::to_rational(exmap & repl) const
 ex power::to_polynomial(exmap & repl) const
 {
 	if (exponent.info(info_flags::posint))
-		return power(basis.to_rational(repl), exponent);
+		return pow(basis.to_rational(repl), exponent);
+	else if (exponent.info(info_flags::negint))
+	{
+		ex basis_pref = collect_common_factors(basis);
+		if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) {
+			// (A*B)^n will be automagically transformed to A^n*B^n
+			ex t = pow(basis_pref, exponent);
+			return t.to_polynomial(repl);
+		}
+		else
+			return pow(replace_with_symbol(pow(basis, _ex_1), repl), -exponent);
+	} 
 	else
 		return replace_with_symbol(*this, repl);
 }
@@ -2353,17 +2516,15 @@ 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;
-	}
+	for (auto & it : seq)
+		s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_rational(repl)));
+
 	ex oc = overall_coeff.to_rational(repl);
 	if (oc.info(info_flags::numeric))
-		return thisexpairseq(s, overall_coeff);
+		return thisexpairseq(std::move(s), overall_coeff);
 	else
-		s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1));
-	return thisexpairseq(s, default_overall_coeff());
+		s.push_back(expair(oc, _ex1));
+	return thisexpairseq(std::move(s), default_overall_coeff());
 }
 
 /** Implementation of ex::to_polynomial() for expairseqs. */
@@ -2371,17 +2532,15 @@ 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_polynomial(repl)));
-		++i;
-	}
+	for (auto & it : seq)
+		s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_polynomial(repl)));
+
 	ex oc = overall_coeff.to_polynomial(repl);
 	if (oc.info(info_flags::numeric))
-		return thisexpairseq(s, overall_coeff);
+		return thisexpairseq(std::move(s), overall_coeff);
 	else
-		s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1));
-	return thisexpairseq(s, default_overall_coeff());
+		s.push_back(expair(oc, _ex1));
+	return thisexpairseq(std::move(s), default_overall_coeff());
 }
 
 
@@ -2400,7 +2559,7 @@ static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
 		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)) {
+			if (is_exactly_a<add>(x) || is_exactly_a<mul>(x) || is_a<power>(x)) {
 				ex f = 1;
 				x = find_common_factor(x, f, repl);
 				x *= f;
@@ -2436,7 +2595,7 @@ static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
 							else
 								v.push_back(t.op(k));
 						}
-						t = (new mul(v))->setflag(status_flags::dynallocated);
+						t = dynallocate<mul>(v);
 						goto term_done;
 					}
 				}
@@ -2446,7 +2605,7 @@ static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
 			t = x;
 term_done:	;
 		}
-		return (new add(terms))->setflag(status_flags::dynallocated);
+		return dynallocate<add>(terms);
 
 	} else if (is_exactly_a<mul>(e)) {
 
@@ -2456,11 +2615,19 @@ term_done:	;
 		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);
+		return dynallocate<mul>(v);
 
 	} else if (is_exactly_a<power>(e)) {
-
-		return e.to_polynomial(repl);
+		const ex e_exp(e.op(1));
+		if (e_exp.info(info_flags::integer)) {
+			ex eb = e.op(0).to_polynomial(repl);
+			ex factor_local(_ex1);
+			ex pre_res = find_common_factor(eb, factor_local, repl);
+			factor *= pow(factor_local, e_exp);
+			return pow(pre_res, e_exp);
+			
+		} else
+			return e.to_polynomial(repl);
 
 	} else
 		return e;
@@ -2471,7 +2638,7 @@ term_done:	;
  *  '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)) {
+	if (is_exactly_a<add>(e) || is_exactly_a<mul>(e) || is_exactly_a<power>(e)) {
 
 		exmap repl;
 		ex factor = 1;