From e7d79ac4ff7654908b7688bc1373624119682f5c Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Wed, 31 Jan 2018 10:04:25 +0100 Subject: [PATCH] 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. --- ginac/normal.cpp | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) 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(p)) { + // Cancel trivial common factor + int ldeg_a = ex_to(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) -- 2.25.4