Avoid unnecessary expansion in sqrfree_yun().

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 60c09d9cc4669c0d4fa859b5c69f09c163fbf64c..f8615499add62a1a8e447a7b8d5f2aaf725b62b0 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1801,31 +1801,51 @@ static epvector sqrfree_yun(const ex &a, const symbol &x)
        ex z = w.diff(x);
        ex g = gcd(w, z);
        if (g.is_zero()) {
-               return epvector{};
+               // manifest zero or hidden zero
+               return {};
        }
        if (g.is_equal(_ex1)) {
-               return epvector{expair(a, _ex1)};
+               // w(x) and w'(x) share no factors: w(x) is square-free
+               return {expair(a, _ex1)};
        }
-       epvector results;
-       ex exponent = _ex0;
+
+       epvector factors;
+       ex i = 0;  // exponent
        do {
                w = quo(w, g, x);
                if (w.is_zero()) {
-                       return results;
+                       // hidden zero
+                       break;
                }
                z = quo(z, g, x) - w.diff(x);
-               exponent = exponent + 1;
+               i += 1;
                if (w.is_equal(x)) {
                        // shortcut for x^n with n ∈ ℕ
-                       exponent += quo(z, w.diff(x), x);
-                       results.push_back(expair(w, exponent));
+                       i += quo(z, w.diff(x), x);
+                       factors.push_back(expair(w, i));
                        break;
                }
                g = gcd(w, z);
                if (!g.is_equal(_ex1)) {
-                       results.push_back(expair(g, exponent));
+                       factors.push_back(expair(g, i));
                }
        } while (!z.is_zero());
+
+       // correct for lost factor
+       // (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.empty() && factors[0].coeff.is_equal(1)) {
+               // multiply factor^1 with lost_factor
+               factors[0].rest *= lost_factor;
+               return factors;
+       }
+       // 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;
 }
 
@@ -1891,10 +1911,14 @@ ex sqrfree(const ex &a, const lst &l)
 
        // convert the argument from something in Q[X] to something in Z[X]
        const numeric lcm = lcm_of_coefficients_denominators(a);
-       const ex tmp = multiply_lcm(a,lcm);
+       const ex tmp = multiply_lcm(a, lcm);
 
        // find the factors
        epvector factors = sqrfree_yun(tmp, x);
+       if (factors.empty()) {
+               // the polynomial was a hidden zero
+               return _ex0;
+       }
 
        // remove symbol x and proceed recursively with the remaining symbols
        args.remove_first();
@@ -1906,21 +1930,10 @@ ex sqrfree(const ex &a, const lst &l)
        }
 
        // Done with recursion, now construct the final result
-       ex result = _ex1;
-       for (auto & it : factors)
-               result *= pow(it.rest, it.coeff);
-
-       // Yun's algorithm does not account for constant factors.  (For univariate
-       // polynomials it works only in the monic case.)  We can correct this by
-       // inserting what has been lost back into the result.  For completeness
-       // we'll also have to recurse down that factor in the remaining variables.
-       if (args.nops()>0)
-               result *= sqrfree(quo(tmp, result, x), args);
-       else
-               result *= quo(tmp, result, x);
+       ex result = mul(factors);
 
        // Put in the rational overall factor again and return
-       return result * lcm.inverse();
+       return result * lcm.inverse();
 }
 
 
@@ -1944,10 +1957,6 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 
        // Factorize denominator and compute cofactors
        epvector yun = sqrfree_yun(denom, x);
-       // Yun's algorithm does not account for constant factors.  (For univariate
-       // polynomials it works only in the monic case.)  We can correct this by
-       // inserting what has been lost back into the result.
-       ex lead_coeff  = quo(denom, mul(yun), x);
        size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
        exvector factor, cofac;
        for (size_t i=0; i<yun.size(); i++) {
@@ -1989,7 +1998,7 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
        // Sum up decomposed fractions
        ex sum = 0;
        for (size_t i=0; i<num_factors; i++)
-               sum += sol(i, 0) / factor[i] / lead_coeff;
+               sum += sol(i, 0) / factor[i];
 
        return red_poly + sum;
 }