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

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index fa09e429f1c563c6ad77298597c88bd8a88aed40..1808bcb90bc9722d5f78121b9c1f613896355e51 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1344,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;
@@ -1361,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
@@ -1370,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)) {
@@ -1382,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