[GiNaC-devel] [PATCH 1/3] divide: try to avoid expanding partially factored expressions.

[Alexei Sheplyakov]

Dear all,

this ginsh code snippet:

collect_common_factors((x*y*a+x*y*b)^(3) + (x*z + x*y)^(2))

gives

x^2*(x*(3*a^2*y^3*b+3*a*y^3*b^2+y^3*b^3+a^3*y^3)+2*z*y+z^2+y^2)

Thus instead of collecting common factors from the terms of sums
(as manual says) collect_common_factors [sometimes] expands term[s]!


---
 ginac/normal.cpp |   66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 files changed, 66 insertions(+), 0 deletions(-)

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 7b90030..49f0bd7 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -614,6 +614,72 @@ #endif
 	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 = (new mul(resv))->setflag(status_flags::dynallocated);
+			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*power(ab, a_exp - 1);
+			return true;
+		}
+		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()) {
-- 
All science is either physics or stamp collecting.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-devel/attachments/20061006/bb94913f/attachment.pgp