Improve gcd(a, b) where one argument is a power of a symbol.

The already implemented recursion
  gcd(x^n, x*p(x)) -> x*gcd(x^(n-1), p(x))
is not ambitious enough: If p(x) has a factor of x, it just goes through
the same step again, and if p(x) has a factor which is a power of x, this
is reapeted many times.

Example:
  gcd(x^n, expand(x^n*(1+x)))
used to go recurse through the gcd routine n times, which could
easily lead to a stack overflow for n about 10^5.

To improve the situation, introduce a special case for gcd where one
argument is a power of a symbol and just cancel the common factor.

This turned out to be the root cause of segfaults in matrix elimination
algorithms, as reported by Patrick Schulz and Vitaly Magerya:
Cf. <https://www.ginac.de/pipermail/ginac-list/2018-January/thread.html>

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 747783c69d39cac29d4bcd5c4bd39f75c4a99afd..436b358235e4e55a7eda1e62ab6e5b330e1f5527 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1470,7 +1470,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
        }
 
        // Some trivial cases
        }
 
        // Some trivial cases

-       ex aex = a.expand(), bex = b.expand();
+       ex aex = a.expand();

        if (aex.is_zero()) {
                if (ca)
                        *ca = _ex0;
        if (aex.is_zero()) {
                if (ca)
                        *ca = _ex0;


@@ -1478,6 +1478,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
                        *cb = _ex1;
                return b;
        }
                        *cb = _ex1;
                return b;
        }

+       ex bex = b.expand();

        if (bex.is_zero()) {
                if (ca)
                        *ca = _ex1;
        if (bex.is_zero()) {
                if (ca)
                        *ca = _ex1;


@@ -1563,7 +1564,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
                        *cb = b;
                return _ex1;
        }
                        *cb = b;
                return _ex1;
        }

-       // move symbols which contained only in one of the polynomials
+       // move symbols which are contained only in one of the polynomials

        // to the end:
        rotate(sym_stats.begin(), vari, sym_stats.end());
 
        // to the end:
        rotate(sym_stats.begin(), vari, sym_stats.end());
 


@@ -1706,17 +1707,27 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
        if (p.is_equal(b)) {
                // a = p^n, b = p, gcd = p
                if (ca)
        if (p.is_equal(b)) {
                // a = p^n, b = p, gcd = p
                if (ca)

-                       *ca = pow(p, a.op(1) - 1);
+                       *ca = pow(p, exp_a - 1);

                if (cb)
                        *cb = _ex1;
                return p;
                if (cb)
                        *cb = _ex1;
                return p;

-       } 
+       }
+       if (is_a<symbol>(p)) {
+               // Cancel trivial common factor
+               int ldeg_a = ex_to<numeric>(exp_a).to_int();
+               int ldeg_b = b.ldegree(p);
+               int min_ldeg = std::min(ldeg_a, ldeg_b);
+               if (min_ldeg > 0) {
+                       ex common = pow(p, min_ldeg);
+                       return gcd(pow(p, ldeg_a - min_ldeg), (b / common).expand(), ca, cb, false) * common;
+               }
+       }

 
        ex p_co, bpart_co;
        ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
 
 
        ex p_co, bpart_co;
        ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
 

-       // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1

        if (p_gcd.is_equal(_ex1)) {
        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)
                if (ca)
                        *ca = a;
                if (cb)