+/** 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, lst & repl)
+{
+ if (is_exactly_a<add>(e)) {
+
+ unsigned num = e.nops();
+ exvector terms; terms.reserve(num);
+ ex gc;
+
+ // Find the common GCD
+ for (unsigned i=0; i<num; i++) {
+ ex x = e.op(i).to_rational(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 (unsigned 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 (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);
+ 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)) {
+
+ unsigned num = e.nops();
+ exvector v; v.reserve(num);
+
+ for (unsigned 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)) {
+
+ ex x = e.to_rational(repl);
+ if (is_exactly_a<power>(x) && x.op(1).info(info_flags::negative))
+ return replace_with_symbol(x, repl);
+ else
+ return x;
+
+ } 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)) {
+
+ lst repl;
+ ex factor = 1;
+ ex r = find_common_factor(e, factor, repl);
+ return factor.subs(repl) * r.subs(repl);
+
+ } else
+ return e;
+}
+
+