From: Richard Kreckel <kreckel@ginac.de>
Date: Wed, 31 Jan 2018 09:04:25 +0000 (+0100)
Subject: Improve gcd(a, b) where one argument is a power of a symbol.
X-Git-Tag: release_1-7-3~5
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=e7d79ac4ff7654908b7688bc1373624119682f5c

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 747783c6..436b3582 100644
--- 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
-	ex aex = a.expand(), bex = b.expand();
+	ex aex = a.expand();
 	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;
 	}
+	ex bex = b.expand();
 	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;
 	}
-	// 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());
 
@@ -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)
-			*ca = pow(p, a.op(1) - 1);
+			*ca = pow(p, exp_a - 1);
 		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);
 
-	// a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
 	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)