From 9364154851069a44a8cebba05dbe818e9a019475 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Thu, 1 Oct 2026 00:04:57 +0200
Subject: [PATCH 1/4] Reduce memory used by multi-modular Groebner bases over
 QQ

Rework the storage used while lifting Groebner bases over the rationals
(msolve -g 2, export_groebner_qq):

- store the modular images of each polynomial in one contiguous
  prime-major block instead of one malloc-ed array per term;
- grow this block by max(4, nprimes/16) primes instead of 16*nthrds:
  the multi-modular loop is mono-threaded, so with -t N most of the
  allocated images were never used;
- drop the per-term array of CRT images (cf_zz): each coefficient is now
  lifted by multi-CRT into a scratch integer and immediately rationally
  reconstructed;
- store only integer numerators in cf_qq: once lifted, all denominators
  were set to 1 (the common denominator is kept in lm);
- compute num * (lcm/den) instead of (num * lcm) / den so that the
  result is not over-allocated.

Also fix the restart paths of core_groebner_qq (uninitialized dlift->S
passed to gb_modpoly_realloc, double free of already lifted polynomials,
leaked hm/cf_64/lead monomials, stale tracer reused when learning again),
a NULL dereference when F4 detects a bad prime while applying the tracer,
and a leak of trace data in duplicate_tracer.

On Katsura-10 the peak heap goes from 103MB to 73MB with one thread and
the peak RSS from 244MB to 70MB with -t 8; on Katsura-11 the peak RSS
goes from 458MB to 308MB. Outputs are unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017J51JHn4LH2YArUKp1VBjv
---
 src/msolve/duplicate.c  |   4 +-
 src/msolve/lifting-gb.c | 433 +++++++++++++++++++++++-----------------
 2 files changed, 256 insertions(+), 181 deletions(-)

