pgcd(), chinrem_gcd(): use appropriate definition of the degree.

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.

diff --git a/check/pgcd_infinite_loop.cpp b/check/pgcd_infinite_loop.cpp
index 013b746afd39d23525ec0900a0f81efb7a7e5bb6..ae39f1d79a4e6d84e7c942b2f1611e1d87f68f8d 100644 (file)
--- 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 bf30b5ce1e16027488374a00a35956923dbb497f..24da934690237e42e26a97caa35b09080eabc672 100644 (file)
--- 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 fac3735e10333390ea42c91a5e229bfd4e5c11df..a33b02d78a7dfc2fdb9fc20818b9c30271c52fb2 100644 (file)
--- 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)