From 3dfe5d31528291d1be6b1fe1180622efdabfebde Mon Sep 17 00:00:00 2001
From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Thu, 21 Nov 2002 22:57:06 +0000
Subject: [PATCH] added collect_common_factors() (is this a good name?)

---
 ginac/normal.cpp | 90 ++++++++++++++++++++++++++++++++++++++++++++++++
 ginac/normal.h   |  3 ++
 2 files changed, 93 insertions(+)

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 5e6aa470..67091f50 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1912,6 +1912,96 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 }
 
 
+/** 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)
+{
+	if (is_a<add>(e)) {
+
+		unsigned num = e.nops();
+		exvector terms; terms.reserve(num);
+		lst repl;
+		ex gc;
+
+		// Find the common GCD
+		for (unsigned i=0; i<num; i++) {
+			ex x = e.op(i).to_rational(repl);
+			if (is_a<add>(x) || is_a<mul>(x)) {
+				ex f = 1;
+				x = find_common_factor(x, f);
+				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.subs(repl);
+
+		// Now divide all terms by the GCD
+		for (unsigned i=0; i<num; i++) {
+			ex x;
+
+			// Try to avoid divide() because it expands the polynomial
+			ex &t = terms[i];
+			if (is_a<mul>(t)) {
+				for (unsigned j=0; j<t.nops(); j++) {
+					if (t.op(j).is_equal(gc)) {
+						exvector v; v.reserve(t.nops());
+						for (unsigned 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).subs(repl);
+						goto term_done;
+					}
+				}
+			}
+
+			divide(t, gc, x);
+			t = x.subs(repl);
+term_done:	;
+		}
+		return (new add(terms))->setflag(status_flags::dynallocated);
+
+	} else if (is_a<mul>(e)) {
+
+		exvector v;
+		for (unsigned i=0; i<e.nops(); i++)
+			v.push_back(find_common_factor(e.op(i), factor));
+		return (new mul(v))->setflag(status_flags::dynallocated);
+
+	} 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_a<add>(e) || is_a<mul>(e)) {
+
+		ex factor = 1;
+		ex r = find_common_factor(e, factor);
+		return factor * r;
+
+	} else
+		return e;
+}
+
+
 /*
  *  Normal form of rational functions
  */
diff --git a/ginac/normal.h b/ginac/normal.h
index be67d8eb..9196f965 100644
--- a/ginac/normal.h
+++ b/ginac/normal.h
@@ -63,6 +63,9 @@ extern ex sqrfree(const ex &a, const lst &l = lst());
 // Square-free partial fraction decomposition of a rational function a(x)
 extern ex sqrfree_parfrac(const ex & a, const symbol & x);
 
+// Collect common factors in sums.
+extern ex collect_common_factors(const ex & e);
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_NORMAL_H__
-- 
2.44.0