diff --git a/src/msolve/duplicate.c b/src/msolve/duplicate.c
index 0bddb533..202de6ed 100644
--- a/src/msolve/duplicate.c
+++ b/src/msolve/duplicate.c
@@ -72,7 +72,9 @@ static inline void duplicate_tracer(
             memcpy(btrace[i]->lm, btrace[0]->lm,
                     (unsigned long)btrace[0]->lml * sizeof(sdm_t));
 
-            /* array of trace data per F4 round */
+            /* array of trace data per F4 round
+             * (replaces the one allocated by initialize_trace) */
+            free(btrace[i]->td);
             btrace[i]->td  = calloc((unsigned long)btrace[0]->ltd, sizeof(td_t));
 
             for(len_t l = 0; l < btrace[0]->ltd; ++l){
diff --git a/src/msolve/lifting-gb.c b/src/msolve/lifting-gb.c
index bda35a8b..de90084a 100644
--- a/src/msolve/lifting-gb.c
+++ b/src/msolve/lifting-gb.c
@@ -31,22 +31,26 @@
 
 typedef struct{
   uint32_t len; /* length of the encoded polynomial */
-  uint32_t **cf_32; /* array of arrays of coefficients
-                       * modulo several primes
-                       */
-  /* Warning: this is free-ed dynamically */
-  mpz_t *cf_zz; /* array which stores CRT lifting of
-                    the coefficients */
-  /* Warning: this is free-ed dynamically */
-  mpz_t *cf_qq; /* array which stores rational coefficients
-                  being lifted, numerators and denominators
-                  are given given as mpz_t
-               */
+  /* Coefficients modulo several primes, stored as one contiguous block of
+   * len * alloc entries in prime-major order: the coefficient of index j
+   * modulo the k-th prime is cf_32[k * len + j] (use MODCF below).
+   * Adding primes only grows the block at its end, and each new prime
+   * fills one contiguous row.
+   * Warning: this is free-ed dynamically (and set to NULL) as soon as the
+   * polynomial has been lifted and checked. */
+  uint32_t *cf_32;
+  /* len integer numerators of the lifted coefficients; they all share the
+   * common denominator lm (hence the polynomial is lm * LM + sum cf_qq[j]*m_j
+   * up to the order of the terms). */
+  mpz_t *cf_qq;
   mpz_t lm; /* stores the leading coefficient (hence all numerators of cf_qq should be divided by lm) */
 } modpolys_struct;
 
 typedef modpolys_struct modpolys_t[1];
 
+/* coefficient of index j of polynomial p modulo the k-th prime */
+#define MODCF(p, j, k) ((p)->cf_32[(size_t)(k) * (p)->len + (j)])
+
 typedef struct {
   uint32_t alloc; /* alloc -> max number of primes */
   uint32_t nprimes; /* number of primes */
@@ -204,14 +208,16 @@ static inline void gb_modpoly_init(gb_modpoly_t modgbs,
   }
   for(uint32_t i = 0; i < ld; i++){
     modgbs->modpolys[i]->len = lens[i];
-    modgbs->modpolys[i]->cf_32 = (uint32_t **)malloc(sizeof(uint32_t *)*lens[i]);
-    modgbs->modpolys[i]->cf_zz = (mpz_t *)malloc(sizeof(mpz_t)*lens[i]);
-    modgbs->modpolys[i]->cf_qq = (mpz_t *)malloc(sizeof(mpz_t)*2*lens[i]);
-    for(uint32_t j = 0; j < lens[i]; j++){
-      modgbs->modpolys[i]->cf_32[j] = calloc(alloc, sizeof(uint32_t));
-      mpz_init(modgbs->modpolys[i]->cf_zz[j]);
+    if(lens[i] > 0){
+      modgbs->modpolys[i]->cf_32 = (uint32_t *)calloc((size_t)lens[i] * alloc,
+                                                      sizeof(uint32_t));
+      modgbs->modpolys[i]->cf_qq = (mpz_t *)malloc(sizeof(mpz_t)*lens[i]);
+    }
+    else{
+      modgbs->modpolys[i]->cf_32 = NULL;
+      modgbs->modpolys[i]->cf_qq = NULL;
     }
-    for(uint32_t j = 0; j < 2 * lens[i]; j++){
+    for(uint32_t j = 0; j < lens[i]; j++){
       mpz_init(modgbs->modpolys[i]->cf_qq[j]);
     }
     mpz_init(modgbs->modpolys[i]->lm);
@@ -220,9 +226,10 @@ static inline void gb_modpoly_init(gb_modpoly_t modgbs,
 }
 
 
+/* Makes room for newalloc more primes. Polynomials whose modular data have
+ * already been free-ed (cf_32 == NULL) are skipped. */
 static inline void gb_modpoly_realloc(gb_modpoly_t modgbs,
-                                      uint32_t newalloc,
-                                      int32_t start){
+                                      uint32_t newalloc){
 
   uint32_t oldalloc = modgbs->alloc;
   modgbs->alloc += newalloc;
@@ -249,18 +256,21 @@ static inline void gb_modpoly_realloc(gb_modpoly_t modgbs,
   for(uint32_t i = oldalloc; i < modgbs->alloc; i++){
     modgbs->cf_64[i] = 0;
   }
-  for(uint32_t i = start; i < modgbs->ld; i++){
-    for(uint32_t j = 0; j < modgbs->modpolys[i]->len; j++){
-      uint32_t *newcfs_pol = (uint32_t *)realloc(modgbs->modpolys[i]->cf_32[j],
-                                                 modgbs->alloc * sizeof(uint32_t));
-      if(newcfs_pol == NULL){
-        fprintf(ERRSTREAM, "Problem when reallocating modgbs (cfs_pol)\n");
-      }
-      modgbs->modpolys[i]->cf_32[j] = newcfs_pol;
-      for(uint32_t k = oldalloc; k < modgbs->alloc; k++){
-        modgbs->modpolys[i]->cf_32[j][k] = 0;
-      }
+  for(uint32_t i = 0; i < modgbs->ld; i++){
+    const size_t len = modgbs->modpolys[i]->len;
+    if(modgbs->modpolys[i]->cf_32 == NULL || len == 0){
+      continue;
+    }
+    uint32_t *newcfs_pol = (uint32_t *)realloc(modgbs->modpolys[i]->cf_32,
+                                               len * modgbs->alloc * sizeof(uint32_t));
+    if(newcfs_pol == NULL){
+      fprintf(ERRSTREAM, "Problem when reallocating modgbs (cfs_pol)\n");
+      exit(1);
     }
+    /* prime-major storage: the new primes are appended at the end */
+    memset(newcfs_pol + len * oldalloc, 0,
+           len * (modgbs->alloc - oldalloc) * sizeof(uint32_t));
+    modgbs->modpolys[i]->cf_32 = newcfs_pol;
   }
 }
 
@@ -279,16 +289,20 @@ static inline void display_gbmodpoly_cf_32(FILE *file,
   for(uint32_t i = 0; i < modgbs->ld; i++){
     uint32_t len = modgbs->modpolys[i]->len;
     fprintf(file, "[%d, ", len);
+    if(modgbs->modpolys[i]->cf_32 == NULL){
+      fprintf(file, "(lifted)],\n");
+      continue;
+    }
     for(uint32_t j = 0; j < len; j++){
       fprintf(file, "[");
       for(uint32_t k = 0; k < modgbs->alloc-1; k++){
-        fprintf(file, "%d, ", modgbs->modpolys[i]->cf_32[j][k]);
+        fprintf(file, "%d, ", MODCF(modgbs->modpolys[i], j, k));
       }
       if(j < len - 1){
-        fprintf(file, "%d], ", modgbs->modpolys[i]->cf_32[j][modgbs->alloc-1]);
+        fprintf(file, "%d], ", MODCF(modgbs->modpolys[i], j, modgbs->alloc-1));
       }
       else{
-        fprintf(file, "%d]\n", modgbs->modpolys[i]->cf_32[j][modgbs->alloc-1]);
+        fprintf(file, "%d]\n", MODCF(modgbs->modpolys[i], j, modgbs->alloc-1));
       }
     }
     fprintf(file, "],\n");
@@ -311,10 +325,11 @@ static inline void display_modpoly(FILE *file,
   }
   display_monomial_single(file, gens, pos, &modgbs->ldm);
 
-  len_t idx, i, k;
+  len_t i, k;
   hm_t *hm    = NULL;
-  idx = modgbs->lmps[pos];
-  hm  = modgbs->hm[idx]+OFFSET;
+  /* modgbs->hm[pos] is the copy of the monomials of the pos-th element
+   * (see gb_modpoly_init), hence it is indexed by pos, not by lmps[pos] */
+  hm  = modgbs->hm[pos]+OFFSET;
 
   ht_t *ht = modgbs->bht;
   const len_t nv  = ht->nv;
@@ -335,25 +350,24 @@ static inline void display_modpoly(FILE *file,
     }
   }
   len_t len = modgbs->modpolys[pos]->len + 1;
+  /* once lifted, all coefficients are integers (common denominator is
+   * stored in lm), cf_qq[i] is the coefficient of index i */
+  mpz_t *cf = modgbs->modpolys[pos]->cf_qq;
   for(i = modgbs->modpolys[pos]->len-1; i > 0 ; i--){
-    if(mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[2*i], 0)>0){
+    if(mpz_cmp_ui(cf[i], 0)>0){
       fprintf(file, "+");
     }
     //coef is -1
-    if(mpz_cmp_si(modgbs->modpolys[pos]->cf_qq[2*i], -1) == 0 && mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[2*i + 1], 1) == 0){
+    if(mpz_cmp_si(cf[i], -1) == 0){
         fprintf(file, "-");
     }
     int star = 0;
     //coef is neither 1 nor -1
-    if((mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[2*i], 1) != 0 || mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[2*i + 1], 1) != 0) && (mpz_cmp_si(modgbs->modpolys[pos]->cf_qq[2*i], -1) != 0 || mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[2*i + 1], 1) != 0)){
-        mpz_out_str(file, 10, modgbs->modpolys[pos]->cf_qq[2*i]);
-        if(mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[2*i + 1], 1) != 0){
-            fprintf(file, "/");
-            mpz_out_str(file, 10, modgbs->modpolys[pos]->cf_qq[2*i + 1]);
-        }
+    if(mpz_cmp_ui(cf[i], 1) != 0 && mpz_cmp_si(cf[i], -1) != 0){
+        mpz_out_str(file, 10, cf[i]);
         fprintf(file, "*");
     }
-    if(mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[2*i], 0) != 0){
+    if(mpz_cmp_ui(cf[i], 0) != 0){
       for (k = 0; k < nv; ++k) {
          if (ht->ev[hm[len-i-1]][evi[k]] == 1) {
             if(star == 1){
@@ -373,7 +387,7 @@ static inline void display_modpoly(FILE *file,
     }
     fflush(file);
   }
-  if(mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[0], 0) > 0){
+  if(mpz_cmp_ui(cf[0], 0) > 0){
     fprintf(file, "+");
   }
   int deg = 0;
@@ -381,28 +395,20 @@ static inline void display_modpoly(FILE *file,
      deg += ht->ev[hm[len-1]][evi[k]];
   }
   if(deg == 0){
-    mpz_out_str(file, 10, modgbs->modpolys[pos]->cf_qq[0]);
-    if(mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[1], 1) != 0){
-      fprintf(file, "/");
-      mpz_out_str(file, 10, modgbs->modpolys[pos]->cf_qq[1]);
-    }
+    mpz_out_str(file, 10, cf[0]);
   }
   else{
-    if(mpz_cmp_si(modgbs->modpolys[pos]->cf_qq[0], -1) == 0 && mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[1], 1) == 0){
+    if(mpz_cmp_si(cf[0], -1) == 0){
         fprintf(file, "-");
     }
     int star = 0;
     //coef is neither 1 nor -1
-    if((mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[0], 1) != 0 || mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[1], 1) != 0) && (mpz_cmp_si(modgbs->modpolys[pos]->cf_qq[0], -1) != 0 || mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[1], 1) != 0)){
-        mpz_out_str(file, 10, modgbs->modpolys[pos]->cf_qq[0]);
-        if(mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[1], 1) != 0){
-            fprintf(file, "/");
-            mpz_out_str(file, 10, modgbs->modpolys[pos]->cf_qq[1]);
-        }
+    if(mpz_cmp_ui(cf[0], 1) != 0 && mpz_cmp_si(cf[0], -1) != 0){
+        mpz_out_str(file, 10, cf[0]);
         fprintf(file, "*");
     }
     fflush(file);
-    if(mpz_cmp_ui(modgbs->modpolys[pos]->cf_qq[0], 0) != 0){
+    if(mpz_cmp_ui(cf[0], 0) != 0){
       for (k = 0; k < nv; ++k) {
          if (ht->ev[hm[len-1]][evi[k]] == 1) {
             if(star == 1){
@@ -461,45 +467,37 @@ static inline void display_lm_gbmodpoly_cf_qq(FILE *file,
   fprintf(file, "]:\n");
 }
 
+/* frees everything but the hash table modgbs->bht */
 static inline void gb_modpoly_without_hash_table_clear(gb_modpoly_t modgbs){
   free(modgbs->primes);
+  modgbs->primes = NULL;
+  free(modgbs->cf_64);
+  modgbs->cf_64 = NULL;
   free(modgbs->ldm);
+  modgbs->ldm = NULL;
   for(uint32_t i = 0; i < modgbs->ld; i++){
+    /* cf_32 is NULL if already free-ed after lifting */
+    free(modgbs->modpolys[i]->cf_32);
+    modgbs->modpolys[i]->cf_32 = NULL;
     for(uint32_t j = 0; j < modgbs->modpolys[i]->len; j++){
-      free(modgbs->modpolys[i]->cf_32[j]);
-      mpz_clear(modgbs->modpolys[i]->cf_zz[j]);
-    }
-    for(uint32_t j = 0; j < 2 * modgbs->modpolys[i]->len; j++){
       mpz_clear(modgbs->modpolys[i]->cf_qq[j]);
     }
     mpz_clear(modgbs->modpolys[i]->lm);
-    free(modgbs->modpolys[i]->cf_32);
-    free(modgbs->modpolys[i]->cf_zz);
     free(modgbs->modpolys[i]->cf_qq);
+    modgbs->modpolys[i]->cf_qq = NULL;
+    free(modgbs->hm[i]);
   }
   free(modgbs->lmps);
+  modgbs->lmps = NULL;
   free(modgbs->hm);
+  modgbs->hm = NULL;
   free(modgbs->modpolys);
+  modgbs->modpolys = NULL;
 }
 
 static inline void gb_modpoly_clear(gb_modpoly_t modgbs){
-  free(modgbs->primes);
-  free(modgbs->ldm);
-
-  for(uint32_t i = 0; i < modgbs->ld; i++){
-    for(uint32_t j = 0; j < 2 * modgbs->modpolys[i]->len; j++){
-      mpz_clear(modgbs->modpolys[i]->cf_qq[j]);
-    }
-    mpz_clear(modgbs->modpolys[i]->lm);
-    free(modgbs->modpolys[i]->cf_qq);
-    free(modgbs->hm[i]);
-  }
-
-  free(modgbs->lmps);
-  free(modgbs->hm);
-  free(modgbs->cf_64);
+  gb_modpoly_without_hash_table_clear(modgbs);
   full_free_hash_table(&modgbs->bht);
-  free(modgbs->modpolys);
 }
 
 static inline int32_t degree(int32_t *mon, int nv){
@@ -630,10 +628,14 @@ static inline int modpgbs_set(gb_modpoly_t modgbs,
       /* hm  = bs->hm[idx]+OFFSET; */
       len = bs->hm[idx][LENGTH];
     }
+    if(modgbs->modpolys[i]->cf_32 == NULL){
+      /* already lifted, modular data are no longer needed */
+      continue;
+    }
     int32_t bc = modgbs->modpolys[i]->len - 1;
+    const uint32_t * const cfs = bs->cf_32[bs->hm[idx][COEFFS]];
     for (j = 1; j < len; ++j) {
-      uint32_t c = bs->cf_32[bs->hm[idx][COEFFS]][j];
-      modgbs->modpolys[i]->cf_32[bc][modgbs->nprimes] = c;
+      MODCF(modgbs->modpolys[i], bc, modgbs->nprimes) = cfs[j];
       bc--;
    }
   }
@@ -734,6 +736,8 @@ static int32_t gb_modular_trace_learning(gb_modpoly_t modgbs,
 
     /* Leading monomials from Grobner basis */
     int32_t *bexp_lm = get_lm_from_bs(bs, bht);
+    /* may be non NULL if the computation has been restarted */
+    free(leadmons[0]);
     leadmons[0] = bexp_lm;
 
     int32_t len = bs->lml;
@@ -784,6 +788,9 @@ static int32_t gb_modular_trace_learning(gb_modpoly_t modgbs,
             }
         }
         if(is_empty){
+            if(bexp_lm2 != NULL){
+              free(bexp_lm);
+            }
             free_basis_without_hash_table(&(bs));
             return is_empty;
         }
@@ -791,6 +798,10 @@ static int32_t gb_modular_trace_learning(gb_modpoly_t modgbs,
 
     /**************************************************/
 
+    if(bexp_lm2 != NULL){
+      /* leadmons[0] is bexp_lm2, bexp_lm is no longer needed */
+      free(bexp_lm);
+    }
     free_basis_without_hash_table(&(bs));
     return is_empty;
 }
@@ -827,6 +838,15 @@ static void gb_modular_trace_application(gb_modpoly_t modgbs,
   int32_t error = 0;
   bs = core_gba(bs_qq, st, &error, lp->p[0]);
   *stf4 = realtime()-rt;
+  /* core_gba returns NULL (and has already free-ed its data) when the
+   * prime has been detected as bad while applying the tracer */
+  if (bs == NULL || error) {
+      bad_primes[0] = 1;
+      if (bs != NULL) {
+          free_basis_and_only_local_hash_table_data(&bs);
+      }
+      return;
+  }
   ht_t **bht = &(bs->ht);
   for(len_t i = 0; i < modgbs->ld; i++){
         if(modgbs->modpolys[i]->len < bs->hm[bs->lmps[i]][LENGTH]-1){
@@ -840,10 +860,6 @@ static void gb_modular_trace_application(gb_modpoly_t modgbs,
             return;
         }
   }
-  if (bs == NULL) {
-      bad_primes[0] = 1;
-      return;
-  }
   int32_t lml = bs->lml;
   if (st->nev > 0) {
       int32_t j = 0;
@@ -893,7 +909,7 @@ static inline void choose_coef_to_lift(gb_modpoly_t modgbs, data_lift_t dlift){
     len_t d = 0;
     len_t len = modgbs->modpolys[i]->len;
     while(d < len - 1 && len > 0){
-      if(modgbs->modpolys[i]->cf_32[d][0]){
+      if(MODCF(modgbs->modpolys[i], d, 0)){
         dlift->coef[i] = d;
         break;
       }
@@ -915,7 +931,7 @@ static inline void incremental_dlift_crt(gb_modpoly_t modgbs, data_lift_t dlift,
   mpz_mul_ui(prod_p[0], mod_p[0], modgbs->primes[modgbs->nprimes - 1 ]);
   for(int32_t k = dlift->lstart; k <= dlift->lend; k++){
       if(modgbs->modpolys[k]->len){
-        uint32_t c = modgbs->modpolys[k]->cf_32[coef[k]][modgbs->nprimes  - 1 ];
+        uint32_t c = MODCF(modgbs->modpolys[k], coef[k], modgbs->nprimes  - 1);
         mpz_CRT_ui(dlift->crt[k], dlift->crt[k], mod_p[0],
                    c, modgbs->primes[modgbs->nprimes - 1 ],
                    prod_p[0], dlift->tmp, 1);
@@ -936,7 +952,7 @@ static inline void incremental_dlift_crt_full(gb_modpoly_t modgbs, data_lift_t d
   mpz_mul_ui(prod_p, mod_p, (uint32_t)newprime);
   for(int32_t k = dl->S; k < dl->end; k++){
       if(modgbs->modpolys[k]->len){
-        uint64_t c = modgbs->modpolys[k]->cf_32[coef[k]][modgbs->nprimes  - 1 ];
+        uint64_t c = MODCF(modgbs->modpolys[k], coef[k], modgbs->nprimes  - 1);
         mpz_CRT_ui(dl->crt[k], dl->crt[k], mod_p,
                    c, newprime, prod_p, dl->tmp, 1);
       }
@@ -945,39 +961,6 @@ static inline void incremental_dlift_crt_full(gb_modpoly_t modgbs, data_lift_t d
 }
 
 
-static inline void crt_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dlift,
-                                   int32_t start, int32_t end){
-  /* Data needed by multi CRT functions */
-  fmpz_comb_t comb;
-  fmpz_comb_temp_t comb_temp;
-
-  fmpz_comb_init(comb, modgbs->primes, modgbs->nprimes);
-  fmpz_comb_temp_init(comb_temp, comb);
-  fmpz_t y;
-  fmpz_init(y);
-
-  modpolys_t *polys = modgbs->modpolys;
-
-  for(int32_t k = start; k < end; k++){
-    if(dlift->check1[k] >= 1){
-      for(int32_t l = 0; l < polys[k]->len; l++){
-        for(uint32_t i = 0; i < modgbs->nprimes-1; i++){
-          modgbs->cf_64[i] = polys[k]->cf_32[l][i];
-        }
-        fmpz_multi_CRT_ui(y, modgbs->cf_64,
-                          comb, comb_temp, 1);
-        fmpz_get_mpz(polys[k]->cf_zz[l], y);
-      }
-    }
-  }
-
-  fmpz_clear(y);
-  fmpz_comb_temp_clear(comb_temp);
-  fmpz_comb_clear(comb);
-
-}
-
-
 static inline void set_recdata(data_lift_t dl, rrec_data_t rd1, rrec_data_t rd2, mpz_t mod_p){
   mpz_fdiv_q_2exp(rd1->N, mod_p, 1);
 
@@ -1003,60 +986,126 @@ static inline void set_recdata(data_lift_t dl, rrec_data_t rd1, rrec_data_t rd2,
 }
 
 /*
+  For the polynomials of index in [start, end) whose witness coefficient has
+  been lifted and checked (check1 > 1), lifts all coefficients by CRT
+  followed by rational reconstruction, and stores the integer numerators in
+  cf_qq and the common denominator in lm.
+
+  The CRT images are only temporaries: each coefficient is lifted with
+  multi-CRT into a single scratch integer which is immediately rationally
+  reconstructed, so that no array of CRT images is kept alive in memory.
+  Denominators are also kept in a scratch array, the size of the longest
+  polynomial handled here.
+
   returns the index of the poly if some of its coeff could not be lifted
-  else returns -1
+  else returns end
  */
-static inline int ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
-                                       int32_t start, int32_t end,
-                                       mpz_t mod_p, rrec_data_t rd1, rrec_data_t rd2){
+static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
+                                           int32_t start, int32_t end,
+                                           mpz_t mod_p, rrec_data_t rd1, rrec_data_t rd2,
+                                           double *st_crt, double *st_rrec){
+
+  modpolys_t *polys = modgbs->modpolys;
 
-  mpz_t rnum, rden;
-  mpz_init(rnum);
+  /* scratch data */
+  uint32_t maxlen = 0;
+  for(int32_t k = start; k < end; k++){
+    if(dl->check1[k] > 1 && polys[k]->len > maxlen){
+      maxlen = polys[k]->len;
+    }
+  }
+  mpz_t rden, zz;
   mpz_init(rden);
+  mpz_init(zz);
+  mpz_t *dens = NULL;
+  if(maxlen > 0){
+    dens = (mpz_t *)malloc(sizeof(mpz_t) * maxlen);
+    for(uint32_t l = 0; l < maxlen; l++){
+      mpz_init(dens[l]);
+    }
+  }
+
+  /* Data needed by multi CRT functions */
+  fmpz_comb_t comb;
+  fmpz_comb_temp_t comb_temp;
+  fmpz_t y;
+  if(maxlen > 0){
+    fmpz_comb_init(comb, modgbs->primes, modgbs->nprimes);
+    fmpz_comb_temp_init(comb_temp, comb);
+    fmpz_init(y);
+  }
 
-  modpolys_t *polys = modgbs->modpolys;
   set_recdata(dl, rd1, rd2, mod_p);
 
+  int32_t ret = end;
+  double rt;
+
   for(int32_t k = start; k < end; k++){
     if(dl->check1[k] > 1 && polys[k]->len > 0){
       mpz_set_ui(dl->tmp, 1);
+      int ok = 1;
       for(int32_t l = 0; l < polys[k]->len; l++){
-        if(ratreconwden(rnum, rden, polys[k]->cf_zz[l], mod_p, dl->den[k], rd1)){
-          mpz_set(polys[k]->cf_qq[2*l], rnum);
-          mpz_set(polys[k]->cf_qq[2*l + 1], rden);
-          mpz_lcm(dl->tmp, dl->tmp, rden);
+        rt = realtime();
+        for(uint32_t i = 0; i < modgbs->nprimes-1; i++){
+          modgbs->cf_64[i] = MODCF(polys[k], l, i);
+        }
+        fmpz_multi_CRT_ui(y, modgbs->cf_64,
+                          comb, comb_temp, 1);
+        fmpz_get_mpz(zz, y);
+        *st_crt += realtime() - rt;
+        rt = realtime();
+        if(ratreconwden(polys[k]->cf_qq[l], rden, zz, mod_p, dl->den[k], rd1)){
+          mpz_swap(dens[l], rden);
+          mpz_lcm(dl->tmp, dl->tmp, dens[l]);
+          *st_rrec += realtime() - rt;
         }
         else{
+          *st_rrec += realtime() - rt;
           mpz_set_ui(dl->gden, 1);
           for(int32_t i = k; i < modgbs->ld; i++){
             dl->check2[i] = 0;
           }
-          mpz_clear(rnum);
-          mpz_clear(rden);
-          return k;
+          ok = 0;
+          break;
         }
       }
+      if(!ok){
+        ret = k;
+        break;
+      }
       mpz_set(polys[k]->lm, dl->den[k]);
       for(int32_t l = 0; l < polys[k]->len; l++){
-        mpz_mul(polys[k]->cf_qq[2*l], polys[k]->cf_qq[2*l], dl->tmp);
-        mpz_divexact(polys[k]->cf_qq[2*l], polys[k]->cf_qq[2*l], polys[k]->cf_qq[2*l+1]);
-        mpz_set_ui(polys[k]->cf_qq[2*l+1], 1);
+        /* num * (lcm / den) rather than (num * lcm) / den: same value, but
+         * the limb array of the (long-lived) result is not over-allocated
+         * by the size of lcm (mpz_divexact does not shrink its output) */
+        mpz_divexact(dens[l], dl->tmp, dens[l]);
+        mpz_mul(polys[k]->cf_qq[l], polys[k]->cf_qq[l], dens[l]);
       }
       mpz_mul(polys[k]->lm, polys[k]->lm, dl->tmp);
       dl->check2[k] += 1;
     }
     if(dl->check1[k] < 1 && polys[k]->len > 0){
-      mpz_clear(rnum);
-      mpz_clear(rden);
-      return k;
+      ret = k;
+      break;
     }
     if(polys[k]->len == 0){
         dl->check2[k] += 1;
     }
   }
-  mpz_clear(rnum);
+
+  if(maxlen > 0){
+    fmpz_clear(y);
+    fmpz_comb_temp_clear(comb_temp);
+    fmpz_comb_clear(comb);
+    for(uint32_t l = 0; l < maxlen; l++){
+      mpz_clear(dens[l]);
+    }
+    free(dens);
+  }
+  mpz_clear(zz);
   mpz_clear(rden);
-  return end;
+
+  return ret;
 }
 
 
@@ -1075,46 +1124,42 @@ static inline int verif_coef(mpz_t num, mpz_t den, uint32_t prime, uint32_t coef
 
 static inline int verif_lifted_basis(gb_modpoly_t modgbs, data_lift_t dl,
                                      int thrds){
-  int b = 1;
-  mpz_t den;
-  mpz_init(den);
   for(int32_t k = 0; k < modgbs->ld; k++){
     if(dl->check1[k]>=1 && dl->check2[k]>0 && dl->check2[k] < NBCHECK ){
-      if(modgbs->modpolys[k]->len == 0){
+      modpolys_t *pol = modgbs->modpolys + k;
+      if(pol[0]->len == 0){
         dl->check2[k]++;
       }
       else{
         for(int i = 0; i < thrds; i++){
-          mp_limb_t prime = modgbs->primes[modgbs->nprimes - (thrds - i) ];
-          for(int32_t c = 0; c < modgbs->modpolys[k]->len; c++){
-            mpz_mul(den, modgbs->modpolys[k]->lm, modgbs->modpolys[k]->cf_qq[2*c+1]);
-            uint32_t coef = modgbs->modpolys[k]->cf_32[c][modgbs->nprimes   - (thrds - i) ];
-            int b = verif_coef(modgbs->modpolys[k]->cf_qq[2*c], den, prime, coef);
-            if(!b){
+          const uint32_t np = modgbs->nprimes - (thrds - i);
+          const mp_limb_t prime = modgbs->primes[np];
+          /* all coefficients share the denominator lm */
+          uint64_t ilc = mpz_fdiv_ui(pol[0]->lm, prime);
+          ilc = mod_p_inverse_32(ilc, prime);
+          for(int32_t c = 0; c < pol[0]->len; c++){
+            uint64_t v = mpz_fdiv_ui(pol[0]->cf_qq[c], prime);
+            v = (v * ilc) % prime;
+            if(v != MODCF(pol[0], c, np)){
               dl->check2[k] = 0;
               mpz_set_ui(dl->gden, 1);
-              b = 0;
-              mpz_clear(den);
-              return b;
+              return 0;
             }
           }
           dl->check2[k]++;
           if(dl->check2[k] >= NBCHECK){
             dl->S = k+1;
             dl->lstart = dl->S;
-            for(uint32_t j = 0; j < modgbs->modpolys[k]->len; j++){
-              free(modgbs->modpolys[k]->cf_32[j]);
-              mpz_clear(modgbs->modpolys[k]->cf_zz[j]);
-            }
-            free(modgbs->modpolys[k]->cf_32);
-            free(modgbs->modpolys[k]->cf_zz);
+            /* modular images are no longer needed */
+            free(pol[0]->cf_32);
+            pol[0]->cf_32 = NULL;
+            break;
           }
         }
       }
     }
   }
-  mpz_clear(den);
-  return b;
+  return 1;
 }
 
 static inline int verif_lifted_rational_wcoef(gb_modpoly_t modgbs, data_lift_t dl,
@@ -1132,7 +1177,7 @@ static inline int verif_lifted_rational_wcoef(gb_modpoly_t modgbs, data_lift_t d
       for(int i = 0; i < thrds; i++){
 
         mp_limb_t prime = modgbs->primes[modgbs->nprimes - (thrds - i) ];
-        uint32_t coef = modgbs->modpolys[k]->cf_32[dl->coef[k]][modgbs->nprimes  - (thrds - i) ];
+        uint32_t coef = MODCF(modgbs->modpolys[k], dl->coef[k], modgbs->nprimes  - (thrds - i));
         int boo = verif_coef(dl->num[k], dl->den[k], prime, coef);
 
         if(!boo){
@@ -1204,11 +1249,9 @@ static void ratrecon_gb(gb_modpoly_t modgbs, data_lift_t dl,
   }
   *st_rrec += realtime()-st;
   if(dl->lstart != dl->start){
-    crt_lift_modgbs(modgbs, dl, dl->lstart, dl->start);
-    *st_crt += realtime() - st;
-    st = realtime();
-    dl->start = ratrecon_lift_modgbs(modgbs, dl, dl->lstart, dl->start, mod_p, recdata1, recdata2);
-    *st_rrec += realtime()-st;
+    dl->start = crt_ratrecon_lift_modgbs(modgbs, dl, dl->lstart, dl->start,
+                                         mod_p, recdata1, recdata2,
+                                         st_crt, st_rrec);
   }
   dl->lstart = dl->start;
 
@@ -1256,8 +1299,7 @@ long max_bit_size_gb(gb_modpoly_t modgbs){
   long nb = 0;
   for(uint32_t i = 0; i < modgbs->ld; i++){
     for(uint32_t j = 0; j < modgbs->modpolys[i]->len; j++){
-      nb = MAX(nb, mpz_sizeinbase(modgbs->modpolys[i]->cf_qq[2*j], 2));
-      nb = MAX(nb, mpz_sizeinbase(modgbs->modpolys[i]->cf_qq[2*j + 1], 2));
+      nb = MAX(nb, mpz_sizeinbase(modgbs->modpolys[i]->cf_qq[j], 2));
     }
     nb = MAX(nb, mpz_sizeinbase(modgbs->modpolys[i]->lm, 2));
   }
@@ -1276,6 +1318,28 @@ long max_bit_size_gb(gb_modpoly_t modgbs){
   - returns -4 if there are too many bad primes
 */
 
+/* Frees the data which are (re)built by a learning phase, so that a new
+ * learning phase can start from scratch after a failed attempt. */
+static void reset_learning_data(mstrace_t msd, md_t *st){
+  for(int i = 0; i < st->nthrds; i++){
+    free(msd->leadmons_current[i]);
+    msd->leadmons_current[i] = NULL;
+    if(i > 0){
+      free(msd->leadmons_ori[i]);
+      msd->leadmons_ori[i] = NULL;
+      if(msd->btrace[i] != NULL){
+        free_trace(&(msd->btrace[i]));
+      }
+    }
+  }
+  /* the tracer learnt with the previous prime must not be applied when
+   * learning again */
+  if(st->tr != NULL){
+    free_trace(&(st->tr));
+  }
+  st->trace_level = NO_TRACER;
+}
+
 gb_modpoly_t *core_groebner_qq(
         gb_modpoly_t *modgbsp,
         bs_t *bs,
@@ -1372,7 +1436,7 @@ gb_modpoly_t *core_groebner_qq(
 
     apply = 1;
 
-    gb_modpoly_realloc((*modgbsp), 1, dlift->S);
+    gb_modpoly_realloc((*modgbsp), 1);
 
 #ifdef DEBUGGBLIFT
     display_gbmodpoly_cf_32(ERRSTREAM, (*modgbsp));
@@ -1411,6 +1475,7 @@ gb_modpoly_t *core_groebner_qq(
       }
 
       gb_modpoly_without_hash_table_clear((*modgbsp));
+      reset_learning_data(msd, st);
 
       free_rrec_data(recdata1);
       free_rrec_data(recdata2);
@@ -1463,7 +1528,14 @@ gb_modpoly_t *core_groebner_qq(
       prime = msd->lp->p[nthrds /* st->nthrds */ - 1];
 
       if((*modgbsp)->alloc <= nprimes + 2){
-        gb_modpoly_realloc((*modgbsp), 16*st->nthrds, dlift->S);
+        /* Room for the images modulo the next primes. The multi-modular
+         * computation is mono-threaded (one prime per step), so this must
+         * not scale with st->nthrds. The increment is kept small (at most
+         * ~1/16 of the current number of primes, but at least 4) since the
+         * unused slots are allocated for every term of the basis which is
+         * not lifted yet, i.e. for most of the basis at the end of the
+         * computation, when memory consumption peaks. */
+        gb_modpoly_realloc((*modgbsp), MAX(4, (uint32_t)nprimes / 16));
       }
 
       gb_modular_trace_application((*modgbsp), msd->mgb,
@@ -1547,6 +1619,7 @@ gb_modpoly_t *core_groebner_qq(
         }
 
         gb_modpoly_without_hash_table_clear((*modgbsp));
+        reset_learning_data(msd, st);
         free_rrec_data(recdata1);
         free_rrec_data(recdata2);
         /* reset info_level to previous setting */
@@ -1655,7 +1728,7 @@ uint64_t export_results_from_groebner_qq(
             int32_t len = gb->modpolys[i]->len;
 
             for(int32_t j = len-1; j >= 0; j--){
-                if (mpz_cmp_ui(gb->modpolys[i]->cf_qq[2 * j], 0) != 0) {
+                if (mpz_cmp_ui(gb->modpolys[i]->cf_qq[j], 0) != 0) {
                     nlen++;
                 }
             }
@@ -1709,11 +1782,11 @@ uint64_t export_results_from_groebner_qq(
         term++;
         if (print_gb == 2) {
             for(int32_t i = l-1; i >= 0; i--) {
-                if(mpz_cmp_ui(gb->modpolys[p]->cf_qq[2*i], 0) != 0) {
+                if(mpz_cmp_ui(gb->modpolys[p]->cf_qq[i], 0) != 0) {
                   for (int32_t n = 0; n < nve; n++) {
                     exp[term * nve + n] = ht->ev[hm[l - i]][evi[n]];
                   }
-                  mpz_set(cf_qq[term], gb->modpolys[p]->cf_qq[2 * i]);
+                  mpz_set(cf_qq[term], gb->modpolys[p]->cf_qq[i]);
 
                   term++;
                 }

From 716720023dc6629d80f2d1ccc4bda998ad4ab737 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Thu, 1 Oct 2026 00:04:57 +0200
Subject: [PATCH 2/4] Fix per-prime leaks in F4 and avoid touching unused hash
 table memory

- free the reducer rows created by symbolic preprocessing in the final
  reduction step (they were leaked at every F4 run with reduce_gb set,
  i.e. once per prime in multi-modular computations);
- when a bad prime is detected while applying the tracer, release the
  local meta data, the matrix and its new rows (finalize_f4 is not
  called on this path);
- copy only the used part (eld entries) of the hash table in
  copy_hash_table, which is called once per prime in the application
  phase.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017J51JHn4LH2YArUKp1VBjv
---
 src/neogb/f4.c   | 59 ++++++++++++++++++++++++++++++++++++++++++++++++
 src/neogb/hash.c |  7 ++++--
 2 files changed, 64 insertions(+), 2 deletions(-)

diff --git a/src/neogb/f4.c b/src/neogb/f4.c
index 7f63b716..f382af94 100644
--- a/src/neogb/f4.c
+++ b/src/neogb/f4.c
@@ -415,6 +415,44 @@ static int32_t initialize_f4(
     return done;
 }
 
+/* frees the mat->np new rows of mat->tr and their coefficient arrays
+ * (used when they are not moved to the basis) */
+static void free_new_matrix_rows(
+    mat_t *mat,
+    const md_t * const md
+    )
+{
+    len_t i;
+    if (mat->tr == NULL) {
+        return;
+    }
+    for (i = 0; i < mat->np; ++i) {
+        hm_t *row = mat->tr[i];
+        if (row == NULL) {
+            continue;
+        }
+        switch (md->ff_bits) {
+            case 8:
+                free(mat->cf_8[row[COEFFS]]);
+                mat->cf_8[row[COEFFS]] = NULL;
+                break;
+            case 16:
+                free(mat->cf_16[row[COEFFS]]);
+                mat->cf_16[row[COEFFS]] = NULL;
+                break;
+            case 32:
+                free(mat->cf_32[row[COEFFS]]);
+                mat->cf_32[row[COEFFS]] = NULL;
+                break;
+            default:
+                break;
+        }
+        free(row);
+        mat->tr[i] = NULL;
+    }
+    mat->np = 0;
+}
+
 static int32_t compute_new_elements(
     mat_t *mat,
     bs_t *bs,
@@ -437,6 +475,10 @@ static int32_t compute_new_elements(
             if (md->info_level > 0) {
                 fprintf(ERRSTREAM, "Wrong number of new elements, bad prime.");
             }
+            /* new rows (and their coefficient arrays) are owned by the
+             * matrix at this point, release them before bailing out */
+            free_new_matrix_rows(mat, md);
+            clear_matrix(mat);
             *errp = 1;
             return 1;
         }
@@ -615,6 +657,20 @@ static void reduce_final_basis(
 
         exact_linear_algebra(mat, bs, bs, md);
 
+        /* The known pivots used as reducers are multiples of basis
+         * elements: their coefficients belong to the basis, but the rows
+         * themselves were allocated by symbolic preprocessing. In this
+         * step, the linear algebra keeps them in its pivot array at the
+         * column of their lead term and does not free them (it only frees
+         * them when detecting a bad prime, in which case mat->np = 0).
+         * Release them here, otherwise they leak at each call. */
+        if (md->ff_bits != 0 && mat->np > 0) {
+            for (i = 0; i < mat->nru; ++i) {
+                free(mat->rr[i]);
+                mat->rr[i] = NULL;
+            }
+        }
+
         free_basis_elements(bs);
 
         convert_sparse_matrix_rows_to_basis_elements(
@@ -716,6 +772,9 @@ bs_t *core_f4(
     }
     if (*errp > 0) {
         free_basis_and_only_local_hash_table_data(&bs);
+        /* also release the local meta data (symbolic hash table, pair
+         * set, ...) and the matrix: finalize_f4 is not called here */
+        free_local_data(&mat, &md);
     } else {
         print_round_information_footer(VERBSTREAM, md);
 
diff --git a/src/neogb/hash.c b/src/neogb/hash.c
index 68f2c4f8..3d445cb4 100644
--- a/src/neogb/hash.c
+++ b/src/neogb/hash.c
@@ -157,7 +157,10 @@ ht_t *copy_hash_table(
     /* keep first entry empty for faster divisibility checks */
     ht->hd  = (hd_t *)calloc(ht->esz, sizeof(hd_t));
 
-    memcpy(ht->hd, bht->hd, (unsigned long)ht->esz * sizeof(hd_t));
+    /* only the first eld entries are in use: copying the whole allocated
+     * size would touch (and thus make resident) memory which is never
+     * read before being overwritten by an insertion */
+    memcpy(ht->hd, bht->hd, (unsigned long)bht->eld * sizeof(hd_t));
     ht->ev  = (exp_t **)malloc(ht->esz * sizeof(exp_t *));
     if (ht->ev == NULL) {
         fprintf(ERRSTREAM, "Computation needs too much memory on this machine,\n");
@@ -171,7 +174,7 @@ ht_t *copy_hash_table(
         fprintf(ERRSTREAM, "initialization failed, esz = %lu,\n", (unsigned long)ht->esz);
         fprintf(ERRSTREAM, "segmentation fault will follow.\n");
     }
-    memcpy(tmp, bht->ev[0], (unsigned long)ht->evl * ht->esz * sizeof(exp_t));
+    memcpy(tmp, bht->ev[0], (unsigned long)ht->evl * bht->eld * sizeof(exp_t));
     ht->eld = bht->eld;
     const hl_t esz  = ht->esz;
     for (j = 0; j < esz; ++j) {

From e26b8bbe34238c64dff8a6865eedce6979fefbc0 Mon Sep 17 00:00:00 2001
From: mohabsafey <mohab.safey@lip6.fr>
Date: Sat, 24 Oct 2026 16:44:21 +0200
Subject: [PATCH 3/4] Speed up the lifting of Groebner bases over QQ

Rational reconstruction now goes through FLINT (fmpq_reconstruct_fmpz_2,
half-gcd based) instead of the quadratic Euclidean algorithm whenever the
bounds allow it, falling back to ratrecon/ratreconwden otherwise.

The number of primes needed is also reduced and made reproducible:
- three witness coefficients are reconstructed and checked per polynomial
  (the last nonzero one and the ones found scanning downwards from 2*len/3
  and len/2) instead of one, so that a spurious witness reconstruction is
  detected before a full lift is attempted;
- the denominator used as a guess is taken from the previous polynomial
  (its lifted denominator, or the one of its witnesses once verified) and
  combined with the first witness denominator for the next witnesses; it
  is no longer reset to 1 after a failed lift or a failed verification;
- the full lift of a polynomial reuses the denominator accumulated so far
  as a guess, so that only the missing factor is reconstructed for each
  coefficient;
- reconstruction bounds keep 16 bits of slack, which makes spurious
  reconstructions unlikely;
- witnesses are reconstructed every rr primes, rr being the largest power
  of 2 <= nprimes/128, instead of a period depending on timings, so that
  the number of primes no longer depends on the machine load.

On 31 benchmark systems (5 random seeds each, single thread) the outputs
are unchanged and the number of primes decreases on 27 of them (up to
26%), while being much more stable across seeds. The average speedup
(geometric mean over all runs) is 1.29, and the speedup on the total time
of the benchmark suite is 1.25. The largest speedups are reached on the
MinRank (minrank) family, up to 1.91, on the critical point (cp_d)
family, up to 1.79, and on the dense random (random2) family, up to 1.71.
The absolute time saved grows with the problem size, and on the eco and
cp_d families so does the speedup.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/msolve/lifting-gb.c | 247 +++++++++++++++++++++++++++++++++++++---
 1 file changed, 229 insertions(+), 18 deletions(-)

diff --git a/src/msolve/lifting-gb.c b/src/msolve/lifting-gb.c
index de90084a..32d0c711 100644
--- a/src/msolve/lifting-gb.c
+++ b/src/msolve/lifting-gb.c
@@ -19,9 +19,67 @@
  * Mohab Safey El Din */
 
 #include<flint/fmpz.h>
+#include<flint/fmpq.h>
 #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
+/* witnesses are reconstructed every rr primes, rr being the largest
+   power of 2 <= nprimes / RR_PERIOD_DIV */
+#define RR_PERIOD_DIV 128
 #ifndef MIN
 #define MIN(x, y) ((x) > (y) ? (y) : (x))
 #endif
@@ -98,6 +156,16 @@ typedef struct{
   int *check1; /* tells whether lifted data are ok with one more prime */
   int *check2; /* tells whether lifted data are ok with two more primes */
   int32_t S;
+  int32_t *coef2; /* second witness coefficient */
+  mpz_t *crt2;
+  mpz_t *num2;
+  mpz_t *den2;
+  mpz_t *cden; /* lcm of the witness denominators */
+  mpz_t pgden; /* denominator guessed from the previous polynomial */
+  int32_t *coef3; /* third witness coefficient */
+  mpz_t *crt3;
+  mpz_t *num3;
+  mpz_t *den3;
 } data_lift_struct;
 
 typedef data_lift_struct data_lift_t[1];
@@ -111,6 +179,7 @@ static inline void data_lift_init(data_lift_t dl,
   dl->lstart = 0;
   dl->nsteps = nsteps;
   dl->S = 0;
+  mpz_init_set_ui(dl->pgden, 1);
   int32_t i;
 
   dl->steps = (int32_t *)calloc(nsteps, sizeof(int32_t));
@@ -141,6 +210,24 @@ static inline void data_lift_init(data_lift_t dl,
   dl->end = npol;
   dl->check1 = calloc(npol, sizeof(int));
   dl->check2 = calloc(npol, sizeof(int));
+  dl->coef2 = calloc(npol, sizeof(int32_t));
+  dl->crt2 = malloc(sizeof(mpz_t) * npol);
+  dl->num2 = malloc(sizeof(mpz_t) * npol);
+  dl->den2 = malloc(sizeof(mpz_t) * npol);
+  dl->cden = malloc(sizeof(mpz_t) * npol);
+  dl->coef3 = calloc(npol, sizeof(int32_t));
+  dl->crt3 = malloc(sizeof(mpz_t) * npol);
+  dl->num3 = malloc(sizeof(mpz_t) * npol);
+  dl->den3 = malloc(sizeof(mpz_t) * npol);
+  for(i = 0; i < npol; i++){
+    mpz_init(dl->crt2[i]);
+    mpz_init(dl->num2[i]);
+    mpz_init(dl->den2[i]);
+    mpz_init(dl->cden[i]);
+    mpz_init(dl->crt3[i]);
+    mpz_init(dl->num3[i]);
+    mpz_init(dl->den3[i]);
+  }
 
 }
 
@@ -165,8 +252,27 @@ static inline void data_lift_clear(data_lift_t dl){
 
   mpz_clear(dl->gden);
   mpz_clear(dl->tmp);
+  mpz_clear(dl->pgden);
   free(dl->check1);
   free(dl->check2);
+  for(int32_t i = 0; i < dl->npol; i++){
+    mpz_clear(dl->crt2[i]);
+    mpz_clear(dl->num2[i]);
+    mpz_clear(dl->den2[i]);
+    mpz_clear(dl->cden[i]);
+    mpz_clear(dl->crt3[i]);
+    mpz_clear(dl->num3[i]);
+    mpz_clear(dl->den3[i]);
+  }
+  free(dl->crt2);
+  free(dl->num2);
+  free(dl->den2);
+  free(dl->cden);
+  free(dl->coef2);
+  free(dl->crt3);
+  free(dl->num3);
+  free(dl->den3);
+  free(dl->coef3);
 
 }
 
@@ -903,20 +1009,33 @@ static void gb_modular_trace_application(gb_modpoly_t modgbs,
 
 }
 
+/* three witness coefficients are chosen per polynomial: the last nonzero
+   one and the nonzero ones found scanning downwards from 2*len/3 and
+   len/2 */
 static inline void choose_coef_to_lift(gb_modpoly_t modgbs, data_lift_t dlift){
   uint32_t ld = modgbs->ld;
   for(int32_t i = 0; i < ld; i++){
-    len_t d = 0;
     len_t len = modgbs->modpolys[i]->len;
-    while(d < len - 1 && len > 0){
+    len_t d = len - 1;
+    while(d > 0 && len > 0){
       if(MODCF(modgbs->modpolys[i], d, 0)){
         dlift->coef[i] = d;
         break;
       }
       else{
-        d++;
+        d--;
       }
     }
+    len_t d2 = 2 * len / 3;
+    while(d2 > 0 && !MODCF(modgbs->modpolys[i], d2, 0)){
+      d2--;
+    }
+    dlift->coef2[i] = d2;
+    len_t d3 = len / 2;
+    while(d3 > 0 && !MODCF(modgbs->modpolys[i], d3, 0)){
+      d3--;
+    }
+    dlift->coef3[i] = d3;
   }
 }
 
@@ -955,6 +1074,12 @@ static inline void incremental_dlift_crt_full(gb_modpoly_t modgbs, data_lift_t d
         uint64_t c = MODCF(modgbs->modpolys[k], coef[k], modgbs->nprimes  - 1);
         mpz_CRT_ui(dl->crt[k], dl->crt[k], mod_p,
                    c, newprime, prod_p, dl->tmp, 1);
+        c = MODCF(modgbs->modpolys[k], dl->coef2[k], modgbs->nprimes  - 1);
+        mpz_CRT_ui(dl->crt2[k], dl->crt2[k], mod_p,
+                   c, newprime, prod_p, dl->tmp, 1);
+        c = MODCF(modgbs->modpolys[k], dl->coef3[k], modgbs->nprimes  - 1);
+        mpz_CRT_ui(dl->crt3[k], dl->crt3[k], mod_p,
+                   c, newprime, prod_p, dl->tmp, 1);
       }
   }
   mpz_set(mod_p, prod_p);
@@ -962,7 +1087,9 @@ static inline void incremental_dlift_crt_full(gb_modpoly_t modgbs, data_lift_t d
 
 
 static inline void set_recdata(data_lift_t dl, rrec_data_t rd1, rrec_data_t rd2, mpz_t mod_p){
-  mpz_fdiv_q_2exp(rd1->N, mod_p, 1);
+  /* keep RR_MARGIN bits of slack (N*D <= mod_p / 2^(1+RR_MARGIN)) so that
+     spurious reconstructions become unlikely */
+  mpz_fdiv_q_2exp(rd1->N, mod_p, 1 + RR_MARGIN);
 
   if(dl->lstart){
 
@@ -1042,7 +1169,10 @@ static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
 
   for(int32_t k = start; k < end; k++){
     if(dl->check1[k] > 1 && polys[k]->len > 0){
-      mpz_set_ui(dl->tmp, 1);
+      /* dl->tmp is the lcm of the denominators found so far, starting
+         from the witness ones; it is used as a guess for the next
+         coefficient, so that only the missing factor is reconstructed */
+      mpz_set(dl->tmp, dl->cden[k]);
       int ok = 1;
       for(int32_t l = 0; l < polys[k]->len; l++){
         rt = realtime();
@@ -1054,14 +1184,16 @@ static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
         fmpz_get_mpz(zz, y);
         *st_crt += realtime() - rt;
         rt = realtime();
-        if(ratreconwden(polys[k]->cf_qq[l], rden, zz, mod_p, dl->den[k], rd1)){
-          mpz_swap(dens[l], rden);
-          mpz_lcm(dl->tmp, dl->tmp, dens[l]);
+        if(lift_ratreconwden(polys[k]->cf_qq[l], rden, zz, mod_p, dl->tmp, rd1)){
+          mpz_mul(dens[l], rden, dl->tmp);
+          mpz_mul(dl->tmp, dl->tmp, rden);
           *st_rrec += realtime() - rt;
         }
         else{
           *st_rrec += realtime() - rt;
-          mpz_set_ui(dl->gden, 1);
+          /* the witnesses of this poly are still valid: their common
+             denominator remains the best guess for reconstructing them */
+          mpz_set(dl->gden, dl->cden[k]);
           for(int32_t i = k; i < modgbs->ld; i++){
             dl->check2[i] = 0;
           }
@@ -1073,7 +1205,7 @@ static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
         ret = k;
         break;
       }
-      mpz_set(polys[k]->lm, dl->den[k]);
+      mpz_set(polys[k]->lm, dl->tmp);
       for(int32_t l = 0; l < polys[k]->len; l++){
         /* num * (lcm / den) rather than (num * lcm) / den: same value, but
          * the limb array of the (long-lived) result is not over-allocated
@@ -1081,7 +1213,6 @@ static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
         mpz_divexact(dens[l], dl->tmp, dens[l]);
         mpz_mul(polys[k]->cf_qq[l], polys[k]->cf_qq[l], dens[l]);
       }
-      mpz_mul(polys[k]->lm, polys[k]->lm, dl->tmp);
       dl->check2[k] += 1;
     }
     if(dl->check1[k] < 1 && polys[k]->len > 0){
@@ -1142,7 +1273,6 @@ static inline int verif_lifted_basis(gb_modpoly_t modgbs, data_lift_t dl,
             v = (v * ilc) % prime;
             if(v != MODCF(pol[0], c, np)){
               dl->check2[k] = 0;
-              mpz_set_ui(dl->gden, 1);
               return 0;
             }
           }
@@ -1179,6 +1309,14 @@ static inline int verif_lifted_rational_wcoef(gb_modpoly_t modgbs, data_lift_t d
         mp_limb_t prime = modgbs->primes[modgbs->nprimes - (thrds - i) ];
         uint32_t coef = MODCF(modgbs->modpolys[k], dl->coef[k], modgbs->nprimes  - (thrds - i));
         int boo = verif_coef(dl->num[k], dl->den[k], prime, coef);
+        if(boo){
+          coef = MODCF(modgbs->modpolys[k], dl->coef2[k], modgbs->nprimes  - (thrds - i));
+          boo = verif_coef(dl->num2[k], dl->den2[k], prime, coef);
+        }
+        if(boo){
+          coef = MODCF(modgbs->modpolys[k], dl->coef3[k], modgbs->nprimes  - (thrds - i));
+          boo = verif_coef(dl->num3[k], dl->den3[k], prime, coef);
+        }
 
         if(!boo){
           for(int32_t kk = k; kk < dl->end; kk++){
@@ -1199,8 +1337,8 @@ static inline int verif_lifted_rational_wcoef(gb_modpoly_t modgbs, data_lift_t d
 static inline int reconstructcoeff(data_lift_t dl, int32_t i, mpz_t mod_p,
                                    rrec_data_t recdata1, rrec_data_t recdata2){
 
-  int b = ratreconwden(dl->num[i], dl->den[i],
-                       dl->crt[i], mod_p, dl->gden, recdata1);
+  int b = lift_ratreconwden(dl->num[i], dl->den[i],
+                            dl->crt[i], mod_p, dl->gden, recdata1);
   if(b){
     mpz_mul(dl->den[i], dl->den[i], dl->gden);
     mpz_gcd(dl->tmp, dl->den[i], dl->num[i]);
@@ -1212,8 +1350,8 @@ static inline int reconstructcoeff(data_lift_t dl, int32_t i, mpz_t mod_p,
   }
   else{
 
-    b = ratrecon(dl->num[i], dl->den[i],
-                 dl->crt[i], mod_p, recdata2);
+    b = lift_ratrecon(dl->num[i], dl->den[i],
+                      dl->crt[i], mod_p, recdata2);
     if(b){
       mpz_set(dl->gden, dl->den[i]);
     }
@@ -1221,6 +1359,54 @@ static inline int reconstructcoeff(data_lift_t dl, int32_t i, mpz_t mod_p,
   return b;
 }
 
+/* reconstructs the second witness using as a guess the lcm of the
+   denominator of the first one and of the one guessed from the previous
+   polynomial; cden[i] is set to the lcm of both witness denominators */
+static inline int reconstructcoeff2(data_lift_t dl, int32_t i, mpz_t mod_p,
+                                    rrec_data_t recdata1, rrec_data_t recdata2){
+  mpz_lcm(dl->cden[i], dl->den[i], dl->pgden);
+  if(!lift_ratreconwden(dl->num2[i], dl->den2[i],
+                        dl->crt2[i], mod_p, dl->cden[i], recdata1)){
+    /* same fallback as for the first witness: balanced bounds, no guess */
+    if(!lift_ratrecon(dl->num2[i], dl->den2[i], dl->crt2[i], mod_p, recdata2)){
+      return 0;
+    }
+    mpz_lcm(dl->cden[i], dl->den[i], dl->den2[i]);
+    mpz_set(dl->gden, dl->cden[i]);
+    return 1;
+  }
+  mpz_mul(dl->den2[i], dl->den2[i], dl->cden[i]);
+  mpz_gcd(dl->tmp, dl->den2[i], dl->num2[i]);
+  mpz_divexact(dl->num2[i], dl->num2[i], dl->tmp);
+  mpz_divexact(dl->den2[i], dl->den2[i], dl->tmp);
+  mpz_lcm(dl->cden[i], dl->den[i], dl->den2[i]);
+  mpz_set(dl->gden, dl->cden[i]);
+  return 1;
+}
+
+/* reconstructs the third witness using the lcm of the denominators of the
+   first two witnesses as a guess; cden[i] is updated accordingly */
+static inline int reconstructcoeff3(data_lift_t dl, int32_t i, mpz_t mod_p,
+                                    rrec_data_t recdata1, rrec_data_t recdata2){
+  if(!lift_ratreconwden(dl->num3[i], dl->den3[i],
+                        dl->crt3[i], mod_p, dl->cden[i], recdata1)){
+    /* same fallback as for the first witness: balanced bounds, no guess */
+    if(!lift_ratrecon(dl->num3[i], dl->den3[i], dl->crt3[i], mod_p, recdata2)){
+      return 0;
+    }
+    mpz_lcm(dl->cden[i], dl->cden[i], dl->den3[i]);
+    mpz_set(dl->gden, dl->cden[i]);
+    return 1;
+  }
+  mpz_mul(dl->den3[i], dl->den3[i], dl->cden[i]);
+  mpz_gcd(dl->tmp, dl->den3[i], dl->num3[i]);
+  mpz_divexact(dl->num3[i], dl->num3[i], dl->tmp);
+  mpz_divexact(dl->den3[i], dl->den3[i], dl->tmp);
+  mpz_lcm(dl->cden[i], dl->cden[i], dl->den3[i]);
+  mpz_set(dl->gden, dl->cden[i]);
+  return 1;
+}
+
 static void ratrecon_gb(gb_modpoly_t modgbs, data_lift_t dl,
                         mpz_t mod_p, mpz_t prod_p,
                         rrec_data_t recdata1, rrec_data_t recdata2,
@@ -1278,8 +1464,27 @@ static void ratrecon_gb(gb_modpoly_t modgbs, data_lift_t dl,
           dl->check1[i]++;
       }
       else{
+        /* guess the denominator from the previous polynomial (its lifted
+           denominator, or the one of its witnesses once they are
+           verified), so that a spurious reconstruction never becomes the
+           guess */
+        if(i > 0){
+          if(dl->check2[i-1] > 0 && modgbs->modpolys[i-1]->len > 0){
+            mpz_set(dl->gden, modgbs->modpolys[i-1]->lm);
+          }
+          else if(dl->check1[i-1] > 1 && modgbs->modpolys[i-1]->len > 0){
+            mpz_set(dl->gden, dl->cden[i-1]);
+          }
+        }
+        mpz_set(dl->pgden, dl->gden);
         int b = reconstructcoeff(dl, i, mod_p,
                                recdata1, recdata2);
+        if(b){
+          b = reconstructcoeff2(dl, i, mod_p, recdata1, recdata2);
+        }
+        if(b){
+          b = reconstructcoeff3(dl, i, mod_p, recdata1, recdata2);
+        }
         if(!b){
           dl->check1[i] = 0;
           break;
@@ -1636,8 +1841,14 @@ gb_modpoly_t *core_groebner_qq(
         ratrecon_gb((*modgbsp), dlift, msd->mod_p, msd->prod_p, recdata1, recdata2,
                     nthrds/* st->nthrds */, &st_crt, &st_rrec, &st_wit);
       }
-      if(st_wit > 2*dlift->rr * stf4){
-        dlift->rr = 2*dlift->rr;
+      /* 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;
+      }
+      if(rr != dlift->rr){
+        dlift->rr = rr;
         if(info_level){
           fprintf(VERBSTREAM, "(->%d)", dlift->rr);
 	        fflush(VERBSTREAM);

From 76c931c7e0da822cbbae18317b7d2808d7f8f2a8 Mon Sep 17 00:00:00 2001
From: mohabsafey <mohab.safey@lip6.fr>
Date: Sun, 25 Oct 2026 10:32:28 +0100
Subject: [PATCH 4/4] Choose rational reconstruction bounds from the guessed
 denominator

When a denominator g is guessed, the value reconstructed is u * g, whose
numerator is about g times the denominator still missing. The bounds
were so far either balanced (N = D = sqrt(M), with M = mod_p / 2^(1 +
RR_MARGIN)) or, once lstart > 0, fixed to D = M^(1/3) and N = M^(2/3)
whatever the size of g. They are now set to D = sqrt(M / g) and
N = M / D, and recomputed whenever the guess changes: for each witness
coefficient, and during the full lift of a polynomial each time the
accumulated denominator grows. Balanced bounds are kept for the
reconstructions done without a guess.

On 38 benchmark systems over QQ (5 random seeds each, single thread)
the outputs are unchanged and the number of primes decreases on 31 of
them (up to 33%), the other 7 being unchanged. Compared with the
previous commit, the average speedup (geometric mean over all runs) is
1.24 and the speedup on the total time of the benchmark suite is 1.28;
compared with v0.11.0 these are 1.76 and 1.78. The largest speedups
over the previous commit are reached on the dense random families
(random2 up to 1.57, random3 up to 1.51), on the MinRank family (up to
1.54) and on the critical point families (cp_d and crit, up to 1.48).
Compared with v0.11.0, the speedups are between 1.95 and 3.00 on these
families (random2 up to 3.00, and 5.48 on the small random2_6_6;
MinRank up to 2.55, cp_d 2.54, crit 2.52, random3 2.41), and between
1.0 and 1.43 on the other ones (kat up to 1.43, chandra 1.39, eco
1.27, bruno 1.27, zimmermann 1.22, henrion 1.17, vor1 1.07, cyclic
1.05, noon 1.07, KLYZ22 1.04). Peak memory decreases on most systems,
and is up to 4 times smaller than with v0.11.0.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/msolve/lifting-gb.c | 57 ++++++++++++++++++++++-------------------
 1 file changed, 31 insertions(+), 26 deletions(-)

diff --git a/src/msolve/lifting-gb.c b/src/msolve/lifting-gb.c
index 32d0c711..14f664a1 100644
--- a/src/msolve/lifting-gb.c
+++ b/src/msolve/lifting-gb.c
@@ -1086,30 +1086,28 @@ static inline void incremental_dlift_crt_full(gb_modpoly_t modgbs, data_lift_t d
 }
 
 
-static inline void set_recdata(data_lift_t dl, rrec_data_t rd1, rrec_data_t rd2, mpz_t mod_p){
-  /* keep RR_MARGIN bits of slack (N*D <= mod_p / 2^(1+RR_MARGIN)) so that
-     spurious reconstructions become unlikely */
-  mpz_fdiv_q_2exp(rd1->N, mod_p, 1 + RR_MARGIN);
-
-  if(dl->lstart){
-
-    mpz_set(rd2->N, rd1->N);
-
-    mpz_sqrt(rd2->N, rd2->N);
-    mpz_set(rd2->D, rd2->N);
-
-    mpz_root(rd1->D, rd1->N, 3);
-    mpz_fdiv_q(rd1->N, rd1->N, rd1->D);
-
-  }
-  else{
-    mpz_sqrt(rd1->N, rd1->N);
-    mpz_set(rd1->D, rd1->N);
-
-    mpz_set(rd2->N, rd1->N);
-    mpz_set(rd2->D, rd1->D);
+/* balanced bounds, used when no denominator is guessed; RR_MARGIN bits of
+   slack are kept (N*D <= mod_p / 2^(1+RR_MARGIN)) so that spurious
+   reconstructions become unlikely */
+static inline void set_recdata(rrec_data_t rd, mpz_t mod_p){
+  mpz_fdiv_q_2exp(rd->N, mod_p, 1 + RR_MARGIN);
+  mpz_sqrt(rd->N, rd->N);
+  mpz_set(rd->D, rd->N);
+}
 
-  }
+/* bounds used to reconstruct u * gden, where gden is a guessed denominator:
+   the numerator is then expected to be about gden times the remaining
+   denominator, so that with M = mod_p / 2^(1+RR_MARGIN), the bounds are
+   D = sqrt(M / gden) and N = M / D */
+static inline void set_recdata_guess(rrec_data_t rd, const mpz_t mod_p,
+                                     const mpz_t gden){
+  mpz_fdiv_q_2exp(rd->N, mod_p, 1 + RR_MARGIN);
+  mpz_fdiv_q(rd->D, rd->N, gden);
+  mpz_sqrt(rd->D, rd->D);
+  if(mpz_sgn(rd->D) == 0){
+    mpz_set_ui(rd->D, 1);
+  }
+  mpz_fdiv_q(rd->N, rd->N, rd->D);
 }
 
 /*
@@ -1162,7 +1160,7 @@ static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
     fmpz_init(y);
   }
 
-  set_recdata(dl, rd1, rd2, mod_p);
+  set_recdata(rd2, mod_p);
 
   int32_t ret = end;
   double rt;
@@ -1173,6 +1171,7 @@ static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
          from the witness ones; it is used as a guess for the next
          coefficient, so that only the missing factor is reconstructed */
       mpz_set(dl->tmp, dl->cden[k]);
+      set_recdata_guess(rd1, mod_p, dl->tmp);
       int ok = 1;
       for(int32_t l = 0; l < polys[k]->len; l++){
         rt = realtime();
@@ -1186,7 +1185,10 @@ static inline int crt_ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
         rt = realtime();
         if(lift_ratreconwden(polys[k]->cf_qq[l], rden, zz, mod_p, dl->tmp, rd1)){
           mpz_mul(dens[l], rden, dl->tmp);
-          mpz_mul(dl->tmp, dl->tmp, rden);
+          if(mpz_cmp_ui(rden, 1)){
+            mpz_mul(dl->tmp, dl->tmp, rden);
+            set_recdata_guess(rd1, mod_p, dl->tmp);
+          }
           *st_rrec += realtime() - rt;
         }
         else{
@@ -1337,6 +1339,7 @@ static inline int verif_lifted_rational_wcoef(gb_modpoly_t modgbs, data_lift_t d
 static inline int reconstructcoeff(data_lift_t dl, int32_t i, mpz_t mod_p,
                                    rrec_data_t recdata1, rrec_data_t recdata2){
 
+  set_recdata_guess(recdata1, mod_p, dl->gden);
   int b = lift_ratreconwden(dl->num[i], dl->den[i],
                             dl->crt[i], mod_p, dl->gden, recdata1);
   if(b){
@@ -1365,6 +1368,7 @@ static inline int reconstructcoeff(data_lift_t dl, int32_t i, mpz_t mod_p,
 static inline int reconstructcoeff2(data_lift_t dl, int32_t i, mpz_t mod_p,
                                     rrec_data_t recdata1, rrec_data_t recdata2){
   mpz_lcm(dl->cden[i], dl->den[i], dl->pgden);
+  set_recdata_guess(recdata1, mod_p, dl->cden[i]);
   if(!lift_ratreconwden(dl->num2[i], dl->den2[i],
                         dl->crt2[i], mod_p, dl->cden[i], recdata1)){
     /* same fallback as for the first witness: balanced bounds, no guess */
@@ -1388,6 +1392,7 @@ static inline int reconstructcoeff2(data_lift_t dl, int32_t i, mpz_t mod_p,
    first two witnesses as a guess; cden[i] is updated accordingly */
 static inline int reconstructcoeff3(data_lift_t dl, int32_t i, mpz_t mod_p,
                                     rrec_data_t recdata1, rrec_data_t recdata2){
+  set_recdata_guess(recdata1, mod_p, dl->cden[i]);
   if(!lift_ratreconwden(dl->num3[i], dl->den3[i],
                         dl->crt3[i], mod_p, dl->cden[i], recdata1)){
     /* same fallback as for the first witness: balanced bounds, no guess */
@@ -1455,7 +1460,7 @@ static void ratrecon_gb(gb_modpoly_t modgbs, data_lift_t dl,
     mpz_set_ui(dl->gden, 1);
   }
 
-  set_recdata(dl, recdata1, recdata2, mod_p);
+  set_recdata(recdata2, mod_p);
 
   st = realtime();
   if(modgbs->nprimes % dl->rr == 0){