diff --git a/AsterX/interface.ccl b/AsterX/interface.ccl
index e520ca38..c729cd03 100644
--- a/AsterX/interface.ccl
+++ b/AsterX/interface.ccl
@@ -7,7 +7,7 @@ INHERITS: CarpetXRegrid ADMBaseX HydroBaseX TmunuBaseX AsterMasks
 USES INCLUDE HEADER: fixmath.hxx
 USES INCLUDE HEADER: loop_device.hxx
 USES INCLUDE HEADER: setup_eos.hxx
-USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DPalenzuela.hxx c2p_1DEntropy.hxx
+USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DRePrimAnd.hxx c2p_1DPalenzuela.hxx c2p_1DEntropy.hxx
 USES INCLUDE HEADER: reconstruct.hxx
 USES INCLUDE HEADER: aster_fd.hxx aster_interp.hxx aster_utils.hxx
 
diff --git a/AsterX/param.ccl b/AsterX/param.ccl
index 30fb5f8a..11dfb461 100644
--- a/AsterX/param.ccl
+++ b/AsterX/param.ccl
@@ -174,6 +174,8 @@ USES CCTK_REAL vw_lim
 USES BOOLEAN Ye_lenient
 USES CCTK_INT max_iter
 USES CCTK_REAL c2p_tol
+USES BOOLEAN soft_root_convergence
+USES CCTK_REAL soft_root_width_factor
 USES BOOLEAN use_press_atmo
 USES CCTK_REAL r_atmo
 USES CCTK_REAL n_rho_atmo
diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx
index 1999e484..886e55be 100644
--- a/AsterX/src/con2prim.cxx
+++ b/AsterX/src/con2prim.cxx
@@ -8,6 +8,7 @@
 #include "c2p.hxx"
 #include "c2p_1DEntropy.hxx"
 #include "c2p_1DPalenzuela.hxx"
+#include "c2p_1DRePrimAnd.hxx"
 #include "c2p_2DNoble.hxx"
 
 #include "aster_utils.hxx"
@@ -21,8 +22,8 @@ using namespace Con2PrimFactory;
 using namespace AsterUtils;
 
 enum class eos_3param { IdealGas, Hybrid, Tabulated };
-enum class c2p_first_t { None, Noble, Palenzuela, Entropy };
-enum class c2p_second_t { None, Noble, Palenzuela, Entropy };
+enum class c2p_first_t { None, Noble, RePrimAnd, Palenzuela, Entropy };
+enum class c2p_second_t { None, Noble, RePrimAnd, Palenzuela, Entropy };
 
 enum C2PFlag : CCTK_INT {
   C2P_INIT = 0,    // initial value
@@ -45,6 +46,8 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
 
   if (CCTK_EQUALS(c2p_prime, "Noble")) {
     c2p_fir = c2p_first_t::Noble;
+  } else if (CCTK_EQUALS(c2p_prime, "RePrimAnd")) {
+    c2p_fir = c2p_first_t::RePrimAnd;
   } else if (CCTK_EQUALS(c2p_prime, "Palenzuela")) {
     c2p_fir = c2p_first_t::Palenzuela;
   } else if (CCTK_EQUALS(c2p_prime, "Entropy")) {
@@ -57,6 +60,8 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
 
   if (CCTK_EQUALS(c2p_second, "Noble")) {
     c2p_sec = c2p_second_t::Noble;
+  } else if (CCTK_EQUALS(c2p_second, "RePrimAnd")) {
+    c2p_sec = c2p_second_t::RePrimAnd;
   } else if (CCTK_EQUALS(c2p_second, "Palenzuela")) {
     c2p_sec = c2p_second_t::Palenzuela;
   } else if (CCTK_EQUALS(c2p_second, "Entropy")) {
@@ -141,19 +146,29 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
     c2p_2DNoble c2p_Noble(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, vw_lim,
                           B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
                           inv_beta_max, Ye_lenient, use_z, use_temperature,
-                          use_press_atmo);
+                          use_press_atmo, soft_root_convergence,
+                          soft_root_width_factor);
+
+    // Construct RePrimAnd c2p object:
+    c2p_1DRePrimAnd c2p_RPA(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, vw_lim,
+                            B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
+                            inv_beta_max, Ye_lenient, use_z, use_temperature,
+                            use_press_atmo, soft_root_convergence,
+                            soft_root_width_factor);
 
     // Construct Palenzuela c2p object:
     c2p_1DPalenzuela c2p_Pal(eos_3p, atmo, max_iter, c2p_tol, alp_thresh,
                              vw_lim, B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
                              inv_beta_max, Ye_lenient, use_z, use_temperature,
-                             use_press_atmo);
+                             use_press_atmo, soft_root_convergence,
+                             soft_root_width_factor);
 
     // Construct Entropy c2p object:
     c2p_1DEntropy c2p_Ent(eos_3p, atmo, max_iter, c2p_tol, alp_thresh, vw_lim,
                           B_lim, rho_BH, eps_BH, vwlim_BH, sigma_max,
                           inv_beta_max, Ye_lenient, use_z, use_temperature,
-                          use_press_atmo);
+                          use_press_atmo, soft_root_convergence,
+                          soft_root_width_factor);
 
     // ----------
 
@@ -270,6 +285,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
                         rep_first);
         break;
       }
+      case c2p_first_t::RePrimAnd: {
+        c2p_RPA.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_first);
+        break;
+      }
       case c2p_first_t::Palenzuela: {
         c2p_Pal.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_first);
         break;
@@ -302,6 +321,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p,
                           rep_second);
           break;
         }
+        case c2p_second_t::RePrimAnd: {
+          c2p_RPA.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_second);
+          break;
+        }
         case c2p_second_t::Palenzuela: {
           c2p_Pal.solve(eos_3p, pv, cv, alp_avg, beta_avg, glo, rep_second);
           break;
@@ -487,6 +510,18 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTS_AsterX_Con2Prim;
   DECLARE_CCTK_PARAMETERS;
 
+  if (CCTK_EQUALS(evolution_eos, "Hybrid") && thermal_eos_atmo) {
+    CCTK_ERROR("Hybrid EOS does not implement *_from_rho_temp_ye; set "
+               "Con2PrimFactory::thermal_eos_atmo = no.");
+  }
+  if (CCTK_EQUALS(evolution_eos, "Tabulated3d") && !use_temperature) {
+    CCTK_ERROR("Tabulated3d requires Con2PrimFactory::use_temperature = yes.");
+  }
+  if (CCTK_EQUALS(evolution_eos, "Tabulated3d") && use_press_atmo) {
+    CCTK_ERROR("Tabulated3d does not support eps_from_rho_press_ye; set "
+               "Con2PrimFactory::use_press_atmo = no.");
+  }
+
   // defining EOS objects
   eos_3param eos_3p_type;
 
@@ -571,21 +606,50 @@ extern "C" void AsterX_Con2Prim_Interpolate_Failed(CCTK_ARGUMENTS) {
           const vec<CCTK_REAL, 6> vely_nbs = get_neighbors(vely, p);
           const vec<CCTK_REAL, 6> velz_nbs = get_neighbors(velz, p);
           const vec<CCTK_REAL, 6> eps_nbs = get_neighbors(eps, p);
+          const vec<CCTK_REAL, 6> press_nbs = get_neighbors(press, p);
           const vec<CCTK_REAL, 6> saved_rho_nbs = get_neighbors(saved_rho, p);
           const vec<CCTK_REAL, 6> saved_velx_nbs = get_neighbors(saved_velx, p);
           const vec<CCTK_REAL, 6> saved_vely_nbs = get_neighbors(saved_vely, p);
           const vec<CCTK_REAL, 6> saved_velz_nbs = get_neighbors(saved_velz, p);
           const vec<CCTK_REAL, 6> saved_eps_nbs = get_neighbors(saved_eps, p);
 
-          CCTK_REAL sum_nbs =
-              sum<6>([&](int i) ARITH_INLINE { return flag_nbs(i); });
-          assert(sum_nbs > 0);
-          rho(p.I) = calc_avg_neighbors(flag_nbs, rho_nbs, saved_rho_nbs);
-          velx(p.I) = calc_avg_neighbors(flag_nbs, velx_nbs, saved_velx_nbs);
-          vely(p.I) = calc_avg_neighbors(flag_nbs, vely_nbs, saved_vely_nbs);
-          velz(p.I) = calc_avg_neighbors(flag_nbs, velz_nbs, saved_velz_nbs);
-          eps(p.I) = calc_avg_neighbors(flag_nbs, eps_nbs, saved_eps_nbs);
-          press(p.I) = (gl_gamma - 1) * eps(p.I) * rho(p.I);
+          const auto good_nb = [&](int i) ARITH_INLINE -> CCTK_REAL {
+            const CCTK_INT flag_i = CCTK_INT(flag_nbs(i));
+            return ((flag_i != C2P_FAIL) && (flag_i != C2P_INIT))
+                       ? CCTK_REAL(1)
+                       : CCTK_REAL(0);
+          };
+
+          const CCTK_REAL sum_nbs =
+              sum<6>([&](int i) ARITH_INLINE { return good_nb(i); });
+          if (sum_nbs <= CCTK_REAL(0)) {
+            return;
+          }
+
+          rho(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                       return good_nb(i) * rho_nbs(i);
+                     }) /
+                     sum_nbs;
+          velx(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                        return good_nb(i) * velx_nbs(i);
+                      }) /
+                      sum_nbs;
+          vely(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                        return good_nb(i) * vely_nbs(i);
+                      }) /
+                      sum_nbs;
+          velz(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                        return good_nb(i) * velz_nbs(i);
+                      }) /
+                      sum_nbs;
+          eps(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                       return good_nb(i) * eps_nbs(i);
+                     }) /
+                     sum_nbs;
+          press(p.I) = sum<6>([&](int i) ARITH_INLINE {
+                         return good_nb(i) * press_nbs(i);
+                       }) /
+                       sum_nbs;
 
           /* reset flag */
           con2prim_flag(p.I) = C2P_AVG;
diff --git a/AsterX/src/estimate_error.cxx b/AsterX/src/estimate_error.cxx
index d8e3d49f..968097e7 100644
--- a/AsterX/src/estimate_error.cxx
+++ b/AsterX/src/estimate_error.cxx
@@ -21,7 +21,7 @@ enum class regrid_t {
 };
 regrid_t regridMethod;
 vector<int> regridVarsI;
