From 49f3345f01b6157b39839500782836f64a916466 Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Sun, 25 Oct 2026 14:10:37 +0100 Subject: [PATCH 1/9] Use FLINT rational reconstruction when lifting parametrizations Rational reconstruction in the lifting of rational parametrizations over QQ (linear forms, multiplication matrix witnesses, trace and determinant, eliminating polynomial and coordinates) now goes through FLINT's half-gcd based fmpq_reconstruct_fmpz_2 instead of the quadratic Euclidean algorithm, falling back to ratrecon/ratreconwden when the bounds do not satisfy 2ND < m. The wrappers introduced for the lifting of Groebner bases are moved from lifting-gb.c to lifting.c so that both liftings share them. Outputs are unchanged; the time spent in CRT and rational reconstruction decreases by 25 to 40% (kat11: 1.22s to 0.70s). Co-Authored-By: Claude Opus 5.5 --- src/msolve/lifting-gb.c | 53 --------------------------------- src/msolve/lifting.c | 65 +++++++++++++++++++++++++++++++++++++---- src/msolve/msolve.c | 4 +-- 3 files changed, 62 insertions(+), 60 deletions(-) diff --git a/src/msolve/lifting-gb.c b/src/msolve/lifting-gb.c index 14f664a1..7f24b640 100644 --- a/src/msolve/lifting-gb.c +++ b/src/msolve/lifting-gb.c @@ -19,61 +19,8 @@ * Mohab Safey El Din */ #include -#include #include "streams.h" -/* rational reconstruction of u mod m with bounds N, D through FLINT - (half-gcd based) instead of the quadratic Euclidean algorithm; - u is assumed to be reduced modulo m */ -static int flint_rrec(mpz_t n, mpz_t d, const mpz_t u, const mpz_t m, - const mpz_t N, const mpz_t D){ - fmpz_t fu, fm, fN, fD; - fmpq_t q; - fmpz_init(fu); fmpz_init(fm); fmpz_init(fN); fmpz_init(fD); - fmpq_init(q); - fmpz_set_mpz(fu, u); - fmpz_set_mpz(fm, m); - fmpz_set_mpz(fN, N); - fmpz_set_mpz(fD, D); - int b = fmpq_reconstruct_fmpz_2(q, fu, fm, fN, fD); - if(b){ - fmpz_get_mpz(n, fmpq_numref(q)); - fmpz_get_mpz(d, fmpq_denref(q)); - } - fmpq_clear(q); - fmpz_clear(fu); fmpz_clear(fm); fmpz_clear(fN); fmpz_clear(fD); - return b; -} - -/* FLINT requires 2ND < m */ -static int flint_rrec_bounds_ok(const mpz_t m, rrec_data_t rd){ - mpz_mul(rd->tmp, rd->N, rd->D); - mpz_mul_2exp(rd->tmp, rd->tmp, 1); - return mpz_cmp(rd->tmp, m) < 0; -} - -/* same as ratrecon, through FLINT when the bounds allow it */ -static int lift_ratrecon(mpz_t n, mpz_t d, mpz_t u, const mpz_t mod, - rrec_data_t rd){ - if(!flint_rrec_bounds_ok(mod, rd)){ - return ratrecon(n, d, u, mod, rd); - } - mpz_mod(u, u, mod); - return flint_rrec(n, d, u, mod, rd->N, rd->D); -} - -/* same as ratreconwden, through FLINT when the bounds allow it */ -static int lift_ratreconwden(mpz_t n, mpz_t d, mpz_t u, const mpz_t mod, - const mpz_t gden, rrec_data_t rd){ - if(!flint_rrec_bounds_ok(mod, rd)){ - return ratreconwden(n, d, u, mod, gden, rd); - } - mpz_mod(u, u, mod); - mpz_mul(rd->r1, u, gden); - mpz_mod(rd->r1, rd->r1, mod); - return flint_rrec(n, d, rd->r1, mod, rd->N, rd->D); -} - #define NBCHECK 2 /* number of bits of slack kept in rational reconstruction bounds */ #define RR_MARGIN 16 diff --git a/src/msolve/lifting.c b/src/msolve/lifting.c index 9eb84c1c..d276c3ca 100644 --- a/src/msolve/lifting.c +++ b/src/msolve/lifting.c @@ -23,8 +23,63 @@ #ifndef MAX #define MAX(a, b) (((a) > (b)) ? (a) : (b)) #endif -/* Initialization of sparse fglm matrix for crt */ #include "msolve-data.h" +#include +#include + +/* rational reconstruction of u mod m with bounds N, D through FLINT + (half-gcd based) instead of the quadratic Euclidean algorithm; + u is assumed to be reduced modulo m */ +static int flint_rrec(mpz_t n, mpz_t d, const mpz_t u, const mpz_t m, + const mpz_t N, const mpz_t D){ + fmpz_t fu, fm, fN, fD; + fmpq_t q; + fmpz_init(fu); fmpz_init(fm); fmpz_init(fN); fmpz_init(fD); + fmpq_init(q); + fmpz_set_mpz(fu, u); + fmpz_set_mpz(fm, m); + fmpz_set_mpz(fN, N); + fmpz_set_mpz(fD, D); + int b = fmpq_reconstruct_fmpz_2(q, fu, fm, fN, fD); + if(b){ + fmpz_get_mpz(n, fmpq_numref(q)); + fmpz_get_mpz(d, fmpq_denref(q)); + } + fmpq_clear(q); + fmpz_clear(fu); fmpz_clear(fm); fmpz_clear(fN); fmpz_clear(fD); + return b; +} + +/* FLINT requires 2ND < m */ +static int flint_rrec_bounds_ok(const mpz_t m, rrec_data_t rd){ + mpz_mul(rd->tmp, rd->N, rd->D); + mpz_mul_2exp(rd->tmp, rd->tmp, 1); + return mpz_cmp(rd->tmp, m) < 0; +} + +/* same as ratrecon, through FLINT when the bounds allow it */ +static int lift_ratrecon(mpz_t n, mpz_t d, mpz_t u, const mpz_t mod, + rrec_data_t rd){ + if(!flint_rrec_bounds_ok(mod, rd)){ + return ratrecon(n, d, u, mod, rd); + } + mpz_mod(u, u, mod); + return flint_rrec(n, d, u, mod, rd->N, rd->D); +} + +/* same as ratreconwden, through FLINT when the bounds allow it */ +static int lift_ratreconwden(mpz_t n, mpz_t d, mpz_t u, const mpz_t mod, + const mpz_t gden, rrec_data_t rd){ + if(!flint_rrec_bounds_ok(mod, rd)){ + return ratreconwden(n, d, u, mod, gden, rd); + } + mpz_mod(u, u, mod); + mpz_mul(rd->r1, u, gden); + mpz_mod(rd->r1, rd->r1, mod); + return flint_rrec(n, d, rd->r1, mod, rd->N, rd->D); +} + +/* Initialization of sparse fglm matrix for crt */ static inline void crt_mpz_matfglm_initset(crt_mpz_matfglm_t crt_mat, sp_matfglm_t *mod_mat){ crt_mat->ncols = mod_mat->ncols; @@ -527,7 +582,7 @@ static inline int old_rat_recon_trace_det(trace_det_fglm_mat_t trace_det, static inline int rat_recon_array(mpz_t *res, mpz_t *crt, int32_t sz, mpz_t modulus, rrec_data_t rdata) { for (int i = 0; i < sz; i++) { - int b = ratrecon(res[2 * i], res[2 * i + 1], crt[i], modulus, rdata); + int b = lift_ratrecon(res[2 * i], res[2 * i + 1], crt[i], modulus, rdata); if (b == 0) return 0; } @@ -572,7 +627,7 @@ static inline int rat_recon_trace_det(trace_det_fglm_mat_t trace_det, mpz_init(gcd); int nlifted = trace_det->w_checked; for(int32_t i = trace_det->w_checked; i < trace_det->nrows; i++){ - int b = ratreconwden(rnum, rden, trace_det->matmul_wcrt[i], modulus, gden, recdata); + int b = lift_ratreconwden(rnum, rden, trace_det->matmul_wcrt[i], modulus, gden, recdata); if(b == 1){ mpz_set(trace_det->matmul_wqq[2*i], rnum); mpz_set(trace_det->matmul_wqq[2*i + 1], rden); @@ -596,7 +651,7 @@ static inline int rat_recon_trace_det(trace_det_fglm_mat_t trace_det, int trace_rec = 0; if(trace_det->done_trace < 2){ int b = - ratreconwden(rnum, rden, trace_det->trace_crt, modulus, gden, recdata); + lift_ratreconwden(rnum, rden, trace_det->trace_crt, modulus, gden, recdata); if (b == 1) { mpz_set(trace_det->trace_num, rnum); mpz_set(trace_det->trace_den, rden); @@ -617,7 +672,7 @@ static inline int rat_recon_trace_det(trace_det_fglm_mat_t trace_det, if(trace_det->done_det < 2 && trace_det->done_trace){ if(trace_rec){ } - int b = ratreconwden(rnum, rden, trace_det->det_crt, modulus, gden, recdata); + int b = lift_ratreconwden(rnum, rden, trace_det->det_crt, modulus, gden, recdata); if (b == 1) { mpz_set(trace_det->det_num, rnum); diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index cb82a156..a93a9157 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -1047,7 +1047,7 @@ static inline int new_rational_reconstruction_mpz_ptr_with_denom_param( //for (deg_t i = 0; i < len; i++) { for (deg_t i = len - 1; i >= 0; i--) { - b = ratreconwden(tmp_num[i], rden, pol[i], modulus, lcm, + b = lift_ratreconwden(tmp_num[i], rden, pol[i], modulus, lcm, rdata); if (b == 0) { @@ -1099,7 +1099,7 @@ static inline int rational_reconstruction_mpz_ptr_with_denom_elim( int b = 0; for(deg_t i = len - 1; i >= 0; i--){ if(check_lift[i] <= 1){ - b = ratreconwden(tmp_num[i], rden, pol[i], modulus, lcm, rdata); + b = lift_ratreconwden(tmp_num[i], rden, pol[i], modulus, lcm, rdata); if(b == 0){ check_lift[i] = 0; return 0; From bed9b068d9b2daef011f425f65f2204ac44ef056 Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Sun, 25 Oct 2026 14:23:17 +0100 Subject: [PATCH 2/9] Make the reconstruction period of parametrizations depend on primes only When lifting rational parametrizations over QQ, rational reconstruction was called every nbdoit iterations, nbdoit being doubled whenever the reconstruction took more than 20% of the modular computations; the number of primes used thus depended on timings and on the machine load. As for the lifting of Groebner bases, the period is now rr primes, rr being the largest power of 2 <= nprimes/128, i.e. max(1, rr/nthrds) iterations, and the reconstruction is attempted once per iteration, after the last good prime, instead of after each prime of the iteration. Outputs are unchanged and the number of primes is now reproducible. Co-Authored-By: Claude Opus 5.5 --- src/msolve/msolve.c | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index a93a9157..b4394938 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -2681,7 +2681,14 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, } /* scrr measures time spent in ratrecon for modular images */ double crr = 0, scrr = 0; - /* CRT + rational reconstruction */ + /* CRT + rational reconstruction; reconstruction is attempted only + * once per iteration, after the last good prime has been used */ + int32_t lastgood = -1; + for (len_t i = 0; i < st->nthrds; i++) { + if (bad_primes[i] == 0) { + lastgood = i; + } + } for (len_t i = 0; i < st->nthrds; i++) { if (bad_primes[i] == 0) { if (rerun == 0) { @@ -2700,7 +2707,8 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, denom, check_lift, modulus, prod_crt, lp->p[i], &result, rnum, rden, recdata, guessed_num, guessed_den, &maxrec, &matrec, &oldmatrec_checked, &matrec_checked, is_lifted, - &mat_lifted, &lin_lifted, doit, nbdoit, st->nthrds, info_level); + &mat_lifted, &lin_lifted, doit && (int32_t)i == lastgood, nbdoit, + st->nthrds, info_level); if (br == 1) { rerun = 0; @@ -2731,13 +2739,21 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, } strat += scrr; - double t = ((double)nbdoit) * ca1; - if ((t == 0) || (scrr >= 0.2 * t && br == 0)) { - nbdoit = 2 * nbdoit; + /* rational reconstruction is called every nbdoit iterations, i.e. + * about every rr primes, rr being the largest power of 2 + * <= nprimes / RR_PERIOD_DIV; the period depends only on the number of + * primes (not on timings) so that the number of primes used is + * reproducible */ + int32_t rr = 1; + while (2 * rr <= nprimes / RR_PERIOD_DIV) { + rr *= 2; + } + rr = MAX(1, rr / st->nthrds); + if (rr != nbdoit) { + nbdoit = rr; lpow2 = 2 * nprimes; - doit = 0; if (info_level) { - fprintf(VERBSTREAM, "\n", nbdoit, scrr, t); + fprintf(VERBSTREAM, "\n", nbdoit); fflush(VERBSTREAM); } } From 3d3d0bd75eb875f636b4e9c36271891b76c75934 Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Sun, 25 Oct 2026 15:38:51 +0100 Subject: [PATCH 3/9] Apply CRT to parametrizations only when reconstructing them When lifting rational parametrizations over QQ, every coefficient of the eliminating polynomial and of the coordinates was lifted by incremental CRT at each prime, which is quadratic in the number of primes, although these CRT images are only used once the trace and determinant witnesses are reconstructed and verified. The modular images of the parametrization are now stored for each prime (in chunks of 64 primes which are never moved), and FLINT's multi-CRT is applied to them only when a rational reconstruction is attempted, one coefficient at a time, right before this coefficient is reconstructed; no CRT is thus applied beyond the first failed reconstruction, and no array of CRT images is kept. The coordinates are lifted one at a time, trying in turn the three bounds which were tried before on all of them. The witnesses are still lifted by incremental CRT. Outputs and numbers of primes are unchanged on 38 benchmark systems (eco, kat, noon, cyclic, henrion, chandra, random, minrank, cp_d, crit and others). The time spent in CRT and rational reconstruction decreases by up to 4.2x (random2_9_9: 20.4s to 4.8s; chandra-11: 50.7s to 14.5s), the total time by up to 17% with one thread and by up to 32% with four, and peak memory is unchanged or slightly lower. Co-Authored-By: Claude Opus 5.5 --- src/msolve/msolve.c | 344 ++++++++++++++++++++++++-------------------- 1 file changed, 191 insertions(+), 153 deletions(-) diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index b4394938..0b29b0fe 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -970,62 +970,128 @@ static inline void normalize_nmod_param(param_t *nmod_param) { } } -/* -On copie nmod_param dans mpz_param -On suppose que les allocations ont ete faites via initialize_mpz_param avant. - */ -static inline void set_mpz_param_nmod(mpz_param_t mpz_param, - param_t *nmod_param) { - mpz_param->elim->length = nmod_param->elim->length; - for (long i = 0; i <= mpz_param->nsols; i++) { - mpz_set_ui(mpz_param->elim->coeffs[i], (nmod_param)->elim->coeffs[i]); - } - mpz_param->elim->length = nmod_param->elim->length; +/* modular images of the parametrization, stored for each prime in one + contiguous block (elim, then the coordinates), so that CRT is only + applied when a rational reconstruction is attempted; blocks are grouped + by chunks of PARAM_IMAGES_CHUNK primes, which are never moved */ +#define PARAM_IMAGES_CHUNK 64 +typedef struct{ + uint32_t **chunks; /* images, prime-major */ + mp_limb_t *primes; /* primes of the stored images */ + int32_t nprimes; + int32_t alloc; /* number of primes which fit in primes */ + deg_t elen; /* number of coefficients of elim */ + deg_t clen; /* number of coefficients of each coordinate */ + nvars_t nc; /* number of coordinates */ + int64_t sz; /* number of coefficients stored per prime */ +} param_images_struct; + +typedef param_images_struct param_images_t[1]; + +static inline void param_images_init(param_images_t im, + const mpz_param_t mpz_param){ + im->elen = mpz_param->elim->length; + im->nc = mpz_param->nvars - 1; + im->clen = im->nc > 0 ? mpz_param->coords[0]->length : 0; + im->sz = im->elen + ((int64_t)im->nc) * im->clen; + im->nprimes = 0; + im->alloc = 0; + im->chunks = NULL; + im->primes = NULL; +} - for (long i = 0; i < nmod_param->denom->length; i++) { - mpz_set_ui(mpz_param->denom->coeffs[i], nmod_param->denom->coeffs[i]); - } - mpz_param->denom->length = nmod_param->denom->length; +static inline void param_images_clear(param_images_t im){ + int32_t nchunks = (im->nprimes + PARAM_IMAGES_CHUNK - 1) / PARAM_IMAGES_CHUNK; + for (int32_t i = 0; i < nchunks; i++) { + free(im->chunks[i]); + } + free(im->chunks); + flint_free(im->primes); + im->chunks = NULL; + im->primes = NULL; + im->nprimes = 0; + im->alloc = 0; +} - for (int j = 0; j < mpz_param->nvars - 1; j++) { +/* coefficients beyond the length of a modular polynomial are zero */ +static inline void copy_nmod_coeffs(uint32_t *dst, const nmod_poly_t pol, + const deg_t len){ + deg_t l = MIN(len, (deg_t)pol->length); + for (deg_t i = 0; i < l; i++) { + dst[i] = pol->coeffs[i]; + } + for (deg_t i = l; i < len; i++) { + dst[i] = 0; + } +} - for (long i = 0; i < nmod_param->coords[j]->length; i++) { - mpz_set_ui(mpz_param->coords[j]->coeffs[i], - nmod_param->coords[j]->coeffs[i]); - } - for (long i = nmod_param->coords[j]->length; - i < nmod_param->elim->length - 1; i++) { - mpz_set_ui(mpz_param->coords[j]->coeffs[i], 0); +static inline void param_images_add(param_images_t im, + const param_t *nmod_param, + const uint32_t prime){ + if (im->nprimes == im->alloc) { + im->alloc += PARAM_IMAGES_CHUNK; + im->primes = flint_realloc(im->primes, im->alloc * sizeof(mp_limb_t)); + im->chunks = realloc(im->chunks, + (im->alloc / PARAM_IMAGES_CHUNK) * sizeof(uint32_t *)); + uint32_t *chunk = malloc(PARAM_IMAGES_CHUNK * im->sz * sizeof(uint32_t)); + if (im->chunks == NULL || chunk == NULL) { + fprintf(ERRSTREAM, "Problem when allocating modular images (%ld primes)\n", + (long)im->alloc); + exit(1); } - mpz_param->coords[j]->length = nmod_param->elim->length - 1; + im->chunks[im->alloc / PARAM_IMAGES_CHUNK - 1] = chunk; + } + uint32_t *dst = im->chunks[im->nprimes / PARAM_IMAGES_CHUNK] + + (im->nprimes % PARAM_IMAGES_CHUNK) * im->sz; + copy_nmod_coeffs(dst, nmod_param->elim, im->elen); + for (nvars_t j = 0; j < im->nc; j++) { + copy_nmod_coeffs(dst + im->elen + j * im->clen, nmod_param->coords[j], + im->clen); } + im->primes[im->nprimes] = prime; + im->nprimes++; } -static inline void crt_lift_mpz_upoly(mpz_upoly_t pol, nmod_poly_t nmod_pol, - mpz_t modulus, int32_t prime, mpz_t prod, - mpz_t tmp, int nthrds) { - len_t i; - /* #pragma omp parallel for num_threads(nthrds) \ */ - /* private(i) schedule(static) */ - for (i = 0; i < pol->length; i++) { - mpz_CRT_ui(pol->coeffs[i], pol->coeffs[i], modulus, nmod_pol->coeffs[i], - prime, prod, tmp, 0); - } +/* data used to apply CRT to the stored images, built once per attempt */ +typedef struct{ + fmpz_comb_t comb; + fmpz_comb_temp_t comb_temp; + mp_limb_t *residues; + fmpz_t y; + mpz_t u; /* CRT of the current coefficient */ +} param_crt_struct; + +typedef param_crt_struct param_crt_t[1]; + +static inline void param_crt_init(param_crt_t pc, const param_images_t im){ + fmpz_comb_init(pc->comb, im->primes, im->nprimes); + fmpz_comb_temp_init(pc->comb_temp, pc->comb); + pc->residues = flint_malloc(im->nprimes * sizeof(mp_limb_t)); + fmpz_init(pc->y); + mpz_init(pc->u); } -/* assumes that all degrees are the same */ -static inline void crt_lift_mpz_param(mpz_param_t mpz_param, - param_t *nmod_param, mpz_t modulus, - mpz_t prod_crt, const int32_t prime, - mpz_t tmp, const int nthrds) { +static inline void param_crt_clear(param_crt_t pc){ + mpz_clear(pc->u); + fmpz_clear(pc->y); + flint_free(pc->residues); + fmpz_comb_temp_clear(pc->comb_temp); + fmpz_comb_clear(pc->comb); +} - /*assumes prod_crt = modulus * prime */ - crt_lift_mpz_upoly(mpz_param->elim, nmod_param->elim, modulus, prime, - prod_crt, tmp, nthrds); - for (len_t i = 0; i < mpz_param->nvars - 1; i++) { - crt_lift_mpz_upoly(mpz_param->coords[i], nmod_param->coords[i], modulus, - prime, prod_crt, tmp, nthrds); - } +/* sets pc->u to the CRT of the stored images of the coefficient at + position pos, in [0, modulus); coefficients are lifted one at a time, + right before being reconstructed, so that no array of CRT images is kept + and no CRT is applied beyond the first failed reconstruction */ +static inline void param_images_crt_coeff(const param_images_t im, + const int64_t pos, + param_crt_struct *pc){ + for (int32_t k = 0; k < im->nprimes; k++) { + pc->residues[k] = im->chunks[k / PARAM_IMAGES_CHUNK] + [(k % PARAM_IMAGES_CHUNK) * im->sz + pos]; + } + fmpz_multi_CRT_ui(pc->y, pc->residues, pc->comb, pc->comb_temp, 0); + fmpz_get_mpz(pc->u, pc->y); } /** @@ -1037,7 +1103,8 @@ static inline void crt_lift_mpz_param(mpz_param_t mpz_param, static inline int new_rational_reconstruction_mpz_ptr_with_denom_param( - mpz_t *recons, mpz_t denominator, mpz_t *pol, deg_t len, mpz_t modulus, + mpz_t *recons, mpz_t denominator, const param_images_t im, int64_t start, + param_crt_struct *pc, deg_t len, mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_t *tmp_num, mpz_t *tmp_den, mpz_t lcm, mpz_t gnum, mpz_t guessed_den, rrec_data_t rdata) { @@ -1047,7 +1114,8 @@ static inline int new_rational_reconstruction_mpz_ptr_with_denom_param( //for (deg_t i = 0; i < len; i++) { for (deg_t i = len - 1; i >= 0; i--) { - b = lift_ratreconwden(tmp_num[i], rden, pol[i], modulus, lcm, + param_images_crt_coeff(im, start + i, pc); + b = lift_ratreconwden(tmp_num[i], rden, pc->u, modulus, lcm, rdata); if (b == 0) { @@ -1078,7 +1146,8 @@ static inline int new_rational_reconstruction_mpz_ptr_with_denom_param( static inline int rational_reconstruction_mpz_ptr_with_denom_elim( - mpz_t *recons, mpz_t denominator, mpz_t *pol, deg_t len, mpz_t modulus, + mpz_t *recons, mpz_t denominator, const param_images_t im, int64_t start, + param_crt_struct *pc, deg_t len, mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_t *tmp_num, mpz_t *tmp_den, int16_t *check_lift, mpz_t lcm, mpz_t gnum, mpz_t guessed_den, rrec_data_t rdata) { @@ -1099,7 +1168,8 @@ static inline int rational_reconstruction_mpz_ptr_with_denom_elim( int b = 0; for(deg_t i = len - 1; i >= 0; i--){ if(check_lift[i] <= 1){ - b = lift_ratreconwden(tmp_num[i], rden, pol[i], modulus, lcm, rdata); + param_images_crt_coeff(im, start + i, pc); + b = lift_ratreconwden(tmp_num[i], rden, pc->u, modulus, lcm, rdata); if(b == 0){ check_lift[i] = 0; return 0; @@ -1134,34 +1204,37 @@ static inline int rational_reconstruction_mpz_ptr_with_denom_elim( static inline int rational_reconstruction_upoly_with_denom_elim( - mpz_upoly_t recons, mpz_t denominator, mpz_upoly_t pol, long len, + mpz_upoly_t recons, mpz_t denominator, const param_images_t im, + int64_t start, param_crt_struct *pc, long len, mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_upoly_t tmp_num, mpz_upoly_t tmp_den, int16_t *check_lift, mpz_t lcm, mpz_t guessed_num, mpz_t guessed_den, rrec_data_t rdata) { return rational_reconstruction_mpz_ptr_with_denom_elim( - recons->coeffs, denominator, pol->coeffs, len, modulus, maxrec, coef, + recons->coeffs, denominator, im, start, pc, len, modulus, maxrec, coef, rnum, rden, tmp_num->coeffs, tmp_den->coeffs, check_lift, lcm, guessed_num, guessed_den, rdata); } static inline int rational_reconstruction_upoly_with_denom_param( - mpz_upoly_t recons, mpz_t denominator, mpz_upoly_t pol, long len, + mpz_upoly_t recons, mpz_t denominator, const param_images_t im, + int64_t start, param_crt_struct *pc, long len, mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_upoly_t tmp_num, mpz_upoly_t tmp_den, mpz_t lcm, mpz_t guessed_num, mpz_t guessed_den, rrec_data_t rdata) { return new_rational_reconstruction_mpz_ptr_with_denom_param( - recons->coeffs, denominator, pol->coeffs, len, modulus, maxrec, coef, + recons->coeffs, denominator, im, start, pc, len, modulus, maxrec, coef, rnum, rden, tmp_num->coeffs, tmp_den->coeffs, lcm, guessed_num, guessed_den, rdata); } -static inline void lift_coordinate(mpz_param_t mpz_param, mpz_param_t tmp_mpz_param, +static inline void lift_coordinate(mpz_param_t mpz_param, + const param_images_t images, param_crt_struct *pc, param_t *nmod_param, mpz_upoly_t numer, mpz_upoly_t denom, int16_t *check_lift, mpz_t modulus, mpz_t prod_crt, int32_t prime, mpq_t *coef, mpz_t rnum, mpz_t rden, rrec_data_t recdata, mpz_t guessed_num, mpz_t guessed_den, deg_t *maxrec, @@ -1178,8 +1251,10 @@ static inline void lift_coordinate(mpz_param_t mpz_param, mpz_param_t tmp_mpz_pa if (is_lifted[0] > 0 && is_lifted[idx + 1] == 0) { int b = rational_reconstruction_upoly_with_denom_param( - mpz_param->coords[idx], denominator, tmp_mpz_param->coords[idx], - nmod_param->coords[idx]->length, modulus, maxrec, coef, rnum, rden, + mpz_param->coords[idx], denominator, images, + images->elen + idx * images->clen, pc, + MIN((deg_t)nmod_param->coords[idx]->length, images->clen), + modulus, maxrec, coef, rnum, rden, numer, denom, lcm, guessed_num, guessed_den, recdata); if(b){ is_lifted[idx + 1] = 1; @@ -1198,7 +1273,8 @@ static inline void lift_coordinate(mpz_param_t mpz_param, mpz_param_t tmp_mpz_pa } -static inline int lift_parametrization(mpz_param_t mpz_param, mpz_param_t tmp_mpz_param, +static inline int lift_parametrization(mpz_param_t mpz_param, + const param_images_t images, param_t *nmod_param, mpz_upoly_t numer, mpz_upoly_t denom, int16_t *check_lift, mpz_t modulus, mpz_t prod_crt, int32_t prime, mpq_t *coef, mpz_t rnum, mpz_t rden, rrec_data_t recdata, mpz_t guessed_num, mpz_t guessed_den, deg_t *maxrec, @@ -1206,6 +1282,11 @@ static inline int lift_parametrization(mpz_param_t mpz_param, mpz_param_t tmp_mp mpz_set_ui(lcm, 1); mpz_set_ui(denominator, 1); + /* CRT is applied to the stored images of each coefficient right before + its rational reconstruction */ + param_crt_t pc; + param_crt_init(pc, images); + int b = 0; if (is_lifted[0] == 0) { @@ -1214,11 +1295,12 @@ static inline int lift_parametrization(mpz_param_t mpz_param, mpz_param_t tmp_mp mpz_swap(recdata->N1, recdata->N2); mpz_swap(recdata->D1, recdata->D2); b = rational_reconstruction_upoly_with_denom_elim( - mpz_param->elim, denominator, tmp_mpz_param->elim, + mpz_param->elim, denominator, images, 0, pc, nmod_param->elim->length, modulus, maxrec, coef, rnum, rden, numer, denom, check_lift, lcm, guessed_num, guessed_den, recdata); if (b == 0) { is_lifted[0] = 0; + param_crt_clear(pc); return 0; } is_lifted[0] = 1; @@ -1240,92 +1322,46 @@ static inline int lift_parametrization(mpz_param_t mpz_param, mpz_param_t tmp_mp mpz_set(guessed_den, lc); - mpz_sqrt(recdata->D2, modulus); - mpz_fdiv_q(recdata->N2, modulus, recdata->D2); - mpz_fdiv_q_2exp(recdata->D2, recdata->D2, 1); - mpz_root(recdata->D1, modulus, 16); - mpz_fdiv_q(recdata->N1, modulus, recdata->D1); - mpz_fdiv_q_2exp(recdata->D1, recdata->D1, 1); - mpz_root(recdata->D, modulus, 32); - mpz_fdiv_q(recdata->N, modulus, recdata->D); - mpz_fdiv_q_2exp(recdata->D, recdata->D, 1); - if(mpz_sizeinbase(guessed_den, 2) <= mpz_sizeinbase(modulus, 2)){ - mpz_set(lcm, guessed_den); - } - else{ - mpz_set_ui(lcm, 1); - } - for(int i = 0; i < nc; i++){ - lift_coordinate(mpz_param, tmp_mpz_param, - nmod_param, numer, denom, check_lift, - modulus, prod_crt, prime, coef, rnum, rden, - recdata, guessed_num, guessed_den, maxrec, - is_lifted, denominator, lcm, c, i, info_level); - if(mpz_sizeinbase(guessed_den, 2) <= mpz_sizeinbase(modulus, 2)){ - mpz_set(lcm, guessed_den); - } - else{ - mpz_set_ui(lcm, 1); - } - } - mpz_swap(recdata->N, recdata->N2); - mpz_swap(recdata->D, recdata->D2); - for(int i = 0; i < nc; i++){ - lift_coordinate(mpz_param, tmp_mpz_param, - nmod_param, numer, denom, check_lift, - modulus, prod_crt, prime, coef, rnum, rden, - recdata, guessed_num, guessed_den, maxrec, - is_lifted, denominator, lcm, c, i, info_level); - if(mpz_sizeinbase(guessed_den, 2) <= mpz_sizeinbase(modulus, 2)){ - mpz_set(lcm, guessed_den); - } - else{ - mpz_set_ui(lcm, 1); - } - } - b = 1; - for(int i = 0; i < nc; i++){ - if(is_lifted[i + 1] == 0){ - b = 0; - } - } - if(b == 1){ - mpz_clear(lc); - return 1; - } + /* bounds tried in turn for each coordinate: D = m^(1/32), then balanced + bounds, then N = m^(1/16) */ + mpz_t bN[3], bD[3]; + for(int k = 0; k < 3; k++){ + mpz_init(bN[k]); + mpz_init(bD[k]); + } + mpz_root(bD[0], modulus, 32); + mpz_fdiv_q(bN[0], modulus, bD[0]); + mpz_fdiv_q_2exp(bD[0], bD[0], 1); + mpz_sqrt(bD[1], modulus); + mpz_fdiv_q(bN[1], modulus, bD[1]); + mpz_fdiv_q_2exp(bD[1], bD[1], 1); + mpz_root(bN[2], modulus, 16); + mpz_fdiv_q(bD[2], modulus, bN[2]); + mpz_fdiv_q_2exp(bN[2], bN[2], 1); - mpz_swap(recdata->N, recdata->N2); - mpz_swap(recdata->D, recdata->D2); - mpz_swap(recdata->N, recdata->N1); - mpz_swap(recdata->D, recdata->D1); - mpz_swap(recdata->N, recdata->D); for(int i = 0; i < nc; i++){ - lift_coordinate(mpz_param, tmp_mpz_param, - nmod_param, numer, denom, check_lift, - modulus, prod_crt, prime, coef, rnum, rden, - recdata, guessed_num, guessed_den, maxrec, - is_lifted, denominator, lcm, c, i, info_level); - if(mpz_sizeinbase(guessed_den, 2) <= mpz_sizeinbase(modulus, 2)){ - mpz_set(lcm, guessed_den); - } - else{ - mpz_set_ui(lcm, 1); + if(is_lifted[i + 1]){ + continue; } - } - b = 1; - for(int i = 0; i < nc; i++){ - if(is_lifted[i + 1] == 0){ - b = 0; + for(int k = 0; k < 3 && is_lifted[i + 1] == 0; k++){ + mpz_set(recdata->N, bN[k]); + mpz_set(recdata->D, bD[k]); + lift_coordinate(mpz_param, images, pc, + nmod_param, numer, denom, check_lift, + modulus, prod_crt, prime, coef, rnum, rden, + recdata, guessed_num, guessed_den, maxrec, + is_lifted, denominator, lcm, c, i, info_level); } } - if(b == 1){ - mpz_clear(lc); - return 1; - } - else{ - mpz_clear(lc); - return 1; + + for(int k = 0; k < 3; k++){ + mpz_clear(bN[k]); + mpz_clear(bD[k]); } + mpq_clear(c); + mpz_clear(lc); + param_crt_clear(pc); + return 1; } /** @@ -1335,7 +1371,8 @@ returns 0 if rational reconstruction failed **/ static inline int rational_reconstruction_param( - mpz_param_t mpz_param, mpz_param_t tmp_mpz_param, param_t *nmod_param, + mpz_param_t mpz_param, param_images_t images, + param_t *nmod_param, nvars_t nlins, nvars_t *linvars, uint32_t *lineqs, trace_det_fglm_mat_t trace_det, sp_matfglm_t *mat, mpz_upoly_t numer, mpz_upoly_t denom, int16_t *check_lift, mpz_t modulus, mpz_t prod_crt, int32_t prime, @@ -1365,8 +1402,9 @@ static inline int rational_reconstruction_param( } /** CRT PART **/ - crt_lift_mpz_param(tmp_mpz_param, nmod_param, modulus, prod_crt, prime, - trace_det->tmp, nthrds); + /* the images of the parametrization are only stored; CRT is applied when + a rational reconstruction is attempted (see lift_parametrization) */ + param_images_add(images, nmod_param, prime); crt_lift_trace_det(trace_det, trace_mod, det_mod, mat, lineqs, nmod_param, @@ -1423,7 +1461,7 @@ static inline int rational_reconstruction_param( mpq_t c; mpq_init(c); - int b = lift_parametrization(mpz_param, tmp_mpz_param, + int b = lift_parametrization(mpz_param, images, nmod_param, numer, denom, check_lift, modulus, prod_crt, prime, coef, rnum, rden, recdata, guessed_num, guessed_den, maxrec, @@ -1982,12 +2020,11 @@ static inline int is_lucky_matmul_prime_ui(uint32_t prime, } -static void free_msolve_trace_qq_secondary_data(mpz_param_t tmp_mpz_param, +static void free_msolve_trace_qq_secondary_data( trace_det_fglm_mat_t trace_det, mpz_t modulus, mpz_t prod_crt, mpq_t result, mpq_t test, mpz_t rnum, mpz_t rden, mpz_upoly_t numer, mpz_upoly_t denom, int16_t *check_lift, mpz_t guessed_num, mpz_t guessed_den, int *is_lifted, rrec_data_t recdata){ - mpz_param_clear(tmp_mpz_param); trace_det_clear(trace_det); mpz_clear(modulus); mpz_clear(prod_crt); @@ -2465,11 +2502,8 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, bsquvars); normalize_nmod_param(nmod_params[0]); - mpz_param_t tmp_mpz_param; - mpz_param_init(tmp_mpz_param); initialize_mpz_param(*mpz_paramp, nmod_params[0]); - initialize_mpz_param(tmp_mpz_param, nmod_params[0]); // attention les longueurs des mpz_param sont fixees par nmod_params[0] // dans des cas exceptionnels, ca peut augmenter avec un autre premier. @@ -2505,9 +2539,11 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, mpz_set_ui(rnum, 0); mpz_init(rden); mpz_set_ui(rden, 1); - set_mpz_param_nmod(tmp_mpz_param, nmod_params[0]); + param_images_t images; + param_images_init(images, *mpz_paramp); + param_images_add(images, nmod_params[0], primeinit); - deg_t nsols = tmp_mpz_param->nsols; + deg_t nsols = (*mpz_paramp)->nsols; mpz_upoly_t numer; mpz_upoly_init2(numer, (nsols + 1), 32 * (nsols + 1)); @@ -2701,7 +2737,7 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, crr = realtime(); if (mcheck == 1) { br = rational_reconstruction_param( - *mpz_paramp, tmp_mpz_param, nmod_params[i], + *mpz_paramp, images, nmod_params[i], bnlins[i], blinvars[i], lineqs_ptr[i], trace_det, bmatrix[i], numer, denom, check_lift, modulus, prod_crt, lp->p[i], &result, rnum, rden, recdata, @@ -2730,7 +2766,8 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, blens_extra_nf, bexps_extra_nf, bcfs_extra_nf, bdata_fglm, bdata_bms, num_gb, leadmons_ori, leadmons_current, bnlins, blinvars, linvars, lineqs_ptr, bsquvars, squvars, lmb_ori, field_char); - free_msolve_trace_qq_secondary_data(tmp_mpz_param, trace_det, modulus, + param_images_clear(images); + free_msolve_trace_qq_secondary_data(trace_det, modulus, prod_crt, result, test, rnum, rden, numer, denom, check_lift, guessed_num, guessed_den, is_lifted, recdata); return -4; @@ -2799,7 +2836,8 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, blens_extra_nf, bexps_extra_nf, bcfs_extra_nf, bdata_fglm, bdata_bms, num_gb, leadmons_ori, leadmons_current, bnlins, blinvars, linvars, lineqs_ptr, bsquvars, squvars, lmb_ori, field_char); - free_msolve_trace_qq_secondary_data(tmp_mpz_param, trace_det, modulus, + param_images_clear(images); + free_msolve_trace_qq_secondary_data(trace_det, modulus, prod_crt, result, test, rnum, rden, numer, denom, check_lift, guessed_num, guessed_den, is_lifted, recdata); return 0; From b3bc41771ab76d1d70b4ebe3860c5b2760bb0f1a Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Sun, 25 Oct 2026 15:47:30 +0100 Subject: [PATCH 4/9] Keep a 16-bit margin in the reconstruction bounds of parametrizations As for the lifting of Groebner bases, the rational reconstruction bounds used when lifting parametrizations over QQ (witnesses, eliminating polynomial and coordinates) are now computed from modulus / 2^RR_MARGIN instead of the modulus, so that 2ND <= modulus / 2^16. This makes spurious reconstructions unlikely: on cyclic7, two spurious reconstructions of the eliminating polynomial, each rejected by the next prime, are avoided. Outputs are unchanged on 35 benchmark systems; the number of primes increases by 0 to 4 on all of them but random2_9_9 (2627 to 2643, the reconstruction period being 16 primes at that point) and cyclic7 (194 to 193). Co-Authored-By: Claude Opus 5.5 --- src/msolve/msolve.c | 40 +++++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index 0b29b0fe..32732d38 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -1322,22 +1322,26 @@ static inline int lift_parametrization(mpz_param_t mpz_param, mpz_set(guessed_den, lc); - /* bounds tried in turn for each coordinate: D = m^(1/32), then balanced - bounds, then N = m^(1/16) */ - mpz_t bN[3], bD[3]; + /* bounds tried in turn for each coordinate: D = M^(1/32), then balanced + bounds, then N = M^(1/16), with M = modulus / 2^RR_MARGIN so that + RR_MARGIN bits of slack are kept */ + mpz_t bN[3], bD[3], mm; for(int k = 0; k < 3; k++){ mpz_init(bN[k]); mpz_init(bD[k]); } - mpz_root(bD[0], modulus, 32); - mpz_fdiv_q(bN[0], modulus, bD[0]); + mpz_init(mm); + mpz_fdiv_q_2exp(mm, modulus, RR_MARGIN); + mpz_root(bD[0], mm, 32); + mpz_fdiv_q(bN[0], mm, bD[0]); mpz_fdiv_q_2exp(bD[0], bD[0], 1); - mpz_sqrt(bD[1], modulus); - mpz_fdiv_q(bN[1], modulus, bD[1]); + mpz_sqrt(bD[1], mm); + mpz_fdiv_q(bN[1], mm, bD[1]); mpz_fdiv_q_2exp(bD[1], bD[1], 1); - mpz_root(bN[2], modulus, 16); - mpz_fdiv_q(bD[2], modulus, bN[2]); + mpz_root(bN[2], mm, 16); + mpz_fdiv_q(bD[2], mm, bN[2]); mpz_fdiv_q_2exp(bN[2], bN[2], 1); + mpz_clear(mm); for(int i = 0; i < nc; i++){ if(is_lifted[i + 1]){ @@ -1421,15 +1425,21 @@ static inline int rational_reconstruction_param( /** CRT DONE **/ /** RATIONAL RECONSTRUCTIONS **/ - mpz_sqrt(recdata->N, modulus); - mpz_fdiv_q(recdata->D, modulus, recdata->N); + /* keep RR_MARGIN bits of slack (2ND <= modulus / 2^RR_MARGIN) so that + spurious reconstructions become unlikely */ + mpz_t mm; + mpz_init(mm); + mpz_fdiv_q_2exp(mm, modulus, RR_MARGIN); + mpz_sqrt(recdata->N, mm); + mpz_fdiv_q(recdata->D, mm, recdata->N); mpz_fdiv_q_2exp(recdata->D, recdata->D, 1); - mpz_root(recdata->D1, modulus, 4); - mpz_fdiv_q(recdata->N1, modulus, recdata->D1); + mpz_root(recdata->D1, mm, 4); + mpz_fdiv_q(recdata->N1, mm, recdata->D1); mpz_fdiv_q_2exp(recdata->D1, recdata->D1, 1); - mpz_root(recdata->D2, modulus, 16); - mpz_fdiv_q(recdata->N2, modulus, recdata->D2); + mpz_root(recdata->D2, mm, 16); + mpz_fdiv_q(recdata->N2, mm, recdata->D2); mpz_fdiv_q_2exp(recdata->D2, recdata->D2, 1); + mpz_clear(mm); mpz_set_ui(guessed_den, 1); mpz_set_ui(rnum, 0); From 0cda8a0721623cbdf2b94a8fbb048c865ade228a Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Sun, 25 Oct 2026 16:24:56 +0100 Subject: [PATCH 5/9] Fall back to bounds from the guessed denominator for parametrizations When lifting parametrizations over QQ, each coefficient of the eliminating polynomial and of the coordinates is now reconstructed by trying in turn the fixed bounds used so far (D = M^(1/16) for the eliminating polynomial; D = M^(1/32), balanced bounds and N = M^(1/16) for the coordinates, with M = modulus / 2^RR_MARGIN), then the bounds set from the denominator accumulated so far, g: D = sqrt(M / g) and N = M / D, as for the lifting of Groebner bases. CRT is applied only once per coefficient. Previously, each fixed bound had to succeed on a whole coordinate. Using only the bounds set from the guessed denominator increases the number of primes by up to 1.9x on several families (eco, noon, henrion, cyclic), whose numerators are much larger than their denominators; the fixed bounds are thus tried first. Outputs are unchanged on 35 benchmark systems and the number of primes never increases. It decreases on KLYZ22 (279 to 225, 18% faster), whose denominators are large, and on the smallest Katsura systems, where some coordinates were lifted later than the eliminating polynomial because no single fixed bound suited all their coefficients (kat8: 22 to 16, kat9: 40 to 31, about 20% faster). From kat10 to kat12 all coordinates are lifted together with the eliminating polynomial, which is the bottleneck, so the number of primes and timings are unchanged (kat12: 323 primes, 325s vs 321s); the same holds for eco8 to eco12, where the fixed bounds always succeed and the fallback is never used. Co-Authored-By: Claude Opus 5.5 --- src/msolve/msolve.c | 129 +++++++++++++++++++++++++------------------- 1 file changed, 75 insertions(+), 54 deletions(-) diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index 32732d38..6308f619 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -1101,22 +1101,44 @@ static inline void param_images_crt_coeff(const param_images_t im, **/ #define RR 1 +/* rational reconstruction of u * lcm modulo modulus, trying in turn the + fixed bounds (bN[k], bD[k]) for k < nb, then the bounds set from the + guessed denominator lcm (D = sqrt(M / lcm), N = M / D); the fixed bounds + suit numerators much larger than denominators, the latter suit large + denominators */ +static inline int param_ratrecon_coeff(mpz_t num, mpz_t den, mpz_t u, + mpz_t modulus, mpz_t lcm, + rrec_data_t rd, mpz_t *bN, mpz_t *bD, + const int nb){ + for (int k = 0; k < nb; k++) { + mpz_set(rd->N, bN[k]); + mpz_set(rd->D, bD[k]); + if (lift_ratreconwden(num, den, u, modulus, lcm, rd)) { + return 1; + } + } + set_recdata_guess(rd, modulus, lcm); + return lift_ratreconwden(num, den, u, modulus, lcm, rd); +} + static inline int new_rational_reconstruction_mpz_ptr_with_denom_param( mpz_t *recons, mpz_t denominator, const param_images_t im, int64_t start, param_crt_struct *pc, deg_t len, mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_t *tmp_num, - mpz_t *tmp_den, mpz_t lcm, mpz_t gnum, mpz_t guessed_den, rrec_data_t rdata) { + mpz_t *tmp_den, mpz_t lcm, mpz_t gnum, mpz_t guessed_den, rrec_data_t rdata, + mpz_t *bN, mpz_t *bD, const int nb) { mpz_set_ui(gnum, 1); int b = 0; - //for (deg_t i = 0; i < len; i++) { + /* lcm is the denominator accumulated so far, used as a guess for the + next coefficient */ for (deg_t i = len - 1; i >= 0; i--) { param_images_crt_coeff(im, start + i, pc); - b = lift_ratreconwden(tmp_num[i], rden, pc->u, modulus, lcm, - rdata); + b = param_ratrecon_coeff(tmp_num[i], rden, pc->u, modulus, lcm, + rdata, bN, bD, nb); if (b == 0) { return 0; @@ -1150,7 +1172,7 @@ static inline int rational_reconstruction_mpz_ptr_with_denom_elim( param_crt_struct *pc, deg_t len, mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_t *tmp_num, mpz_t *tmp_den, int16_t *check_lift, mpz_t lcm, mpz_t gnum, - mpz_t guessed_den, rrec_data_t rdata) { + mpz_t guessed_den, rrec_data_t rdata, mpz_t *bN, mpz_t *bD, const int nb) { mpz_set_ui(gnum, 1); @@ -1169,7 +1191,8 @@ static inline int rational_reconstruction_mpz_ptr_with_denom_elim( for(deg_t i = len - 1; i >= 0; i--){ if(check_lift[i] <= 1){ param_images_crt_coeff(im, start + i, pc); - b = lift_ratreconwden(tmp_num[i], rden, pc->u, modulus, lcm, rdata); + b = param_ratrecon_coeff(tmp_num[i], rden, pc->u, modulus, lcm, + rdata, bN, bD, nb); if(b == 0){ check_lift[i] = 0; return 0; @@ -1209,12 +1232,12 @@ static inline int rational_reconstruction_upoly_with_denom_elim( mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_upoly_t tmp_num, mpz_upoly_t tmp_den, int16_t *check_lift, mpz_t lcm, mpz_t guessed_num, - mpz_t guessed_den, rrec_data_t rdata) { + mpz_t guessed_den, rrec_data_t rdata, mpz_t *bN, mpz_t *bD, const int nb) { return rational_reconstruction_mpz_ptr_with_denom_elim( recons->coeffs, denominator, im, start, pc, len, modulus, maxrec, coef, rnum, rden, tmp_num->coeffs, tmp_den->coeffs, check_lift, lcm, guessed_num, - guessed_den, rdata); + guessed_den, rdata, bN, bD, nb); } static inline int rational_reconstruction_upoly_with_denom_param( @@ -1222,13 +1245,13 @@ static inline int rational_reconstruction_upoly_with_denom_param( int64_t start, param_crt_struct *pc, long len, mpz_t modulus, deg_t *maxrec, mpq_t *coef, mpz_t rnum, mpz_t rden, mpz_upoly_t tmp_num, mpz_upoly_t tmp_den, mpz_t lcm, mpz_t guessed_num, - mpz_t guessed_den, rrec_data_t rdata) { + mpz_t guessed_den, rrec_data_t rdata, mpz_t *bN, mpz_t *bD, const int nb) { return new_rational_reconstruction_mpz_ptr_with_denom_param( recons->coeffs, denominator, im, start, pc, len, modulus, maxrec, coef, rnum, rden, tmp_num->coeffs, tmp_den->coeffs, lcm, guessed_num, - guessed_den, rdata); + guessed_den, rdata, bN, bD, nb); } @@ -1239,7 +1262,7 @@ static inline void lift_coordinate(mpz_param_t mpz_param, mpz_t modulus, mpz_t prod_crt, int32_t prime, mpq_t *coef, mpz_t rnum, mpz_t rden, rrec_data_t recdata, mpz_t guessed_num, mpz_t guessed_den, deg_t *maxrec, int *is_lifted, mpz_t denominator, mpz_t lcm, mpq_t c, const int idx, - const int info_level){ + mpz_t *bN, mpz_t *bD, const int nb, const int info_level){ if(mpz_sizeinbase(guessed_den, 2) <= mpz_sizeinbase(modulus, 2)){ mpz_set(lcm, guessed_den); @@ -1255,7 +1278,7 @@ static inline void lift_coordinate(mpz_param_t mpz_param, images->elen + idx * images->clen, pc, MIN((deg_t)nmod_param->coords[idx]->length, images->clen), modulus, maxrec, coef, rnum, rden, - numer, denom, lcm, guessed_num, guessed_den, recdata); + numer, denom, lcm, guessed_num, guessed_den, recdata, bN, bD, nb); if(b){ is_lifted[idx + 1] = 1; mpz_set(mpq_denref(c), denominator); @@ -1287,19 +1310,47 @@ static inline int lift_parametrization(mpz_param_t mpz_param, param_crt_t pc; param_crt_init(pc, images); + /* fixed bounds tried first, with M = modulus / 2^RR_MARGIN so that + RR_MARGIN bits of slack are kept: D = M^(1/16) for elim; D = M^(1/32), + balanced bounds, then N = M^(1/16) for the coordinates; the bounds + set from the guessed denominator are tried last */ + mpz_t bN[4], bD[4], mm; + for(int k = 0; k < 4; k++){ + mpz_init(bN[k]); + mpz_init(bD[k]); + } + mpz_init(mm); + mpz_fdiv_q_2exp(mm, modulus, RR_MARGIN); + /* elim */ + mpz_root(bD[0], mm, 16); + mpz_fdiv_q(bN[0], mm, bD[0]); + mpz_fdiv_q_2exp(bD[0], bD[0], 1); + /* coordinates */ + mpz_root(bD[1], mm, 32); + mpz_fdiv_q(bN[1], mm, bD[1]); + mpz_fdiv_q_2exp(bD[1], bD[1], 1); + mpz_sqrt(bD[2], mm); + mpz_fdiv_q(bN[2], mm, bD[2]); + mpz_fdiv_q_2exp(bD[2], bD[2], 1); + mpz_root(bN[3], mm, 16); + mpz_fdiv_q(bD[3], mm, bN[3]); + mpz_fdiv_q_2exp(bN[3], bN[3], 1); + mpz_clear(mm); + int b = 0; if (is_lifted[0] == 0) { - mpz_swap(recdata->N, recdata->N2); - mpz_swap(recdata->D, recdata->D2); - mpz_swap(recdata->N1, recdata->N2); - mpz_swap(recdata->D1, recdata->D2); b = rational_reconstruction_upoly_with_denom_elim( mpz_param->elim, denominator, images, 0, pc, nmod_param->elim->length, modulus, maxrec, coef, rnum, rden, numer, - denom, check_lift, lcm, guessed_num, guessed_den, recdata); + denom, check_lift, lcm, guessed_num, guessed_den, recdata, + bN, bD, 1); if (b == 0) { is_lifted[0] = 0; + for(int k = 0; k < 4; k++){ + mpz_clear(bN[k]); + mpz_clear(bD[k]); + } param_crt_clear(pc); return 0; } @@ -1322,43 +1373,19 @@ static inline int lift_parametrization(mpz_param_t mpz_param, mpz_set(guessed_den, lc); - /* bounds tried in turn for each coordinate: D = M^(1/32), then balanced - bounds, then N = M^(1/16), with M = modulus / 2^RR_MARGIN so that - RR_MARGIN bits of slack are kept */ - mpz_t bN[3], bD[3], mm; - for(int k = 0; k < 3; k++){ - mpz_init(bN[k]); - mpz_init(bD[k]); - } - mpz_init(mm); - mpz_fdiv_q_2exp(mm, modulus, RR_MARGIN); - mpz_root(bD[0], mm, 32); - mpz_fdiv_q(bN[0], mm, bD[0]); - mpz_fdiv_q_2exp(bD[0], bD[0], 1); - mpz_sqrt(bD[1], mm); - mpz_fdiv_q(bN[1], mm, bD[1]); - mpz_fdiv_q_2exp(bD[1], bD[1], 1); - mpz_root(bN[2], mm, 16); - mpz_fdiv_q(bD[2], mm, bN[2]); - mpz_fdiv_q_2exp(bN[2], bN[2], 1); - mpz_clear(mm); - for(int i = 0; i < nc; i++){ if(is_lifted[i + 1]){ continue; } - for(int k = 0; k < 3 && is_lifted[i + 1] == 0; k++){ - mpz_set(recdata->N, bN[k]); - mpz_set(recdata->D, bD[k]); - lift_coordinate(mpz_param, images, pc, - nmod_param, numer, denom, check_lift, - modulus, prod_crt, prime, coef, rnum, rden, - recdata, guessed_num, guessed_den, maxrec, - is_lifted, denominator, lcm, c, i, info_level); - } + lift_coordinate(mpz_param, images, pc, + nmod_param, numer, denom, check_lift, + modulus, prod_crt, prime, coef, rnum, rden, + recdata, guessed_num, guessed_den, maxrec, + is_lifted, denominator, lcm, c, i, bN + 1, bD + 1, 3, + info_level); } - for(int k = 0; k < 3; k++){ + for(int k = 0; k < 4; k++){ mpz_clear(bN[k]); mpz_clear(bD[k]); } @@ -1433,12 +1460,6 @@ static inline int rational_reconstruction_param( mpz_sqrt(recdata->N, mm); mpz_fdiv_q(recdata->D, mm, recdata->N); mpz_fdiv_q_2exp(recdata->D, recdata->D, 1); - mpz_root(recdata->D1, mm, 4); - mpz_fdiv_q(recdata->N1, mm, recdata->D1); - mpz_fdiv_q_2exp(recdata->D1, recdata->D1, 1); - mpz_root(recdata->D2, mm, 16); - mpz_fdiv_q(recdata->N2, mm, recdata->D2); - mpz_fdiv_q_2exp(recdata->D2, recdata->D2, 1); mpz_clear(mm); mpz_set_ui(guessed_den, 1); From 41ac4220358162fbc07000deaef2eb8b72b4337e Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Sun, 25 Oct 2026 17:28:00 +0100 Subject: [PATCH 6/9] Fix a leak in mulmat_reconstruct and remove dead code from lifting.c mulmat_reconstruct allocated ncols integers for the row denominators but cleared only nrows of them, leaking ncols - nrows integers at each call when multiplication matrices are lifted (-L 1). Also remove twelve static functions of lifting.c which are no longer called: the former lifting of multiplication matrices (crt_mpz_matfglm_*, mpz_matfglm_*, mpq_matfglm_*, crt_lift_mat, rat_recon_matfglm, check_matrix_and_linear_forms), replaced by mulmat_reconstruct, as well as old_rat_recon_trace_det and check_linear_forms. The corresponding types are kept in msolve-data.h, which is installed. Outputs are unchanged, with and without -L 1. Co-Authored-By: Claude Opus 5.5 --- src/msolve/lifting.c | 339 +------------------------------------------ 1 file changed, 1 insertion(+), 338 deletions(-) diff --git a/src/msolve/lifting.c b/src/msolve/lifting.c index d276c3ca..11f8fcfe 100644 --- a/src/msolve/lifting.c +++ b/src/msolve/lifting.c @@ -79,175 +79,6 @@ static int lift_ratreconwden(mpz_t n, mpz_t d, mpz_t u, const mpz_t mod, return flint_rrec(n, d, rd->r1, mod, rd->N, rd->D); } -/* Initialization of sparse fglm matrix for crt */ -static inline void crt_mpz_matfglm_initset(crt_mpz_matfglm_t crt_mat, - sp_matfglm_t *mod_mat){ - crt_mat->ncols = mod_mat->ncols; - crt_mat->nrows = mod_mat->nrows; - uint64_t sz = crt_mat->nrows * crt_mat->ncols; - crt_mat->dense_mat = (mpz_t *)malloc(sz * sizeof(mpz_t)); - for (uint64_t i = 0; i < sz; i++){ - mpz_init_set_ui(crt_mat->dense_mat[i], mod_mat->dense_mat[i]); - } - long diff = crt_mat->ncols - crt_mat->nrows; - crt_mat->triv_idx = malloc(diff * sizeof(uint32_t)); - crt_mat->triv_pos = malloc(diff * sizeof(uint32_t)); - crt_mat->dense_idx = malloc(crt_mat->nrows * sizeof(uint32_t)); - crt_mat->dst = malloc(crt_mat->nrows * sizeof(uint32_t)); - - for (long i = 0; i < diff; i++){ - crt_mat->triv_idx[i] = mod_mat->triv_idx[i]; - crt_mat->triv_pos[i] = mod_mat->triv_pos[i]; - } - for (long i = 0; i < crt_mat->nrows; i++){ - crt_mat->dense_idx[i] = mod_mat->dense_idx[i]; - crt_mat->dst[i] = mod_mat->dst[i]; - } -} - -static inline void crt_mpz_matfglm_clear(crt_mpz_matfglm_t crt_mat){ - uint64_t sz = crt_mat->nrows * crt_mat->ncols; - - for (uint64_t i = 0; i < sz; i++){ - mpz_clear(crt_mat->dense_mat[i]); - } - free(crt_mat->dense_mat); - - free(crt_mat->triv_idx); - free(crt_mat->triv_pos); - free(crt_mat->dense_idx); - free(crt_mat->dst); -} - -/* Initialization of sparse mpz fglm matrix */ -static inline void mpz_matfglm_initset(mpz_matfglm_t mpz_mat, - sp_matfglm_t *mod_mat) { - mpz_mat->ncols = mod_mat->ncols; - mpz_mat->nrows = mod_mat->nrows; - mpz_mat->dense_mat = malloc(mpz_mat->ncols * mpz_mat->nrows * sizeof(mpz_t)); - uint64_t sz = mpz_mat->nrows * mpz_mat->ncols; - for (uint64_t i = 0; i < sz; i++) { - mpz_init(mpz_mat->dense_mat[i]); - } - mpz_mat->denoms = malloc(mpz_mat->nrows * sizeof(mpz_t)); - for (uint64_t i = 0; i < mpz_mat->nrows; i++) { - mpz_init(mpz_mat->denoms[i]); - } - long diff = mpz_mat->ncols - mpz_mat->nrows; - mpz_mat->triv_idx = calloc(diff, sizeof(uint32_t)); - mpz_mat->triv_pos = calloc(diff, sizeof(uint32_t)); - mpz_mat->dense_idx = calloc(mpz_mat->nrows, sizeof(uint32_t)); - mpz_mat->dst = calloc(mpz_mat->nrows, sizeof(uint32_t)); - - for (long i = 0; i < diff; i++) { - mpz_mat->triv_idx[i] = mod_mat->triv_idx[i]; - mpz_mat->triv_pos[i] = mod_mat->triv_pos[i]; - } - for (long i = 0; i < mpz_mat->nrows; i++) { - mpz_mat->dense_idx[i] = mod_mat->dense_idx[i]; - mpz_mat->dst[i] = mod_mat->dst[i]; - } -} - -static inline void mpz_matfglm_clear(mpz_matfglm_t mpz_mat) { - - uint64_t sz = mpz_mat->nrows * mpz_mat->ncols; - for (uint64_t i = 0; i < sz; i++) { - mpz_clear(mpz_mat->dense_mat[i]); - } - free(mpz_mat->dense_mat); - for (uint64_t i = 0; i < mpz_mat->nrows; i++) { - mpz_clear(mpz_mat->denoms[i]); - } - free(mpz_mat->denoms); - - free(mpz_mat->triv_idx); - free(mpz_mat->triv_pos); - free(mpz_mat->dense_idx); - free(mpz_mat->dst); -} - -/* Initialization of sparse fglm matrix for rational reconstruction */ -static inline void mpq_matfglm_initset(mpq_matfglm_t mpq_mat, - sp_matfglm_t *mod_mat) { - mpq_mat->ncols = mod_mat->ncols; - mpq_mat->nrows = mod_mat->nrows; - mpq_mat->dense_mat = - malloc(2 * mpq_mat->ncols * mpq_mat->nrows * sizeof(mpz_t)); - mpq_mat->denoms = malloc(mpq_mat->nrows * sizeof(mpz_t)); - for (uint32_t i = 0; i < mpq_mat->nrows; i++) { - mpz_init(mpq_mat->denoms[i]); - } - uint64_t nc = 2 * mpq_mat->ncols; - - for (uint32_t i = 0; i < mpq_mat->nrows; i++) { - uint64_t c = 2 * i * mpq_mat->ncols; - for (uint32_t j = 0; j < nc; j++) { - mpz_init(mpq_mat->dense_mat[c + j]); - j++; - mpz_init(mpq_mat->dense_mat[c + j]); - } - } - long diff = mpq_mat->ncols - mpq_mat->nrows; - mpq_mat->triv_idx = calloc(diff, sizeof(uint32_t)); - mpq_mat->triv_pos = calloc(diff, sizeof(uint32_t)); - mpq_mat->dense_idx = calloc(mpq_mat->nrows, sizeof(uint32_t)); - mpq_mat->dst = calloc(mpq_mat->nrows, sizeof(uint32_t)); - - for (long i = 0; i < diff; i++) { - mpq_mat->triv_idx[i] = mod_mat->triv_idx[i]; - mpq_mat->triv_pos[i] = mod_mat->triv_pos[i]; - } - for (long i = 0; i < mpq_mat->nrows; i++) { - mpq_mat->dense_idx[i] = mod_mat->dense_idx[i]; - mpq_mat->dst[i] = mod_mat->dst[i]; - } -} - -static inline void mpq_matfglm_clear(mpq_matfglm_t mpq_mat) { - /* entries of mpq_mat->denoms are cleared previously when building mpz_mat */ - for (uint32_t i = 0; i < mpq_mat->nrows; i++) { - mpz_clear(mpq_mat->denoms[i]); - } - free(mpq_mat->denoms); - - uint64_t nc = 2 * mpq_mat->ncols; - for (uint32_t i = 0; i < mpq_mat->nrows; i++) { - uint64_t c = 2 * i * mpq_mat->ncols; - for (uint32_t j = 0; j < nc; j++) { - mpz_clear(mpq_mat->dense_mat[c + j]); - } - } - free(mpq_mat->dense_mat); - - free(mpq_mat->triv_idx); - free(mpq_mat->triv_pos); - free(mpq_mat->dense_idx); - free(mpq_mat->dst); -} - -static inline void mpq_matfglm_partial_clear(mpq_matfglm_t mpq_mat) { - /* entries of mpq_mat->denoms are cleared previously when building mpz_mat */ - for (uint32_t i = 0; i < mpq_mat->nrows; i++) { - mpz_clear(mpq_mat->denoms[i]); - } - free(mpq_mat->denoms); - - uint64_t nc = mpq_mat->ncols; - for (uint32_t i = 0; i < mpq_mat->nrows; i++) { - uint64_t c = 2 * i * mpq_mat->ncols; - for (uint32_t j = 0; j < nc; j++) { - mpz_clear(mpq_mat->dense_mat[c + 2 * j]); - } - } - free(mpq_mat->dense_mat); - - free(mpq_mat->triv_idx); - free(mpq_mat->triv_pos); - free(mpq_mat->dense_idx); - free(mpq_mat->dst); -} - static inline void copy_modular_matrix(trace_det_fglm_mat_t trace_det, sp_matfglm_t **mod_mat, uint32_t newalloc, mp_limb_t prime){ deg_t nrows = (*mod_mat)->nrows; @@ -503,16 +334,6 @@ static inline void crt_lift_trace_det(trace_det_fglm_mat_t trace_det, } } -static inline void crt_lift_mat(crt_mpz_matfglm_t mat, sp_matfglm_t *mod_mat, - mpz_t modulus, mpz_t prod_crt, - const int32_t prime, mpz_t tmp, - const int32_t nrows, const int nthrds) { - /*assumes prod_crt = modulus * prime */ - const uint64_t sz = mat->nrows * mat->ncols; - crt_lift_dense_rows(mat->dense_mat, mod_mat->dense_mat, nrows * mat->ncols, - sz, modulus, prod_crt, prime, tmp, nthrds); -} - static inline void build_linear_forms(mpz_t *mpz_linear_forms, mpz_t *mpq_linear_forms, uint32_t sz, int nlins) { @@ -552,33 +373,6 @@ static inline void build_linear_forms(mpz_t *mpz_linear_forms, } -static inline int old_rat_recon_trace_det(trace_det_fglm_mat_t trace_det, - rrec_data_t recdata, mpz_t modulus, - mpz_t rnum, mpz_t rden, mpz_t gden){ - if(trace_det->done_trace < 2){ - int b = ratrecon(rnum, rden, trace_det->trace_crt, modulus, recdata); - if(b == 1){ - mpz_set(trace_det->trace_num, rnum); - mpz_set(trace_det->trace_den, rden); - } - else - { - return 0; - } - } - if(trace_det->done_det < 2){ - int b = ratrecon(rnum, rden, trace_det->det_crt, modulus, recdata); - if(b == 1){ - mpz_set(trace_det->det_num, rnum); - mpz_set(trace_det->det_den, rden); - } - else{ - return 0; - } - } - return 1; -} - static inline int rat_recon_array(mpz_t *res, mpz_t *crt, int32_t sz, mpz_t modulus, rrec_data_t rdata) { for (int i = 0; i < sz; i++) { @@ -800,26 +594,6 @@ static inline int check_trace_det_data(trace_det_fglm_mat_t trace_det, return 0; } -static inline int check_linear_forms(trace_det_fglm_mat_t trace_det, - uint32_t *modular_linear_forms, - int32_t prime){ - int32_t nlins = trace_det->nlins; - int32_t nv = trace_det ->nv; - nvars_t sz = nv + 1; - for(int32_t i = 0; i < nlins; i++){ - nvars_t n = i *sz; - nvars_t n2 = i * (sz + 1); - for(int32_t j = 0; j < sz; j++){ - int b = check_lifted_coeff(trace_det->mpz_linear_forms[n2 + j], - trace_det->mpz_linear_forms[n2+sz], - modular_linear_forms[n], - prime); - if(!b) return 0; - } - } - return 1; -} - static inline void check_matrix(trace_det_fglm_mat_t trace_det, sp_matfglm_t *mod_mat, int32_t prime){ @@ -904,7 +678,7 @@ static inline void mulmat_reconstruct(trace_det_fglm_mat_t trace_det, } } - for(int32_t i = 0; i < nrows; i++){ + for(int32_t i = 0; i < ncols; i++){ mpz_clear(rowdenoms[i]); } free(rowdenoms); @@ -921,117 +695,6 @@ static inline void mulmat_reconstruct(trace_det_fglm_mat_t trace_det, } } -static inline void -rat_recon_matfglm(mpq_matfglm_t mpq_mat, crt_mpz_matfglm_t crt_mat, - mpz_t modulus, rrec_data_t rdata, - mpz_t rnum, mpz_t rden, mpz_t *guessed_den, deg_t *matrec, int *mat_lifted) { - const uint32_t nrows = crt_mat->nrows; - const uint32_t ncols = crt_mat->ncols; - - mpz_t lcm, coef; - mpz_init(lcm); - mpz_init(coef); - mpz_t gden; - mpz_init(gden); - if((0==1) && (*matrec)>1){ - mpz_set(gden, *guessed_den); - } - else{ - mpz_set_ui(gden, 1); - } - for (uint32_t i = *matrec; i < nrows; i++) { - uint64_t c = i * ncols; - mpz_set_ui(lcm, 1); - - for (uint32_t j = 0; j < ncols; j++) { - int b = 0; - mpz_mul(coef, crt_mat->dense_mat[c + j], lcm); - mpz_mod(coef, coef, modulus); -// b = ratrecon(rnum, rden, coef, modulus, rdata); - b = ratreconwden(rnum, rden, coef, modulus, gden, rdata); - if (b == 1) { - mpz_set(mpq_mat->dense_mat[2 * c + 2 * j], rnum); - mpz_set(mpq_mat->dense_mat[2 * c + 2 * j + 1], rden); - mpz_mul(mpq_mat->dense_mat[2 * c + 2 * j + 1], - mpq_mat->dense_mat[2 * c + 2 * j + 1], lcm); - mpz_mul(mpq_mat->dense_mat[2 * c + 2 * j + 1], - mpq_mat->dense_mat[2 * c + 2 * j + 1], gden); - mpz_lcm(lcm, lcm, rden); - } else { - mpz_clear(gden); - mpz_clear(lcm); - mpz_clear(coef); - return ; - } - } - mpz_set_ui(mpq_mat->denoms[i], 1); - for (uint32_t j = 0; j < ncols; j++) { - mpz_lcm(mpq_mat->denoms[i], mpq_mat->denoms[i], - mpq_mat->dense_mat[2 * c + 2 * j + 1]); - } - - for(uint32_t j = 0; j < ncols; j++){ - mpz_mul(mpq_mat->dense_mat[2 * c+ 2 * j], mpq_mat->dense_mat[2 * c+ 2 * j], - mpq_mat->denoms[i]); - mpz_divexact(mpq_mat->dense_mat[2 * c + 2 * j], mpq_mat->dense_mat[2 * c + 2 * j], - mpq_mat->dense_mat[2 * c + 2 * j + 1]); - } - (*matrec)++; - } - mpz_clear(coef); - mpz_clear(lcm); - mpz_clear(gden); - - *mat_lifted = 1; -} - - -static inline void check_matrix_and_linear_forms(mpq_matfglm_t mpq_mat, sp_matfglm_t *mod_mat, - mpz_t *mpz_lin, uint32_t *mod_lin, const int nlins, const nvars_t nv, - int *mat_lifted, int* lin_lifted, deg_t *oldmatrec_checked, deg_t *matrec_checked, const uint32_t prime){ - *oldmatrec_checked = *matrec_checked; - /* checks lifted linear forms */ - nvars_t sz = nv + 1; - uint32_t cand; - if(*lin_lifted < 2){ - for(int i = 0; i < nlins; i++){ - nvars_t n = i * sz; - nvars_t n2 = i * (sz + 1); - uint64_t inv = mpz_fdiv_ui(mpz_lin[n2 + sz], prime); - inv = mod_p_inverse_32(inv, prime); - for(nvars_t j = 0; j < sz; j++){ - cand = (((uint64_t)mpz_fdiv_ui(mpz_lin[n2+j], prime)) * inv) % prime; - if(cand != mod_lin[n + j]){ - *lin_lifted = 0; - *matrec_checked = 0; - return; - } - } - } - *lin_lifted = 2; - } - - /*checks multiplication matrix */ - mod_mat->charac = prime; - - for(int32_t i = (*oldmatrec_checked); i < mpq_mat->nrows; i++){ - uint64_t lc = mpz_fdiv_ui(mpq_mat->denoms[i], prime); - lc = mod_p_inverse_32(lc, prime); - uint32_t nc = i*mpq_mat->ncols; - for(uint32_t j = 0 ; j < mpq_mat->ncols; j++){ - uint32_t mod = mpz_fdiv_ui(mpq_mat->dense_mat[2*(nc+j)], prime); - if(mod_mat->dense_mat[nc+j] != ( ((uint64_t)mod) * lc )% prime){ - *mat_lifted = 0; - *matrec_checked = MAX(0, i); - return; - } - } - *matrec_checked = i + 1; - } - *mat_lifted = 2; - -} - void initialize_rrec_data(rrec_data_t recdata) { mpz_init(recdata->r0); mpz_set_ui(recdata->r0, 0); From 6986d4ad4e667ed6b73726ce7fd2ab75c74483ad Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Sun, 25 Oct 2026 17:14:49 +0100 Subject: [PATCH 7/9] Reconstruct witnesses of eliminating polynomials with a common denominator The coefficients of the monic eliminating polynomial share a common denominator (the leading coefficient of the primitive integer polynomial). On dense random, MinRank and critical point systems, this denominator is about as large as the numerators, so reconstructing a single witness coefficient (trace or determinant) needs about twice the number of bits of these numerators, and the eliminating polynomial cannot be attempted before. Four witness coefficients (the trace and determinant ones, and the first nonzero ones from degrees nsols/2 and nsols/4) are now also reconstructed simultaneously: a short vector (d, n_1, .., n_4) with n_j = d u_j mod m is computed by LLL (FLINT's fmpz_lll), which needs about bits(n) + bits(d) / 4 bits of modulus. The candidate is accepted when its entries are below 2^(4 bits(m) / 5 - RR_MARGIN) and when it matches the images of the next prime; d is then used as a guess for the eliminating polynomial, which can then be reconstructed with its usual bounds. As LLL is more expensive than rational reconstruction, it is tried only while no single witness is known, at 16 primes and then each time the number of primes grows by 5%. Outputs are unchanged on 35 benchmark systems. On the 18 dense random, MinRank, critical point (cp_d, crit) systems, the number of primes decreases by 32 to 38%, and the total time by 20 to 39% on those taking more than 0.3s (geometric mean speedup 1.44; random2_9_9: 2643 to 1651 primes, 81.7s to 56.1s; random3_6_6: 763 to 495 primes, 33.5s to 21.7s). On the other families the number of primes is unchanged and timings are within noise (+-3%), the LLL attempts made before a single witness is known costing a few hundredths of a second. Co-Authored-By: Claude Opus 5.5 --- src/msolve/lifting.c | 2 + src/msolve/msolve.c | 172 ++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 171 insertions(+), 3 deletions(-) diff --git a/src/msolve/lifting.c b/src/msolve/lifting.c index 11f8fcfe..c570c7d8 100644 --- a/src/msolve/lifting.c +++ b/src/msolve/lifting.c @@ -26,6 +26,8 @@ #include "msolve-data.h" #include #include +#include +#include /* rational reconstruction of u mod m with bounds N, D through FLINT (half-gcd based) instead of the quadratic Euclidean algorithm; diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index 6308f619..e93f825d 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -1094,6 +1094,140 @@ static inline void param_images_crt_coeff(const param_images_t im, fmpz_get_mpz(pc->u, pc->y); } +/* Witnesses of the eliminating polynomial reconstructed simultaneously: + its coefficients share a common denominator (the leading coefficient of + the primitive integer polynomial), which can be much larger than their + numerators. Reconstructing k coefficients u_1, .., u_k with a common + denominator d, i.e. finding a short vector (d, n_1, .., n_k) with + n_j = d u_j mod m by LLL, needs about bits(n) + bits(d) / k bits of + modulus instead of bits(n) + bits(d) for a single coefficient. The + denominator found is used as a guess for the eliminating polynomial. + As LLL is more expensive than rational reconstruction, it is tried at + a geometric schedule: each time the number of primes grows by + VEC_WIT_GROWTH percent. */ +#define VEC_WIT_K 4 +#define VEC_WIT_START 16 +#define VEC_WIT_GROWTH 5 + +typedef struct{ + int k; /* number of witnesses, 0 if disabled */ + deg_t idx[VEC_WIT_K]; /* indices of the witness coefficients in elim */ + int state; /* 0: nothing, 1: candidate to be checked with + the next prime, 2: candidate checked */ + int32_t next; /* number of primes of the next attempt */ + mpz_t den; /* common denominator */ + mpz_t num[VEC_WIT_K]; /* numerators */ +} vec_wit_struct; + +typedef vec_wit_struct vec_wit_t[1]; + +/* witnesses: the trace and det witnesses, then the first nonzero + coefficients of elim from degrees nsols/2 and nsols/4 upwards */ +static inline void vec_wit_init(vec_wit_t vw, const nmod_poly_t elim, + const deg_t tridx, const deg_t detidx){ + deg_t nsols = elim->length - 1; + mpz_init(vw->den); + for (int j = 0; j < VEC_WIT_K; j++) { + mpz_init(vw->num[j]); + } + vw->state = 0; + vw->next = VEC_WIT_START; + vw->k = 0; + if (nsols < 4 * VEC_WIT_K) { + return; + } + deg_t cand[VEC_WIT_K] = {tridx, detidx, nsols / 2, nsols / 4}; + for (int j = 0; j < VEC_WIT_K; j++) { + deg_t i = cand[j]; + int again = 1; + while (again && i < nsols) { + again = (elim->coeffs[i] == 0); + for (int l = 0; l < vw->k; l++) { + again = again || (vw->idx[l] == i); + } + if (again) { + i++; + } + } + if (i < nsols) { + vw->idx[vw->k] = i; + vw->k++; + } + } + if (vw->k < 2) { + vw->k = 0; + } +} + +static inline void vec_wit_clear(vec_wit_t vw){ + mpz_clear(vw->den); + for (int j = 0; j < VEC_WIT_K; j++) { + mpz_clear(vw->num[j]); + } +} + +/* a candidate is accepted when the max-norm of (d, n_1, .., n_k) is at most + 2^(k bits(m) / (k + 1) - RR_MARGIN), the analogue of the bounds with + RR_MARGIN bits of slack used for rational reconstruction */ +static inline void vec_wit_try(vec_wit_t vw, const param_images_t im, + const mpz_t modulus){ + const int k = vw->k; + param_crt_t pc; + param_crt_init(pc, im); + fmpz_mat_t B; + fmpz_mat_init(B, k + 1, k + 1); + fmpz_one(fmpz_mat_entry(B, 0, 0)); + for (int j = 0; j < k; j++) { + param_images_crt_coeff(im, vw->idx[j], pc); + fmpz_set_mpz(fmpz_mat_entry(B, 0, j + 1), pc->u); + fmpz_set_mpz(fmpz_mat_entry(B, j + 1, j + 1), modulus); + } + param_crt_clear(pc); + fmpz_lll_t fl; + fmpz_lll_context_init_default(fl); + fmpz_lll(B, NULL, fl); + + slong bits = (mpz_sizeinbase(modulus, 2) * k) / (k + 1) - RR_MARGIN; + int ok = !fmpz_is_zero(fmpz_mat_entry(B, 0, 0)) && bits > 0; + for (int j = 0; j <= k && ok; j++) { + ok = (fmpz_bits(fmpz_mat_entry(B, 0, j)) <= bits); + } + if (ok) { + int neg = fmpz_sgn(fmpz_mat_entry(B, 0, 0)) < 0; + fmpz_get_mpz(vw->den, fmpz_mat_entry(B, 0, 0)); + for (int j = 0; j < k; j++) { + fmpz_get_mpz(vw->num[j], fmpz_mat_entry(B, 0, j + 1)); + } + if (neg) { + mpz_neg(vw->den, vw->den); + for (int j = 0; j < k; j++) { + mpz_neg(vw->num[j], vw->num[j]); + } + } + vw->state = 1; + } + fmpz_mat_clear(B); +} + +/* checks the candidate with the modular images of elim for a new prime */ +static inline void vec_wit_check(vec_wit_t vw, const nmod_poly_t elim, + const uint32_t prime){ + uint64_t d = mpz_fdiv_ui(vw->den, prime); + if (d == 0) { + vw->state = 0; + return; + } + uint64_t dinv = mod_p_inverse_32(d, prime); + for (int j = 0; j < vw->k; j++) { + uint64_t c = mpz_fdiv_ui(vw->num[j], prime); + if ((c * dinv) % prime != elim->coeffs[vw->idx[j]]) { + vw->state = 0; + return; + } + } + vw->state = 2; +} + /** la sortie est recons / denominator @@ -1402,7 +1536,7 @@ returns 0 if rational reconstruction failed **/ static inline int rational_reconstruction_param( - mpz_param_t mpz_param, param_images_t images, + mpz_param_t mpz_param, param_images_t images, vec_wit_t vec_wit, param_t *nmod_param, nvars_t nlins, nvars_t *linvars, uint32_t *lineqs, trace_det_fglm_mat_t trace_det, sp_matfglm_t *mat, mpz_upoly_t numer, @@ -1432,6 +1566,14 @@ static inline int rational_reconstruction_param( mulmat_reconstruct(trace_det, prod_crt); } + /* checks the candidate common denominator of the vector witnesses */ + if (vec_wit->state == 1) { + vec_wit_check(vec_wit, nmod_param->elim, prime); + if (vec_wit->state == 2 && info_level) { + fprintf(VERBSTREAM, "[v]"); + } + } + /** CRT PART **/ /* the images of the parametrization are only stored; CRT is applied when a rational reconstruction is attempted (see lift_parametrization) */ @@ -1445,6 +1587,16 @@ static inline int rational_reconstruction_param( mpz_mul_ui(modulus, modulus, prime); + /* simultaneous reconstruction of the vector witnesses, while the + eliminating polynomial is not lifted and no single witness is known */ + if (vec_wit->k && vec_wit->state == 0 && is_lifted[0] == 0 + && trace_det->done_trace < 2 && trace_det->done_det < 2 + && images->nprimes >= vec_wit->next) { + vec_wit_try(vec_wit, images, modulus); + vec_wit->next = MAX(images->nprimes + 1, + images->nprimes + images->nprimes * VEC_WIT_GROWTH / 100); + } + if (doit == 0) { return 0; } @@ -1483,7 +1635,12 @@ static inline int rational_reconstruction_param( } } - if (trace_det->done_trace > 1 || trace_det->done_det > 1) { + if (vec_wit->state == 2) { + mpz_lcm(guessed_den, guessed_den, vec_wit->den); + } + + if (trace_det->done_trace > 1 || trace_det->done_det > 1 + || vec_wit->state == 2) { mpz_t denominator; mpz_init_set_ui(denominator, 1); @@ -2573,6 +2730,8 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, param_images_t images; param_images_init(images, *mpz_paramp); param_images_add(images, nmod_params[0], primeinit); + vec_wit_t vec_wit; + vec_wit_init(vec_wit, nmod_params[0]->elim, tridx, detidx); deg_t nsols = (*mpz_paramp)->nsols; @@ -2761,6 +2920,11 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, if (rerun == 0) { mcheck = check_param_modular(*mpz_paramp, nmod_params[i], lp->p[i], is_lifted, check_lift, trace_det, info_level); + if (is_lifted[0] == 0) { + /* elim is wrong: the common denominator found by the vector + witnesses is discarded */ + vec_wit->state = 0; + } } else{ check_elim_poly(nmod_params[i]->elim, lp->p[i], numer->coeffs, denom->coeffs, check_lift); @@ -2768,7 +2932,7 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, crr = realtime(); if (mcheck == 1) { br = rational_reconstruction_param( - *mpz_paramp, images, nmod_params[i], + *mpz_paramp, images, vec_wit, nmod_params[i], bnlins[i], blinvars[i], lineqs_ptr[i], trace_det, bmatrix[i], numer, denom, check_lift, modulus, prod_crt, lp->p[i], &result, rnum, rden, recdata, @@ -2797,6 +2961,7 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, blens_extra_nf, bexps_extra_nf, bcfs_extra_nf, bdata_fglm, bdata_bms, num_gb, leadmons_ori, leadmons_current, bnlins, blinvars, linvars, lineqs_ptr, bsquvars, squvars, lmb_ori, field_char); + vec_wit_clear(vec_wit); param_images_clear(images); free_msolve_trace_qq_secondary_data(trace_det, modulus, prod_crt, result, test, rnum, rden, numer, denom, check_lift, @@ -2867,6 +3032,7 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, blens_extra_nf, bexps_extra_nf, bcfs_extra_nf, bdata_fglm, bdata_bms, num_gb, leadmons_ori, leadmons_current, bnlins, blinvars, linvars, lineqs_ptr, bsquvars, squvars, lmb_ori, field_char); + vec_wit_clear(vec_wit); param_images_clear(images); free_msolve_trace_qq_secondary_data(trace_det, modulus, prod_crt, result, test, rnum, rden, numer, denom, check_lift, From 9d7596ad8479013c22536836087eab572c5d684b Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Mon, 26 Oct 2026 12:20:32 +0100 Subject: [PATCH 8/9] Try a more skewed reconstruction bound for eliminating polynomials The eliminating polynomial is reconstructed from its highest degree down, carrying the lcm of the denominators already found, so most coefficients only bring a few new denominator bits and are bound by the size of their numerators. The bound D = M^(1/64) is now tried before D = M^(1/16), which lowers the number of primes needed by 2 to 5% on katsura and eco systems (kat12: 323 -> 313, eco13: 238 -> 231) with identical outputs. Random and minrank systems are unaffected. Co-Authored-By: Claude Opus 5.5 --- src/msolve/msolve.c | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index e93f825d..e1b6909c 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -1445,30 +1445,34 @@ static inline int lift_parametrization(mpz_param_t mpz_param, param_crt_init(pc, images); /* fixed bounds tried first, with M = modulus / 2^RR_MARGIN so that - RR_MARGIN bits of slack are kept: D = M^(1/16) for elim; D = M^(1/32), + RR_MARGIN bits of slack are kept: D = M^(1/64), then D = M^(1/16) for + elim; D = M^(1/32), balanced bounds, then N = M^(1/16) for the coordinates; the bounds set from the guessed denominator are tried last */ - mpz_t bN[4], bD[4], mm; - for(int k = 0; k < 4; k++){ + mpz_t bN[5], bD[5], mm; + for(int k = 0; k < 5; k++){ mpz_init(bN[k]); mpz_init(bD[k]); } mpz_init(mm); mpz_fdiv_q_2exp(mm, modulus, RR_MARGIN); /* elim */ - mpz_root(bD[0], mm, 16); + mpz_root(bD[0], mm, 64); mpz_fdiv_q(bN[0], mm, bD[0]); mpz_fdiv_q_2exp(bD[0], bD[0], 1); - /* coordinates */ - mpz_root(bD[1], mm, 32); + mpz_root(bD[1], mm, 16); mpz_fdiv_q(bN[1], mm, bD[1]); mpz_fdiv_q_2exp(bD[1], bD[1], 1); - mpz_sqrt(bD[2], mm); + /* coordinates */ + mpz_root(bD[2], mm, 32); mpz_fdiv_q(bN[2], mm, bD[2]); mpz_fdiv_q_2exp(bD[2], bD[2], 1); - mpz_root(bN[3], mm, 16); - mpz_fdiv_q(bD[3], mm, bN[3]); - mpz_fdiv_q_2exp(bN[3], bN[3], 1); + mpz_sqrt(bD[3], mm); + mpz_fdiv_q(bN[3], mm, bD[3]); + mpz_fdiv_q_2exp(bD[3], bD[3], 1); + mpz_root(bN[4], mm, 16); + mpz_fdiv_q(bD[4], mm, bN[4]); + mpz_fdiv_q_2exp(bN[4], bN[4], 1); mpz_clear(mm); int b = 0; @@ -1478,10 +1482,10 @@ static inline int lift_parametrization(mpz_param_t mpz_param, mpz_param->elim, denominator, images, 0, pc, nmod_param->elim->length, modulus, maxrec, coef, rnum, rden, numer, denom, check_lift, lcm, guessed_num, guessed_den, recdata, - bN, bD, 1); + bN, bD, 2); if (b == 0) { is_lifted[0] = 0; - for(int k = 0; k < 4; k++){ + for(int k = 0; k < 5; k++){ mpz_clear(bN[k]); mpz_clear(bD[k]); } @@ -1515,11 +1519,11 @@ static inline int lift_parametrization(mpz_param_t mpz_param, nmod_param, numer, denom, check_lift, modulus, prod_crt, prime, coef, rnum, rden, recdata, guessed_num, guessed_den, maxrec, - is_lifted, denominator, lcm, c, i, bN + 1, bD + 1, 3, + is_lifted, denominator, lcm, c, i, bN + 2, bD + 2, 3, info_level); } - for(int k = 0; k < 4; k++){ + for(int k = 0; k < 5; k++){ mpz_clear(bN[k]); mpz_clear(bD[k]); } From 29eb86b016f51b719742528214d893c0bcd9ccd0 Mon Sep 17 00:00:00 2001 From: mohabsafey Date: Wed, 28 Oct 2026 17:36:27 +0100 Subject: [PATCH 9/9] Fix small memory leaks in the setup of msolve_trace_qq AddressSanitizer reports a few constant-size leaks in the multi-modular setup, independent of the number of primes. Free invalid_gens once the input is imported, free the per-thread copies of the linear-form data and of the current leading monomials in free_msolve_trace_qq_initial_data (zero-initializing the pointer arrays), free the monomial basis when initial_modular_step fails, and free the FGLM data of the first call before the second call in the non-shape-position case. Also free the parametrization on the failure paths of nmod_fglm_compute_trace_data, drop a name buffer that was immediately overwritten when adding a linear form, and count the terms of the added linear form in gens->nterms so their coefficients are cleared. Co-Authored-By: Claude Opus 5.5 --- src/fglm/fglm_core.c | 2 ++ src/msolve/msolve.c | 29 +++++++++++++++++++++-------- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/src/fglm/fglm_core.c b/src/fglm/fglm_core.c index 85e8c66e..f11c0fe5 100644 --- a/src/fglm/fglm_core.c +++ b/src/fglm/fglm_core.c @@ -1480,6 +1480,7 @@ param_t *nmod_fglm_compute_trace_data(sp_matfglm_t *matrix, mod_t prime, nvars) == 0){ fprintf(ERRSTREAM, "Matrix is not invertible (there should be a bug)\n"); + free_fglm_param(param); return NULL; } @@ -1526,6 +1527,7 @@ param_t *nmod_fglm_compute_trace_data(sp_matfglm_t *matrix, mod_t prime, fprintf(ERRSTREAM, "Matrix is not invertible (there should be a bug)\n"); } *success = 0; + free_fglm_param(param); return NULL; } if (right_param == 1) { diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index e1b6909c..5fac84b1 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -585,7 +585,6 @@ static int add_linear_form_to_input_system(data_gens_ff_t *gens, * the end */ gens->vnames = realloc(gens->vnames, (unsigned long)gens->nvars * sizeof(char *)); - gens->vnames[gens->nvars - 1] = (char *)malloc(2 * sizeof(char)); gens->vnames[gens->nvars - 1] = extra_var; int32_t *old_exps = gens->exps; gens->exps = @@ -610,6 +609,7 @@ static int add_linear_form_to_input_system(data_gens_ff_t *gens, } else { gens->mpz_cfs = realloc(gens->mpz_cfs, (unsigned long)2 * len_new * (sizeof(mpz_t *))); + gens->nterms += nvars_new; for (i = 2 * len_old; i < 2 * len_new; i += 2) { gens->mpz_cfs[i] = (mpz_t *)malloc(sizeof(mpz_t)); mpz_init(*(gens->mpz_cfs[i])); @@ -686,7 +686,6 @@ static int add_random_linear_form_to_input_system(data_gens_ff_t *gens, * the end */ gens->vnames = realloc(gens->vnames, (unsigned long)gens->nvars * sizeof(char *)); - gens->vnames[gens->nvars - 1] = (char *)malloc(2 * sizeof(char)); gens->vnames[gens->nvars - 1] = extra_var; int32_t *old_exps = gens->exps; gens->exps = @@ -711,6 +710,7 @@ static int add_random_linear_form_to_input_system(data_gens_ff_t *gens, } else { gens->mpz_cfs = realloc(gens->mpz_cfs, (unsigned long)2 * len_new * (sizeof(mpz_t *))); + gens->nterms += nvars_new; for (i = 2 * len_old; i < 2 * len_new; i += 2) { gens->mpz_cfs[i] = (mpz_t *)malloc(sizeof(mpz_t)); mpz_init(*(gens->mpz_cfs[i])); @@ -1958,6 +1958,7 @@ static int32_t *initial_modular_step( if(md->info_level > 1){ fprintf (VERBSTREAM,"------------------------------------------------------------------------------------------------------\n"); } + free(lmb); free_basis_without_hash_table(&(bs)); free(bs); return NULL; @@ -1976,6 +1977,7 @@ static int32_t *initial_modular_step( *success = 0; *dim = 0; *dquot_ori = dquot; + free(lmb); free_basis_without_hash_table(&(bs)); free(bs); return NULL; @@ -1987,6 +1989,11 @@ static int32_t *initial_modular_step( md->fglm_rtime = realtime(); param_t **test_nmod_param = (param_t **)malloc(sizeof(param_t *)); + /* the second call allocates its own fglm data */ + free_fglm_data(*bdata_fglm); + *bdata_fglm = NULL; + free_fglm_bms_data(*bdata_bms); + *bdata_bms = NULL; *test_nmod_param = nmod_fglm_compute_trace_data(*bmatrix, fc, bs->ht->nv, *bsz, *nlins_ptr, linvars, @@ -2270,7 +2277,12 @@ static void free_msolve_trace_qq_initial_data(int *invalid_gens, md_t *st, prime free_fglm_param(nmod_params[i]); } free(leadmons_ori[i]); - //free(leadmons_current[i]); + free(leadmons_current[i]); + free(lineqs_ptr[i]); + if (i > 0) { + free(blinvars[i]); + free(bsquvars[i]); + } if (bs[i] != NULL) { free_basis(&(bs[i])); free(bs[i]); @@ -2427,6 +2439,7 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, prime_start, nr_primes, pbm_file, 0 /*truncate_lifting */, info_level)) { free(st); + free(invalid_gens); return -3; } @@ -2440,7 +2453,7 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, /* read in ideal, move coefficients to integers */ import_input_data(bs_qq, st, 0, st->ngens_input, lens, exps, cfs, invalid_gens); - //free(invalid_gens); + free(invalid_gens); invalid_gens = NULL; print_initial_statistics(VERBSTREAM, st); @@ -2534,15 +2547,15 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, /* data for linear forms */ nvars_t nlins = 0; /*number of linear forms*/ nvars_t *bnlins = (nvars_t *)calloc(st->nthrds, sizeof(nvars_t)); - nvars_t **blinvars = (nvars_t **)malloc( - st->nthrds * sizeof(nvars_t *)); /*indicates which variables are linear*/ + nvars_t **blinvars = (nvars_t **)calloc( + st->nthrds, sizeof(nvars_t *)); /*indicates which variables are linear*/ nvars_t *linvars = calloc(bs_qq->ht->nv, sizeof(nvars_t)); blinvars[0] = linvars; uint32_t **lineqs_ptr = - malloc(st->nthrds * sizeof(uint32_t *)); /*coeffs of linear forms*/ + calloc(st->nthrds, sizeof(uint32_t *)); /*coeffs of linear forms*/ /*data for squared variables*/ - nvars_t **bsquvars = (nvars_t **)malloc(st->nthrds * sizeof(nvars_t *)); + nvars_t **bsquvars = (nvars_t **)calloc(st->nthrds, sizeof(nvars_t *)); nvars_t *squvars = calloc(nr_vars - 1, sizeof(nvars_t)); bsquvars[0] = squvars;