From e65cbada7281938475aa0cf31a3d56405ac91d4e Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Fri, 21 Aug 2020 17:47:19 +0200 Subject: [PATCH] Avoid unnecessary expansion in sqrfree_yun(). Don't make the polynomial primitive inside Yun's algorithm. Computing the primitive part will incur an extra expansion of the polynomial and this may distroy a pre-factored structure in the other variables. In the multi- variate case, working with primitive polynomials isn't sufficient anyways to avoid a final division to discover lost factors. --- ginac/normal.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 6ad61137..f8615499 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1797,8 +1797,7 @@ ex lcm(const ex &a, const ex &b, bool check_args) * @return vector of expairs (factor, exponent), sorted by exponents */ static epvector sqrfree_yun(const ex &a, const symbol &x) { - // make polynomial primitive w.r.t. x - ex w = a.primpart(x); + ex w = a; ex z = w.diff(x); ex g = gcd(w, z); if (g.is_zero()) { @@ -1814,6 +1813,10 @@ static epvector sqrfree_yun(const ex &a, const symbol &x) ex i = 0; // exponent do { w = quo(w, g, x); + if (w.is_zero()) { + // hidden zero + break; + } z = quo(z, g, x) - w.diff(x); i += 1; if (w.is_equal(x)) { @@ -1829,20 +1832,18 @@ static epvector sqrfree_yun(const ex &a, const symbol &x) } while (!z.is_zero()); // correct for lost factor - // (NB: This factor is not necessarily the unit and content lost when - // we made the polynomial primitive since Yun's algorithm only finds - // factors up to a unit.) + // (being based on GCDs, Yun's algorithm only finds factors up to a unit) const ex lost_factor = quo(a, mul{factors}, x); if (lost_factor.is_equal(_ex1)) { // trivial lost factor return factors; } - if (factors[0].coeff.is_equal(1)) { + if (!factors.empty() && factors[0].coeff.is_equal(1)) { // multiply factor^1 with lost_factor factors[0].rest *= lost_factor; return factors; } - // only factors with higher exponent: prepend lost_factor^1 to the results + // no factor^1: prepend lost_factor^1 to the results epvector results = {expair(lost_factor, 1)}; std::move(factors.begin(), factors.end(), std::back_inserter(results)); return results; -- 2.31.1