-vector<array<int, Loop::dim>> indextypeRegridVars;
+vector<array<int, Loop::dim> > indextypeRegridVars;
 
 extern "C" void AsterX_EstimateError_Setup(CCTK_ARGUMENTS) {
   DECLARE_CCTK_PARAMETERS;
diff --git a/AsterX/src/fluxes.cxx b/AsterX/src/fluxes.cxx
index 113b703a..b337bd4a 100644
--- a/AsterX/src/fluxes.cxx
+++ b/AsterX/src/fluxes.cxx
@@ -113,19 +113,20 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType *eos_3p, const rec_var_t rec_var,
   // Prebind the velocity slice for this direction once
   const auto gf_vel_dir_i = gf_vels(dir_i);
 
-  const auto reconstruct_pt =
-      [=] CCTK_DEVICE(const GF3D2<const CCTK_REAL> &var, const PointDesc &p,
-                      bool gf_is_rho, bool gf_is_press) {
-        return reconstruct<vec<CCTK_REAL, 2>>(var, p, reconstruction, dir_i,
-                                              gf_is_rho, gf_is_press, press,
-                                              gf_vel_dir_i, reconstruct_params);
-      };
+  const auto reconstruct_pt = [=] CCTK_DEVICE(const GF3D2<const CCTK_REAL> &var,
+                                              const PointDesc &p,
+                                              bool gf_is_rho,
+                                              bool gf_is_press) {
+    return reconstruct<vec<CCTK_REAL, 2> >(var, p, reconstruction, dir_i,
+                                           gf_is_rho, gf_is_press, press,
+                                           gf_vel_dir_i, reconstruct_params);
+  };
   const auto reconstruct_loworder =
       [=] CCTK_DEVICE(const GF3D2<const CCTK_REAL> &var, const PointDesc &p,
                       bool gf_is_rho, bool gf_is_press) {
-        return reconstruct<vec<CCTK_REAL, 2>>(var, p, reconstruction_LO, dir_i,
-                                              gf_is_rho, gf_is_press, press,
-                                              gf_vel_dir_i, reconstruct_params);
+        return reconstruct<vec<CCTK_REAL, 2> >(
+            var, p, reconstruction_LO, dir_i, gf_is_rho, gf_is_press, press,
+            gf_vel_dir_i, reconstruct_params);
       };
 
   const auto calcflux =
diff --git a/Con2PrimFactory/interface.ccl b/Con2PrimFactory/interface.ccl
index 1ae03a77..35e41db4 100644
--- a/Con2PrimFactory/interface.ccl
+++ b/Con2PrimFactory/interface.ccl
@@ -8,6 +8,7 @@ INCLUDES HEADER: cons.hxx IN cons.hxx
 INCLUDES HEADER: c2p_2DNoble.hxx IN c2p_2DNoble.hxx 
 INCLUDES HEADER: c2p_1DPalenzuela.hxx IN c2p_1DPalenzuela.hxx 
 INCLUDES HEADER: c2p_1DEntropy.hxx IN c2p_1DEntropy.hxx 
+INCLUDES HEADER: c2p_1DRePrimAnd.hxx IN c2p_1DRePrimAnd.hxx
 USES INCLUDE HEADER: setup_eos.hxx 
 USES INCLUDE HEADER: roots.hxx
 USES INCLUDE HEADER: aster_utils.hxx
diff --git a/Con2PrimFactory/param.ccl b/Con2PrimFactory/param.ccl
index acfb7609..cf35d95a 100644
--- a/Con2PrimFactory/param.ccl
+++ b/Con2PrimFactory/param.ccl
@@ -6,6 +6,7 @@ no
 restricted :
 
 KEYWORD c2p_prime "Name of the main con2prim scheme" STEERABLE=ALWAYS {
+  "RePrimAnd" :: "RePrimAnd"
   "Noble" ::"Noble"
   "Palenzuela" ::"Palenzuela"
   "Entropy" ::"Entropy"
@@ -14,6 +15,7 @@ KEYWORD c2p_prime "Name of the main con2prim scheme" STEERABLE=ALWAYS {
 "Noble"
 
 KEYWORD c2p_second "Name of the backup con2prim scheme" STEERABLE=ALWAYS {
+  "RePrimAnd" :: "RePrimAnd"
   "Noble" ::"Noble"
   "Palenzuela" ::"Palenzuela"
   "Entropy" ::"Entropy"
@@ -52,6 +54,7 @@ CCTK_REAL B_lim "Upper limit for the value of magnetization" STEERABLE=ALWAYS
   0:            :: "Must be positive"
 }  100.0
 
+# Safe default for C2P robustness
 CCTK_REAL vw_lim "Upper limit for the value of velocity * lorentz_factor " STEERABLE=ALWAYS
 {
   0:            :: "Must be positive"
@@ -71,6 +74,19 @@ CCTK_REAL c2p_tol "c2p torelance for root finding" STEERABLE=ALWAYS
   0:* :: "Larger than zero"
 } 1e-10
 
+BOOLEAN soft_root_convergence
+"Allow near-converged root solutions to be accepted as success (warning-only)."
+STEERABLE=ALWAYS
+{
+} "no"
+
+CCTK_REAL soft_root_width_factor
+"Relative bracket-width factor (x c2p_tol) used for near-converged root acceptance. Ignored unless soft_root_convergence=yes."
+STEERABLE=ALWAYS
+{
+  1.0:* :: "Must be >= 1"
+} 1.0
+
 CCTK_REAL tauFluid_atmo "Minimum conserved internal energy of fluid" STEERABLE=ALWAYS
 {
   [0.0:* :: "Larger than zero"
diff --git a/Con2PrimFactory/src/c2p.hxx b/Con2PrimFactory/src/c2p.hxx
index c14e9e7d..e92a4081 100644
--- a/Con2PrimFactory/src/c2p.hxx
+++ b/Con2PrimFactory/src/c2p.hxx
@@ -51,6 +51,8 @@ protected:
   bool use_zprim;
   bool use_temp;
   bool use_press_atmo;
+  bool soft_root_convergence{false};
+  CCTK_REAL soft_root_width_factor{1.0};
 
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
   get_Ssq_Exact(const vec<CCTK_REAL, 3> &mom,
@@ -447,9 +449,12 @@ c2p::cons_floors_and_ceilings(const EOSType *eos_3p, cons_vars &cv,
   // Compute Bsq
   vec<CCTK_REAL, 3> B_low = calc_contraction(glo, cv.dBvec);
   const CCTK_REAL BsqL = calc_contraction(B_low, cv.dBvec);
+  //const CCTK_REAL tauF_atmo =
+  //    std::max(cv.dens * atmo.eps_atmo, sqrt_detg * tauFluid_atmo);
   const CCTK_REAL tau_lim = 0.5 * BsqL / sqrt_detg;
 
   if (cv.tau <= tau_lim) {
+    //cv.tau = tau_lim + tauF_atmo;
     cv.tau = tau_lim + sqrt_detg * tauFluid_atmo;
   }
 
diff --git a/Con2PrimFactory/src/c2p_1DEntropy.hxx b/Con2PrimFactory/src/c2p_1DEntropy.hxx
index 1a2b46b3..558d15dd 100644
--- a/Con2PrimFactory/src/c2p_1DEntropy.hxx
+++ b/Con2PrimFactory/src/c2p_1DEntropy.hxx
@@ -17,7 +17,8 @@ public:
       CCTK_REAL tol, CCTK_REAL alp_thresh_in, CCTK_REAL vwlim, CCTK_REAL B_lim,
       CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in,
       CCTK_REAL sigma_max_in, CCTK_REAL inv_beta_max_in, bool ye_len,
-      bool use_z, bool use_temperature, bool use_pressure_atmo);
+      bool use_z, bool use_temperature, bool use_pressure_atmo,
+      bool soft_root_conv, CCTK_REAL soft_root_width_factor_in);
 
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
   get_Ssq_Exact(const vec<CCTK_REAL, 3> &mom,
@@ -66,7 +67,8 @@ CCTK_HOST CCTK_DEVICE
         CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in,
         CCTK_REAL vwlim_BH_in, CCTK_REAL sigma_max_in,
         CCTK_REAL inv_beta_max_in, bool ye_len, bool use_z,
-        bool use_temperature, bool use_pressure_atmo) {
+        bool use_temperature, bool use_pressure_atmo, bool soft_root_conv,
+        CCTK_REAL soft_root_width_factor_in) {
 
   // Base
   atmo = atm;
@@ -86,6 +88,8 @@ CCTK_HOST CCTK_DEVICE
   use_zprim = use_z;
   use_temp = use_temperature;
   use_press_atmo = use_pressure_atmo;
+  soft_root_convergence = soft_root_conv;
+  soft_root_width_factor = fmax(CCTK_REAL(1.0), soft_root_width_factor_in);
 
   // Derived
   GammaIdealFluid = eos_3p->gamma;
@@ -270,10 +274,11 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
                      const CCTK_REAL alp, const vec<CCTK_REAL, 3> &beta,
                      const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const {
 
-  // ROOTSTAT status = ROOTSTAT::SUCCESS;
+  ROOTSTAT status = ROOTSTAT::SUCCESS;
   rep.iters = 0;
   rep.adjust_cons = false;
   rep.set_atmo = false;
+  rep.soft_root_conv = false;
   rep.status = c2p_report::SUCCESS;
 
   /* Check validity of the 3-metric and compute its inverse */
@@ -361,6 +366,12 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
   // const CCTK_INT minbits = std::numeric_limits<CCTK_REAL>::digits - 4;
   const CCTK_INT maxiters = maxIterations;
 
+  const CCTK_REAL f_a0 = fn(a);
+  const CCTK_REAL f_b0 = fn(b);
+  if ((!isfinite(f_a0)) || (!isfinite(f_b0)) || (f_a0 * f_b0 > 0.0)) {
+    status = ROOTSTAT::NOT_BRACKETED;
+  }
+
   auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters);
 
   CCTK_REAL xEntropy_Sol = 0.5 * (result.first + result.second);
@@ -389,14 +400,32 @@ c2p_1DEntropy::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
     return;
   }
 
-  if (abs(result.first - result.second) >
-      tolerance_0 * min(abs(result.first), abs(result.second))) {
-
-    // set status to root not converged
-    rep.set_root_conv();
-    cv = cv_const;
-    // status = ROOTSTAT::NOT_CONVERGED;
-    return;
+  const CCTK_REAL root_width = abs(result.first - result.second);
+  const CCTK_REAL strict_width_tol =
+      tolerance_0 * min(abs(result.first), abs(result.second));
+  if (root_width > strict_width_tol) {
+    bool accept_soft = false;
+    if (soft_root_convergence) {
+      const CCTK_REAL scale =
+          fmax(CCTK_REAL(0.0), fmax(abs(result.first), abs(result.second)));
+      const CCTK_REAL soft_width_tol =
+          soft_root_width_factor * tolerance * scale;
+      accept_soft = std::isfinite(root_width) && std::isfinite(soft_width_tol) &&
+                    (root_width <= soft_width_tol);
+    }
+    if (!accept_soft) {
+      // set status to root not converged / not bracketed
+      if (status == ROOTSTAT::NOT_BRACKETED) {
+        rep.set_root_bracket();
+      } else {
+        rep.set_root_conv();
+        status = ROOTSTAT::NOT_CONVERGED;
+      }
+      cv = cv_const;
+      (void)status;
+      return;
+    }
+    rep.set_soft_root_conv();
   }
 
   // set to atmo if computed rho is below floor density
diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx
index 66d64f10..f1cd37bd 100644
--- a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx
+++ b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx
@@ -18,7 +18,8 @@ public:
       CCTK_REAL tol, CCTK_REAL alp_thresh_in, CCTK_REAL vwlim, CCTK_REAL B_lim,
       CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in,
       CCTK_REAL sigma_max_in, CCTK_REAL inv_beta_max_in, bool ye_len,
-      bool use_z, bool use_temperature, bool use_pressure_atmo);
+      bool use_z, bool use_temperature, bool use_pressure_atmo,
+      bool soft_root_conv, CCTK_REAL soft_root_width_factor_in);
 
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
   get_Ssq_Exact(const vec<CCTK_REAL, 3> &mom,
@@ -68,7 +69,8 @@ CCTK_HOST CCTK_DEVICE
         CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in,
         CCTK_REAL vwlim_BH_in, CCTK_REAL sigma_max_in,
         CCTK_REAL inv_beta_max_in, bool ye_len, bool use_z,
-        bool use_temperature, bool use_pressure_atmo) {
+        bool use_temperature, bool use_pressure_atmo, bool soft_root_conv,
+        CCTK_REAL soft_root_width_factor_in) {
 
   // Base
   atmo = atm;
@@ -88,6 +90,8 @@ CCTK_HOST CCTK_DEVICE
   use_zprim = use_z;
   use_temp = use_temperature;
   use_press_atmo = use_pressure_atmo;
+  soft_root_convergence = soft_root_conv;
+  soft_root_width_factor = fmax(CCTK_REAL(1.0), soft_root_width_factor_in);
 
   // Derived
   GammaIdealFluid = eos_3p->gamma;
@@ -311,10 +315,11 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
                         const CCTK_REAL alp, const vec<CCTK_REAL, 3> &beta,
                         const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const {
 
-  // ROOTSTAT status = ROOTSTAT::SUCCESS;
+  ROOTSTAT status = ROOTSTAT::SUCCESS;
   rep.iters = 0;
   rep.adjust_cons = false;
   rep.set_atmo = false;
+  rep.soft_root_conv = false;
   rep.status = c2p_report::SUCCESS;
 
   /* Check validity of the 3-metric and compute its inverse */
@@ -379,6 +384,13 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
     return;
   }
 
+  // Clamp Ye in the local working conservatives so the root solve and prim
+  // reconstruction both use EOS-valid Ye.
+  CCTK_REAL Ye_raw = cv.DYe / cv.dens;
+  const CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max);
+  const bool ye_clipped = (Ye_raw < eos_3p->rgye.min) || (Ye_raw > eos_3p->rgye.max);
+  cv.DYe = cv.dens * Ye;
+
   // Find x, this is the recovery process
   // const CCTK_INT minbits = std::numeric_limits<CCTK_REAL>::digits - 4;
   // const CCTK_INT maxiters = maxIterations;
@@ -420,34 +432,30 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
     //}
   }
 
-  auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters);
-
-  // hybrid: prefer endpoint with smaller |f|, else midpoint
-  CCTK_REAL a_root = result.first;
-  CCTK_REAL b_root = result.second;
-  CCTK_REAL fa = fn(a_root);
-  CCTK_REAL fb = fn(b_root);
-
-  CCTK_REAL xPalenzuela_Sol;
-  if (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) {
-    // exact root or smaller residual at b
-    xPalenzuela_Sol = b_root;
-  } else if (std::abs(fa) < std::abs(fb)) {
-    // smaller residual at a
-    xPalenzuela_Sol = a_root;
-  } else {
-    // fall back to midpoint
-    xPalenzuela_Sol = CCTK_REAL(0.5) * (a_root + b_root);
+  const CCTK_REAL f_a0 = fn(a);
+  const CCTK_REAL f_b0 = fn(b);
+  if ((!isfinite(f_a0)) || (!isfinite(f_b0)) || (f_a0 * f_b0 > 0.0)) {
+    status = ROOTSTAT::NOT_BRACKETED;
   }
 
-  // for now, we do not use the above hybrid scheme
-  xPalenzuela_Sol = 0.5 * (result.first + result.second);
+  auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters);
 
-  // if (abs(fn(result.first)) < abs(fn(result.second))) {
-  //  xPalenzuela_Sol = result.first;
-  //} else {
-  //  xPalenzuela_Sol = result.second;
-  //}
+  // Legacy endpoint-preference selector kept for reference; below we use the
+  // midpoint rule for xPalenzuela_Sol.
+  // CCTK_REAL a_root = result.first;
+  // CCTK_REAL b_root = result.second;
+  // CCTK_REAL fa = fn(a_root);
+  // CCTK_REAL fb = fn(b_root);
+  // CCTK_REAL xPalenzuela_Sol;
+  // if (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) {
+  //   xPalenzuela_Sol = b_root;
+  // } else if (std::abs(fa) < std::abs(fb)) {
+  //   xPalenzuela_Sol = a_root;
+  // } else {
+  //   xPalenzuela_Sol = CCTK_REAL(0.5) * (a_root + b_root);
+  // }
+
+  CCTK_REAL xPalenzuela_Sol = CCTK_REAL(0.5) * (result.first + result.second);
 
   xPalenzuelaToPrim(xPalenzuela_Sol, Ssq, Bsq, BiSi, eos_3p, pv, cv, gup, glo);
 
@@ -467,6 +475,10 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
     return;
   }
 
+  if (ye_clipped) {
+    rep.adjust_cons = true;
+  }
+
   // General comment:
   // One could think of expressing the following condition
   // in a way that is "safe" against NaNs and infs. First, this only makes
@@ -483,14 +495,47 @@ c2p_1DPalenzuela::solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
 
   // TODO: have an explicit check on max_iters, e.g.:
   // if (rep.iters >= maxiters || abs(fn(xPalenzuela_Sol)) > tolerance) {
-  if (abs(result.first - result.second) >
-      tolerance_0 * min(abs(result.first), abs(result.second))) {
+  const CCTK_REAL root_width = abs(result.first - result.second);
+  const CCTK_REAL strict_width_tol =
+      tolerance_0 * min(abs(result.first), abs(result.second));
+  if (root_width > strict_width_tol) {
+    bool accept_soft = false;
+    if (soft_root_convergence) {
+      const CCTK_REAL scale =
+          fmax(CCTK_REAL(0.0), fmax(abs(result.first), abs(result.second)));
+      const CCTK_REAL soft_width_tol =
+          soft_root_width_factor * tolerance * scale;
+      accept_soft = std::isfinite(root_width) && std::isfinite(soft_width_tol) &&
+                    (root_width <= soft_width_tol);
+    }
+    if (!accept_soft) {
+      // set status to root not converged / not bracketed
+      if (status == ROOTSTAT::NOT_BRACKETED) {
+        rep.set_root_bracket();
+      } else {
+        rep.set_root_conv();
+        status = ROOTSTAT::NOT_CONVERGED;
+      }
+      cv = cv_const;
+      (void)status;
+      return;
+    }
+    rep.set_soft_root_conv();
+  }
 
-    // set status to root not converged
-    rep.set_root_conv();
-    cv = cv_const;
-    // status = ROOTSTAT::NOT_CONVERGED;
-    return;
+  const vec<CCTK_REAL, 3> v_low_pre = calc_contraction(glo, pv.vel);
+  const CCTK_REAL vsq_pre = calc_contraction(v_low_pre, pv.vel);
+  const CCTK_REAL sol_v = std::sqrt(std::max(CCTK_REAL(0), vsq_pre));
+  if (sol_v > v_lim) {
+    pv.rho = cv.dens / w_lim;
+    pv.vel *= v_lim / (sol_v + CCTK_REAL(1e-300));
+    pv.w_lor = w_lim;
+    const auto rgeps_lim = eos_3p->range_eps_from_rho_ye(pv.rho, pv.Ye);
+    pv.eps = fmin(fmax(rgeps_lim.min, pv.eps), rgeps_lim.max);
+    pv.press = eos_3p->press_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+    pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+    pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+    rep.adjust_cons = true;
   }
 
   // set to atmo if computed rho is below floor density
diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx
new file mode 100644
index 00000000..d791b7a3
--- /dev/null
+++ b/Con2PrimFactory/src/c2p_1DRePrimAnd.hxx
@@ -0,0 +1,318 @@
+#ifndef C2P_1DREPRIMAND_HXX
+#define C2P_1DREPRIMAND_HXX
+
+#include <algorithm>
+#include <cmath>
+
+#include "c2p.hxx"
+#include "c2p_report.hxx"
+#include "prims.hxx"
+#include "cons.hxx"
+#include "roots.hxx"
+
+#include "c2p_1DRePrimAnd_intervals.hxx"
+#include "c2p_1DRePrimAnd_rootfinder.hxx"
+#include "c2p_1DRePrimAnd_internals.hxx"
+
+namespace Con2PrimFactory {
+
+class c2p_1DRePrimAnd : public c2p {
+public:
+  CCTK_REAL GammaIdealFluid;
+
+  template <typename EOSType>
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DRePrimAnd(
+      const EOSType *eos_3p, const atmosphere &atm_in, CCTK_INT maxIter,
+      CCTK_REAL tol, CCTK_REAL alp_thresh_in, CCTK_REAL vwlim,
+      CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in,
+      CCTK_REAL vwlim_BH_in, CCTK_REAL sigma_max_in,
+      CCTK_REAL inv_beta_max_in, bool ye_len, bool use_z,
+      bool use_temperature, bool use_pressure_atmo, bool soft_root_conv,
+      CCTK_REAL soft_root_width_factor_in) {
+
+    atmo = atm_in;
+    maxIterations = maxIter;
+    tolerance = tol;
+    alp_thresh = alp_thresh_in;
+    vw_lim = vwlim;
+    w_lim = sqrt(1.0 + vw_lim * vw_lim);
+    v_lim = vw_lim / w_lim;
+    Bsq_lim = B_lim * B_lim;
+    rho_BH = rho_BH_in;
+    eps_BH = eps_BH_in;
+    vwlim_BH = vwlim_BH_in;
+    ye_lenient = ye_len;
+    use_zprim = use_z;
+    use_temp = use_temperature;
+    use_press_atmo = use_pressure_atmo;
+    soft_root_convergence = soft_root_conv;
+    soft_root_width_factor = fmax(CCTK_REAL(1.0), soft_root_width_factor_in);
+    sigma_max = sigma_max_in;
+    inv_beta_max = inv_beta_max_in;
+
+    GammaIdealFluid = eos_3p->gamma;
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
+  set_to_nan(prim_vars &pv, cons_vars &cv) const {
+    pv.set_to_nan();
+    cv.set_to_nan();
+  }
+
+  template <typename EOSType>
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
+  solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
+        const CCTK_REAL alp, const vec<CCTK_REAL, 3> &beta,
+        const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const {
+    rep.iters = 0;
+    rep.adjust_cons = false;
+    rep.set_atmo = false;
+    rep.soft_root_conv = false;
+    rep.status = c2p_report::SUCCESS;
+
+    const CCTK_REAL spatial_detg = calc_det(glo);
+    const CCTK_REAL sqrt_detg = sqrt(spatial_detg);
+    const bool minor1{glo(X, X) > 0.0};
+    const bool minor2{glo(X, X) * glo(Y, Y) - glo(X, Y) * glo(X, Y) > 0.0};
+    const bool minor3{spatial_detg > 0.0};
+    if (!(minor1 && minor2 && minor3)) {
+      rep.set_invalid_detg(sqrt_detg);
+      set_to_nan(pv, cv);
+      return;
+    }
+
+    const smat<CCTK_REAL, 3> gup = calc_inv(glo, spatial_detg);
+
+    const cons_vars cv_const = cv;
+
+    cv.dens /= sqrt_detg;
+    cv.tau /= sqrt_detg;
+    cv.mom /= sqrt_detg;
+    cv.dBvec /= sqrt_detg;
+    cv.DYe /= sqrt_detg;
+    cv.DEnt /= sqrt_detg;
+
+    const CCTK_REAL Ssq =
+        calc_contraction(calc_contraction(gup, cv.mom), cv.mom);
+    const CCTK_REAL Bsq =
+        calc_contraction(calc_contraction(glo, cv.dBvec), cv.dBvec);
+    const CCTK_REAL BiSi = cv.dBvec(X) * cv.mom(X) + cv.dBvec(Y) * cv.mom(Y) +
+                           cv.dBvec(Z) * cv.mom(Z);
+
+    if ((!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) ||
+        (!isfinite(BiSi)) || (!isfinite(cv.DYe)) || (!isfinite(cv.DEnt))) {
+      rep.set_nans_in_cons(cv.dens, Ssq, Bsq, BiSi, cv.DYe);
+      set_to_nan(pv, cv);
+      return;
+    }
+
+    if (Bsq > Bsq_lim) {
+      rep.set_B_limit(Bsq);
+      set_to_nan(pv, cv);
+      return;
+    }
+
+    // Ye needs to be a non-const lvalue for EOS calls
+    CCTK_REAL Ye_raw = cv.DYe / cv.dens;
+    const CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max);
+    const bool ye_clipped =
+        (Ye_raw < eos_3p->rgye.min) || (Ye_raw > eos_3p->rgye.max);
+
+    typename RePrimAnd::froot<EOSType>::cache cache{};
+    RePrimAnd::froot<EOSType> f(
+        eos_3p, Ye, cv.dens, cv.tau / cv.dens, Ssq / (cv.dens * cv.dens),
+        (BiSi * BiSi) / (cv.dens * cv.dens * cv.dens), Bsq / cv.dens, cache);
+
+    ROOTSTAT prep_stat{};
+    interval<CCTK_REAL> mu_br = f.initial_bracket(prep_stat);
+    if (prep_stat != ROOTSTAT::SUCCESS) {
+      if (prep_stat == ROOTSTAT::NOT_BRACKETED)
+        rep.set_prep_root_bracket();
+      else
+        rep.set_prep_root_conv();
+      cv = cv_const;
+      return;
+    }
+
+    interval<CCTK_REAL> rgrho_iv{eos_3p->rgrho.min, eos_3p->rgrho.max};
+    RePrimAnd::rarecase<EOSType> rc(mu_br, rgrho_iv, f);
+    if (rc.rho_too_big) {
+      rep.set_range_rho(cv.dens, cv.dens);
+      cv = cv_const;
+      return;
+    }
+    if (rc.rho_too_small) {
+      rep.set_atmo_set();
+      pv.Bvec = cv.dBvec;
+      atmo.set(pv, cv, glo);
+      return;
+    }
+    mu_br = rc.bracket;
+
+    const CCTK_INT maxiters = maxIterations;
+
+    ROOTSTAT root_stat{};
+    interval<CCTK_REAL> result =
+        findroot_no_deriv(f, mu_br, tolerance, (unsigned int)maxiters, root_stat);
+    rep.iters = (CCTK_INT)cache.calls;
+
+    if (root_stat == ROOTSTAT::NOT_CONVERGED) {
+      bool accept_soft = false;
+      if (soft_root_convergence) {
+        const CCTK_REAL a = fmin(result.min(), result.max());
+        const CCTK_REAL b = fmax(result.min(), result.max());
+        if (std::isfinite(a) && std::isfinite(b) && b >= a) {
+          const CCTK_REAL width = b - a;
+          const CCTK_REAL scale =
+              fmax(CCTK_REAL(0.0), fmax(std::abs(a), std::abs(b)));
+          const CCTK_REAL width_tol =
+              soft_root_width_factor * tolerance * scale;
+          const CCTK_REAL mu_ref =
+              fmax(CCTK_REAL(0.5) * (a + b), CCTK_REAL(1.0e-300));
+          const bool stopif_ok = f.stopif(mu_ref, width, tolerance);
+          accept_soft = (width <= width_tol) || stopif_ok;
+        }
+      }
+      if (!accept_soft) {
+        rep.set_root_conv();
+        cv = cv_const;
+        return;
+      }
+      rep.set_soft_root_conv();
+      root_stat = ROOTSTAT::SUCCESS;
+    }
+
+    if (root_stat != ROOTSTAT::SUCCESS) {
+      if (rc.rho_big) {
+        rep.set_range_rho(cv.dens, cv.dens);
+        cv = cv_const;
+        return;
+      }
+      if (rc.rho_small) {
+        rep.set_atmo_set();
+        pv.Bvec = cv.dBvec;
+        atmo.set(pv, cv, glo);
+        return;
+      }
+      rep.set_root_bracket();
+      cv = cv_const;
+      return;
+    }
+
+    CCTK_REAL mu = cache.lmu;
+    const CCTK_REAL mu_lo = fmin(result.min(), result.max());
+    const CCTK_REAL mu_hi = fmax(result.min(), result.max());
+
+    if (!std::isfinite(mu) || mu < mu_lo || mu > mu_hi) {
+      auto fn = [&](CCTK_REAL x) { return f(x); };
+      const CCTK_REAL a_root = result.min();
+      const CCTK_REAL b_root = result.max();
+      const CCTK_REAL fa = fn(a_root);
+      const CCTK_REAL fb = fn(b_root);
+      mu = (fb == (CCTK_REAL)0 || std::abs(fb) < std::abs(fa)) ? b_root
+           : (std::abs(fa) < std::abs(fb))                     ? a_root
+                                     : CCTK_REAL(0.5) * (a_root + b_root);
+
+      // Ensure cache corresponds to the final chosen mu.
+      (void)f(mu);
+    }
+
+    // ------------------------------------------------------------------
+    // Use the EOS-consistent RePrimAnd cached primitives
+    // ------------------------------------------------------------------
+    pv.rho = cache.rho;
+    pv.Ye = cache.ye;
+    pv.eps = cache.eps;
+    pv.press = cache.press;
+    pv.w_lor = cache.w;
+
+    // If root-cache thermodynamics were clipped to EOS bounds, force
+    // conservative recomputation for consistency with adjusted primitives.
+    if (ye_clipped) {
+      rep.adjust_cons = true;
+    }
+    {
+      const auto rgeps = eos_3p->range_eps_from_rho_ye(pv.rho, pv.Ye);
+      if (cache.eps_raw < rgeps.min || cache.eps_raw > rgeps.max) {
+        rep.adjust_cons = true;
+      }
+    }
+
+    pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+    pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+
+    // ------------------------------------------------------------------
+    // Velocity: RePrimAnd magnetic-aware reconstruction
+    //   v^i = mu * x * ( r^i + (r.b) * mu * b^i )
+    // where r^i = S^i / D, b^i = B^i / sqrt(D), (r.b) = (B^i S_i)/(D^(3/2)).
+    // ------------------------------------------------------------------
+
+    const vec<CCTK_REAL, 3> mom_up = calc_contraction(gup, cv.mom);
+
+    const CCTK_REAL D = cv.dens;
+    const CCTK_REAL sqD = std::sqrt(D);
+
+    const vec<CCTK_REAL, 3> r_u = mom_up * (1.0 / (D + 1e-300));
+    const vec<CCTK_REAL, 3> b_u = cv.dBvec * (1.0 / (sqD + 1e-300));
+    const CCTK_REAL rb = BiSi / (D * (sqD + 1e-300));
+
+    pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u);
+
+    const CCTK_REAL sol_v = std::sqrt(std::max(CCTK_REAL(0), cache.vsqr));
+    if (sol_v > v_lim) {
+      pv.rho = cv.dens / w_lim;
+      pv.vel *= v_lim / (sol_v + CCTK_REAL(1e-300));
+      pv.w_lor = w_lim;
+      const auto rgeps_lim = eos_3p->range_eps_from_rho_ye(pv.rho, pv.Ye);
+      pv.eps = fmin(fmax(rgeps_lim.min, pv.eps), rgeps_lim.max);
+      pv.press = eos_3p->press_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+      pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+      pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+      rep.adjust_cons = true;
+    } else {
+      pv.w_lor = cache.w;
+    }
+
+    pv.Bvec = cv.dBvec;
+    const vec<CCTK_REAL, 3> Elow = calc_cross_product(pv.Bvec, pv.vel);
+    pv.E = calc_contraction(gup, Elow);
+
+    if (pv.rho < atmo.rho_cut) {
+      rep.set_atmo_set();
+      atmo.set(pv, cv, glo);
+      return;
+    }
+
+    c2p::prims_floors_and_ceilings(eos_3p, pv, cv, alp, beta, glo, rep);
+
+    if (rep.adjust_cons) {
+      cv.from_prim(pv, glo);
+      cv.dBvec = cv_const.dBvec;
+    } else {
+      cv = cv_const;
+      cv.DEnt = cv.dens * pv.entropy;
+    }
+  }
+
+  template <typename EOSType>
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
+  solve(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
+        const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const {
+    // Backward-compatible overload for legacy call sites.
+    const CCTK_REAL alp = 1.0;
+    const vec<CCTK_REAL, 3> beta{0.0, 0.0, 0.0};
+    solve(eos_3p, pv, cv, alp, beta, glo, rep);
+  }
+
+  template <typename EOSType>
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
+  operator()(const EOSType *eos_3p, prim_vars &pv, cons_vars &cv,
+             const smat<CCTK_REAL, 3> & /*gup_unused*/,
+             const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const {
+    solve(eos_3p, pv, cv, glo, rep);
+  }
+};
+
+} // namespace Con2PrimFactory
+
+#endif
diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx
new file mode 100644
index 00000000..a00e0f9d
--- /dev/null
+++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_internals.hxx
@@ -0,0 +1,309 @@
+#ifndef C2P_1DREPRIMAND_INTERNALS_HXX
+#define C2P_1DREPRIMAND_INTERNALS_HXX
+
+#include <cmath>
+#include <limits>
+
+#include "prims.hxx"
+#include "cons.hxx"
+#include "c2p_1DRePrimAnd_intervals.hxx"
+#include "c2p_1DRePrimAnd_rootfinder.hxx"
+
+namespace Con2PrimFactory {
+namespace RePrimAnd {
+
+// --------------------------- f_upper ---------------------------------
+
+class f_upper {
+public:
+  using value_t = CCTK_REAL;
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE
+  f_upper(CCTK_REAL h0_, CCTK_REAL rsqr_, CCTK_REAL rbsqr_, CCTK_REAL bsqr_)
+      : h0(h0_), h0sqr(h0_ * h0_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_) {}
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
+  x_from_mu(CCTK_REAL mu) const {
+    return 1.0 / (1.0 + mu * bsqr);
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
+  rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const {
+    return x * (rsqr * x + mu * (x + 1.0) * rbsqr);
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
+  new_h0w_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const {
+    return std::sqrt(h0sqr + rfsqr_from_mu_x(mu, x));
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE auto
+  operator()(CCTK_REAL mu) const -> std::pair<CCTK_REAL, CCTK_REAL> {
+    const CCTK_REAL x = x_from_mu(mu);
+    const CCTK_REAL xsq = x * x;
+    const CCTK_REAL hw = new_h0w_from_mu_x(mu, x);
+    const CCTK_REAL b = x * (xsq * rsqr + mu * (1.0 + x + xsq) * rbsqr);
+    const CCTK_REAL f = mu * hw - 1.0;
+    const CCTK_REAL df = (h0sqr + b) / hw;
+    return {f, df};
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval<CCTK_REAL>
+  initial_bracket() const {
+    CCTK_REAL mu_min = 1.0 / std::sqrt(h0sqr + rsqr);
+    const CCTK_REAL mu0 = 1.0 / h0;
+    const CCTK_REAL rfs_min = rfsqr_from_mu_x(mu0, x_from_mu(mu0));
+    CCTK_REAL mu_max = 1.0 / std::sqrt(h0sqr + rfs_min);
+    const CCTK_REAL margin = 10 * std::numeric_limits<CCTK_REAL>::epsilon();
+    mu_max *= 1.0 + margin;
+    mu_min *= 1.0 - margin;
+    if (mu_max <= mu_min) {
+      mu_min = 0.0;
+      mu_max = mu0 * (1.0 + margin);
+    }
+    return {mu_min, mu_max};
+  }
+
+  CCTK_REAL h0, h0sqr, rsqr, rbsqr, bsqr;
+};
+
+// --------------------------- froot -----------------------------------
+
+template <typename EOSType> class froot {
+public:
+  using value_t = CCTK_REAL;
+
+  struct cache {
+    CCTK_REAL ye{0}, lmu{0}, x{0};
+    CCTK_REAL rho{0}, rho_raw{0};
+    CCTK_REAL eps{0}, eps_raw{0};
+    CCTK_REAL etot{0}, etot_raw{0};
+    CCTK_REAL press{0};
+    CCTK_REAL vsqr{0};
+    CCTK_REAL w{1};
+    unsigned int calls{0};
+  };
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE
+  froot(const EOSType *eos_, CCTK_REAL valid_ye, CCTK_REAL d_, CCTK_REAL qtot_,
+        CCTK_REAL rsqr_, CCTK_REAL rbsqr_, CCTK_REAL bsqr_, cache &last_)
+      : eos(eos_), d(d_), qtot(qtot_), rsqr(rsqr_), rbsqr(rbsqr_), bsqr(bsqr_),
+        brosqr(rsqr_ * bsqr_ - rbsqr_), last(last_) {
+    last.ye = valid_ye;
+    last.calls = 0;
+
+    // Build h0 from EOS min bounds
+    const CCTK_REAL rhomin = eos->rgrho.min;
+    const auto rgeps0 = eos->range_eps_from_rho_ye(rhomin, valid_ye);
+    CCTK_REAL epsmin = rgeps0.min;
+    const CCTK_REAL pmin =
+        eos->press_from_rho_eps_ye(rhomin, epsmin, valid_ye);
+    h0 = 1.0 + epsmin + pmin / rhomin;
+    h0 = fmax(h0, CCTK_REAL(1.0) + std::numeric_limits<CCTK_REAL>::epsilon());
+
+    const CCTK_REAL zsqrinf = rsqr / (h0 * h0);
+    const CCTK_REAL wsqrinf = 1.0 + zsqrinf;
+    winf = std::sqrt(wsqrinf);
+    vsqrinf = zsqrinf / wsqrinf;
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
+  operator()(CCTK_REAL mu) {
+    cache &c = last;
+    c.lmu = mu;
+    c.x = x_from_mu(mu);
+
+    const CCTK_REAL rfsqr = rfsqr_from_mu_x(mu, c.x);
+    const CCTK_REAL qf = qf_from_mu_x(mu, c.x);
+
+    c.vsqr = rfsqr * mu * mu;
+    if (c.vsqr >= vsqrinf) {
+      c.vsqr = vsqrinf;
+      c.w = winf;
+    } else {
+      c.w = 1.0 / std::sqrt(1.0 - c.vsqr);
+    }
+
+    c.rho_raw = d / c.w;
+    c.rho = fmin(fmax(eos->rgrho.min, c.rho_raw), eos->rgrho.max);
+
+    c.eps_raw = get_eps_raw(mu, qf, rfsqr, c.w);
+    const auto rgeps = eos->range_eps_from_rho_ye(c.rho, c.ye);
+    c.eps = fmin(fmax(rgeps.min, c.eps_raw), rgeps.max);
+
+    c.etot_raw = c.rho_raw * (1.0 + c.eps_raw);
+    c.etot = c.rho * (1.0 + c.eps);
+
+    c.press = eos->press_from_rho_eps_ye(c.rho, c.eps, c.ye);
+    ++c.calls;
+
+    const CCTK_REAL inv_etot =
+        CCTK_REAL(1.0) /
+        fmax(c.etot, std::numeric_limits<CCTK_REAL>::min());
+    const CCTK_REAL one_plus_a = (c.etot + c.press) * inv_etot;
+    const CCTK_REAL h = (c.etot + c.press) / c.rho;
+
+    const CCTK_REAL hbw_raw = one_plus_a * (1.0 + qf - mu * rfsqr);
+    const CCTK_REAL hbw = fmax(hbw_raw, h / c.w);
+    const CCTK_REAL newmu = 1.0 / (hbw + rfsqr * mu);
+
+    return mu - newmu;
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval<CCTK_REAL>
+  initial_bracket(ROOTSTAT &status) const {
+    status = ROOTSTAT::SUCCESS;
+    CCTK_REAL mu_max = CCTK_REAL(1) / h0;
+
+    if (rsqr >= h0 * h0) {
+      const int ndigits2 = 36;
+      const CCTK_REAL margin = std::ldexp(CCTK_REAL(1), 3 - ndigits2);
+      RePrimAnd::f_upper g(h0, rsqr, rbsqr, bsqr);
+      const CCTK_REAL mu_try =
+          findroot_using_deriv(g, status, ndigits2, ndigits2 + 4);
+
+      if (status != ROOTSTAT::SUCCESS || !std::isfinite(mu_try)) {
+        if (status == ROOTSTAT::SUCCESS)
+          status = ROOTSTAT::NOT_CONVERGED;
+        return {CCTK_REAL(0), CCTK_REAL(1) / h0};
+      }
+
+      mu_max = mu_try * (CCTK_REAL(1) + margin);
+    }
+
+    return {CCTK_REAL(0), mu_max};
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval<CCTK_REAL>
+  initial_bracket() const {
+    ROOTSTAT status{};
+    return initial_bracket(status);
+  }
+
+  // RePrimAnd stop condition used by the TOMS748-style bracket solver.
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE bool
+  stopif(CCTK_REAL mu, CCTK_REAL dmu, CCTK_REAL acc) const {
+    return std::abs(dmu) * last.w * last.w < mu * acc;
+  }
+
+  // Public parameters accessed by rarecase/f_rare
+  CCTK_REAL h0{1}, d{0}, qtot{0}, rsqr{0}, rbsqr{0}, bsqr{0}, brosqr{0};
+  CCTK_REAL winf{1}, vsqrinf{0};
+
+private:
+  const EOSType *eos;
+  cache &last;
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
+  x_from_mu(CCTK_REAL mu) const {
+    return 1.0 / (1.0 + mu * bsqr);
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
+  rfsqr_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const {
+    return x * (rsqr * x + mu * (x + 1.0) * rbsqr);
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
+  qf_from_mu_x(CCTK_REAL mu, CCTK_REAL x) const {
+    const CCTK_REAL mux = mu * x;
+    return qtot - (bsqr + mux * mux * brosqr) * CCTK_REAL(0.5);
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE static CCTK_REAL
+  get_eps_raw(CCTK_REAL mu, CCTK_REAL qf, CCTK_REAL rfsqr, CCTK_REAL w) {
+    return w * (qf - mu * rfsqr * (1.0 - mu * w / (1.0 + w)));
+  }
+};
+
+// -------------------- f_rare and rarecase -------------------
+
+template <typename EOSType> class f_rare {
+  const CCTK_REAL v2targ;
+  const froot<EOSType> &f;
+
+public:
+  using value_t = CCTK_REAL;
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE
+  f_rare(CCTK_REAL wtarg, const froot<EOSType> &f_)
+      : v2targ(1.0 - 1.0 / (wtarg * wtarg)), f(f_) {}
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE auto
+  operator()(CCTK_REAL mu) const -> std::pair<CCTK_REAL, CCTK_REAL> {
+    const CCTK_REAL x = 1.0 / (1.0 + mu * f.bsqr);
+    const CCTK_REAL xsq = x * x;
+    const CCTK_REAL rfsq = x * (f.rsqr * x + mu * (x + 1.0) * f.rbsqr);
+    const CCTK_REAL vsq = mu * mu * rfsq;
+    const CCTK_REAL y = vsq - v2targ;
+    const CCTK_REAL dy =
+        CCTK_REAL(2) * mu * x * (xsq * f.rsqr + mu * (xsq + x + 1.0) * f.rbsqr);
+    return {y, dy};
+  }
+};
+
+template <typename EOSType> class rarecase {
+public:
+  interval<CCTK_REAL> bracket;
+  bool rho_too_big{false}, rho_big{false}, rho_too_small{false},
+      rho_small{false};
+
+  CCTK_HOST CCTK_DEVICE rarecase(const interval<CCTK_REAL> ibracket,
+                                 const interval<CCTK_REAL> rgrho,
+                                 const froot<EOSType> &f) {
+    CCTK_REAL muc0 = ibracket.min();
+    CCTK_REAL muc1 = ibracket.max();
+    const int ndigits = 30;
+
+    if (f.d > rgrho.max()) {
+      const CCTK_REAL wc = f.d / rgrho.max();
+      if (wc > f.winf) {
+        rho_too_big = true;
+      } else {
+        f_rare<EOSType> g(wc, f);
+        if (g(muc1).first <= 0.0) {
+          rho_too_big = true;
+        } else if (g(muc0).first < 0.0) {
+          ROOTSTAT status{};
+          CCTK_REAL mucc =
+              findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2);
+          if (status == ROOTSTAT::SUCCESS) {
+            muc0 = fmax(muc0, mucc);
+            rho_big = true;
+          } else {
+            rho_too_big = true;
+          }
+        }
+      }
+    }
+
+    if (f.d < f.winf * rgrho.min()) {
+      const CCTK_REAL wc = f.d / rgrho.min();
+      if (wc < 1.0) {
+        rho_too_small = true;
+      } else {
+        f_rare<EOSType> g(wc, f);
+        if (g(muc0).first >= 0.0) {
+          rho_too_small = true;
+        } else if (g(muc1).first > 0.0) {
+          ROOTSTAT status{};
+          CCTK_REAL mucc =
+              findroot_using_deriv(g, ibracket, status, ndigits, ndigits + 2);
+          if (status == ROOTSTAT::SUCCESS) {
+            muc1 = fmin(muc1, mucc);
+            rho_small = true;
+          } else {
+            rho_too_small = true;
+          }
+        }
+      }
+    }
+
+    bracket = interval<CCTK_REAL>{muc0, muc1};
+  }
+};
+
+} // namespace RePrimAnd
+} // namespace Con2PrimFactory
+
+#endif
diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx
new file mode 100644
index 00000000..ee516a59
--- /dev/null
+++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_intervals.hxx
@@ -0,0 +1,40 @@
+#ifndef C2P_1DREPRIMAND_INTERVALS_HXX
+#define C2P_1DREPRIMAND_INTERVALS_HXX
+
+namespace Con2PrimFactory {
+
+template <class T> class interval {
+  T lo_, hi_;
+
+public:
+  using value_t = T;
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval()
+      : lo_(T(0)), hi_(T(0)) {}
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE interval(T lo, T hi)
+      : lo_(lo), hi_(hi) {}
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T min() const {
+    return lo_;
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T max() const {
+    return hi_;
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T &min() { return lo_; }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T &max() { return hi_; }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T width() const {
+    return hi_ - lo_;
+  }
+
+  CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE T mid() const {
+    return (lo_ + hi_) * T(0.5);
+  }
+};
+
+} // namespace Con2PrimFactory
+#endif
diff --git a/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx
new file mode 100644
index 00000000..01728641
--- /dev/null
+++ b/Con2PrimFactory/src/c2p_1DRePrimAnd_rootfinder.hxx
@@ -0,0 +1,242 @@
+#ifndef C2P_1DREPRIMAND_ROOTFINDER_HXX
+#define C2P_1DREPRIMAND_ROOTFINDER_HXX
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <type_traits>
+#include <utility>
+
+#include "roots.hxx" // Algo::brent
+#include "c2p_1DRePrimAnd_intervals.hxx"
+#include "c2p_toms748.hxx"
+#include "c2p_utils.hxx"
+
+namespace Con2PrimFactory {
+
+namespace RePrimAnd_rootdetail {
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool same_sign(T a,
+                                                                        T b) {
+  return (a > T(0) && b > T(0)) || (a < T(0) && b < T(0));
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool
+in_open_interval(T x, T a, T b) {
+  return (x > a) && (x < b);
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T midpoint(T a,
+                                                                     T b) {
+  return (a + b) * T(0.5);
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T secant(T a, T fa,
+                                                                   T b, T fb) {
+  const T denom = (fb - fa);
+  if (denom == T(0))
+    return std::numeric_limits<T>::quiet_NaN();
+  return (a * fb - b * fa) / denom;
+}
+
+} // namespace RePrimAnd_rootdetail
+
+// Derivative-assisted Newton with bracketing, Brent fallback on pathology.
+template <class F, class T = typename F::value_t>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T
+findroot_using_deriv(F &f, interval<T> bracket, ROOTSTAT &errs, int digits,
+                     unsigned int max_calls = 256) {
+  using namespace RePrimAnd_rootdetail;
+
+  if (max_calls < 4u) {
+    errs = ROOTSTAT::NOT_CONVERGED;
+    return std::numeric_limits<T>::quiet_NaN();
+  }
+
+  T a = bracket.min();
+  T b = bracket.max();
+  if (b < a) {
+    const T tmp = a;
+    a = b;
+    b = tmp;
+  }
+
+  auto eval_pair = [&](T x) -> std::pair<T, T> {
+    const auto val = f(x);
+    return {T(val.first), T(val.second)};
+  };
+
+  auto va = eval_pair(a);
+  auto vb = eval_pair(b);
+  T fa = va.first;
+  T fb = vb.first;
+
+  // Match RePrimAnd endpoint precedence: right endpoint first.
+  if (fb == T(0)) {
+    errs = ROOTSTAT::SUCCESS;
+    return b;
+  }
+  if (fa == T(0)) {
+    errs = ROOTSTAT::SUCCESS;
+    return a;
+  }
+  if (!std::isfinite(fa) || !std::isfinite(fb) || same_sign(fa, fb)) {
+    errs = ROOTSTAT::NOT_BRACKETED;
+    return std::numeric_limits<T>::quiet_NaN();
+  }
+
+  const int minbits =
+      (digits > 0) ? digits : (std::numeric_limits<T>::digits - 4);
+  const T tol = std::ldexp(T(1), -minbits);
+  const unsigned int newton_calls = max_calls - 2u;
+
+  // Same initial guess as RePrimAnd reference.
+  T x = secant(a, fa, b, fb);
+  if (!std::isfinite(x) || !in_open_interval(x, a, b))
+    x = midpoint(a, b);
+
+  int stalled = 0;
+  for (unsigned int calls = 0; calls < newton_calls; ++calls) {
+    const auto vx = eval_pair(x);
+    const T fx = vx.first;
+    const T dfx = vx.second;
+
+    if (!std::isfinite(fx) || !std::isfinite(dfx))
+      break;
+    if (fx == T(0)) {
+      errs = ROOTSTAT::SUCCESS;
+      return x;
+    }
+
+    if (same_sign(fa, fx)) {
+      a = x;
+      fa = fx;
+    } else {
+      b = x;
+      fb = fx;
+    }
+
+    const T width = b - a;
+    const T scale = std::max(T(1), std::max(std::abs(a), std::abs(b)));
+    if (std::abs(width) <= tol * scale) {
+      errs = ROOTSTAT::SUCCESS;
+      return x;
+    }
+
+    T x_new = std::numeric_limits<T>::quiet_NaN();
+    const T dfx_floor = std::sqrt(std::numeric_limits<T>::epsilon());
+    if (std::abs(dfx) > dfx_floor)
+      x_new = x - fx / dfx;
+    if (!std::isfinite(x_new) || !in_open_interval(x_new, a, b))
+      x_new = secant(a, fa, b, fb);
+    if (!std::isfinite(x_new) || !in_open_interval(x_new, a, b))
+      x_new = midpoint(a, b);
+
+    const T step_scale = std::max(T(1), std::abs(x));
+    if (std::abs(x_new - x) <= std::numeric_limits<T>::epsilon() * step_scale)
+      ++stalled;
+    else
+      stalled = 0;
+    x = x_new;
+
+    if (stalled >= 2)
+      break;
+  }
+
+  // Fallback to robust bracketing in pathological derivative/stall cases.
+  auto fn = [&](T x0) -> T { return eval_pair(x0).first; };
+  if (!std::isfinite(a) || !std::isfinite(b) || !std::isfinite(fa) ||
+      !std::isfinite(fb) || same_sign(fa, fb)) {
+    errs = ROOTSTAT::NOT_BRACKETED;
+    return std::numeric_limits<T>::quiet_NaN();
+  }
+
+  int iters = 0;
+  auto ab = Algo::brent(fn, a, b, minbits, int(max_calls), iters);
+  errs = (iters >= int(max_calls)) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS;
+  if (errs != ROOTSTAT::SUCCESS)
+    return std::numeric_limits<T>::quiet_NaN();
+
+  const T a_root = ab.first, b_root = ab.second;
+  const T fa2 = fn(a_root), fb2 = fn(b_root);
+  if (!std::isfinite(fa2) || !std::isfinite(fb2)) {
+    errs = ROOTSTAT::NOT_CONVERGED;
+    return std::numeric_limits<T>::quiet_NaN();
+  }
+
+  if (fb2 == T(0) || std::abs(fb2) < std::abs(fa2))
+    return b_root;
+  if (std::abs(fa2) < std::abs(fb2))
+    return a_root;
+  return (a_root + b_root) * T(0.5);
+}
+
+// Derivative-free TOMS748 adapter, Brent fallback on pathological stall.
+template <class F, class T = typename F::value_t>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval<T>
+findroot_no_deriv(F &f, interval<T> bracket, T acc, unsigned int max_calls,
+                  ROOTSTAT &errs) {
+  if (max_calls < 10u) {
+    errs = ROOTSTAT::NOT_CONVERGED;
+    T a = bracket.min();
+    T b = bracket.max();
+    if (b < a) {
+      const T tmp = a;
+      a = b;
+      b = tmp;
+    }
+    return {a, b};
+  }
+
+  interval<T> res = toms748_solve(f, bracket, acc, max_calls, errs);
+  if (errs == ROOTSTAT::SUCCESS || errs == ROOTSTAT::NOT_BRACKETED)
+    return res;
+
+  T a = res.min();
+  T b = res.max();
+  if (b < a) {
+    const T tmp = a;
+    a = b;
+    b = tmp;
+  }
+  if (!std::isfinite(a) || !std::isfinite(b))
+    return res;
+
+  auto fn = [&](T x) -> T { return f(x); };
+  const T fa = fn(a), fb = fn(b);
+  if (fb == T(0))
+    return {b, b};
+  if (fa == T(0))
+    return {a, a};
+  if (!std::isfinite(fa) || !std::isfinite(fb) || fa * fb > T(0)) {
+    errs = ROOTSTAT::NOT_BRACKETED;
+    return {std::numeric_limits<T>::lowest(), std::numeric_limits<T>::max()};
+  }
+
+  const int minbits = std::numeric_limits<T>::digits - 4;
+  int iters = 0;
+  auto ab = Algo::brent(fn, a, b, minbits, int(max_calls), iters);
+  errs = (iters >= int(max_calls)) ? ROOTSTAT::NOT_CONVERGED : ROOTSTAT::SUCCESS;
+
+  return {ab.first, ab.second};
+}
+
+template <class F, class T = typename F::value_t>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T
+findroot_using_deriv(F &f, ROOTSTAT &errs, int digits,
+                     unsigned int max_calls = 256) {
+  return findroot_using_deriv(f, f.initial_bracket(), errs, digits, max_calls);
+}
+
+template <class F, class T = typename F::value_t>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval<T>
+findroot_no_deriv(F &f, T tol, unsigned int max_calls, ROOTSTAT &errs) {
+  return findroot_no_deriv(f, f.initial_bracket(), tol, max_calls, errs);
+}
+
+} // namespace Con2PrimFactory
+#endif
diff --git a/Con2PrimFactory/src/c2p_2DNoble.hxx b/Con2PrimFactory/src/c2p_2DNoble.hxx
index 7de97925..d3397dc2 100644
--- a/Con2PrimFactory/src/c2p_2DNoble.hxx
+++ b/Con2PrimFactory/src/c2p_2DNoble.hxx
@@ -20,7 +20,8 @@ public:
       CCTK_REAL tol, CCTK_REAL alp_thresh_in, CCTK_REAL vwlim, CCTK_REAL B_lim,
       CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in,
       CCTK_REAL sigma_max_in, CCTK_REAL inv_beta_max_in, bool ye_len,
-      bool use_z, bool use_temperature, bool use_pressure_atmo);
+      bool use_z, bool use_temperature, bool use_pressure_atmo,
+      bool soft_root_conv, CCTK_REAL soft_root_width_factor_in);
 
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
   get_Ssq_Exact(const vec<CCTK_REAL, 3> &mom,
@@ -78,7 +79,8 @@ CCTK_HOST
         CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in,
         CCTK_REAL vwlim_BH_in, CCTK_REAL sigma_max_in,
         CCTK_REAL inv_beta_max_in, bool ye_len, bool use_z,
-        bool use_temperature, bool use_pressure_atmo) {
+        bool use_temperature, bool use_pressure_atmo, bool soft_root_conv,
+        CCTK_REAL soft_root_width_factor_in) {
 
   // Base
   atmo = atm;
@@ -98,6 +100,8 @@ CCTK_HOST
   use_zprim = use_z;
   use_temp = use_temperature;
   use_press_atmo = use_pressure_atmo;
+  soft_root_convergence = soft_root_conv;
+  soft_root_width_factor = fmax(CCTK_REAL(1.0), soft_root_width_factor_in);
 
   // Derived
   GammaIdealFluid = eos_3p->gamma;
@@ -286,6 +290,7 @@ c2p_2DNoble::solve(const EOSType *eos_3p, prim_vars &pv, prim_vars &pv_seeds,
   rep.iters = 0;
   rep.adjust_cons = false;
   rep.set_atmo = false;
+  rep.soft_root_conv = false;
   rep.status = c2p_report::SUCCESS;
 
   /* Check validity of the 3-metric and compute its inverse */
@@ -522,18 +527,26 @@ c2p_2DNoble::solve(const EOSType *eos_3p, prim_vars &pv, prim_vars &pv_seeds,
   //}
 
   if (fabs(errx) > tolerance) {
-    // set status to root not converged
-    rep.set_root_conv();
-    // status = ROOTSTAT::NOT_CONVERGED;
-    cv = cv_const;
-    return;
+    bool accept_soft = false;
+    if (soft_root_convergence && std::isfinite(errx)) {
+      const CCTK_REAL soft_tol = soft_root_width_factor * tolerance;
+      accept_soft = (fabs(errx) <= soft_tol);
+    }
+    if (!accept_soft) {
+      // set status to root not converged
+      rep.set_root_conv();
+      // status = ROOTSTAT::NOT_CONVERGED;
+      cv = cv_const;
+      return;
+    }
+    rep.set_soft_root_conv();
   }
 
   // Check for bad untrapped divergences
   if ((!isfinite(f)) || (!isfinite(df))) {
     rep.set_root_bracket();
-    // TODO: currently, divergences in f and df are marked as failed bracketing.
-    // Change this.
+    cv = cv_const;
+    return;
   }
 
   /* Calculate primitives from Z and W */
diff --git a/Con2PrimFactory/src/c2p_report.hxx b/Con2PrimFactory/src/c2p_report.hxx
index f612015b..6c6824b5 100644
--- a/Con2PrimFactory/src/c2p_report.hxx
+++ b/Con2PrimFactory/src/c2p_report.hxx
@@ -65,6 +65,9 @@ public:
   /// Whether the artificial atmosphere was enforced.
   bool set_atmo{false};
 
+  /// Root solve accepted via near-convergence policy (warning-only).
+  bool soft_root_conv{false};
+
   /// Number of calls to the EOS needed for the root finding.
   CCTK_INT iters;
 
@@ -205,6 +208,12 @@ public:
     adjust_cons = true;
   }
 
+  CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
+  set_soft_root_conv() {
+    status = SUCCESS;
+    soft_root_conv = true;
+  }
+
   CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
   debug_message() const {
     switch (status) {
@@ -216,6 +225,9 @@ public:
       if (adjust_cons) {
         printf("Conserved variables have been changed. \n");
       }
+      if (soft_root_conv) {
+        printf("Root solve accepted via near-convergence policy. \n");
+      }
       break;
     case INVALID_DETG:
       printf("Invalid metric determinant: detg = %16.8e \n", detg);
@@ -233,8 +245,7 @@ public:
       printf("Density out of range, dens = %16.8e, rho = %16.8e \n", dens, rho);
       break;
     case RANGE_EPS:
-      printf("Specific energy was out of range! eps = %16.8e \n",
-             eps);
+      printf("Specific energy was out of range! eps = %16.8e \n", eps);
       break;
     case SPEED_LIMIT:
       printf("Speed limit exceeded, vx, vy, vz = %16.8e, %16.8e, %16.8e \n",
@@ -259,7 +270,8 @@ public:
       printf("Preparatory root finding failed (faulty bracketing) \n");
       break;
     case ERR_CODE_NOT_SET:
-      printf("Error code not set! Report was initialized but never changed. \n");
+      printf(
+          "Error code not set! Report was initialized but never changed. \n");
       break;
     default:
       assert(false);
diff --git a/Con2PrimFactory/src/c2p_toms748.hxx b/Con2PrimFactory/src/c2p_toms748.hxx
new file mode 100644
index 00000000..114b4d6f
--- /dev/null
+++ b/Con2PrimFactory/src/c2p_toms748.hxx
@@ -0,0 +1,271 @@
+#ifndef C2P_TOMS748_HXX
+#define C2P_TOMS748_HXX
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <type_traits>
+#include <utility>
+
+#include "c2p_1DRePrimAnd_intervals.hxx"
+#include "c2p_utils.hxx"
+
+namespace Con2PrimFactory {
+
+namespace toms748_detail {
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool same_sign(T a,
+                                                                        T b) {
+  return (a > T(0) && b > T(0)) || (a < T(0) && b < T(0));
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool
+in_open_interval(T x, T a, T b) {
+  return (x > a) && (x < b);
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T midpoint(T a,
+                                                                     T b) {
+  return (a + b) * T(0.5);
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T
+secant_step(T a, T fa, T b, T fb) {
+  const T denom = fb - fa;
+  if (denom == T(0))
+    return std::numeric_limits<T>::quiet_NaN();
+  return (a * fb - b * fa) / denom;
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline T
+iqi_step(T a, T fa, T b, T fb, T c, T fc) {
+  const T dab = (fa - fb);
+  const T dac = (fa - fc);
+  const T dba = (fb - fa);
+  const T dbc = (fb - fc);
+  const T dca = (fc - fa);
+  const T dcb = (fc - fb);
+
+  if (dab == T(0) || dac == T(0) || dba == T(0) || dbc == T(0) ||
+      dca == T(0) || dcb == T(0))
+    return std::numeric_limits<T>::quiet_NaN();
+
+  return a * fb * fc / (dab * dac) + b * fa * fc / (dba * dbc) +
+         c * fa * fb / (dca * dcb);
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval<T>
+invalid_interval() {
+  return {std::numeric_limits<T>::lowest(), std::numeric_limits<T>::max()};
+}
+
+template <class F, class T> struct has_stopif {
+private:
+  template <class U>
+  static auto test(int)
+      -> decltype(std::declval<U &>().stopif(std::declval<T>(),
+                                             std::declval<T>(),
+                                             std::declval<T>()),
+                  std::true_type{});
+  template <class> static std::false_type test(...);
+
+public:
+  static constexpr bool value = decltype(test<F>(0))::value;
+};
+
+template <class F, class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool
+call_stopif_impl(F &f, T left, T right, T tol, std::true_type) {
+  // Preserve RePrimAnd contract: stopif(left, right-left, tol)
+  return f.stopif(left, right - left, tol);
+}
+
+template <class F, class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool
+call_stopif_impl(F &, T left, T right, T tol, std::false_type) {
+  const T width = std::abs(right - left);
+  const T scale = std::max(T(1), std::max(std::abs(left), std::abs(right)));
+  return width <= tol * scale;
+}
+
+template <class F, class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline bool
+call_stopif(F &f, T left, T right, T tol) {
+  using has_stopif_t = std::integral_constant<bool, has_stopif<F, T>::value>;
+  return call_stopif_impl(f, left, right, tol, has_stopif_t{});
+}
+
+template <class T>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
+update_bracket(T x, T fx, T &a, T &fa, T &b, T &fb) {
+  if (same_sign(fa, fx)) {
+    a = x;
+    fa = fx;
+  } else {
+    b = x;
+    fb = fx;
+  }
+}
+
+} // namespace toms748_detail
+
+template <class F, class T = typename F::value_t>
+CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline interval<T>
+toms748_solve(F &f, interval<T> bracket, T tol, unsigned int max_calls,
+              ROOTSTAT &errs) {
+  using namespace toms748_detail;
+
+  T a = bracket.min();
+  T b = bracket.max();
+  if (b < a) {
+    const T tmp = a;
+    a = b;
+    b = tmp;
+  }
+
+  if (!std::isfinite(a) || !std::isfinite(b)) {
+    errs = ROOTSTAT::NOT_BRACKETED;
+    return invalid_interval<T>();
+  }
+
+  T fa = f(a);
+  T fb = f(b);
+  unsigned int calls = 2u;
+
+  // Match RePrimAnd endpoint precedence: right endpoint first.
+  if (fb == T(0)) {
+    errs = ROOTSTAT::SUCCESS;
+    return {b, b};
+  }
+  if (fa == T(0)) {
+    errs = ROOTSTAT::SUCCESS;
+    return {a, a};
+  }
+
+  if (!std::isfinite(fa) || !std::isfinite(fb) || same_sign(fa, fb)) {
+    errs = ROOTSTAT::NOT_BRACKETED;
+    return invalid_interval<T>();
+  }
+
+  if (max_calls <= 2u) {
+    errs = ROOTSTAT::NOT_CONVERGED;
+    return {a, b};
+  }
+
+  if (call_stopif(f, a, b, tol)) {
+    errs = ROOTSTAT::SUCCESS;
+    return {a, b};
+  }
+
+  T prev_width = b - a;
+
+  while (calls < max_calls) {
+    const T a0 = a, b0 = b;
+    const T fa0 = fa, fb0 = fb;
+    const T width0 = b - a;
+
+    T x1 = secant_step(a, fa, b, fb);
+    if (!std::isfinite(x1) || !in_open_interval(x1, a, b))
+      x1 = midpoint(a, b);
+
+    T fx1 = f(x1);
+    ++calls;
+
+    if (!std::isfinite(fx1)) {
+      if (calls >= max_calls)
+        break;
+      x1 = midpoint(a, b);
+      fx1 = f(x1);
+      ++calls;
+      if (!std::isfinite(fx1)) {
+        errs = ROOTSTAT::NOT_CONVERGED;
+        return {a, b};
+      }
+    }
+
+    if (fx1 == T(0)) {
+      errs = ROOTSTAT::SUCCESS;
+      return {x1, x1};
+    }
+
+    update_bracket(x1, fx1, a, fa, b, fb);
+    if (call_stopif(f, a, b, tol)) {
+      errs = ROOTSTAT::SUCCESS;
+      return {a, b};
+    }
+    if (calls >= max_calls)
+      break;
+
+    T x2 = iqi_step(a0, fa0, b0, fb0, x1, fx1);
+    if (!std::isfinite(x2) || !in_open_interval(x2, a, b))
+      x2 = midpoint(a, b);
+
+    const T width1 = b - a;
+    if ((x2 - a) < T(0.125) * width1 || (b - x2) < T(0.125) * width1)
+      x2 = midpoint(a, b);
+
+    T fx2 = f(x2);
+    ++calls;
+
+    if (!std::isfinite(fx2)) {
+      if (calls >= max_calls)
+        break;
+      x2 = midpoint(a, b);
+      fx2 = f(x2);
+      ++calls;
+      if (!std::isfinite(fx2)) {
+        errs = ROOTSTAT::NOT_CONVERGED;
+        return {a, b};
+      }
+    }
+
+    if (fx2 == T(0)) {
+      errs = ROOTSTAT::SUCCESS;
+      return {x2, x2};
+    }
+
+    update_bracket(x2, fx2, a, fa, b, fb);
+    if (call_stopif(f, a, b, tol)) {
+      errs = ROOTSTAT::SUCCESS;
+      return {a, b};
+    }
+
+    const T new_width = b - a;
+    if (calls < max_calls &&
+        (new_width > T(0.5) * width0 || new_width > T(0.9) * prev_width)) {
+      const T xm = midpoint(a, b);
+      const T fxm = f(xm);
+      ++calls;
+
+      if (!std::isfinite(fxm)) {
+        errs = ROOTSTAT::NOT_CONVERGED;
+        return {a, b};
+      }
+      if (fxm == T(0)) {
+        errs = ROOTSTAT::SUCCESS;
+        return {xm, xm};
+      }
+
+      update_bracket(xm, fxm, a, fa, b, fb);
+      if (call_stopif(f, a, b, tol)) {
+        errs = ROOTSTAT::SUCCESS;
+        return {a, b};
+      }
+    }
+
+    prev_width = b - a;
+  }
+
+  errs = ROOTSTAT::NOT_CONVERGED;
+  return {a, b};
+}
+
+} // namespace Con2PrimFactory
+
+#endif
diff --git a/Con2PrimFactory/src/test.cxx b/Con2PrimFactory/src/test.cxx
index b8f87764..59cb61c6 100644
--- a/Con2PrimFactory/src/test.cxx
+++ b/Con2PrimFactory/src/test.cxx
@@ -6,6 +6,7 @@
 #include "c2p.hxx"
 #include "c2p_1DEntropy.hxx"
 #include "c2p_1DPalenzuela.hxx"
+#include "c2p_1DRePrimAnd.hxx"
 #include "c2p_2DNoble.hxx"
 
 #include "c2p_utils.hxx"
@@ -66,17 +67,21 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) {
   //  Ye_lenient, use_z, use_temperature)
   c2p_2DNoble c2p_Noble(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH,
                         eps_BH, vwlim_BH, sigma_max, inv_beta_max, true, false,
-                        false, false);
+                        false, false, false, 1.0);
   c2p_1DPalenzuela c2p_Pal(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH,
                            eps_BH, vwlim_BH, sigma_max, inv_beta_max, true,
-                           false, false, false);
+                           false, false, false, false, 1.0);
+  c2p_1DRePrimAnd c2p_RPA(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH,
+                          eps_BH, vwlim_BH, sigma_max, inv_beta_max, true,
+                          false, false, false, false, 1.0);
   c2p_1DEntropy c2p_Ent(eos_3p_ig, atmo, 100, 1e-8, alp_thresh, 1, 1, rho_BH,
                         eps_BH, vwlim_BH, sigma_max, inv_beta_max, true, false,
-                        false, false);
+                        false, false, false, 1.0);
 
   // Construct error report object:
   c2p_report rep_Noble;
   c2p_report rep_Pal;
+  c2p_report rep_RPA;
   c2p_report rep_Ent;
 
   // Set primitive seeds
@@ -211,6 +216,45 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) {
   */
   rep_Pal.debug_message();
 
+  // Testing C2P RePrimAnd
+  CCTK_VINFO("Testing C2P RePrimAnd...");
+  c2p_RPA.solve(eos_3p_ig, pv, cv_all, alp, beta, g, rep_RPA);
+
+  printf("pv_seeds, pv: \n"
+         "rho: %f, %f \n"
+         "eps: %f, %f \n"
+         "Ye: %f, %f \n"
+         "press: %f, %f \n"
+         "temperature: %f, %f \n"
+         "entropy: %f, %f \n"
+         "velx: %f, %f \n"
+         "vely: %f, %f \n"
+         "velz: %f, %f \n"
+         "Bx: %f, %f \n"
+         "By: %f, %f \n"
+         "Bz: %f, %f \n",
+         pv_seeds.rho, pv.rho, pv_seeds.eps, pv.eps, pv_seeds.Ye, pv.Ye,
+         pv_seeds.press, pv.press, pv_seeds.temperature, pv.temperature,
+         pv_seeds.entropy, pv.entropy, pv_seeds.vel(0), pv.vel(0),
+         pv_seeds.vel(1), pv.vel(1), pv_seeds.vel(2), pv.vel(2),
+         pv_seeds.Bvec(0), pv.Bvec(0), pv_seeds.Bvec(1), pv.Bvec(1),
+         pv_seeds.Bvec(2), pv.Bvec(2));
+  printf("cv: \n"
+         "dens: %f \n"
+         "tau: %f \n"
+         "momx: %f \n"
+         "momy: %f \n"
+         "momz: %f \n"
+         "DYe: %f \n"
+         "dBx: %f \n"
+         "dBy: %f \n"
+         "dBz: %f \n"
+         "DEnt: %f \n",
+         cv_all.dens, cv_all.tau, cv_all.mom(0), cv_all.mom(1), cv_all.mom(2),
+         cv_all.DYe, cv_all.dBvec(0), cv_all.dBvec(1), cv_all.dBvec(2),
+         cv_all.DEnt);
+  rep_RPA.debug_message();
+
   // Testing C2P Entropy
   CCTK_VINFO("Testing C2P Entropy...");
   c2p_Ent.solve(eos_3p_ig, pv, cv_all, alp, beta, g, rep_Ent);