GCD avoids to produce expanded expressions (Alexei's patch).

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index cb7e35ef1db41c318bcb5d9fdd4339d6f4c8f1ee..1808bcb90bc9722d5f78121b9c1f613896355e51 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-2004 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2005 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
@@ -20,7 +20,7 @@
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
 #include <algorithm>
@@ -1247,11 +1247,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
                        ex dummy;
                        if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
                                g *= gc;
-                               ex lc = g.lcoeff(x);
-                               if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
-                                       return -g;
-                               else
-                                       return g;
+                               return g;
                        }
                }
 
@@ -1348,10 +1344,12 @@ factored_b:
        // Input polynomials of the form poly^n are sometimes also trivial
        if (is_exactly_a<power>(a)) {
                ex p = a.op(0);
+               const ex& exp_a = a.op(1);
                if (is_exactly_a<power>(b)) {
-                       if (p.is_equal(b.op(0))) {
+                       ex pb = b.op(0);
+                       const ex& exp_b = b.op(1);
+                       if (p.is_equal(pb)) {
                                // a = p^n, b = p^m, gcd = p^min(n, m)
-                               ex exp_a = a.op(1), exp_b = b.op(1);
                                if (exp_a < exp_b) {
                                        if (ca)
                                                *ca = _ex1;
@@ -1365,7 +1363,32 @@ factored_b:
                                                *cb = _ex1;
                                        return power(p, exp_b);
                                }
-                       }
+                       } else {
+                               ex p_co, pb_co;
+                               ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args);
+                               if (p_gcd.is_equal(_ex1)) {
+                                       // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==>
+                                       // gcd(a,b) = 1
+                                       if (ca)
+                                               *ca = a;
+                                       if (cb)
+                                               *cb = b;
+                                       return _ex1;
+                                       // XXX: do I need to check for p_gcd = -1?
+                               } else {
+                                       // there are common factors:
+                                       // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
+                                       // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
+                                       if (exp_a < exp_b) {
+                                               return power(p_gcd, exp_a)*
+                                                       gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
+                                       } else {
+                                               return power(p_gcd, exp_b)*
+                                                       gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
+                                       }
+                               } // p_gcd.is_equal(_ex1)
+                       } // p.is_equal(pb)
+
                } else {
                        if (p.is_equal(b)) {
                                // a = p^n, b = p, gcd = p
@@ -1374,8 +1397,24 @@ factored_b:
                                if (cb)
                                        *cb = _ex1;
                                return p;
+                       } 
+
+                       ex p_co, bpart_co;
+                       ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
+
+                       if (p_gcd.is_equal(_ex1)) {
+                               // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
+                               if (ca)
+                                       *ca = a;
+                               if (cb)
+                                       *cb = b;
+                               return _ex1;
+                       } else {
+                               // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+                               return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
                        }
-               }
+               } // is_exactly_a<power>(b)
+
        } else if (is_exactly_a<power>(b)) {
                ex p = b.op(0);
                if (p.is_equal(a)) {
@@ -1386,6 +1425,23 @@ factored_b:
                                *cb = power(p, b.op(1) - 1);
                        return p;
                }
+
+               ex p_co, apart_co;
+               const ex& exp_b(b.op(1));
+               ex p_gcd = gcd(a, p, &apart_co, &p_co, false);
+               if (p_gcd.is_equal(_ex1)) {
+                       // b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1
+                       if (ca)
+                               *ca = a;
+                       if (cb)
+                               *cb = b;
+                       return _ex1;
+               } else {
+                       // there are common factors:
+                       // a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+
+                       return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false);
+               } // p_gcd.is_equal(_ex1)
        }
 #endif