From 4ffb3cbb3ab5f642e461bcbf8fb29743752c5d58 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Wed, 24 Jan 2018 22:37:24 +0100 Subject: [PATCH] Speed up special cases of square-free factorization. Square-free factorization of polynomials containing a factor which is a high power P of a symbol x did scale like O(P) in space and time. This patch introduces a shortcut in Yun's algorithm, such that the computation is only O(1) in space and time. This makes it possible to compute, say sqrfree(x^P + x^(P+1)) => (1+x)*x^P with P=123456789. Found this to be a bottleneck while debugging one of Vitaly Magerya's examples. --- ginac/normal.cpp | 78 +++++++++++++++++++++++++----------------------- 1 file changed, 40 insertions(+), 38 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 231ef7e8..9f8b7b4a 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1785,32 +1785,39 @@ ex lcm(const ex &a, const ex &b, bool check_args) * @param a multivariate polynomial over Z[X], treated here as univariate * polynomial in x (needs not be expanded). * @param x variable to factor in - * @return vector of factors sorted in ascending degree */ -static exvector sqrfree_yun(const ex &a, const symbol &x) + * @return vector of expairs (factor, exponent), sorted by exponents */ +static epvector sqrfree_yun(const ex &a, const symbol &x) { - exvector res; ex w = a; ex z = w.diff(x); ex g = gcd(w, z); if (g.is_zero()) { - return res; + return epvector{}; } if (g.is_equal(_ex1)) { - res.push_back(a); - return res; + return epvector{expair(a, _ex1)}; } - ex y; + epvector results; + ex exponent = _ex0; do { w = quo(w, g, x); if (w.is_zero()) { return res; } - y = quo(z, g, x); - z = y - w.diff(x); + z = quo(z, g, x) - w.diff(x); + exponent = exponent + 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)); + break; + } g = gcd(w, z); - res.push_back(g); + if (!g.is_equal(_ex1)) { + results.push_back(expair(g, exponent)); + } } while (!z.is_zero()); - return res; + return results; } @@ -1878,30 +1885,28 @@ ex sqrfree(const ex &a, const lst &l) const ex tmp = multiply_lcm(a,lcm); // find the factors - exvector factors = sqrfree_yun(tmp, x); + epvector factors = sqrfree_yun(tmp, x); - // construct the next list of symbols with the first element popped - lst newargs = args; - newargs.remove_first(); + // remove symbol x and proceed recursively with the remaining symbols + args.remove_first(); // recurse down the factors in remaining variables - if (newargs.nops()>0) { + if (args.nops()>0) { for (auto & it : factors) - it = sqrfree(it, newargs); + it.rest = sqrfree(it.rest, args); } // Done with recursion, now construct the final result ex result = _ex1; - int p = 1; for (auto & it : factors) - result *= pow(it, p++); + 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 (newargs.nops()>0) - result *= sqrfree(quo(tmp, result, x), newargs); + if (args.nops()>0) + result *= sqrfree(quo(tmp, result, x), args); else result *= quo(tmp, result, x); @@ -1929,24 +1934,21 @@ ex sqrfree_parfrac(const ex & a, const symbol & x) //clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl; // Factorize denominator and compute cofactors - exvector yun = sqrfree_yun(denom, x); -//clog << "yun factors: " << exprseq(yun) << endl; - size_t num_yun = yun.size(); - exvector factor; factor.reserve(num_yun); - exvector cofac; cofac.reserve(num_yun); - for (size_t i=0; i(yun.back().coeff).to_int(); + exvector factor, cofac; + for (size_t i=0; i(yun[i].coeff); + for (size_t j=0; j