-/** Rationalization of non-rational functions.
- * This function converts a general expression to a rational polynomial
- * by replacing all non-rational subexpressions (like non-rational numbers,
- * non-integer powers or functions like sin(), cos() etc.) to temporary
- * symbols. This makes it possible to use functions like gcd() and divide()
- * on non-rational functions by applying to_rational() on the arguments,
- * calling the desired function and re-substituting the temporary symbols
- * in the result. To make the last step possible, all temporary symbols and
- * their associated expressions are collected in the list specified by the
- * repl_lst parameter in the form {symbol == expression}, ready to be passed
- * as an argument to ex::subs().
- *
- * @param repl_lst collects a list of all temporary symbols and their replacements
- * @return rationalized expression */
-ex ex::to_rational(lst &repl_lst) const
+
+/** 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, exmap & repl)
+{
+ if (is_exactly_a<add>(e)) {
+
+ size_t num = e.nops();
+ exvector terms; terms.reserve(num);
+ ex gc;
+
+ // Find the common GCD
+ 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) || is_a<power>(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 (size_t 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 (size_t j=0; j<t.nops(); j++) {
+ if (t.op(j).is_equal(gc)) {
+ exvector v; v.reserve(t.nops());
+ for (size_t 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)) {
+
+ size_t num = e.nops();
+ exvector v; v.reserve(num);
+
+ 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);
+
+ } else if (is_exactly_a<power>(e)) {
+ 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 *= power(factor_local, e_exp);
+ return power(pre_res, e_exp);
+
+ } else
+ return e.to_polynomial(repl);
+
+ } 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) || is_exactly_a<power>(e)) {
+
+ exmap repl;
+ ex factor = 1;
+ ex r = find_common_factor(e, factor, repl);
+ return factor.subs(repl, subs_options::no_pattern) * r.subs(repl, subs_options::no_pattern);
+
+ } else
+ return e;
+}
+
+
+/** Resultant of two expressions e1,e2 with respect to symbol s.
+ * Method: Compute determinant of Sylvester matrix of e1,e2,s. */
+ex resultant(const ex & e1, const ex & e2, const ex & s)