mod_gcd: naive hack to chose a 'good' prime (to speed up gcd computation).

[ginac.git] / ginac / polynomial / mod_gcd.cpp

diff --git a/ginac/polynomial/mod_gcd.cpp b/ginac/polynomial/mod_gcd.cpp
index 8aa84808c1f64a5975086eaa743d78c93f4e4455..82d5f427f397abd6b8e6cce68ac84eb6c03a72e5 100644 (file)
--- a/ginac/polynomial/mod_gcd.cpp
+++ b/ginac/polynomial/mod_gcd.cpp
@@ -104,8 +104,18 @@ void mod_gcd(upoly& result, upoly A, upoly B)
        ring_t q(0);
        upoly H;
 
-       ring_t p(2);
+       int count = 0;
+       const ring_t p_threshold = ring_t(1) << (8*sizeof(void *));
+       ring_t p = isqrt(std::min(max_coeff(A), max_coeff(B)));
        while (true) {
+               if (count >= 8) {
+                       count = 0;
+                       if (p < p_threshold)
+                               p <<= 1;
+                       else
+                               p = p + (p >> 4);
+               } else 
+                       ++count;
                find_next_prime(p, g);
 
                // Map the polynomials onto Z/p[x]