- collect_common_factors() works better with negative powers

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index d88f3dde4925ead551a23a0362749fe68a4b7037..e1d7349fa509b600cf8321c54a5ca1ed0a93180d 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -6,7 +6,7 @@
  *  computation, square-free factorization and rational function normalization. */
 
 /*
- *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -471,14 +471,14 @@ ex decomp_rational(const ex &a, const symbol &x)
 }
 
 
-/** Pseudo-remainder of polynomials a(x) and b(x) in Z[x].
+/** Pseudo-remainder of polynomials a(x) and b(x) in Q[x].
  *
  *  @param a  first polynomial in x (dividend)
  *  @param b  second polynomial in x (divisor)
  *  @param x  a and b are polynomials in x
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
- *  @return pseudo-remainder of a(x) and b(x) in Z[x] */
+ *  @return pseudo-remainder of a(x) and b(x) in Q[x] */
 ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
        if (b.is_zero())
@@ -523,14 +523,14 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
 }
 
 
-/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Z[x].
+/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Q[x].
  *
  *  @param a  first polynomial in x (dividend)
  *  @param b  second polynomial in x (divisor)
  *  @param x  a and b are polynomials in x
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
- *  @return sparse pseudo-remainder of a(x) and b(x) in Z[x] */
+ *  @return sparse pseudo-remainder of a(x) and b(x) in Q[x] */
 ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args)
 {
        if (b.is_zero())
@@ -1741,6 +1741,7 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
        return res;
 }
 
+
 /** Compute a square-free factorization of a multivariate polynomial in Q[X].
  *
  *  @param a  multivariate polynomial over Q[X]
@@ -1842,6 +1843,7 @@ ex sqrfree(const ex &a, const lst &l)
        return result * lcm.inverse();
 }
 
+
 /** Compute square-free partial fraction decomposition of rational function
  *  a(x).
  *
@@ -2413,4 +2415,106 @@ ex expairseq::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, lst & repl)
+{
+       if (is_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_a<add>(x) || is_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_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_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_a<power>(e)) {
+
+               ex x = e.to_rational(repl);
+               if (is_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_a<add>(e) || is_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;
+}
+
+
 } // namespace GiNaC