From 1f14fb79bacd8de080c73a930bf9f8a5d514cbe0 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Mon, 19 Dec 2022 18:14:32 +0100 Subject: [PATCH] [BUGFIX] Fix corner cases of factor_multivariate(). In some rare cases, factor() failed with what looks like an endless loop trying evaluation homomorphisms with ever increasing moduli (or maybe it just took eons to complete). Whether this happened depended on which symbol was left unevaluated for univariate factorization. As a workaround for this problem, we now apply the evaluation homomorphism for each symbol in a round-robin mode until one of them succeeds. This somewhat speeds up factorization of multivariate polynomials even in those cases where the old implementation terminated in time. Thanks to Vitaly Magerya for reporting this. --- check/exam_factor.cpp | 37 +++++++++ ginac/factor.cpp | 172 ++++++++++++++++++++++++------------------ 2 files changed, 137 insertions(+), 72 deletions(-) diff --git a/check/exam_factor.cpp b/check/exam_factor.cpp index d7ba1174..55213579 100644 --- a/check/exam_factor.cpp +++ b/check/exam_factor.cpp @@ -280,6 +280,42 @@ static unsigned exam_factor_wang() return result; } +static unsigned exam_factor_magerya() +{ + // In 2017, Vitaly Magerya reported a class of biviariate polynomials + // where Hensel lifting sometimes failed to terminate. + // https://www.ginac.de/pipermail/ginac-list/2017-December/002162.html + unsigned result = 0; + ex e; + symbol x("x"), y("y"); + + e = (1+2*x+y)*(1+2*x-y)*(2*x-y)*(2*x+y); + result += check_factor_expanded(e); + + e = (7+4*x+y)*(-5+2*x-y)*(-6+6*x+y)*y*(10+2*x+y); + result += check_factor_expanded(e); + + e = (8+6*x-y)*(-5+4*x-y)*(-5+6*x+y)*(-2+2*x-y)*(2+4*x+y); + result += check_factor_expanded(e); + + e = -(-4+4*x+5*y)*(1+4*x+5*y)*(2+3*y)*(1+2*x-y)*(4+2*x+y); + result += check_factor_expanded(e); + + e = (-3+y-2*x)*(4+3*y-4*x)*(3+3*y+2*x)*(-2+3*y+2*x)*(-1+4*y+3*x); + result += check_factor_expanded(e); + + e = (-9+7*x+y)*(-5+6*x+y)*(4+2*x+y)*(5+2*x-y)*(7+9*x-y)*(8+6*x-y); + result += check_factor_expanded(e); + + e = pow(2*x-y,2)*(-1+6*x-y)*(-1+3*x-y)*(-2+4*x-y)*(1+4*x-y)*(4*x-y)*(2+4*x-y); + result += check_factor_expanded(e); + + e = (5+2*y-3*x)*(-4+4*y+3*x)*(-3+4*y-2*x)*(4+5*y-x)*(3*y+2*x)*(-1+3*y+5*x)*(5+3*y+4*x); + result += check_factor_expanded(e); + + return result; +} + unsigned exam_factor() { unsigned result = 0; @@ -291,6 +327,7 @@ unsigned exam_factor() result += exam_factor3(); cout << '.' << flush; result += exam_factor_content(); cout << '.' << flush; result += exam_factor_wang(); cout << '.' << flush; + result += exam_factor_magerya(); cout << '.' << flush; return result; } diff --git a/ginac/factor.cpp b/ginac/factor.cpp index bbefcf3f..f38757c3 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -144,7 +144,8 @@ typedef std::vector umodpoly; typedef std::vector upoly; typedef vector upvec; -// COPY FROM UPOLY.HPP + +// COPY FROM UPOLY.H // CHANGED size_t -> int !!! template static int degree(const T& p) @@ -206,7 +207,7 @@ canonicalize(T& p, const typename T::size_type hint = std::numeric_limits& basis) /** Berlekamp's modular factorization. * - * The implementation follows the algorithm in chapter 8 of [GCL]. + * The implementation follows algorithm 8.4 of [GCL]. * * @param[in] a modular polynomial * @param[out] upv vector containing modular factors. if upv was not empty the @@ -1027,7 +1030,7 @@ static void modsqrfree(const umodpoly& a, upvec& factors, vector& mult) /** Distinct degree factorization (DDF). * - * The implementation follows the algorithm in chapter 8 of [GCL]. + * The implementation follows algorithm 8.8 of [GCL]. * * @param[in] a_ modular polynomial * @param[out] degrees vector containing the degrees of the factors of the @@ -1215,7 +1218,7 @@ static inline cl_I calc_bound(const upoly& a) /** Hensel lifting as used by factor_univariate(). * - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.1 of [GCL]. * * @param[in] a_ primitive univariate polynomials * @param[in] p prime number that does not divide lcoeff(a) @@ -1648,7 +1651,7 @@ static vector multivar_diophant(const vector& a_, const ex& x, const ex& * with deg(s_i) < deg(a_i) * and with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r * - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.3 of [GCL]. * * @param[in] a vector of modular univariate polynomials * @param[in] x symbol @@ -1703,7 +1706,7 @@ static void change_modulus(const cl_modint_ring& R, umodpoly& a) * * Solves s*a + t*b == 1 mod p^k given a,b. * - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.3 of [GCL]. * * @param[in] a polynomial * @param[in] b polynomial @@ -1762,7 +1765,7 @@ static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned * s_1*b_1 + ... + s_r*b_r == x^m mod p^k * with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r * - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.3 of [GCL]. * * @param a vector with univariate polynomials mod p^k * @param x symbol @@ -1846,7 +1849,7 @@ static ex make_modular(const ex& e, const cl_modint_ring& R) * s_1*b_1 + ... + s_r*b_r == c mod * with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r * - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.2 of [GCL]. * * @param a_ vector of multivariate factors mod p^k * @param x symbol (equiv. x_1 in [GCL]) @@ -1956,7 +1959,7 @@ static vector multivar_diophant(const vector& a_, const ex& x, const ex& } /** Multivariate Hensel lifting. - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.4 of [GCL]. * Since we don't have a data type for modular multivariate polynomials, the * respective operations are done in a GiNaC::ex and the function * make_modular() is then called to make the coefficient modular p^l. @@ -2142,25 +2145,24 @@ static bool checkdivisors(const lst& f) * * @param[in] u multivariate polynomial to be factored * @param[in] vn leading coefficient of u in x (x==first symbol in syms) - * @param[in] syms set of symbols that appear in u + * @param[in] x first symbol that appears in u + * @param[in] syms_wox remaining symbols that appear in u * @param[in] f lst containing the factors of the leading coefficient vn * @param[in,out] modulus integer modulus for random number generation (i.e. |a_i| < modulus) * @param[out] u0 returns the evaluated (univariate) polynomial * @param[out] a returns the valid evaluation points. must have initial size equal * number of symbols-1 before calling generate_set */ -static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst& f, +static void generate_set(const ex& u, const ex& vn, const ex& x, const exset& syms_wox, const lst& f, numeric& modulus, ex& u0, vector& a) { - const ex& x = *syms.begin(); while ( true ) { ++modulus; // generate a set of integers ... u0 = u; ex vna = vn; ex vnatry; - exset::const_iterator s = syms.begin(); - ++s; + auto s = syms_wox.begin(); for ( size_t i=0; i(fnum.op(i)) ) { - s = syms.begin(); - ++s; + s = syms_wox.begin(); for ( size_t j=0; j(cont) ) { - return unit * factor_sqrfree(cont) * factor_sqrfree(pp); - } - - // factor leading coefficient - ex vn = pp.collect(x).lcoeff(x); - ex vnlst; - if ( is_a(vn) ) { - vnlst = lst{vn}; - } - else { - ex vnfactors = factor(vn); - vnlst = put_factors_into_lst(vnfactors); - } - - const unsigned int maxtrials = 3; - numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3; - vector a(syms.size()-1, 0); - - // try now to factorize until we are successful - while ( true ) { +struct factorization_ctx { + const ex poly, x; // polynomial, first symbol x... + const exset syms_wox; // ...remaining symbols w/o x + ex unit, cont, pp; // unit * cont * pp == poly + ex vn, vnlst; // leading coeff, factors of leading coeff + numeric modulus; // incremented each time we try + /** returns factors or empty if it did not succeed */ + ex try_next_evaluation_homomorphism() + { + constexpr unsigned maxtrials = 3; + vector a(syms_wox.size(), 0); unsigned int trialcount = 0; unsigned int prime; @@ -2256,7 +2227,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) while ( trialcount < maxtrials ) { // generate a set of valid evaluation points - generate_set(pp, vn, syms, ex_to(vnlst), modulus, u, a); + generate_set(pp, vn, x, syms_wox, ex_to(vnlst), modulus, u, a); ufac = factor_univariate(u, x, prime); ufaclst = put_factors_into_lst(ufac); @@ -2265,7 +2236,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) if ( factor_count <= 1 ) { // irreducible - return poly; + return lst{pp}; } if ( min_factor_count < 0 ) { // first time here @@ -2297,8 +2268,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) vector ftilde(vnlst.nops()-1); for ( size_t i=0; i epv; - auto s = syms.begin(); - ++s; + auto s = syms_wox.begin(); for ( size_t i=0; i ctx_in_x; + for (auto x : syms) { + exset syms_wox; // remaining syms w/o x + copy_if(syms.begin(), syms.end(), + inserter(syms_wox, syms_wox.end()), [x](const ex& y){ return y != x; }); + + factorization_ctx ctx = {.poly = poly, .x = x, + .syms_wox = syms_wox}; + + // make polynomial primitive + poly.unitcontprim(x, ctx.unit, ctx.cont, ctx.pp); + if ( !is_a(ctx.cont) ) { + // content is a polynomial in one or more of remaining syms, let's start over + return ctx.unit * factor_sqrfree(ctx.cont) * factor_sqrfree(ctx.pp); + } + + // find factors of leading coefficient + ctx.vn = ctx.pp.collect(x).lcoeff(x); + ctx.vnlst = put_factors_into_lst(factor(ctx.vn)); + + ctx.modulus = (ctx.vnlst.nops() > 3) ? ctx.vnlst.nops() : 3; + + ctx_in_x.push_back(ctx); + } + + // try an evaluation homomorphism for each context in round-robin + auto ctx = ctx_in_x.begin(); + while ( true ) { + + ex res = ctx->try_next_evaluation_homomorphism(); + if ( res != lst{} ) { - ex result = cont * unit; + // found the factors + ex result = ctx->cont * ctx->unit; for ( size_t i=0; ix, unit, cont, pp); + result *= unit * cont * pp; } return result; } + + // switch context for next symbol + if (++ctx == ctx_in_x.end()) { + ctx = ctx_in_x.begin(); + } } } -- 2.49.0