From 468538650454e524bd017fbe63e4b8c404c5531f Mon Sep 17 00:00:00 2001 From: Nick Lewycky Date: Wed, 21 Oct 2026 01:59:33 -0700 Subject: [PATCH 1/5] Fix integer overflows in the heuristic univariate gcd bound_valuation kept coefficient bit counts in an int, and the evaluation shifted coefficients by degree * n bits through an unsigned. Coefficients of 2^31 or more bits, or a degree * n of 2^32 or more, did not fit: the evaluation point or the shift came out wrong, so the gcd was wrong or reconstruct_polynomial wrote past its coefficient array. Count bits in size_t, and skip the heuristic (falling back to the subresultant gcd) when the shift does not fit in an unsigned. Co-Authored-By: Claude Opus 5.5 --- src/upolynomial/gcd.c | 24 ++++++++++++----- test/poly/CMakeLists.txt | 1 + test/poly/test_upolynomial_gcd.cpp | 41 ++++++++++++++++++++++++++++++ 3 files changed, 59 insertions(+), 7 deletions(-) create mode 100644 test/poly/test_upolynomial_gcd.cpp diff --git a/src/upolynomial/gcd.c b/src/upolynomial/gcd.c index d2d3eac0..96c89e1c 100644 --- a/src/upolynomial/gcd.c +++ b/src/upolynomial/gcd.c @@ -26,6 +26,7 @@ #include "utils/statistics.h" #include +#include STAT_DECLARE(int, upolynomial, gcd_euclid) STAT_DECLARE(int, upolynomial, gcd_euclid_extended) @@ -399,6 +400,8 @@ static lp_upolynomial_t* reconstruct_polynomial(size_t max_size, lp_integer_t* p // Basically compute digits modulo x_value int d = 0; while (integer_sgn(lp_Z, p_value)) { + // With 2^pow from bound_valuation the value has at most max_size digits + assert((size_t) d < max_size); integer_div_rem_pow2_Z(&div, &rem, p_value, pow); if (integer_bits(&rem) + 1 >= pow) { @@ -432,12 +435,12 @@ 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) { +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; @@ -445,7 +448,7 @@ int bound_valuation(const lp_upolynomial_t* A, const lp_upolynomial_t* B, const 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; } @@ -453,7 +456,7 @@ int bound_valuation(const lp_upolynomial_t* A, const lp_upolynomial_t* B, const 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; } @@ -461,7 +464,7 @@ int bound_valuation(const lp_upolynomial_t* A, const lp_upolynomial_t* B, const 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,7 +509,7 @@ 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); @@ -515,6 +518,13 @@ lp_upolynomial_t* upolynomial_gcd_heuristic(const lp_upolynomial_t* A, const lp_ 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) + size_t deg_A = lp_upolynomial_degree(A); + if (n > UINT_MAX || (deg_A > 0 && n > UINT_MAX / deg_A)) { + break; + } + // Evaluate A and B evaluate_polynomial(A, &A_cont, n, &A_v); evaluate_polynomial(B, &B_cont, n, &B_v); diff --git a/test/poly/CMakeLists.txt b/test/poly/CMakeLists.txt index 49bcf872..e054677d 100644 --- a/test/poly/CMakeLists.txt +++ b/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") diff --git a/test/poly/test_upolynomial_gcd.cpp b/test/poly/test_upolynomial_gcd.cpp new file mode 100644 index 00000000..d1b098f1 --- /dev/null +++ b/test/poly/test_upolynomial_gcd.cpp @@ -0,0 +1,41 @@ +#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN +#include "../polyxx/doctest.h" + +#include +#include +#include +#include + +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]); + } +} From f2c634d873684c83e098a50074723451efe04ff2 Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Wed, 21 Oct 2026 16:20:47 +0200 Subject: [PATCH 2/5] Update the PR to avoid unnecessary recomputations. --- src/upolynomial/gcd.c | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/upolynomial/gcd.c b/src/upolynomial/gcd.c index 96c89e1c..4e41ea82 100644 --- a/src/upolynomial/gcd.c +++ b/src/upolynomial/gcd.c @@ -435,6 +435,7 @@ static lp_upolynomial_t* reconstruct_polynomial(size_t max_size, lp_integer_t* p return p; } +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 @@ -516,14 +517,12 @@ lp_upolynomial_t* upolynomial_gcd_heuristic(const lp_upolynomial_t* A, const lp_ 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; - // The evaluation shifts by up to deg(A) * n bits, an unsigned: past - // that, no heuristic (the caller falls back to the subresultant gcd) - size_t deg_A = lp_upolynomial_degree(A); - if (n > UINT_MAX || (deg_A > 0 && n > UINT_MAX / deg_A)) { - break; - } + while (D == 0 && n < max_n && (attempts --)) { // Evaluate A and B evaluate_polynomial(A, &A_cont, n, &A_v); From 1958ca8c36d67c7a36178b1ee1b622a711f28dfa Mon Sep 17 00:00:00 2001 From: Nick Lewycky Date: Wed, 21 Oct 2026 20:59:55 -0700 Subject: [PATCH 3/5] Remove bound_valuation from the Debian symbols file bound_valuation is now static, so libpoly.so.0 no longer exports it, and dpkg-gensymbols would report the listed symbol as missing. Co-Authored-By: Claude Opus 5.5 --- debian/libpoly0.symbols | 1 - 1 file changed, 1 deletion(-) diff --git a/debian/libpoly0.symbols b/debian/libpoly0.symbols index 2beec4b5..611b212f 100644 --- a/debian/libpoly0.symbols +++ b/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 From d390ef29c3bce688737ddace9fcd94ba2cd6662a Mon Sep 17 00:00:00 2001 From: Nick Lewycky Date: Wed, 21 Oct 2026 20:59:32 -0700 Subject: [PATCH 4/5] Fail the reconstruction instead of asserting on too many digits reconstruct_polynomial asserted that the value has at most max_size digits, and release builds compile the assert out, so a violation would write past the dense polynomial. It now frees its temporaries and returns 0, and the heuristic treats that as a failed attempt (the next power, or the subresultant fallback). Co-Authored-By: Claude Opus 5.5 --- src/upolynomial/gcd.c | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/upolynomial/gcd.c b/src/upolynomial/gcd.c index 4e41ea82..4618a5ab 100644 --- a/src/upolynomial/gcd.c +++ b/src/upolynomial/gcd.c @@ -398,10 +398,18 @@ 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; + size_t d = 0; while (integer_sgn(lp_Z, p_value)) { - // With 2^pow from bound_valuation the value has at most max_size digits - assert((size_t) d < max_size); + // With 2^pow from bound_valuation the value has at most max_size digits. + // More digits mean no gcd of this degree: fail (the caller tries again + // or falls back) rather than write past the dense polynomial + if (d >= max_size) { + upolynomial_dense_destruct(&p_d); + integer_destruct(&div); + integer_destruct(&rem); + integer_destruct(&P); + return 0; + } integer_div_rem_pow2_Z(&div, &rem, p_value, pow); if (integer_bits(&rem) + 1 >= pow) { @@ -539,7 +547,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; } From dba332c4dab459e861bee872dd629ff4ecb876b7 Mon Sep 17 00:00:00 2001 From: Thomas Hader Date: Thu, 22 Oct 2026 13:28:19 +0200 Subject: [PATCH 5/5] Only one free/exit in `reconstruct_polynomial` --- src/upolynomial/gcd.c | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/src/upolynomial/gcd.c b/src/upolynomial/gcd.c index 4618a5ab..33c3b88d 100644 --- a/src/upolynomial/gcd.c +++ b/src/upolynomial/gcd.c @@ -399,17 +399,9 @@ static lp_upolynomial_t* reconstruct_polynomial(size_t max_size, lp_integer_t* p // Basically compute digits modulo x_value size_t d = 0; - while (integer_sgn(lp_Z, p_value)) { - // With 2^pow from bound_valuation the value has at most max_size digits. - // More digits mean no gcd of this degree: fail (the caller tries again - // or falls back) rather than write past the dense polynomial - if (d >= max_size) { - upolynomial_dense_destruct(&p_d); - integer_destruct(&div); - integer_destruct(&rem); - integer_destruct(&P); - return 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) { @@ -423,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; + + // check if we were successful + if (integer_sgn(lp_Z, p_value) == 0) { + upolynomial_dense_touch(&p_d, d - 1); - // Now, multiply with cont - upolynomial_dense_mult_c(&p_d, lp_Z, cont); + // We only care about primitive GCDs in the reconstruction + upolynomial_dense_mk_primitive_Z(&p_d, 1); - // Get the sparse representation - lp_upolynomial_t* p = upolynomial_dense_to_upolynomial(&p_d, lp_Z); + // 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);