Effect: fixes (rare) GCD miscalculation. pgcd() treats polynomials Z_p[x_0, ..., x_n] as R[x_0, ..., x_{n - 1}], where R = Z_p[x_n]. Therefore one should use correct definition of the degree (i.e. compute it w.r.t. x_0, ..., x_{n-1}, but NOT w.r.t. x_n!). One should use appropriate definition of degree (that is, w.r.t. x_0, ..., x_n) in chinrem_gcd() too. Thanks to Aless Lasaruk for a bug report. --- check/pgcd_infinite_loop.cpp | 8 ++++++-- ginac/polynomial/mgcd.cpp | 12 ++++++++---- ginac/polynomial/pgcd.cpp | 10 ++++++---- 3 files changed, 20 insertions(+), 10 deletions(-) diff --git a/check/pgcd_infinite_loop.cpp b/check/pgcd_infinite_loop.cpp index 013b746..ae39f1d 100644 --- a/check/pgcd_infinite_loop.cpp +++ b/check/pgcd_infinite_loop.cpp @@ -26,12 +26,16 @@ static const string srep("\ int main(int argc, char** argv) { - cout << "Checking for more pgcd() bugs (infinite loop) ... " << flush; + cout << "Checking for more pgcd() bugs (infinite loop, miscalculation) ... " << flush; parser the_parser; ex e = the_parser(srep); const symbol x = ex_to<symbol>(the_parser.get_syms()["x"]); ex g = gcd(e, e.diff(x)); + ex should_be = the_parser(string("u^4*z^2")); + if (!(g-should_be).expand().is_zero()) { + cout << "GCD was miscomputed. " << flush; + return 1; + } cout << "not found. " << flush; - cout << g << endl << flush; return 0; } diff --git a/ginac/polynomial/mgcd.cpp b/ginac/polynomial/mgcd.cpp index bf30b5c..24da934 100644 --- a/ginac/polynomial/mgcd.cpp +++ b/ginac/polynomial/mgcd.cpp @@ -27,6 +27,7 @@ #include "primes_factory.h" #include "divide_in_z_p.h" #include "poly_cra.h" +#include <numeric> // std::accumulate #include <cln/integer.h> @@ -58,12 +59,15 @@ ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars) const cln::cl_I g_lc = cln::gcd(a_lc, b_lc); const ex& x(vars.back()); - int n = std::min(A.degree(x), B.degree(x)); + exp_vector_t n = std::min(degree_vector(A, vars), degree_vector(B, vars)); + const int nTot = std::accumulate(n.begin(), n.end(), 0); const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient()); const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient()); - const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)* + + const cln::cl_I lcoeff_limit = (cln::cl_I(1) << nTot)*cln::abs(g_lc)* std::min(A_max_coeff, B_max_coeff); + cln::cl_I q = 0; ex H = 0; @@ -83,8 +87,8 @@ ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars) const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars); const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p); Cp = (Cp*numeric(nlc)).expand().smod(pnum); - int cp_deg = Cp.degree(x); - if (cp_deg == 0) + exp_vector_t cp_deg = degree_vector(Cp, vars); + if (zerop(cp_deg)) return numeric(g_lc); if (zerop(q)) { H = Cp; diff --git a/ginac/polynomial/pgcd.cpp b/ginac/polynomial/pgcd.cpp index fac3735..a33b02d 100644 --- a/ginac/polynomial/pgcd.cpp +++ b/ginac/polynomial/pgcd.cpp @@ -80,7 +80,8 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p) const ex lc_gcd = euclid_gcd(AL, BL, mainvar, p); // The estimate of degree of the gcd of Ab and Bb - int gcd_deg = std::min(degree(Aprim, mainvar), degree(Bprim, mainvar)); + exp_vector_t gcd_deg = std::min(degree_vector(Aprim, restvars), + degree_vector(Bprim, restvars)); eval_point_finder::value_type b; eval_point_finder find_eval_point(p); @@ -105,8 +106,11 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p) const cln::cl_I correct_lc = smod(lcb_gcd*recip(Cblc, p), p); Cb = (Cb*numeric(correct_lc)).smod(pn); + const exp_vector_t img_gcd_deg = degree_vector(Cb, restvars); + // Test for relatively prime polynomials + if (zerop(img_gcd_deg)) + return cont_gcd; // Test for unlucky homomorphisms - const int img_gcd_deg = Cb.degree(restvars.back()); if (img_gcd_deg < gcd_deg) { // The degree decreased, previous homomorphisms were // bad, so we have to start it all over. @@ -121,8 +125,6 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p) // evaluation point is bad. Skip it. continue; } - if (img_gcd_deg == 0) - return cont_gcd; // Image has the same degree as the previous one // (or at least not higher than the limit) -- 1.7.1