Fix integer overflows in the heuristic univariate gcd

debian/libpoly0.symbols

@@ -1,5 +1,4 @@
 libpoly.so.0 libpoly0 #MINVER#
- bound_valuation@Base 1.3
  coefficient_add@Base 1.3
  coefficient_add_monomial@Base 1.3
  coefficient_add_mul@Base 1.3

src/upolynomial/gcd.c

@@ -26,6 +26,7 @@
 #include "utils/statistics.h"
 
 #include <assert.h>
+#include <limits.h>
 
 STAT_DECLARE(int, upolynomial, gcd_euclid)
 STAT_DECLARE(int, upolynomial, gcd_euclid_extended)
@@ -397,8 +398,10 @@ static lp_upolynomial_t* reconstruct_polynomial(size_t max_size, lp_integer_t* p
   integer_pow(lp_Z, &P, &P, pow);
 
   // Basically compute digits modulo x_value
-  int d = 0;
-  while (integer_sgn(lp_Z, p_value)) {
+  size_t d = 0;
+  // With 2^pow from bound_valuation the value has at most max_size digits.
+  // More digits mean no gcd of this degree
+  while (integer_sgn(lp_Z, p_value) && d < max_size) {
     integer_div_rem_pow2_Z(&div, &rem, p_value, pow);
 
     if (integer_bits(&rem) + 1 >= pow) {
@@ -412,16 +415,22 @@ static lp_upolynomial_t* reconstruct_polynomial(size_t max_size, lp_integer_t* p
 
     d ++;
   }
-  upolynomial_dense_touch(&p_d, d - 1);
 
-  // We only care about primitive GCDs in the reconstruction
-  upolynomial_dense_mk_primitive_Z(&p_d, 1);
+  lp_upolynomial_t* p = NULL;
 
-  // Now, multiply with cont
-  upolynomial_dense_mult_c(&p_d, lp_Z, cont);
+  // check if we were successful
+  if (integer_sgn(lp_Z, p_value) == 0) {
+    upolynomial_dense_touch(&p_d, d - 1);
 
-  // Get the sparse representation
-  lp_upolynomial_t* p = upolynomial_dense_to_upolynomial(&p_d, lp_Z);
+    // We only care about primitive GCDs in the reconstruction
+    upolynomial_dense_mk_primitive_Z(&p_d, 1);
+
+    // Now, multiply with cont
+    upolynomial_dense_mult_c(&p_d, lp_Z, cont);
+
+    // Get the sparse representation
+    p = upolynomial_dense_to_upolynomial(&p_d, lp_Z);
+  }
 
   // Free temporaries
   upolynomial_dense_destruct(&p_d);
@@ -432,36 +441,37 @@ static lp_upolynomial_t* reconstruct_polynomial(size_t max_size, lp_integer_t* p
   return p;
 }
 
-int bound_valuation(const lp_upolynomial_t* A, const lp_upolynomial_t* B, const lp_integer_t* A_cont, const lp_integer_t* B_cont) {
+static
+size_t bound_valuation(const lp_upolynomial_t* A, const lp_upolynomial_t* B, const lp_integer_t* A_cont, const lp_integer_t* B_cont) {
   // Get the highest power 2^n such that
   //  2^n >= 2*min(|pp(A)|, |pp(B)|) + 2
 
-  int A_max = 0;
-  int B_max = 0;
+  size_t A_max = 0;
+  size_t B_max = 0;
   size_t k;
 
   lp_integer_t tmp;
   integer_construct_from_int(lp_Z, &tmp, 0);
 
   for (k = 0; k < A->size; ++ k) {
     integer_div_Z(&tmp, &A->monomials[k].coefficient, A_cont);
-    int bits = integer_bits(&tmp);
+    size_t bits = integer_bits(&tmp);
     if (bits > A_max) {
       A_max = bits;
     }
   }
 
   for (k = 0; k < B->size; ++ k) {
     integer_div_Z(&tmp, &B->monomials[k].coefficient, B_cont);
-    int bits = integer_bits(&tmp);
+    size_t bits = integer_bits(&tmp);
     if (bits > B_max) {
       B_max = bits;
     }
   }
 
   integer_destruct(&tmp);
 
-  int bits_min = A_max > B_max ? B_max : A_max;
+  size_t bits_min = A_max > B_max ? B_max : A_max;
 
   return bits_min + 2;
 }
@@ -506,14 +516,19 @@ lp_upolynomial_t* upolynomial_gcd_heuristic(const lp_upolynomial_t* A, const lp_
   }
 
   // The number we use for valuation 2^n
-  int n = bound_valuation(A, B, &A_cont, &B_cont);
+  size_t n = bound_valuation(A, B, &A_cont, &B_cont);
 
   lp_integer_t A_v, B_v, D_v;
   integer_construct_from_int(lp_Z, &A_v, 0);
   integer_construct_from_int(lp_Z, &B_v, 0);
   integer_construct_from_int(lp_Z, &D_v, 0);
 
-  while (D == 0 && (attempts --)) {
+  // The evaluation shifts by up to deg(A) * n bits, an unsigned: past
+  // that, no heuristic (the caller falls back to the subresultant gcd)
+  const size_t deg_A = lp_upolynomial_degree(A);
+  const size_t max_n = deg_A > 0 ? UINT_MAX / deg_A : UINT_MAX;
+
+  while (D == 0 && n < max_n && (attempts --)) {
 
     // Evaluate A and B
     evaluate_polynomial(A, &A_cont, n, &A_v);
@@ -530,7 +545,7 @@ lp_upolynomial_t* upolynomial_gcd_heuristic(const lp_upolynomial_t* A, const lp_
     D = reconstruct_polynomial(lp_upolynomial_degree(A) + 1, &D_v, n, &d);
 
     // Check if it divides both (B is smaller degree)
-    if (!lp_upolynomial_divides(D, B) || !lp_upolynomial_divides(D, A)) {
+    if (D && (!lp_upolynomial_divides(D, B) || !lp_upolynomial_divides(D, A))) {
       lp_upolynomial_delete(D);
       D = 0;
     }

test/poly/CMakeLists.txt

@@ -1,5 +1,6 @@
 set(poly_tests
     test_feasible_int_set
+    test_upolynomial_gcd
 )
 
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -Wextra -std=c++11")

test/poly/test_upolynomial_gcd.cpp

@@ -0,0 +1,41 @@
+#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
+#include "../polyxx/doctest.h"
+
+#include <gmp.h>
+#include <poly.h>
+#include <integer.h>
+#include <upolynomial.h>
+
+TEST_CASE("upolynomial::gcd_shift_past_unsigned") {
+  // g = x^64 - 2^(2^26): the heuristic would shift by degree * n > 2^32
+  // bits, past an unsigned, so the gcd comes from the subresultant algorithm
+  lp_integer_t g_coeffs[65];
+  for (int i = 0; i < 65; i ++) {
+    lp_integer_construct_from_int(lp_Z, &g_coeffs[i], 0);
+  }
+  lp_integer_assign_int(lp_Z, &g_coeffs[0], 1);
+  mpz_mul_2exp(&g_coeffs[0], &g_coeffs[0], 1 << 26);
+  mpz_neg(&g_coeffs[0], &g_coeffs[0]);
+  lp_integer_assign_int(lp_Z, &g_coeffs[64], 1);
+  lp_upolynomial_t* g = lp_upolynomial_construct(lp_Z, 64, g_coeffs);
+
+  int p_coeffs[2] = { 1, 1 };
+  int q_coeffs[2] = { 2, 1 };
+  lp_upolynomial_t* p = lp_upolynomial_construct_from_int(lp_Z, 1, p_coeffs);
+  lp_upolynomial_t* q = lp_upolynomial_construct_from_int(lp_Z, 1, q_coeffs);
+  lp_upolynomial_t* gp = lp_upolynomial_mul(g, p);
+  lp_upolynomial_t* gq = lp_upolynomial_mul(g, q);
+
+  lp_upolynomial_t* gcd = lp_upolynomial_gcd(gp, gq);
+  CHECK(lp_upolynomial_cmp(gcd, g) == 0);
+
+  lp_upolynomial_delete(gcd);
+  lp_upolynomial_delete(gq);
+  lp_upolynomial_delete(gp);
+  lp_upolynomial_delete(q);
+  lp_upolynomial_delete(p);
+  lp_upolynomial_delete(g);
+  for (int i = 0; i < 65; i ++) {
+    lp_integer_destruct(&g_coeffs[i]);
+  }
+}