From 54c1eca199508d821795ca08506f6e7788ebee91 Mon Sep 17 00:00:00 2001
From: Jens Vollinga <vollinga@thep.physik.uni-mainz.de>
Date: Thu, 30 Jun 2005 17:58:53 +0000
Subject: [PATCH] collect_common_factors now works with powers [A.Sheplyakov]

---
 ginac/normal.cpp | 59 +++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 56 insertions(+), 3 deletions(-)

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 64c34ac2..db1ca467 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -714,6 +714,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_iterator itrb = b.begin(); itrb != b.end(); ++itrb) {
+			sym_desc_vec sym_stats;
+			get_symbol_stats(a, *itrb, sym_stats);
+			if (!divide_in_z(qbar, *itrb, q, sym_stats.begin()))
+				return false;
+
+			qbar = q;
+		}
+		return true;
+	}
+
 	// Main symbol
 	const ex &x = var->sym;
 
@@ -1478,6 +1503,26 @@ 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;
+		}
+	}
+
 	// Gather symbol statistics
 	sym_desc_vec sym_stats;
 	get_symbol_stats(a, b, sym_stats);
@@ -2459,8 +2504,16 @@ term_done:	;
 		return (new mul(v))->setflag(status_flags::dynallocated);
 
 	} 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::posint)) {
+			ex eb = e.op(0).to_polynomial(repl);
+			ex factor_local(_ex1);
+			ex pre_res = find_common_factor(eb, factor_local, repl);
+			factor *= power(factor_local, e_exp);
+			return power(pre_res, e_exp);
+			
+		} else
+			return e.to_polynomial(repl);
 
 	} else
 		return e;
@@ -2471,7 +2524,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;
-- 
2.44.0