Merge pull request 'Audit hardening (S1+S2): numerics, newton_core refactor, solver diagnostics, validation' (#39) from fix/b1-v3-c1-quick-wins into main

.gitea/workflows/cpp-tests.yml

@@ -41,7 +41,11 @@ jobs:
         run: cmake -S code -B build -DCMAKE_BUILD_TYPE=Release
 
       - name: Build
-        run: nice -n 19 cmake --build build --target conformallab_tests -j$(nproc)
+        # Serial build (-j1): the 800m-limited container OOM-kills cc1plus when
+        # several Eigen+GoogleTest translation units compile in parallel
+        # (`-j$(nproc)`), failing test-fast non-deterministically regardless of
+        # branch content.  Mirrors the Pi-protection already used by test-cgal.
+        run: nice -n 19 cmake --build build --target conformallab_tests -j1
 
       - name: Run tests
         run: >
@@ -160,9 +164,18 @@ jobs:
       - name: shellcheck (scripts/**/*.sh, severity=warning, strict)
         run: bash scripts/quality/shellcheck.sh --strict
 
+      - name: Install lcov (coverage gate)
+        run: apt-get install -y --no-install-recommends lcov
+
+      - name: Coverage gate (fast-test suite, ramp-up mode)
+        # SKIP_COVERAGE_GATE=1: reports numbers but does not fail on threshold.
+        # Remove once I5 is resolved (coverage is wired to the full CGAL suite,
+        # not just the 26 fast tests).  Threshold: 80% line / 70% branch / 90% func.
+        run: SKIP_COVERAGE_GATE=1 bash scripts/quality/coverage.sh
+
       - name: Summary
         if: always()
         run: |
-          echo "QUALITY ▸ all four gates passed."
+          echo "QUALITY ▸ all gates passed."
           echo "  see scripts/quality/README.md for the full catalogue"
-          echo "  (sanitizers, clang-tidy, coverage, etc. are local-only)"
+          echo "  (sanitizers, clang-tidy, coverage report in build-coverage/)"

code/include/CGAL/Discrete_circle_packing.h

@@ -8,7 +8,7 @@
 \ingroup PkgConformalMapRef
 
 User-facing entry for the **face-based** circle-packing functional of
-Bobenko-Pinkall-Springborn 2010.  See `cp_euclidean_functional.hpp`
+Bobenko-Pinkall-Springborn 2015.  See `cp_euclidean_functional.hpp`
 for the underlying algorithm and `doc/architecture/phase-9a-validation.md`
 for the line-by-line mapping to the Java original
 `CPEuclideanFunctional.java`.
@@ -44,7 +44,7 @@ namespace CGAL {
 \ingroup PkgConformalMapConcepts
 \brief Traits class for `discrete_circle_packing_euclidean()` —
 declares the kernel, mesh and property-map types used by the
-BPS-2010 face-based circle-packing functional.
+BPS-2015 face-based circle-packing functional.
 
 Primary template; specialise it for non-`Surface_mesh` triangle meshes.
 */
@@ -96,6 +96,9 @@ struct Default_cp_euclidean_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
 Result of `discrete_circle_packing_euclidean`.  Carries face DOFs
 `ρ_f = log R_f` rather than the vertex DOFs of the classical modes.
 */
+/// Why the Newton solve terminated (re-exported from the solver; I1 audit).
+using Newton_status = ::conformallab::NewtonStatus;
+
 template <typename FT = double>
 struct Circle_packing_result
 {
@@ -108,6 +111,14 @@ struct Circle_packing_result
     FT   gradient_norm = FT(0);
     /// `true` iff `gradient_norm < gradient_tolerance` at exit.
     bool converged     = false;
+
+    /// Why the solve stopped: `Converged` / `MaxIterations` /
+    /// `LinearSolverFailed` / `LineSearchStalled`.
+    Newton_status status = Newton_status::MaxIterations;
+    /// `true` iff any Newton step fell back to SparseQR (gauge singularity).
+    bool sparse_qr_fallback_used = false;
+    /// Smallest `|Dᵢᵢ|` of the last LDLT (near-singularity proxy; 0 if none).
+    FT   min_ldlt_pivot = FT(0);
 };
 
 // ── Entry function ────────────────────────────────────────────────────────────
@@ -115,7 +126,7 @@ struct Circle_packing_result
 /*!
 \ingroup PkgConformalMapRef
 
-Compute the BPS-2010 face-based circle-packing of `mesh`.
+Compute the BPS-2015 face-based circle-packing of `mesh`.
 
 \tparam TriangleMesh  A `CGAL::Surface_mesh<P>`.
 \tparam NamedParameters  Optional CGAL named-parameter pack.
@@ -132,7 +143,7 @@ Compute the BPS-2010 face-based circle-packing of `mesh`.
 \returns A `Circle_packing_result<FT>` with `ρ_f` per face.
 
 \pre `mesh` is a triangle mesh.
-\pre `φ_f` and `θ_e` satisfy the BPS-2010 admissibility conditions
+\pre `φ_f` and `θ_e` satisfy the BPS-2015 admissibility conditions
      (Σ_f φ_f = 2π·χ + Σ_e (π − θ_e), see paper §6).
 */
 template <typename TriangleMesh,
@@ -188,6 +199,9 @@ auto discrete_circle_packing_euclidean(
     result.iterations    = nr.iterations;
     result.gradient_norm = nr.grad_inf_norm;
     result.converged     = nr.converged;
+    result.status                  = nr.status;
+    result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used;
+    result.min_ldlt_pivot          = static_cast<FT>(nr.min_ldlt_pivot);
 
     // ── output_uv_map (Phase 8b-Lite extension) ────────────────────────────
     //
@@ -201,7 +215,7 @@ auto discrete_circle_packing_euclidean(
     //
     // For Phase 8b-Lite we deliberately don't fake it.  If the caller
     // supplies `output_uv_map(pmap)` we throw `std::runtime_error` with a
-    // clear pointer to Phase 9c (BPS-2010 §6 face-based circle-packing
+    // clear pointer to Phase 9c (BPS-2015 §6 face-based circle-packing
     // layout, ~150 lines, on the porting roadmap).  Failing loudly is
     // better than silently writing zeros.
     //

code/include/CGAL/Discrete_conformal_map.h

@@ -78,6 +78,10 @@ namespace CGAL {
 // Result type
 // ════════════════════════════════════════════════════════════════════════════
 
+/// Why the Newton solve terminated (re-exported from the solver; I1 audit).
+/// `Converged` / `MaxIterations` / `LinearSolverFailed` / `LineSearchStalled`.
+using Newton_status = ::conformallab::NewtonStatus;
+
 /*!
 \ingroup PkgConformalMapRef
 
@@ -100,9 +104,18 @@ struct Conformal_map_result
     /// `true` iff `gradient_norm < gradient_tolerance`.
     bool converged = false;
 
+    /// Why the solve stopped (more informative than `converged` alone): one of
+    /// `Converged` / `MaxIterations` / `LinearSolverFailed` / `LineSearchStalled`.
+    Newton_status status = Newton_status::MaxIterations;
+
     /// `true` iff the linear solver used the SparseQR fallback at any
     /// Newton step (gauge mode on closed mesh without pinned vertex).
     bool sparse_qr_fallback_used = false;
+
+    /// Smallest `|Dᵢᵢ|` of the last LDLT factorisation — a cheap
+    /// near-singularity proxy (small ⇒ ill-conditioned solve). `0` if LDLT
+    /// never succeeded.
+    FT min_ldlt_pivot = FT(0);
 };
 
 
@@ -272,6 +285,9 @@ auto discrete_conformal_map_euclidean(
     result.iterations    = nr.iterations;
     result.gradient_norm = nr.grad_inf_norm;
     result.converged     = nr.converged;
+    result.status                  = nr.status;
+    result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used;
+    result.min_ldlt_pivot          = static_cast<FT>(nr.min_ldlt_pivot);
 
     // ── 8. Optional layout step (Phase 8b-Lite extension) ──────────────────
     //
@@ -410,6 +426,9 @@ auto discrete_conformal_map_spherical(
     result.iterations    = nr.iterations;
     result.gradient_norm = nr.grad_inf_norm;
     result.converged     = nr.converged;
+    result.status                  = nr.status;
+    result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used;
+    result.min_ldlt_pivot          = static_cast<FT>(nr.min_ldlt_pivot);
 
     // Optional 3-D layout step (point on S²)
     auto uv_param = parameters::get_parameter(
@@ -459,6 +478,14 @@ struct Hyper_ideal_map_result
     FT   gradient_norm = FT(0);
     /// `true` iff `gradient_norm < gradient_tolerance` at exit.
     bool converged     = false;
+
+    /// Why the solve stopped: `Converged` / `MaxIterations` /
+    /// `LinearSolverFailed` / `LineSearchStalled`.
+    Newton_status status = Newton_status::MaxIterations;
+    /// `true` iff any Newton step fell back to SparseQR (gauge singularity).
+    bool sparse_qr_fallback_used = false;
+    /// Smallest `|Dᵢᵢ|` of the last LDLT (near-singularity proxy; 0 if none).
+    FT   min_ldlt_pivot = FT(0);
 };
 
 /*!
@@ -538,6 +565,9 @@ auto discrete_conformal_map_hyper_ideal(
     result.iterations    = nr.iterations;
     result.gradient_norm = nr.grad_inf_norm;
     result.converged     = nr.converged;
+    result.status                  = nr.status;
+    result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used;
+    result.min_ldlt_pivot          = static_cast<FT>(nr.min_ldlt_pivot);
 
     // Optional Poincaré-disk layout (2-D in the unit disk).
     auto uv_param = parameters::get_parameter(

code/include/CGAL/Discrete_inversive_distance.h

@@ -8,7 +8,7 @@
 \ingroup PkgConformalMapRef
 
 User-facing entry for the **vertex-based** inversive-distance circle-
-packing functional of Luo (2004), with the Bowers-Stephenson (2004)
+packing functional of Luo 2004, with the Bowers-Stephenson (2004)
 initialisation.  See `inversive_distance_functional.hpp` for the
 underlying algorithm and `doc/roadmap/research-track.md` (item 9a.2)
 for the research-track classification — this functional has **no Java
@@ -98,7 +98,7 @@ struct Default_inversive_distance_traits<CGAL::Surface_mesh<typename K::Point_3>
 /*!
 \ingroup PkgConformalMapRef
 
-Compute the Luo-2004 vertex-based inversive-distance circle packing of `mesh`.
+Compute the Luo 2004 vertex-based inversive-distance circle packing of `mesh`.
 
 The per-edge constant `I_ij` is computed once at the start from the input
 3-D geometry via the Bowers-Stephenson identity
@@ -207,6 +207,9 @@ auto discrete_inversive_distance_map(
     result.iterations    = nr.iterations;
     result.gradient_norm = nr.grad_inf_norm;
     result.converged     = nr.converged;
+    result.status                  = nr.status;
+    result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used;
+    result.min_ldlt_pivot          = static_cast<FT>(nr.min_ldlt_pivot);
 
     // ── Optional layout step (Phase 8b-Lite extension) ─────────────────────
     //

code/include/constants.hpp

@@ -17,4 +17,28 @@ constexpr double PI     = 3.14159265358979323846264338328;
 /// 2π (full turn).
 constexpr double TWO_PI = 2.0 * PI;
 
+// ── Numerical thresholds and bounds (N4/N6 audit) ──────────────────────────
+
+/// Edge length floor for log-scale functionals (euclidean_functional, spherical_functional).
+/// exp(-30.0) ≈ 3e-7: edges below this length are treated as degenerate/zero.
+/// Rationale: avoids log(tiny) overflow and enforces a minimum edge scale.
+constexpr double LOG_EDGE_LENGTH_FLOOR = -30.0;
+
+/// Hyper-ideal scale factor lower bound (hyper_ideal_functional:179-181).
+/// If b (scale) becomes negative (infeasible), clamp to 0.01 to keep geometry valid.
+/// Rationale: ensures b > 0 maintains the Lorentzian model's causal structure.
+constexpr double HYPER_IDEAL_SCALE_FLOOR = 0.01;
+
+/// Sharpness β of the optional C¹ smooth-barrier scale floor (N3 audit).
+/// Only used when the hyper-ideal clamp mode is `SmoothBarrier`; larger β
+/// makes the softplus transition tighter around `HYPER_IDEAL_SCALE_FLOOR`
+/// (β·floor ≈ 1, so β ≈ 100 keeps the barrier ≈ identity for b ≳ 0.05 while
+/// staying C¹ across b = 0).  See `clamp_hyper_ideal_scale`.
+constexpr double HYPER_IDEAL_SCALE_SHARPNESS = 100.0;
+
+/// Domain guard for asin in spherical_l (spherical_geometry:34).
+/// Clamps exp(λ/2) to slightly below 1.0 so asin stays in [−π/2, π/2].
+/// Rationale: asin(x) is undefined for |x| > 1; this prevents IEEE Inf/NaN.
+constexpr double ASIN_DOMAIN_GUARD = 1.0 - 1e-15;
+
 } // namespace conformallab

code/include/cp_euclidean_functional.hpp

@@ -19,7 +19,7 @@
 // │  The variable is ρ_f = log R_f.                                           │
 // │  Adjacent face-circles intersect at a prescribed angle θ_e per edge.      │
 // │                                                                           │
-// │  This is the FACE-DUAL of the classical vertex-based Luo (2004)           │
+// │  This is the FACE-DUAL of the classical vertex-based Luo 2004           │
 // │  inversive-distance circle packing implemented in                          │
 // │  inversive_distance_functional.hpp (Phase 9a.2).  The relation            │
 // │     I_ij  =  cos θ_e                                                       │
@@ -33,7 +33,7 @@
 // │    θ_e per edge   intersection angle of the two face-circles              │
 // │    φ_f per face   target sum of corner-angles inside the face             │
 // │                                                                           │
-// │  Energy  (BPS-2010 §6)                                                    │
+// │  Energy  (BPS-2015 §6)                                                    │
 // │    E(ρ) =  Σ_f  φ_f · ρ_f                                                  │
 // │         +  Σ_{(h,f=face(h)):                                              │
 // │              [ if opposite face exists ]                                   │
@@ -52,7 +52,7 @@
 // │    Per interior hf:   −(p + θ*)  added to G[face(h)]                      │
 // │    Per boundary hf:   −2 θ*       added to G[face(h)]                     │
 // │                                                                           │
-// │  Hessian (analytic, BPS-2010 eq. 6.8; Java getHessian lines 127-166)     │
+// │  Hessian (analytic, BPS-2015 eq. 6.8; Java getHessian lines 127-166)     │
 // │    Per interior undirected edge e (connecting faces j and k):             │
 // │      h_jk = sin θ / (cosh Δρ − cos θ)                                     │
 // │      H[j,j] += h_jk,  H[k,k] += h_jk,  H[j,k] −= h_jk,  H[k,j] −= h_jk    │
@@ -90,7 +90,7 @@ using CPEMapD = ConformalMesh::Property_map<Edge_index, double>;
 // ── Persistent map bundle ─────────────────────────────────────────────────────
 
 /// Bundle of the three property maps consumed by the CP-Euclidean
-/// (Bobenko-Pinkall-Springborn 2010) circle-packing functional.
+/// (Bobenko-Pinkall-Springborn 2015) circle-packing functional.
 struct CPEuclideanMaps {
     CPFMapI f_idx;     ///< DOF index per face (−1 = pinned)
     CPEMapD theta_e;   ///< intersection angle per edge (default π/2 = orthogonal)

code/include/euclidean_functional.hpp

@@ -162,7 +162,7 @@ inline void compute_euclidean_lambda0_from_mesh(ConformalMesh& mesh, EuclideanMa
         if (len > 1e-15)
             m.lambda0[e] = 2.0 * std::log(len);
         else
-            m.lambda0[e] = -30.0;   // degenerate edge
+            m.lambda0[e] = LOG_EDGE_LENGTH_FLOOR;   // degenerate edge
     }
 }
 

code/include/euclidean_hessian.hpp

@@ -17,7 +17,8 @@
 // │  side lengths lij = exp(Λ̃ij/2):                                         │
 // │                                                                          │
 // │    t12 = −l12+l23+l31,  t23 = l12−l23+l31,  t31 = l12+l23−l31         │
-// │    l123 = l12+l23+l31,  denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area     │
+// │    l123 = l12+l23+l31,  denom2 = 8·Area  (Area via Kahan's stable      │
+// │                          side-length formula — see euclidean_cot_weights)│
 // │                                                                          │
 // │  Cotangent at vertex k (opposite t_opp, adjacent t_a and t_b):          │
 // │    cot_k = (t_opp · l123 − t_a · t_b) / denom2                         │
@@ -48,6 +49,7 @@
 #include <vector>
 #include <array>
 #include <cmath>
+#include <algorithm>
 #include <stdexcept>
 
 namespace conformallab {
@@ -86,11 +88,25 @@ inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31)
         return {0.0, 0.0, 0.0, false};
 
     const double l123 = l12 + l23 + l31;
-    const double denom2_sq = t12 * t23 * t31 * l123;
-    if (denom2_sq <= 0.0) return {0.0, 0.0, 0.0, false};
 
-    // denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area
+    // Triangle area via Kahan's numerically stable formula.  Sort the side
-    const double denom2 = 2.0 * std::sqrt(denom2_sq);
+    // lengths a ≥ b ≥ c and evaluate with the cancellation-avoiding grouping
+    //   Area = ¼·√[ (a+(b+c))·(c−(a−b))·(c+(a−b))·(a+(b−c)) ].
+    // This replaces the naive 2·√(t12·t23·t31·l123): that product of
+    // near-equal-length differences loses precision for needle/cap triangles,
+    // where it would feed large relative error straight into the cotangent
+    // weights and the linear solve.  The result denom2 = 8·Area is identical
+    // up to rounding for well-shaped triangles, but accurate for slivers.
+    double a = l12, b = l23, c = l31;
+    if (a < b) std::swap(a, b);
+    if (a < c) std::swap(a, c);
+    if (b < c) std::swap(b, c);   // now a ≥ b ≥ c
+
+    const double kahan = (a + (b + c)) * (c - (a - b))
+                       * (c + (a - b)) * (a + (b - c));
+    if (kahan <= 0.0) return {0.0, 0.0, 0.0, false};
+
+    const double denom2 = 8.0 * (0.25 * std::sqrt(kahan));   // = 8·Area
 
     // Formula: cot_k = (t_opp · l123 − t_a · t_b) / denom2
     //   cot1: t_opp=t23, t_a=t12, t_b=t31   (v1 opposite l23)

code/include/hyper_ideal_functional.hpp

@@ -31,6 +31,7 @@
 //   // res.energy, res.gradient
 
 #include "conformal_mesh.hpp"
+#include "constants.hpp"
 #include "hyper_ideal_geometry.hpp"
 #include "hyper_ideal_utility.hpp"
 #include <CGAL/boost/graph/iterator.h>
@@ -157,6 +158,55 @@ static inline std::size_t hidx(Halfedge_index h) { return halfedge_to_index(h);
 // HyperIdealFunctional.java's defensive behaviour (lines 122-127 of the
 // Java original); this keeps the FD perturbation regime well-defined.
 
+// ── Scale-floor clamp mode (N3 audit) ─────────────────────────────────────────
+//
+// The vertex scale b must stay positive for the hyper-ideal geometry to be
+// valid.  Two ways to enforce that are offered, selectable per call:
+//
+//   • HardJava (DEFAULT) — the original `b < 0 → HYPER_IDEAL_SCALE_FLOOR` snap.
+//     Bit-for-bit faithful to HyperIdealFunctional.java, so the Java
+//     golden-oracle parity tests hold exactly.  It is only C⁰ at b = 0: a
+//     Newton step that crosses the feasibility boundary sees a kink, which can
+//     stall convergence (numerical-stability audit N3).  This is the
+//     production default precisely to preserve parity.
+//
+//   • SmoothBarrier — a C¹ softplus floor
+//         b ↦ floor + softplus_β(b − floor),   softplus_β(x) = log(1+e^{βx})/β
+//     which is ≈ b for b well above the floor (to machine precision once
+//     β·(b−floor) ≳ 35) and decays smoothly to `floor` as b → −∞, with a
+//     continuous derivative everywhere.  This removes the N3 kink and is the
+//     mathematically clean choice, but it perturbs values near the boundary
+//     and therefore deviates from the Java oracle — opt in when robustness of
+//     a boundary-crossing solve matters more than strict parity.
+//
+// Both modes share the same floor (`HYPER_IDEAL_SCALE_FLOOR`); SmoothBarrier
+// additionally uses `HYPER_IDEAL_SCALE_SHARPNESS` (β).  The `a < 0 → 0` edge
+// clamp is unaffected (a = 0 is a genuine geometric floor, not flagged by N3).
+enum class HyperIdealScaleClamp {
+    HardJava,       ///< b<0 → floor.  C⁰, Java-parity-faithful (default).
+    SmoothBarrier   ///< C¹ softplus floor; clean but deviates from Java near b=0.
+};
+
+/// Apply the vertex-scale floor to `b` under the chosen clamp `mode`.
+/// `HardJava` reproduces the original snap; `SmoothBarrier` is the C¹
+/// softplus floor (see `HyperIdealScaleClamp`).  `variable` mirrors the
+/// call-site guard (only clamp DOFs that are actually free / variable).
+inline double clamp_hyper_ideal_scale(double b, bool variable,
+                                      HyperIdealScaleClamp mode) noexcept
+{
+    if (!variable) return b;
+    if (mode == HyperIdealScaleClamp::HardJava)
+        return b < 0.0 ? HYPER_IDEAL_SCALE_FLOOR : b;
+
+    // SmoothBarrier: floor + softplus_β(b − floor), evaluated stably.
+    const double beta = HYPER_IDEAL_SCALE_SHARPNESS;
+    const double x    = beta * (b - HYPER_IDEAL_SCALE_FLOOR);
+    // softplus_β(x) = log1p(e^{βx})/β, with the standard large-x guard
+    // (for x ≳ 35, log1p(e^x) == x to double precision) to avoid overflow.
+    const double softplus = (x > 35.0) ? x : std::log1p(std::exp(x));
+    return HYPER_IDEAL_SCALE_FLOOR + softplus / beta;
+}
+
 /// Six per-face angle outputs computed from local DOFs (see
 /// `face_angles_from_local_dofs`).  Used by the block-FD Hessian.
 struct FaceAngleOutputs {
@@ -175,15 +225,16 @@ struct FaceAngleOutputs {
 inline FaceAngleOutputs face_angles_from_local_dofs(
     double b1, double b2, double b3,
     double a12, double a23, double a31,
-    bool   v1b, bool v2b, bool v3b)
+    bool   v1b, bool v2b, bool v3b,
+    HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
 {
     // Same defensive clamps as compute_face_angles.
     if (v1b && v2b && a12 < 0.0) a12 = 0.0;
     if (v2b && v3b && a23 < 0.0) a23 = 0.0;
     if (v3b && v1b && a31 < 0.0) a31 = 0.0;
-    if (v1b && b1  < 0.0) b1 = 0.01;
+    b1 = clamp_hyper_ideal_scale(b1, v1b, clamp);
-    if (v2b && b2  < 0.0) b2 = 0.01;
+    b2 = clamp_hyper_ideal_scale(b2, v2b, clamp);
-    if (v3b && b3  < 0.0) b3 = 0.01;
+    b3 = clamp_hyper_ideal_scale(b3, v3b, clamp);
 
     double l12 = lij(b1, b2, a12, v1b, v2b);
     double l23 = lij(b2, b3, a23, v2b, v3b);
@@ -246,7 +297,8 @@ static FaceAngles compute_face_angles(
     const ConformalMesh&       mesh,
     Face_index                 f,
     const std::vector<double>& x,
-    const HyperIdealMaps&      m)
+    const HyperIdealMaps&      m,
+    HyperIdealScaleClamp       clamp = HyperIdealScaleClamp::HardJava)
 {
     Halfedge_index h0 = mesh.halfedge(f);
     Halfedge_index h1 = mesh.next(h0);
@@ -276,9 +328,9 @@ static FaceAngles compute_face_angles(
     if (fa.v1b && fa.v2b && fa.a12 < 0.0) fa.a12 = 0.0;
     if (fa.v2b && fa.v3b && fa.a23 < 0.0) fa.a23 = 0.0;
     if (fa.v3b && fa.v1b && fa.a31 < 0.0) fa.a31 = 0.0;
-    if (fa.v1b && fa.b1  < 0.0)           fa.b1  = 0.01;
+    fa.b1 = clamp_hyper_ideal_scale(fa.b1, fa.v1b, clamp);
-    if (fa.v2b && fa.b2  < 0.0)           fa.b2  = 0.01;
+    fa.b2 = clamp_hyper_ideal_scale(fa.b2, fa.v2b, clamp);
-    if (fa.v3b && fa.b3  < 0.0)           fa.b3  = 0.01;
+    fa.b3 = clamp_hyper_ideal_scale(fa.b3, fa.v3b, clamp);
 
     double l12 = lij(fa.b1, fa.b2, fa.a12, fa.v1b, fa.v2b);
     double l23 = lij(fa.b2, fa.b3, fa.a23, fa.v2b, fa.v3b);
@@ -390,7 +442,8 @@ inline HyperIdealResult evaluate_hyper_ideal(
     const std::vector<double>& x,
     const HyperIdealMaps&      m,
     bool                       need_energy   = true,
-    bool                       need_gradient = true)
+    bool                       need_gradient = true,
+    HyperIdealScaleClamp       clamp = HyperIdealScaleClamp::HardJava)
 {
     HyperIdealResult res;
 
@@ -406,7 +459,7 @@ inline HyperIdealResult evaluate_hyper_ideal(
         Halfedge_index h1 = mesh.next(h0);
         Halfedge_index h2 = mesh.next(h1);
 
-        FaceAngles fa = compute_face_angles(mesh, f, x, m);
+        FaceAngles fa = compute_face_angles(mesh, f, x, m, clamp);
 
         // Store computed angles into temporary arrays.
         // h_alpha[h] = α for the edge of h in this face.

code/include/hyper_ideal_hessian.hpp

@@ -64,7 +64,8 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
     const HyperIdealMaps&      m,
-    double                     eps = 1e-5)
+    double                     eps   = 1e-5,
+    HyperIdealScaleClamp       clamp = HyperIdealScaleClamp::HardJava)
 {
     const int n = hyper_ideal_dimension(mesh, m);
     std::vector<Eigen::Triplet<double>> trips;
@@ -77,8 +78,8 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian(
         xp[sj] = x[sj] + eps;
         xm[sj] = x[sj] - eps;
 
-        auto Gp = evaluate_hyper_ideal(mesh, xp, m, /*energy=*/false).gradient;
+        auto Gp = evaluate_hyper_ideal(mesh, xp, m, /*energy=*/false, /*grad=*/true, clamp).gradient;
-        auto Gm = evaluate_hyper_ideal(mesh, xm, m, /*energy=*/false).gradient;
+        auto Gm = evaluate_hyper_ideal(mesh, xm, m, /*energy=*/false, /*grad=*/true, clamp).gradient;
 
         xp[sj] = xm[sj] = x[sj];   // restore
 
@@ -101,9 +102,10 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian_sym(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
     const HyperIdealMaps&      m,
-    double                     eps = 1e-5)
+    double                     eps   = 1e-5,
+    HyperIdealScaleClamp       clamp = HyperIdealScaleClamp::HardJava)
 {
-    auto H = hyper_ideal_hessian(mesh, x, m, eps);
+    auto H = hyper_ideal_hessian(mesh, x, m, eps, clamp);
     Eigen::SparseMatrix<double> Ht = H.transpose();
     return (H + Ht) * 0.5;
 }
@@ -142,7 +144,8 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
     const HyperIdealMaps&      m,
-    double                     eps = 1e-5)
+    double                     eps   = 1e-5,
+    HyperIdealScaleClamp       clamp = HyperIdealScaleClamp::HardJava)
 {
     const int n = hyper_ideal_dimension(mesh, m);
     std::vector<Eigen::Triplet<double>> trips;
@@ -187,9 +190,9 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd(
             vm[j] -= eps;
 
             auto Op = face_angles_from_local_dofs(
-                vp[0], vp[1], vp[2], vp[3], vp[4], vp[5], v1b, v2b, v3b);
+                vp[0], vp[1], vp[2], vp[3], vp[4], vp[5], v1b, v2b, v3b, clamp);
             auto Om = face_angles_from_local_dofs(
-                vm[0], vm[1], vm[2], vm[3], vm[4], vm[5], v1b, v2b, v3b);
+                vm[0], vm[1], vm[2], vm[3], vm[4], vm[5], v1b, v2b, v3b, clamp);
 
             const double Gp[6] = {
                 Op.beta1, Op.beta2, Op.beta3,
@@ -221,9 +224,10 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd_sym(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
     const HyperIdealMaps&      m,
-    double                     eps = 1e-5)
+    double                     eps   = 1e-5,
+    HyperIdealScaleClamp       clamp = HyperIdealScaleClamp::HardJava)
 {
-    auto H = hyper_ideal_hessian_block_fd(mesh, x, m, eps);
+    auto H = hyper_ideal_hessian_block_fd(mesh, x, m, eps, clamp);
     Eigen::SparseMatrix<double> Ht = H.transpose();
     return (H + Ht) * 0.5;
 }

code/include/inversive_distance_functional.hpp

@@ -315,6 +315,43 @@ inline std::vector<double> inversive_distance_gradient(
     return G;
 }
 
+/// Per-face contribution to the Inversive-Distance gradient, as a pure
+/// 3→3 kernel of the three local vertex DOFs `(u1,u2,u3)` and the three
+/// edge inversive distances `(I12,I23,I31)`.  No mesh, no property maps —
+/// used by the per-face block-FD Hessian in `inversive_distance_hessian.hpp`.
+///
+/// Returns the three values this face *adds* to the global gradient at
+/// `(v1,v2,v3)`.  Since `G_v = Θ_v − Σ_faces α_v`, the per-face contribution
+/// is the **negative** corner angles `(−α₁,−α₂,−α₃)`.  A non-real circle
+/// configuration (any `ℓ² ≤ 0`) contributes nothing — exactly mirroring the
+/// face-skip (`continue`) in `inversive_distance_gradient`, so the block-FD
+/// Hessian and the full-FD Hessian see the same per-face support.
+struct IDFaceGradContribs {
+    double g1; ///< contribution to G at v₁  (= −α₁)
+    double g2; ///< contribution to G at v₂  (= −α₂)
+    double g3; ///< contribution to G at v₃  (= −α₃)
+};
+
+inline IDFaceGradContribs inversive_distance_face_grad_contribs(
+    double u1, double u2, double u3,
+    double I12, double I23, double I31)
+{
+    double l12sq = id_detail::edge_length_squared(u1, u2, I12);
+    double l23sq = id_detail::edge_length_squared(u2, u3, I23);
+    double l31sq = id_detail::edge_length_squared(u3, u1, I31);
+
+    // Same face-skip as inversive_distance_gradient: a non-real circle
+    // configuration has no limiting angle, so the face contributes 0.
+    if (l12sq <= 0.0 || l23sq <= 0.0 || l31sq <= 0.0)
+        return {0.0, 0.0, 0.0};
+
+    // euclidean_angles deliberately keeps its limiting (valid=false) angles
+    // here — the convex C¹ extension — matching the gradient's choice not to
+    // skip on !fa.valid.  Corner at v_k = fa.alpha_k (see gradient Pass-2 trace).
+    auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq));
+    return {-fa.alpha1, -fa.alpha2, -fa.alpha3};
+}
+
 /// Inversive-Distance energy `E(u) = ∫₀¹ ⟨G(t·u), u⟩ dt`, evaluated
 /// with 10-point Gauss-Legendre (constants shared with `euclidean_energy`).
 inline double inversive_distance_energy(

code/include/inversive_distance_hessian.hpp

@@ -0,0 +1,187 @@
+#pragma once
+// Copyright (c) 2024-2026 Tarik Moussa.
+// SPDX-License-Identifier: MIT
+
+// inversive_distance_hessian.hpp
+//
+// Phase 9a.2 — Hessian of the inversive-distance circle-packing functional
+//              (Luo 2004 / Bowers-Stephenson 2004).
+//
+// ┌──────────────────────────────────────────────────────────────────────────┐
+// │  Implementation strategy                                                  │
+// │                                                                           │
+// │  TWO finite-difference Hessian implementations are provided here,         │
+// │  mirroring the HyperIdeal pair in `hyper_ideal_hessian.hpp`:              │
+// │                                                                           │
+// │    1.  `inversive_distance_hessian` — full finite-difference baseline.    │
+// │        Cost ≈ n × (cost of a full gradient evaluation) = O(n · F).        │
+// │        Kept as the cross-validation reference for the block-FD variant.   │
+// │                                                                           │
+// │    2.  `inversive_distance_hessian_block_fd` — per-face block-FD.         │
+// │        Each face contributes to the gradient through exactly 3 vertex     │
+// │        DOFs (u₁,u₂,u₃); we FD the 3×3 local Jacobian of that face's       │
+// │        gradient contribution and scatter it.  Cost ≈ F × 6 face-angle     │
+// │        evaluations = O(F).  Speed-up factor ≈ n/6 over full-FD.           │
+// │                                                                           │
+// │  Why the block-FD is correct (locality lemma):                            │
+// │     G_v = Θ_v − Σ_{f ∋ v} α_v(f),   and α_v(f) depends ONLY on the 3      │
+// │  vertex DOFs of face f.  Hence  ∂G_x/∂y = Σ_{f: x,y ∈ {v1,v2,v3}(f)}      │
+// │  ∂(−α_x)/∂y at f, so accumulating per-face 3×3 blocks reproduces the      │
+// │  full Hessian (identical to O(ε²)).                                       │
+// │                                                                           │
+// │  An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in        │
+// │  `doc/roadmap/research-track.md` as Phase 9a.2-analytic; it would take     │
+// │  the cost from O(F)·(FD constant) to a single O(F) analytic pass.         │
+// └──────────────────────────────────────────────────────────────────────────┘
+
+#include "inversive_distance_functional.hpp"
+#include <Eigen/Sparse>
+#include <vector>
+#include <cmath>
+
+namespace conformallab {
+
+/// Full finite-difference Inversive-Distance Hessian (baseline).
+/// Cost: `n` full-gradient evaluations ≈ `O(n·F)`.  Use for small meshes
+/// or as a correctness reference for the block-FD variant.
+inline Eigen::SparseMatrix<double> inversive_distance_hessian(
+    const ConformalMesh&         mesh,
+    const std::vector<double>&   x,
+    const InversiveDistanceMaps& m,
+    double                       eps = 1e-5)
+{
+    const int n = inversive_distance_dimension(mesh, m);
+    std::vector<Eigen::Triplet<double>> trips;
+    trips.reserve(static_cast<std::size_t>(n) * 16);
+
+    std::vector<double> xp = x, xm = x;
+    for (int j = 0; j < n; ++j) {
+        const std::size_t sj = static_cast<std::size_t>(j);
+        xp[sj] = x[sj] + eps;
+        xm[sj] = x[sj] - eps;
+
+        auto Gp = inversive_distance_gradient(mesh, xp, m);
+        auto Gm = inversive_distance_gradient(mesh, xm, m);
+
+        xp[sj] = xm[sj] = x[sj];   // restore
+
+        for (int i = 0; i < n; ++i) {
+            double val = (Gp[static_cast<std::size_t>(i)]
+                        - Gm[static_cast<std::size_t>(i)]) / (2.0 * eps);
+            if (std::abs(val) > 1e-15)
+                trips.emplace_back(i, j, val);
+        }
+    }
+
+    Eigen::SparseMatrix<double> H(n, n);
+    H.setFromTriplets(trips.begin(), trips.end());
+    return H;
+}
+
+/// Symmetrised full-FD Inversive-Distance Hessian: `(H + Hᵀ)/2`.
+inline Eigen::SparseMatrix<double> inversive_distance_hessian_sym(
+    const ConformalMesh&         mesh,
+    const std::vector<double>&   x,
+    const InversiveDistanceMaps& m,
+    double                       eps = 1e-5)
+{
+    auto H = inversive_distance_hessian(mesh, x, m, eps);
+    Eigen::SparseMatrix<double> Ht = H.transpose();
+    return (H + Ht) * 0.5;
+}
+
+// ── Block-FD Hessian ──────────────────────────────────────────────────────────
+//
+// The 3 local DOFs of a face f are (u_{v1}, u_{v2}, u_{v3}).  For each free
+// local DOF we recompute the face's gradient contribution (−α₁,−α₂,−α₃) at
+// x ± ε along that axis and read off the 3×3 Jacobian.  The result scatters
+// into the global Hessian via the DOF-index lookup.
+//
+// Cost: F × 6 face-angle evaluations (3 DOFs × 2 directions) vs n×F for
+// full-FD — a speed-up of ≈ n/6, i.e. ~hundreds× on large closed meshes.
+/// Per-face block-FD Inversive-Distance Hessian.  Uses the locality lemma
+/// `∂G_x/∂y = Σ_{f: x,y ∈ {v1,v2,v3}(f)} ∂(−α_x)/∂y` to perturb only the 3
+/// face-local DOFs at a time, giving an `F·6` face-evaluation budget vs `n·F`
+/// for full-FD.  Mathematically equivalent to `inversive_distance_hessian`
+/// up to O(ε²) FD rounding.
+inline Eigen::SparseMatrix<double> inversive_distance_hessian_block_fd(
+    const ConformalMesh&         mesh,
+    const std::vector<double>&   x,
+    const InversiveDistanceMaps& m,
+    double                       eps = 1e-5)
+{
+    const int n = inversive_distance_dimension(mesh, m);
+    std::vector<Eigen::Triplet<double>> trips;
+    trips.reserve(9 * mesh.number_of_faces());
+
+    for (auto f : mesh.faces()) {
+        Halfedge_index h0 = mesh.halfedge(f);
+        Halfedge_index h1 = mesh.next(h0);
+        Halfedge_index h2 = mesh.next(h1);
+
+        Vertex_index v1 = mesh.source(h0);
+        Vertex_index v2 = mesh.source(h1);
+        Vertex_index v3 = mesh.source(h2);
+        Edge_index   e12 = mesh.edge(h0);
+        Edge_index   e23 = mesh.edge(h1);
+        Edge_index   e31 = mesh.edge(h2);
+
+        // Local DOF indices: (u1, u2, u3).  Pinned slots = -1.
+        const int idx[3] = { m.v_idx[v1], m.v_idx[v2], m.v_idx[v3] };
+
+        const double I12 = m.I_e[e12];
+        const double I23 = m.I_e[e23];
+        const double I31 = m.I_e[e31];
+
+        // Local DOF values (0 for pinned).
+        const double vals[3] = {
+            id_detail::dof_val(idx[0], x),
+            id_detail::dof_val(idx[1], x),
+            id_detail::dof_val(idx[2], x)
+        };
+
+        for (int j = 0; j < 3; ++j) {
+            if (idx[j] < 0) continue;   // never perturb a pinned DOF
+
+            double vp[3], vm[3];
+            for (int k = 0; k < 3; ++k) { vp[k] = vm[k] = vals[k]; }
+            vp[j] += eps;
+            vm[j] -= eps;
+
+            auto Cp = inversive_distance_face_grad_contribs(
+                vp[0], vp[1], vp[2], I12, I23, I31);
+            auto Cm = inversive_distance_face_grad_contribs(
+                vm[0], vm[1], vm[2], I12, I23, I31);
+
+            const double Gp[3] = { Cp.g1, Cp.g2, Cp.g3 };
+            const double Gm[3] = { Cm.g1, Cm.g2, Cm.g3 };
+
+            for (int i = 0; i < 3; ++i) {
+                if (idx[i] < 0) continue;   // pinned: contributes nothing
+                const double val = (Gp[i] - Gm[i]) / (2.0 * eps);
+                if (std::abs(val) > 1e-15)
+                    trips.emplace_back(idx[i], idx[j], val);
+            }
+        }
+    }
+
+    Eigen::SparseMatrix<double> H(n, n);
+    H.setFromTriplets(trips.begin(), trips.end());
+    return H;
+}
+
+/// Symmetrised block-FD Inversive-Distance Hessian: `(H + Hᵀ)/2` of
+/// `inversive_distance_hessian_block_fd(...)` for solvers requiring strict
+/// symmetry.
+inline Eigen::SparseMatrix<double> inversive_distance_hessian_block_fd_sym(
+    const ConformalMesh&         mesh,
+    const std::vector<double>&   x,
+    const InversiveDistanceMaps& m,
+    double                       eps = 1e-5)
+{
+    auto H = inversive_distance_hessian_block_fd(mesh, x, m, eps);
+    Eigen::SparseMatrix<double> Ht = H.transpose();
+    return (H + Ht) * 0.5;
+}
+
+} // namespace conformallab

code/include/mesh_io.hpp

@@ -29,6 +29,7 @@
 #include <CGAL/boost/graph/helpers.h>
 #include <string>
 #include <stdexcept>
+#include <cmath>
 
 namespace conformallab {
 
@@ -62,6 +63,12 @@ inline ConformalMesh load_mesh(const std::string& filename)
         throw std::runtime_error(
             "conformallab: mesh from " + filename +
             " is not triangulated (functionals require triangle faces)");
+    for (auto v : mesh.vertices()) {
+        const auto& p = mesh.point(v);
+        if (!std::isfinite(p.x()) || !std::isfinite(p.y()) || !std::isfinite(p.z()))
+            throw std::runtime_error(
+                "conformallab: non-finite vertex coordinate in " + filename);
+    }
     return mesh;
 }
 

code/include/newton_solver.hpp

@@ -36,6 +36,7 @@
 #include "hyper_ideal_hessian.hpp"
 #include "cp_euclidean_functional.hpp"
 #include "inversive_distance_functional.hpp"
+#include "inversive_distance_hessian.hpp"
 #include <Eigen/SparseCholesky>
 #include <Eigen/SparseQR>
 #include <Eigen/OrderingMethods>
@@ -47,12 +48,42 @@ namespace conformallab {
 
 // ── Result ────────────────────────────────────────────────────────────────────
 
-/// Result of `newton_solve(...)` — converged DOF vector + diagnostics.
+/// Why a Newton solve terminated (I1 audit).  Distinguishes the three
+/// non-convergent exits that `converged == false` previously conflated.
+enum class NewtonStatus {
+    Converged,           ///< `grad_inf_norm < tol` — success.
+    MaxIterations,       ///< hit `max_iter` without reaching `tol`.
+    LinearSolverFailed,  ///< both LDLT and SparseQR failed on H·Δx = −G.
+    LineSearchStalled    ///< the line search could not reduce the residual.
+};
+
+/// Human-readable name for a `NewtonStatus` (for logs / test messages).
+inline const char* to_string(NewtonStatus s) noexcept
+{
+    switch (s) {
+        case NewtonStatus::Converged:          return "Converged";
+        case NewtonStatus::MaxIterations:      return "MaxIterations";
+        case NewtonStatus::LinearSolverFailed: return "LinearSolverFailed";
+        case NewtonStatus::LineSearchStalled:  return "LineSearchStalled";
+    }
+    return "Unknown";
+}
+
+/// Result of a Newton solve — converged DOF vector + diagnostics.
 struct NewtonResult {
     std::vector<double> x;            ///< DOF vector at termination.
-    int                 iterations;   ///< Newton steps taken.
+    int                 iterations;   ///< Newton steps actually completed.
     double              grad_inf_norm;///< max |Gᵢ| at termination.
-    bool                converged;    ///< `true` iff `grad_inf_norm < tol`.
+    bool                converged;    ///< `true` iff `status == Converged`.
+
+    // I1: why it stopped (more informative than the `converged` bool alone).
+    NewtonStatus        status = NewtonStatus::MaxIterations;
+
+    // N7: linear-algebra conditioning diagnostics.
+    bool   sparse_qr_fallback_used = false; ///< any iteration fell back to SparseQR.
+    double min_ldlt_pivot          = 0.0;   ///< smallest |Dᵢᵢ| of the last LDLT
+                                            ///< factorisation (≈ near-singularity
+                                            ///< proxy); 0 if LDLT never succeeded.
 };
 
 // ── Internal helpers ──────────────────────────────────────────────────────────
@@ -103,14 +134,18 @@ inline Eigen::VectorXd solve_with_fallback(
 /// Solve `A·x = rhs` with the same SimplicialLDLT → SparseQR fallback
 /// strategy used inside all three Newton solvers.  If `fallback_used`
 /// is non-null, it is set to `true` iff the SparseQR fallback ran.
+/// If `ok` is non-null, it is set to `true` iff at least one solver succeeded.
 /// Returns `Eigen::VectorXd::Zero(rhs.size())` if both solvers fail.
 inline Eigen::VectorXd solve_linear_system(
     const Eigen::SparseMatrix<double>& A,
     const Eigen::VectorXd&             rhs,
-    bool*                              fallback_used = nullptr)
+    bool*                              fallback_used = nullptr,
+    bool*                              ok            = nullptr)
 {
-    bool ok = false;
+    bool ok_local = false;
-    return detail::solve_with_fallback(A, rhs, ok, fallback_used);
+    auto x = detail::solve_with_fallback(A, rhs, ok_local, fallback_used);
+    if (ok) *ok = ok_local;
+    return x;
 }
 
 namespace detail {   // re-open for the remaining helpers
@@ -137,6 +172,14 @@ namespace detail {   // re-open for the remaining helpers
 // If neither phase satisfies Armijo, return the best (smallest-residual) point
 // visited.  If nothing beat ‖G(x)‖, return x unchanged and set *improved=false,
 // so the caller can stop cleanly instead of taking the old divergent full step.
+//
+// B2 (api-performance audit): the gradient computed at the accepted trial point
+// is normally discarded (only its norm is kept), forcing the next Newton
+// iteration to recompute G(x) from scratch.  When `accepted_grad` is non-null it
+// receives the gradient vector at the returned point, so the caller can reuse it
+// as the next iteration's convergence gradient.  It is only written when the
+// returned point differs from `x` (i.e. `*improved == true`); on a pure stall
+// (`*improved == false`) the caller must recompute G(x) itself.
 template <typename GradFn>
 inline std::vector<double> line_search(
     const std::vector<double>& x,
@@ -145,6 +188,7 @@ inline std::vector<double> line_search(
     double                     norm0,
     GradFn&&                   grad_fn,
     bool*                      improved      = nullptr,
+    std::vector<double>*       accepted_grad = nullptr,
     int                        max_halvings  = 20,
     double                     c1            = 1e-4)
 {
@@ -155,14 +199,18 @@ inline std::vector<double> line_search(
     std::vector<double> best_x    = x;
     double              best_norm = norm0;
 
-    // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew`.
+    std::vector<double> trial_grad;   // gradient at the most recent trial point
+    std::vector<double> best_grad;    // gradient at best_x (B2)
+
+    // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew` and the
+    // gradient there in `trial_grad`.
     auto eval = [&](const Eigen::VectorXd& dir, double alpha) -> double {
         for (int i = 0; i < n; ++i)
             xnew[static_cast<std::size_t>(i)] =
                 x[static_cast<std::size_t>(i)] + alpha * dir[i];
-        auto Gnew = grad_fn(xnew);
+        trial_grad = grad_fn(xnew);
         double s = 0.0;
-        for (double v : Gnew) s += v * v;
+        for (double v : trial_grad) s += v * v;
         return std::sqrt(s);
     };
 
@@ -170,9 +218,10 @@ inline std::vector<double> line_search(
     double alpha = 1.0;
     for (int ls = 0; ls < max_halvings; ++ls) {
         double norm_new = eval(dx, alpha);
-        if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; }
+        if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; }
         if (norm_new * norm_new <= (1.0 - 2.0 * c1 * alpha) * norm0_sq) {
             if (improved) *improved = true;
+            if (accepted_grad) *accepted_grad = trial_grad;
             return xnew;
         }
         alpha *= 0.5;
@@ -185,9 +234,10 @@ inline std::vector<double> line_search(
         for (int ls = 0; ls < max_halvings; ++ls) {
             double thresh_sq = norm0_sq - 2.0 * c1 * alpha * dsd_sq;
             double norm_new  = eval(d_sd, alpha);
-            if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; }
+            if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; }
             if (thresh_sq >= 0.0 && norm_new * norm_new <= thresh_sq) {
                 if (improved) *improved = true;
+                if (accepted_grad) *accepted_grad = trial_grad;
                 return xnew;
             }
             alpha *= 0.5;
@@ -196,10 +246,169 @@ inline std::vector<double> line_search(
 
     // ── Both phases failed Armijo — never take the divergent full step. ───────
     //    Return the best point seen; if none improved, stay put and signal stall.
-    if (improved) *improved = (best_norm < norm0);
+    const bool any_improvement = (best_norm < norm0);
+    if (improved) *improved = any_improvement;
+    if (accepted_grad && any_improvement) *accepted_grad = best_grad;
     return best_x;
 }
 
+// ── Cached linear solver (B4) ─────────────────────────────────────────────────
+//
+// solve_with_fallback rebuilds a fresh SimplicialLDLT every call, re-running the
+// symbolic analysis (fill-reducing reordering / COLAMD) each Newton iteration —
+// even though the Hessian sparsity pattern is iteration-invariant for a fixed
+// mesh + DOF layout.  This object keeps one persistent LDLT and runs
+// `analyzePattern` once, then only `factorize` per iteration.
+//
+// FD-built Hessians (hyper-ideal, inversive-distance) prune near-zero triplets,
+// so their pattern can shift slightly between iterations; we detect that via a
+// change in `nonZeros()` and re-analyze, which keeps the cache correct and
+// self-healing.  The SparseQR fallback for singular/rank-deficient systems is
+// preserved verbatim (built fresh on the rare path).  The fallback semantics are
+// identical to `solve_with_fallback`, so behaviour is unchanged.
+struct NewtonLinearSolver {
+    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
+    bool          analyzed = false;
+    Eigen::Index  last_nnz = -1;
+    bool          fallback_used  = false;  ///< true iff the last solve used SparseQR
+    double        last_min_pivot = 0.0;    ///< N7: smallest |Dᵢᵢ| of the last
+                                           ///< successful LDLT (0 if it fell back)
+
+    Eigen::VectorXd solve(const Eigen::SparseMatrix<double>& A,
+                          const Eigen::VectorXd&             rhs,
+                          bool&                              ok)
+    {
+        fallback_used  = false;
+        last_min_pivot = 0.0;
+        if (!analyzed || A.nonZeros() != last_nnz) {
+            ldlt.analyzePattern(A);
+            analyzed = true;
+            last_nnz = A.nonZeros();
+        }
+        ldlt.factorize(A);
+        if (ldlt.info() == Eigen::Success) {
+            Eigen::VectorXd x = ldlt.solve(rhs);
+            if (ldlt.info() == Eigen::Success) {
+                ok = true;
+                // N7: smallest |Dᵢᵢ| in the LDLᵀ diagonal — a cheap
+                // near-singularity proxy (small ⇒ ill-conditioned step).
+                const auto D = ldlt.vectorD();
+                last_min_pivot = (D.size() > 0) ? D.cwiseAbs().minCoeff() : 0.0;
+                return x;
+            }
+        }
+        // Fallback: SparseQR — handles singular/rank-deficient A (gauge modes).
+        fallback_used = true;
+        Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr(A);
+        if (qr.info() == Eigen::Success) {
+            Eigen::VectorXd x = qr.solve(rhs);
+            if (qr.info() == Eigen::Success) { ok = true; return x; }
+        }
+        ok = false;
+        return Eigen::VectorXd::Zero(rhs.size());
+    }
+};
+
+// ── Unified Newton core (H2) ──────────────────────────────────────────────────
+//
+// All five DCE Newton solvers (Euclidean, Spherical, HyperIdeal, CP-Euclidean,
+// Inversive-Distance) ran a near-identical loop that differed only in (a) the
+// gradient function, (b) the Hessian function, and (c) the spherical sign quirk
+// (the concave spherical energy factorises −H).  This template captures that one
+// loop so a fix lands once instead of five times.
+//
+// `grad(x) -> std::vector<double>`        : the energy gradient G(x).
+// `hess(x) -> Eigen::SparseMatrix<double>`: the TRUE Hessian H(x) (un-negated).
+// `concave`                               : if true, solve (−H)·Δx = G (spherical);
+//                                           otherwise H·Δx = −G.  Both are the same
+//                                           Newton system; `concave` only selects
+//                                           the PSD matrix actually factorised.
+//
+// Folds in B2 (reuse the line-search gradient as the next convergence gradient),
+// B4 (cached symbolic factorization), and B5 (no dead solver variable).
+template <typename GradFn, typename HessFn>
+inline NewtonResult newton_core(
+    std::vector<double> x,
+    GradFn&&            grad,
+    HessFn&&            hess,
+    bool                concave,
+    double              tol,
+    int                 max_iter)
+{
+    const int n = static_cast<int>(x.size());
+
+    NewtonResult res;
+    res.converged     = false;
+    res.iterations    = 0;
+    res.grad_inf_norm = 0.0;
+
+    NewtonLinearSolver solver;
+
+    // B2: gradient carried across the loop; the first one is the only "extra"
+    // evaluation — every later iteration reuses the line-search gradient.
+    std::vector<double> G_std = grad(x);
+
+    for (int iter = 0; iter < max_iter; ++iter) {
+        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
+
+        const double inf_norm = G.cwiseAbs().maxCoeff();
+        if (inf_norm < tol) {
+            res.converged     = true;
+            res.status        = NewtonStatus::Converged;   // I1
+            res.grad_inf_norm = inf_norm;
+            res.iterations    = iter;
+            res.x             = x;
+            return res;
+        }
+
+        Eigen::SparseMatrix<double> H = hess(x);
+
+        // Newton system: factor the PSD matrix (H for convex, −H for concave).
+        Eigen::SparseMatrix<double> A = concave ? Eigen::SparseMatrix<double>(-H) : H;
+        Eigen::VectorXd rhs = concave ? Eigen::VectorXd(G) : Eigen::VectorXd(-G);
+
+        bool ok = false;
+        Eigen::VectorXd dx = solver.solve(A, rhs, ok);
+        // N7: accumulate conditioning diagnostics from this solve.
+        res.sparse_qr_fallback_used = res.sparse_qr_fallback_used || solver.fallback_used;
+        res.min_ldlt_pivot          = solver.last_min_pivot;
+        if (!ok) {                       // I1/H1: linear solver failed this step
+            res.iterations = iter;       //        (iter steps actually completed)
+            res.status     = NewtonStatus::LinearSolverFailed;
+            break;
+        }
+
+        // Globalised line search (Armijo + steepest-descent fallback).
+        // d_sd = −(H·G) is the merit-function steepest descent for f = ½‖G‖²,
+        // using the TRUE (un-negated) Hessian for both signs.
+        const double norm0 = G.norm();
+        Eigen::VectorXd d_sd = -(H * G);
+        bool improved = true;
+        std::vector<double> G_next;
+        x = detail::line_search(x, dx, d_sd, norm0, grad, &improved, &G_next);
+        if (!improved) {                 // I1/H1: line search stalled — stop cleanly
+            res.iterations = iter;
+            res.status     = NewtonStatus::LineSearchStalled;
+            break;
+        }
+
+        G_std = std::move(G_next);       // B2: reuse for the next convergence check
+        res.iterations = iter + 1;
+    }
+
+    // If we fell out of the loop without converging or breaking, it was max_iter
+    // (the default `res.status` is already MaxIterations; a break has set its own).
+
+    // Final gradient norm (reached only on max-iter / break — recomputed to be
+    // correct after a stalled or failed step).
+    auto G_final = grad(x);
+    double inf_final = 0.0;
+    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
+    res.grad_inf_norm = inf_final;
+    res.x             = x;
+    return res;
+}
+
 } // namespace detail
 
 // ── Euclidean Newton solver ────────────────────────────────────────────────────
@@ -233,63 +442,21 @@ inline NewtonResult newton_euclidean(
     double                     tol      = 1e-8,
     int                        max_iter = 200)
 {
-    std::vector<double> x = x0;
+    // Layout is loop-invariant — decide the Hessian variant once (B3): cyclic
-    const int n = static_cast<int>(x.size());
+    // layout (edge DOFs present) → analytic cyclic Hessian, which covers the
-
+    // vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian omits.
-    NewtonResult res;
-    res.converged      = false;
-    res.iterations     = 0;
-    res.grad_inf_norm  = 0.0;
-
-    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
-
-    for (int iter = 0; iter < max_iter; ++iter) {
-        // ── Gradient ──────────────────────────────────────────────────────────
-        auto G_std = euclidean_gradient(mesh, x, m);
-        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
-
-        double inf_norm = G.cwiseAbs().maxCoeff();
-        if (inf_norm < tol) {
-            res.converged     = true;
-            res.grad_inf_norm = inf_norm;
-            res.iterations    = iter;
-            res.x             = x;
-            return res;
-        }
-
-        // ── Hessian + solve H·Δx = −G  (SparseQR fallback for singular H) ──
-        // Cyclic layout (edge DOFs present) → analytic cyclic Hessian, covering
-        // the vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian
-        // omits.  Vertex-only layout → analytic cotangent Laplacian.
     bool has_edge_dof = false;
     for (auto e : mesh.edges())
         if (m.e_idx[e] >= 0) { has_edge_dof = true; break; }
-        auto H = has_edge_dof ? euclidean_hessian_analytic(mesh, x, m)
-                              : euclidean_hessian(mesh, x, m);
-        bool ok = false;
-        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
-        if (!ok) break;
 
-        // ── Globalised line search (Armijo + steepest-descent fallback) ───────
+    return detail::newton_core(
-        double norm0 = G.norm();
+        std::move(x0),
-        Eigen::VectorXd d_sd = -(H * G);   // merit-function steepest descent
+        [&](const std::vector<double>& xc) { return euclidean_gradient(mesh, xc, m); },
-        bool improved = true;
+        [&](const std::vector<double>& xc) {
-        x = detail::line_search(x, dx, d_sd, norm0,
+            return has_edge_dof ? euclidean_hessian_analytic(mesh, xc, m)
-            [&](const std::vector<double>& xnew) {
+                                : euclidean_hessian(mesh, xc, m);
-                return euclidean_gradient(mesh, xnew, m);
+        },
-            }, &improved);
+        /*concave=*/false, tol, max_iter);
-        if (!improved) break;   // line search stalled — stop cleanly
-
-        res.iterations = iter + 1;
-    }
-
-    // Report final gradient norm
-    auto G_final = euclidean_gradient(mesh, x, m);
-    double inf_final = 0.0;
-    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
-    res.grad_inf_norm = inf_final;
-    res.x             = x;
-    return res;
 }
 
 // ── Spherical Newton solver ───────────────────────────────────────────────────
@@ -322,55 +489,14 @@ inline NewtonResult newton_spherical(
     double                     tol      = 1e-8,
     int                        max_iter = 200)
 {
-    std::vector<double> x = x0;
+    // Concave energy: the Hessian H is NSD, so newton_core factors −H (PSD) and
-    const int n = static_cast<int>(x.size());
+    // solves (−H)·Δx = G — the same Newton system, just the sign that keeps LDLT
-
+    // well-defined.  The merit steepest-descent inside still uses the true H.
-    NewtonResult res;
+    return detail::newton_core(
-    res.converged      = false;
+        std::move(x0),
-    res.iterations     = 0;
+        [&](const std::vector<double>& xc) { return spherical_gradient(mesh, xc, m); },
-    res.grad_inf_norm  = 0.0;
+        [&](const std::vector<double>& xc) { return spherical_hessian(mesh, xc, m); },
-
+        /*concave=*/true, tol, max_iter);
-    for (int iter = 0; iter < max_iter; ++iter) {
-        // ── Gradient ──────────────────────────────────────────────────────────
-        auto G_std = spherical_gradient(mesh, x, m);
-        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
-
-        double inf_norm = G.cwiseAbs().maxCoeff();
-        if (inf_norm < tol) {
-            res.converged     = true;
-            res.grad_inf_norm = inf_norm;
-            res.iterations    = iter;
-            res.x             = x;
-            return res;
-        }
-
-        // ── Hessian: negate to get PSD; solve (−H)·Δx = G  ──────────────────
-        auto H    = spherical_hessian(mesh, x, m);
-        auto negH = Eigen::SparseMatrix<double>(-H);
-        bool ok   = false;
-        Eigen::VectorXd dx = detail::solve_with_fallback(negH, G, ok);
-        if (!ok) break;
-
-        // ── Globalised line search (Armijo + steepest-descent fallback) ───────
-        // d_sd uses the actual (un-negated) Hessian: ∇f = H·G for f = ½‖G‖².
-        double norm0 = G.norm();
-        Eigen::VectorXd d_sd = -(H * G);
-        bool improved = true;
-        x = detail::line_search(x, dx, d_sd, norm0,
-            [&](const std::vector<double>& xnew) {
-                return spherical_gradient(mesh, xnew, m);
-            }, &improved);
-        if (!improved) break;
-
-        res.iterations = iter + 1;
-    }
-
-    auto G_final = spherical_gradient(mesh, x, m);
-    double inf_final = 0.0;
-    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
-    res.grad_inf_norm = inf_final;
-    res.x             = x;
-    return res;
 }
 
 // ── HyperIdeal Newton solver ──────────────────────────────────────────────────
@@ -395,6 +521,12 @@ inline NewtonResult newton_spherical(
 /// \param max_iter  Maximum Newton iterations.  Default: 200.
 /// \param hess_eps  Finite-difference step for Hessian approximation.  Default: 1e-5.
 ///                  (Phase 9b will replace this with an analytic Hessian.)
+/// \param clamp     Vertex-scale floor mode (N3 audit).  Default `HardJava`
+///                  reproduces the Java oracle's hard `b<0 → floor` snap
+///                  (C⁰, parity-faithful).  `SmoothBarrier` uses the C¹
+///                  softplus floor — smoother near the feasibility boundary
+///                  but deviates from the Java golden values.  See
+///                  `HyperIdealScaleClamp` in hyper_ideal_functional.hpp.
 /// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
 ///
 /// \see Springborn (2020), Theorem 1.3 for the strict convexity proof.
@@ -405,55 +537,20 @@ inline NewtonResult newton_hyper_ideal(
     const HyperIdealMaps&      m,
     double                     tol      = 1e-8,
     int                        max_iter = 200,
-    double                     hess_eps = 1e-5)
+    double                     hess_eps = 1e-5,
+    HyperIdealScaleClamp       clamp    = HyperIdealScaleClamp::HardJava)
 {
-    std::vector<double> x = x0;
+    // block-FD Hessian (~33–1166× faster than full-FD); `clamp` selects the
-    const int n = static_cast<int>(x.size());
+    // vertex-scale floor mode (N3).  Convex energy → factor H directly.
-
+    return detail::newton_core(
-    NewtonResult res;
+        std::move(x0),
-    res.converged      = false;
+        [&](const std::vector<double>& xc) {
-    res.iterations     = 0;
+            return evaluate_hyper_ideal(mesh, xc, m, /*energy=*/false, /*grad=*/true, clamp).gradient;
-    res.grad_inf_norm  = 0.0;
+        },
-
+        [&](const std::vector<double>& xc) {
-    for (int iter = 0; iter < max_iter; ++iter) {
+            return hyper_ideal_hessian_block_fd_sym(mesh, xc, m, hess_eps, clamp);
-        // ── Gradient ──────────────────────────────────────────────────────────
+        },
-        auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false).gradient;
+        /*concave=*/false, tol, max_iter);
-        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
-
-        double inf_norm = G.cwiseAbs().maxCoeff();
-        if (inf_norm < tol) {
-            res.converged     = true;
-            res.grad_inf_norm = inf_norm;
-            res.iterations    = iter;
-            res.x             = x;
-            return res;
-        }
-
-        // ── Hessian (numerical FD) + solve H·Δx = −G ─────────────────────────
-        auto H  = hyper_ideal_hessian_sym(mesh, x, m, hess_eps);
-        bool ok = false;
-        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
-        if (!ok) break;
-
-        // ── Globalised line search (Armijo + steepest-descent fallback) ───────
-        double norm0 = G.norm();
-        Eigen::VectorXd d_sd = -(H * G);
-        bool improved = true;
-        x = detail::line_search(x, dx, d_sd, norm0,
-            [&](const std::vector<double>& xnew) {
-                return evaluate_hyper_ideal(mesh, xnew, m, false).gradient;
-            }, &improved);
-        if (!improved) break;
-
-        res.iterations = iter + 1;
-    }
-
-    auto G_final = evaluate_hyper_ideal(mesh, x, m, false).gradient;
-    double inf_final = 0.0;
-    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
-    res.grad_inf_norm = inf_final;
-    res.x             = x;
-    return res;
 }
 
 // ── CP-Euclidean Newton solver (Phase 9a.1) ───────────────────────────────────
@@ -461,7 +558,7 @@ inline NewtonResult newton_hyper_ideal(
 /// Solve the CP-Euclidean circle-packing problem: find ρ ∈ ℝ^F such that the
 /// per-face angle sums match φ_f at every free face.
 ///
-/// The CP-Euclidean energy (Bobenko-Pinkall-Springborn 2010 §6) is strictly
+/// The CP-Euclidean energy (Bobenko-Pinkall-Springborn 2015 §6) is strictly
 /// convex on its open domain of validity, so the Hessian H is PSD and the
 /// solution is unique up to the gauge mode pinned by `f_idx == −1`.
 /// `cp_euclidean_hessian` provides the analytic 2×2-per-edge formula
@@ -479,7 +576,7 @@ inline NewtonResult newton_hyper_ideal(
 /// \note  Unlike the Euclidean solver, the CP-Euclidean Hessian is exact
 ///        (analytic), so the SparseQR fallback only triggers in genuine
 ///        gauge-singular situations (no pinned face).
-/// \see   doc/architecture/phase-9a-validation.md §1 for the BPS-2010 mapping.
+/// \see   doc/architecture/phase-9a-validation.md §1 for the BPS-2015 mapping.
 inline NewtonResult newton_cp_euclidean(
     ConformalMesh&             mesh,
     std::vector<double>        x0,
@@ -487,50 +584,12 @@ inline NewtonResult newton_cp_euclidean(
     double                     tol      = 1e-8,
     int                        max_iter = 200)
 {
-    std::vector<double> x = x0;
+    // Convex energy with an exact analytic Hessian (BPS-2015 2×2-per-edge).
-    const int n = static_cast<int>(x.size());
+    return detail::newton_core(
-
+        std::move(x0),
-    NewtonResult res;
+        [&](const std::vector<double>& xc) { return cp_euclidean_gradient(mesh, xc, m); },
-    res.converged      = false;
+        [&](const std::vector<double>& xc) { return cp_euclidean_hessian(mesh, xc, m); },
-    res.iterations     = 0;
+        /*concave=*/false, tol, max_iter);
-    res.grad_inf_norm  = 0.0;
-
-    for (int iter = 0; iter < max_iter; ++iter) {
-        auto G_std = cp_euclidean_gradient(mesh, x, m);
-        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
-
-        double inf_norm = G.cwiseAbs().maxCoeff();
-        if (inf_norm < tol) {
-            res.converged     = true;
-            res.grad_inf_norm = inf_norm;
-            res.iterations    = iter;
-            res.x             = x;
-            return res;
-        }
-
-        auto H = cp_euclidean_hessian(mesh, x, m);
-        bool ok = false;
-        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
-        if (!ok) break;
-
-        double norm0 = G.norm();
-        Eigen::VectorXd d_sd = -(H * G);
-        bool improved = true;
-        x = detail::line_search(x, dx, d_sd, norm0,
-            [&](const std::vector<double>& xnew) {
-                return cp_euclidean_gradient(mesh, xnew, m);
-            }, &improved);
-        if (!improved) break;
-
-        res.iterations = iter + 1;
-    }
-
-    auto G_final = cp_euclidean_gradient(mesh, x, m);
-    double inf_final = 0.0;
-    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
-    res.grad_inf_norm = inf_final;
-    res.x             = x;
-    return res;
 }
 
 // ── Inversive-Distance Newton solver (Phase 9a.2) ─────────────────────────────
@@ -542,10 +601,11 @@ inline NewtonResult newton_cp_euclidean(
 /// domain where every triangle satisfies the inequalities.  Luo's 1-form is
 /// closed there, so the path-integral energy is well-defined.
 ///
-/// MVP implementation: the Hessian is computed by **finite differences** of
+/// The Hessian is computed by **per-face block finite differences**
-/// the analytic gradient (same pattern as the Phase 4a HyperIdeal solver).
+/// (`inversive_distance_hessian_block_fd_sym`, same pattern as the HyperIdeal
-/// An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in
+/// solver after Phase 9b) — O(F) face evaluations instead of the O(n·F) of the
-/// `doc/roadmap/research-track.md` as Phase 9a.2-analytic.
+/// full-FD baseline.  An analytic Hessian via Glickenstein 2011 eq. (4.6) is
+/// tracked in `doc/roadmap/research-track.md` as Phase 9a.2-analytic.
 ///
 /// \param mesh      Input triangle mesh.
 /// \param x0        Initial DOF vector (length = number of free vertices).
@@ -568,81 +628,14 @@ inline NewtonResult newton_inversive_distance(
     int                           max_iter = 200,
     double                        hess_eps = 1e-5)
 {
-    std::vector<double> x = x0;
+    // Convex energy; block-FD Hessian (~n/6× faster than full-FD).
-    const int n = static_cast<int>(x.size());
+    return detail::newton_core(
-
+        std::move(x0),
-    NewtonResult res;
+        [&](const std::vector<double>& xc) { return inversive_distance_gradient(mesh, xc, m); },
-    res.converged      = false;
+        [&](const std::vector<double>& xc) {
-    res.iterations     = 0;
+            return inversive_distance_hessian_block_fd_sym(mesh, xc, m, hess_eps);
-    res.grad_inf_norm  = 0.0;
+        },
-
+        /*concave=*/false, tol, max_iter);
-    // Local FD Hessian builder — n × (cost of gradient eval).
-    auto build_hessian = [&](const std::vector<double>& xc) -> Eigen::SparseMatrix<double> {
-        std::vector<Eigen::Triplet<double>> trips;
-        trips.reserve(static_cast<std::size_t>(n) * 16);   // sparse heuristic
-
-        std::vector<double> xp = xc, xm = xc;
-        for (int j = 0; j < n; ++j) {
-            const std::size_t sj = static_cast<std::size_t>(j);
-            xp[sj] = xc[sj] + hess_eps;
-            xm[sj] = xc[sj] - hess_eps;
-
-            auto Gp = inversive_distance_gradient(mesh, xp, m);
-            auto Gm = inversive_distance_gradient(mesh, xm, m);
-
-            xp[sj] = xm[sj] = xc[sj];   // restore
-
-            for (int i = 0; i < n; ++i) {
-                double val = (Gp[static_cast<std::size_t>(i)]
-                            - Gm[static_cast<std::size_t>(i)])
-                             / (2.0 * hess_eps);
-                if (std::abs(val) > 1e-15)
-                    trips.emplace_back(i, j, val);
-            }
-        }
-        Eigen::SparseMatrix<double> H(n, n);
-        H.setFromTriplets(trips.begin(), trips.end());
-        // Symmetrise — FD rounding may introduce tiny asymmetries.
-        Eigen::SparseMatrix<double> Ht = H.transpose();
-        return (H + Ht) * 0.5;
-    };
-
-    for (int iter = 0; iter < max_iter; ++iter) {
-        auto G_std = inversive_distance_gradient(mesh, x, m);
-        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
-
-        double inf_norm = G.cwiseAbs().maxCoeff();
-        if (inf_norm < tol) {
-            res.converged     = true;
-            res.grad_inf_norm = inf_norm;
-            res.iterations    = iter;
-            res.x             = x;
-            return res;
-        }
-
-        auto H = build_hessian(x);
-        bool ok = false;
-        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
-        if (!ok) break;
-
-        double norm0 = G.norm();
-        Eigen::VectorXd d_sd = -(H * G);
-        bool improved = true;
-        x = detail::line_search(x, dx, d_sd, norm0,
-            [&](const std::vector<double>& xnew) {
-                return inversive_distance_gradient(mesh, xnew, m);
-            }, &improved);
-        if (!improved) break;
-
-        res.iterations = iter + 1;
-    }
-
-    auto G_final = inversive_distance_gradient(mesh, x, m);
-    double inf_final = 0.0;
-    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
-    res.grad_inf_norm = inf_final;
-    res.x             = x;
-    return res;
 }
 
 } // namespace conformallab

code/include/serialization.hpp

@@ -93,19 +93,44 @@ inline std::vector<double> load_result_json(
 {
     using json = nlohmann::json;
     std::ifstream ifs(path);
-    if (!ifs) throw std::runtime_error("Cannot open: " + path);
+    if (!ifs) throw std::runtime_error("conformallab: cannot open: " + path);
-    json j; ifs >> j;
 
+    // V1: wrap parse + field extraction so nlohmann exceptions (parse_error,
+    // type_error, out_of_range) surface as std::runtime_error with the path.
+    json j;
+    try {
+        ifs >> j;
+    } catch (const json::exception& e) {
+        throw std::runtime_error(
+            "conformallab: malformed JSON in " + path + ": " + e.what());
+    }
+
+    try {
         if (geom && j.contains("geometry"))
             *geom = j["geometry"].get<std::string>();
 
+        // V4: validate required top-level key before accessing it.
+        if (!j.contains("dof_vector"))
+            throw std::runtime_error(
+                "conformallab: result JSON missing field 'dof_vector' in " + path);
         std::vector<double> x = j.at("dof_vector").get<std::vector<double>>();
 
         if (res) {
+            // V4: validate nested solver keys before accessing.
+            if (!j.contains("solver"))
+                throw std::runtime_error(
+                    "conformallab: result JSON missing field 'solver' in " + path);
+            const auto& s = j.at("solver");
+            for (const char* key : {"converged", "iterations", "grad_inf_norm"}) {
+                if (!s.contains(key))
+                    throw std::runtime_error(
+                        std::string("conformallab: result JSON missing field 'solver.")
+                        + key + "' in " + path);
+            }
             res->x             = x;
-        res->converged     = j["solver"]["converged"].get<bool>();
+            res->converged     = s.at("converged").get<bool>();
-        res->iterations    = j["solver"]["iterations"].get<int>();
+            res->iterations    = s.at("iterations").get<int>();
-        res->grad_inf_norm = j["solver"]["grad_inf_norm"].get<double>();
+            res->grad_inf_norm = s.at("grad_inf_norm").get<double>();
         }
 
         if (layout2d && j.contains("layout") && j["layout"]["dim"] == 2) {
@@ -118,6 +143,10 @@ inline std::vector<double> load_result_json(
         }
 
         return x;
+    } catch (const json::exception& e) {
+        throw std::runtime_error(
+            "conformallab: malformed JSON in " + path + ": " + e.what());
+    }
 }
 
 // ════════════════════════════════════════════════════════════════════════════
@@ -256,8 +285,22 @@ inline std::vector<double> load_result_xml(
         else if (line.find("<Solver") != std::string::npos) {
             if (res) {
                 res->converged = (detail_xml::xml_get_attr(line, "converged") == "true");
-                res->iterations    = std::stoi(detail_xml::xml_get_attr(line, "iterations"));
+                // V2: stoi/stod throw std::invalid_argument on empty or non-numeric
-                res->grad_inf_norm = std::stod(detail_xml::xml_get_attr(line, "grad_inf_norm"));
+                // attribute values; wrap and rethrow as runtime_error with context.
+                try {
+                    auto iter_str = detail_xml::xml_get_attr(line, "iterations");
+                    auto grad_str = detail_xml::xml_get_attr(line, "grad_inf_norm");
+                    if (iter_str.empty())
+                        throw std::runtime_error("missing attribute 'iterations'");
+                    if (grad_str.empty())
+                        throw std::runtime_error("missing attribute 'grad_inf_norm'");
+                    res->iterations    = std::stoi(iter_str);
+                    res->grad_inf_norm = std::stod(grad_str);
+                } catch (const std::exception& e) {
+                    throw std::runtime_error(
+                        "conformallab: malformed XML Solver element in "
+                        + path + ": " + e.what());
+                }
             }
         }
         // DOF vector

code/include/spherical_functional.hpp

@@ -35,6 +35,7 @@
 // └──────────────────────────────────────────────────────────────────────────┘
 
 #include "conformal_mesh.hpp"
+#include "constants.hpp"
 #include "spherical_geometry.hpp"
 #include "gauss_legendre.hpp"
 #include <CGAL/boost/graph/iterator.h>
@@ -176,7 +177,7 @@ inline void compute_spherical_lambda0_from_mesh(ConformalMesh& mesh, SphericalMa
         if (half_sin > 1e-15)
             m.lambda0[e] = 2.0 * std::log(half_sin);
         else
-            m.lambda0[e] = -30.0;  // very short edge: essentially 0
+            m.lambda0[e] = LOG_EDGE_LENGTH_FLOOR;  // very short edge: essentially 0
     }
 }
 

code/include/spherical_geometry.hpp

@@ -31,7 +31,7 @@ constexpr double PI_SPHER = PI;
 inline double spherical_l(double lambda)
 {
     double half = std::exp(lambda * 0.5);
-    if (half >= 1.0) half = 1.0 - 1e-15;
+    if (half >= 1.0) half = ASIN_DOMAIN_GUARD;
     if (half <= 0.0) return 0.0;
     return 2.0 * std::asin(half);
 }

code/tests/cgal/test_cgal_traits_mvp.cpp

@@ -111,6 +111,11 @@ TEST(CGALDiscreteConformalMap, SingleTriangleConverges)
     EXPECT_GE(result.iterations, 0);
     EXPECT_EQ(result.u_per_vertex.size(), num_vertices(mesh));
 
+    // I1/N7 propagated through the CGAL layer: status enum + diagnostics.
+    EXPECT_EQ(result.status, CGAL::Newton_status::Converged);
+    EXPECT_FALSE(result.sparse_qr_fallback_used);   // pinned vertex → no gauge mode
+    EXPECT_GE(result.min_ldlt_pivot, 0.0);
+
     // With the default flat-disc target curvature and the first vertex pinned,
     // the natural-theta equilibrium is at u = 0 — Newton should accept x0=0.
     for (double u : result.u_per_vertex)

code/tests/cgal/test_euclidean_hessian.cpp

@@ -61,6 +61,47 @@ TEST(EuclideanHessian, CotWeights_RightIsoscelesTriangle)
     EXPECT_NEAR(cw.cot3, 1.0, 1e-12);   // 45° at v3
 }
 
+// ════════════════════════════════════════════════════════════════════════════
+// N5: cotangent weights on a SLIVER triangle match a high-precision reference.
+//
+// The area is now computed by Kahan's stable side-length formula instead of
+// the naive 2·√(t12·t23·t31·l123).  On a thin (sliver) triangle the naive
+// product of near-equal-length differences loses precision, biasing every
+// cotangent weight.  Here we cross-check against an independent law-of-cosines
+// reference in long double (80-bit on x86 CI — a genuine high-precision oracle;
+// equal to double on ARM64, where it still serves as a regression cross-check).
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(EuclideanHessian, CotWeights_SliverMatchesHighPrecisionReference)
+{
+    // Thin isosceles: short base l12, two near-unit legs (apex angle ≈ 0.01 rad).
+    const double l12 = 0.01, l23 = 1.0, l31 = 1.0;
+    auto cw = euclidean_cot_weights(l12, l23, l31);
+    ASSERT_TRUE(cw.valid);
+
+    // cot at the vertex opposite side `opp`, with adjacent sides s1, s2:
+    //   cos = (s1² + s2² − opp²) / (2·s1·s2),  sin = √((1−cos)(1+cos)).
+    auto cot_ref = [](long double opp, long double s1, long double s2) {
+        long double cosA = (s1 * s1 + s2 * s2 - opp * opp) / (2.0L * s1 * s2);
+        long double sinA = std::sqrt((1.0L - cosA) * (1.0L + cosA));
+        return static_cast<double>(cosA / sinA);
+    };
+    const double c1 = cot_ref(l23, l12, l31);  // v1 opposite l23
+    const double c2 = cot_ref(l31, l12, l23);  // v2 opposite l31
+    const double c3 = cot_ref(l12, l23, l31);  // v3 opposite l12 (needle angle)
+
+    auto rel = [](double a, double b) {
+        return std::abs(a - b) / std::max(1.0, std::abs(b));
+    };
+    EXPECT_LT(rel(cw.cot1, c1), 1e-9);
+    EXPECT_LT(rel(cw.cot2, c2), 1e-9);
+    EXPECT_LT(rel(cw.cot3, c3), 1e-9);
+
+    EXPECT_TRUE(std::isfinite(cw.cot1) &&
+                std::isfinite(cw.cot2) &&
+                std::isfinite(cw.cot3));
+}
+
 // ════════════════════════════════════════════════════════════════════════════
 // Hessian is symmetric: H[i,j] == H[j,i]
 // ════════════════════════════════════════════════════════════════════════════

code/tests/cgal/test_hyper_ideal_hessian.cpp

@@ -338,3 +338,81 @@ TEST(HyperIdealHessian, BlockFD_FasterThanFullFD)
     EXPECT_GE(ms_full, 3 * ms_block)
         << "Block-FD should be at least 3× faster than full-FD on this mesh";
 }
+
+// ════════════════════════════════════════════════════════════════════════════
+// N3 — Scale-floor clamp modes (HardJava vs SmoothBarrier)
+//
+// The vertex-scale floor can be applied two ways (HyperIdealScaleClamp):
+//   • HardJava (default): b<0 → floor.  Faithful to Java, but only C⁰ — and
+//     in fact value-discontinuous at b=0 (jumps from `floor` to 0).
+//   • SmoothBarrier: floor + softplus_β(b−floor).  C¹ everywhere, ≈ b away
+//     from the floor, and keeps b ≥ floor > 0 smoothly.
+// These tests pin the contract of both modes and the default-preservation.
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(HyperIdealClamp, SmoothBarrierIsContinuousAndC1AcrossZero)
+{
+    const auto Hard   = HyperIdealScaleClamp::HardJava;
+    const auto Smooth = HyperIdealScaleClamp::SmoothBarrier;
+    auto f = [](double b, HyperIdealScaleClamp mode) {
+        return clamp_hyper_ideal_scale(b, /*variable=*/true, mode);
+    };
+
+    const double eps = 1e-6;
+
+    // HardJava jumps in VALUE at b=0 (floor on the left, 0 on the right).
+    EXPECT_GT(std::abs(f(-eps, Hard) - f(+eps, Hard)), 0.5 * HYPER_IDEAL_SCALE_FLOOR);
+
+    // SmoothBarrier is continuous in value across b=0 …
+    EXPECT_NEAR(f(-eps, Smooth), f(+eps, Smooth), 1e-4);
+
+    // … and C¹: the one-sided slopes across b=0 agree (the HardJava slopes,
+    // 0 on the left and 1 on the right, would not).
+    auto slope = [&](double b0, HyperIdealScaleClamp mode) {
+        return (f(b0 + eps, mode) - f(b0 - eps, mode)) / (2.0 * eps);
+    };
+    const double sL = slope(-1e-3, Smooth);
+    const double sR = slope(+1e-3, Smooth);
+    EXPECT_NEAR(sL, sR, 5e-2) << "SmoothBarrier derivative should be continuous";
+}
+
+TEST(HyperIdealClamp, SmoothBarrierStaysAboveFloorAndApproachesIdentity)
+{
+    const auto Smooth = HyperIdealScaleClamp::SmoothBarrier;
+    auto f = [&](double b) { return clamp_hyper_ideal_scale(b, true, Smooth); };
+
+    // At or above the floor for any input (mathematically > floor; for very
+    // negative b the softplus underflows to 0 in double, giving exactly floor).
+    EXPECT_GE(f(-1e6), HYPER_IDEAL_SCALE_FLOOR);
+    EXPECT_GE(f(-1.0), HYPER_IDEAL_SCALE_FLOOR);
+    EXPECT_GT(f(0.0),  0.0);
+    EXPECT_LT(f(-1e6) - HYPER_IDEAL_SCALE_FLOOR, 1e-9);   // converges down to floor
+
+    // ≈ identity well above the floor (the normal operating regime).
+    EXPECT_NEAR(f(1.0),  1.0,  1e-9);
+    EXPECT_NEAR(f(5.0),  5.0,  1e-12);
+
+    // Pinned DOFs are never clamped, regardless of mode.
+    EXPECT_EQ(clamp_hyper_ideal_scale(-3.0, /*variable=*/false, Smooth), -3.0);
+}
+
+// Away from the feasibility boundary (all b = 1 ≫ floor), the two modes must
+// produce the same Hessian — SmoothBarrier only differs near b ≈ floor, so the
+// default-mode parity results are not perturbed in the normal regime.
+TEST(HyperIdealClamp, ModesAgreeAwayFromBoundary)
+{
+    auto mesh = make_tetrahedron();
+    auto m    = setup_hyper_ideal_maps(mesh);
+    const int n = assign_all_dof_indices(mesh, m);
+    auto x = natural_x(mesh, m);   // b = 1, a = 0.5 — all well above the floor
+
+    auto H_hard = hyper_ideal_hessian_block_fd_sym(
+        mesh, x, m, 1e-5, HyperIdealScaleClamp::HardJava);
+    auto H_soft = hyper_ideal_hessian_block_fd_sym(
+        mesh, x, m, 1e-5, HyperIdealScaleClamp::SmoothBarrier);
+
+    Eigen::MatrixXd Dh(H_hard), Ds(H_soft);
+    EXPECT_LT((Dh - Ds).cwiseAbs().maxCoeff(), 1e-9)
+        << "clamp modes must agree when every scale is well above the floor";
+    (void)n;
+}

code/tests/cgal/test_lawson_hyperideal.cpp

@@ -194,6 +194,54 @@ TEST(LawsonHyperIdeal, ConvergenceGoldenVector_JavaXVal)
     }
 }
 
+// ════════════════════════════════════════════════════════════════════════════
+// N3 — SmoothBarrier clamp mode converges to the same golden solution
+//
+// The C¹ smooth-barrier scale floor is an opt-in alternative to the Java hard
+// clamp.  On this well-posed problem the scales stay well above the floor
+// throughout the solve, so the barrier is ≈ identity and the converged vector
+// must match the same golden classes as the HardJava run above — demonstrating
+// the smooth mode is a drop-in that does not move the solution when the clamp
+// region is never entered.
+// ════════════════════════════════════════════════════════════════════════════
+TEST(LawsonHyperIdeal, SmoothBarrierConvergesToGoldenVector)
+{
+    std::set<Edge_index> original;
+    ConformalMesh m = make_lawson_square_tiled(&original);
+
+    HyperIdealMaps maps = setup_hyper_ideal_maps(m);
+    const int n = assign_all_dof_indices(m, maps);
+    ASSERT_EQ(n, 4 + 18);
+
+    for (auto e : m.edges())
+        if (original.count(e)) maps.theta_e[e] = PI / 2.0;
+
+    std::vector<double> x0(static_cast<std::size_t>(n), 1.0);
+    auto res = newton_hyper_ideal(m, x0, maps, /*tol=*/1e-10, /*max_iter=*/200,
+                                  /*hess_eps=*/1e-5,
+                                  HyperIdealScaleClamp::SmoothBarrier);
+    ASSERT_TRUE(res.converged)
+        << "SmoothBarrier HyperIdeal Newton did not converge; ||G||="
+        << res.grad_inf_norm;
+
+    constexpr double b_gold      = 1.1462158341786262;
+    constexpr double a_orig_gold = 1.7627471737467797;
+    constexpr double a_aux_gold  = 2.633915794495759;
+    const double tol = 1e-5;
+
+    for (auto v : m.vertices()) {
+        const int iv = maps.v_idx[v];
+        ASSERT_GE(iv, 0);
+        EXPECT_NEAR(res.x[static_cast<std::size_t>(iv)], b_gold, tol);
+    }
+    for (auto e : m.edges()) {
+        const int ie = maps.e_idx[e];
+        ASSERT_GE(ie, 0);
+        const double gold = original.count(e) ? a_orig_gold : a_aux_gold;
+        EXPECT_NEAR(res.x[static_cast<std::size_t>(ie)], gold, tol);
+    }
+}
+
 // ════════════════════════════════════════════════════════════════════════════
 // Branch-points variant — Java HyperIdealConvergenceTest...WithBranchPoints
 //

code/tests/cgal/test_layout.cpp

@@ -24,6 +24,7 @@
 #include "serialization.hpp"
 #include <gtest/gtest.h>
 #include <cmath>
+#include <fstream>
 #include <vector>
 #include <filesystem>
 
@@ -370,3 +371,67 @@ TEST(Serialization, XML_RoundTrip)
 
     std::filesystem::remove(path);
 }
+
+// ════════════════════════════════════════════════════════════════════════════
+// I2: serialization throws on malformed / schema-violating input
+//
+// These tests cover the throw paths in load_result_json / load_result_xml
+// that were previously untested (I2 in the test-coverage audit).
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(Serialization, LoadResultJson_ThrowsOnMissingFile)
+{
+    EXPECT_THROW(load_result_json("/tmp/does_not_exist_conflab.json"),
+                 std::runtime_error);
+}
+
+TEST(Serialization, LoadResultJson_ThrowsOnMalformedJson)
+{
+    const std::string path = "/tmp/conflab_bad.json";
+    { std::ofstream ofs(path); ofs << "{not valid json!!!"; }
+    EXPECT_THROW(load_result_json(path), std::runtime_error);
+    std::filesystem::remove(path);
+}
+
+TEST(Serialization, LoadResultJson_ThrowsOnMissingDofVector)
+{
+    // V4: a valid JSON object but without the required "dof_vector" key.
+    const std::string path = "/tmp/conflab_nodof.json";
+    { std::ofstream ofs(path); ofs << R"({"geometry":"euclidean"})"; }
+    EXPECT_THROW(load_result_json(path), std::runtime_error);
+    std::filesystem::remove(path);
+}
+
+TEST(Serialization, LoadResultJson_ThrowsOnMissingSolverField)
+{
+    // V4: has dof_vector but solver block is absent when res != nullptr.
+    const std::string path = "/tmp/conflab_nosolver.json";
+    { std::ofstream ofs(path); ofs << R"({"dof_vector":[0.1,0.2]})"; }
+    NewtonResult res;
+    EXPECT_THROW(load_result_json(path, &res), std::runtime_error);
+    std::filesystem::remove(path);
+}
+
+TEST(Serialization, LoadResultXml_ThrowsOnMissingFile)
+{
+    EXPECT_THROW(load_result_xml("/tmp/does_not_exist_conflab.xml"),
+                 std::runtime_error);
+}
+
+TEST(Serialization, LoadResultXml_ThrowsOnMalformedSolverAttribute)
+{
+    // V2: non-numeric attribute causes stoi/stod error that must surface
+    // as runtime_error, not std::invalid_argument.
+    const std::string path = "/tmp/conflab_badxml.xml";
+    {
+        std::ofstream ofs(path);
+        ofs << R"(<?xml version="1.0"?>
+<ConformalResult geometry="euclidean" vertices="4" faces="4">
+  <Solver converged="true" iterations="not_a_number" grad_inf_norm="1e-9"/>
+  <DOFVector n="1">0.0</DOFVector>
+</ConformalResult>)";
+    }
+    NewtonResult res;
+    EXPECT_THROW(load_result_xml(path, &res), std::runtime_error);
+    std::filesystem::remove(path);
+}

code/tests/cgal/test_mesh_io.cpp

@@ -18,6 +18,7 @@
 #include <cmath>
 #include <cstdio>
 #include <filesystem>
+#include <fstream>
 #include <stdexcept>
 
 using namespace conformallab;
@@ -111,6 +112,30 @@ TEST(MeshIO, LoadMeshThrowsOnMissingFile)
     EXPECT_THROW(load_mesh("/tmp/does_not_exist_conflab.off"), std::runtime_error);
 }
 
+// ════════════════════════════════════════════════════════════════════════════
+// I3: load_mesh throws on non-triangulated mesh
+//
+// The quad-strip mesh has 4-vertex faces; load_mesh must reject it at the
+// I/O boundary rather than letting the quad faces flow silently into the
+// triangle-only functionals.
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(MeshIO, LoadMeshThrowsOnNonTriangulatedMesh)
+{
+    // Write a minimal OFF quad-face mesh (2 quads, not triangles).
+    auto path = tmp_path("quad.off");
+    {
+        std::ofstream ofs(path);
+        ofs << "OFF\n6 2 0\n"
+            << "0 0 0\n1 0 0\n1 1 0\n0 1 0\n"
+            << "2 0 0\n2 1 0\n"
+            << "4 0 1 2 3\n"
+            << "4 1 4 5 2\n";
+    }
+    EXPECT_THROW(load_mesh(path), std::runtime_error);
+    rm(path);
+}
+
 // ════════════════════════════════════════════════════════════════════════════
 // Vertex positions survive a round-trip (OFF)
 // ════════════════════════════════════════════════════════════════════════════
@@ -171,3 +196,42 @@ TEST(MeshIO, SaveLoadConvenienceWrappers)
 
     rm(path);
 }
+
+// ════════════════════════════════════════════════════════════════════════════
+// V3: load_mesh throws on non-finite vertex coordinate (NaN / Inf)
+//
+// A mesh file containing NaN or Inf vertex coordinates must be rejected at
+// the I/O boundary — before the coordinates can silently poison the solver.
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(MeshIO, LoadMeshThrowsOnNaNCoordinate)
+{
+    // Write a minimal OFF triangle with a NaN in the first vertex.
+    auto path = tmp_path("nan.off");
+    {
+        std::ofstream ofs(path);
+        ofs << "OFF\n3 1 0\n"
+            << "nan 0.0 0.0\n"
+            << "1.0 0.0 0.0\n"
+            << "0.0 1.0 0.0\n"
+            << "3 0 1 2\n";
+    }
+    EXPECT_THROW(load_mesh(path), std::runtime_error);
+    rm(path);
+}
+
+TEST(MeshIO, LoadMeshThrowsOnInfCoordinate)
+{
+    // Write a minimal OFF triangle with Inf in the y-coordinate.
+    auto path = tmp_path("inf.off");
+    {
+        std::ofstream ofs(path);
+        ofs << "OFF\n3 1 0\n"
+            << "0.0 0.0 0.0\n"
+            << "1.0 inf 0.0\n"
+            << "0.0 1.0 0.0\n"
+            << "3 0 1 2\n";
+    }
+    EXPECT_THROW(load_mesh(path), std::runtime_error);
+    rm(path);
+}

code/tests/cgal/test_newton_phase9a.cpp

@@ -16,6 +16,7 @@
 #include "newton_solver.hpp"
 #include "cp_euclidean_functional.hpp"
 #include "inversive_distance_functional.hpp"
+#include "inversive_distance_hessian.hpp"
 #include "mesh_builder.hpp"
 #include "conformal_mesh.hpp"
 
@@ -232,6 +233,49 @@ TEST(NewtonPhase9a, InversiveDistance_PerturbedTetrahedron_Converges)
     EXPECT_LT(res.grad_inf_norm, 1e-8);
 }
 
+// ════════════════════════════════════════════════════════════════════════════
+// 6b. Inversive-Distance block-FD Hessian == full-FD Hessian (B1 port)
+//
+// The solver was switched from the O(n·F) full-FD Hessian to the O(F)
+// per-face block-FD Hessian.  The two must agree to FD rounding on any
+// configuration — including a perturbed (non-equilibrium) one, where the
+// off-diagonal coupling is non-trivial.  This is the locality-lemma
+// cross-validation that licenses the faster path.
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(NewtonPhase9a, InversiveDistance_BlockFDHessianMatchesFullFD)
+{
+    auto mesh = make_tetrahedron();
+    auto m    = setup_inversive_distance_maps(mesh);
+    compute_inversive_distance_init_from_mesh(mesh, m);
+
+    // Pin one vertex; index the rest.
+    auto vit = mesh.vertices().begin();
+    m.v_idx[*vit++] = -1;
+    int n = 0;
+    for (; vit != mesh.vertices().end(); ++vit) m.v_idx[*vit] = n++;
+
+    // Evaluate the two Hessians away from equilibrium so off-diagonals are live.
+    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
+    for (int i = 0; i < n; ++i) x[static_cast<std::size_t>(i)] = 0.07 * (i + 1);
+
+    auto H_full  = inversive_distance_hessian_sym(mesh, x, m);
+    auto H_block = inversive_distance_hessian_block_fd_sym(mesh, x, m);
+
+    ASSERT_EQ(H_full.rows(),  H_block.rows());
+    ASSERT_EQ(H_full.cols(),  H_block.cols());
+
+    Eigen::MatrixXd Df = Eigen::MatrixXd(H_full);
+    Eigen::MatrixXd Db = Eigen::MatrixXd(H_block);
+    double max_abs_diff = (Df - Db).cwiseAbs().maxCoeff();
+    EXPECT_LT(max_abs_diff, 1e-7)
+        << "block-FD and full-FD Hessians must agree to FD rounding;\n"
+        << "max |Δ| = " << max_abs_diff;
+
+    // Sanity: the matrices are non-trivial (not both accidentally zero).
+    EXPECT_GT(Df.cwiseAbs().maxCoeff(), 1e-3);
+}
+
 // ════════════════════════════════════════════════════════════════════════════
 // 7. CP-Euclidean Newton — uses analytic Hessian (NOT FD)
 //

code/tests/cgal/test_newton_solver.cpp

@@ -500,4 +500,240 @@ TEST(SparseQRFallback, Euclidean_ClosedMeshNoPinConverges)
         << "Euclidean Newton on closed tetrahedron (no pin) must converge via SparseQR; "
            "grad_inf_norm = " << res.grad_inf_norm;
     EXPECT_LT(res.grad_inf_norm, 1e-8);
+    EXPECT_EQ(res.status, NewtonStatus::Converged);
+    // N7: the closed mesh with no pinned vertex has a gauge-singular Hessian.
+    // That near-singularity must be observable in the diagnostics — either the
+    // SparseQR fallback fired, or (when LDLT "succeeds" with a ~0 pivot, which is
+    // exactly the silent case N7 targets) the smallest LDLT pivot is tiny.
+    EXPECT_TRUE(res.sparse_qr_fallback_used || res.min_ldlt_pivot < 1e-6)
+        << "gauge singularity should surface via fallback flag or min_ldlt_pivot; "
+        << "fallback=" << res.sparse_qr_fallback_used
+        << " min_pivot=" << res.min_ldlt_pivot;
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// C1: solve_linear_system exposes the ok flag via the new out-parameter
+//
+// The C1 audit finding was that solve_linear_system silently dropped the
+// internal ok flag, making double-solver failure invisible to callers.
+// These tests verify that the new optional bool* ok parameter is correctly
+// wired: ok=true for successful solves, and the parameter is optional (null
+// is safe).
+//
+// Note: constructing a matrix where Eigen's SparseQR hard-fails is not
+// practical — SparseQR reports Eigen::Success even for rank-0 inputs,
+// returning a zero min-norm solution.  The observable contract for callers is
+// therefore: ok=true whenever a solver reports success; ok=false only when
+// both report failure (an edge case Eigen essentially never produces).  The
+// fix in C1 ensures that if that edge case ever fires it propagates to the
+// caller — previously it was silently swallowed.
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(SparseQRFallback, OkFlag_TrueOnFullRankSystem)
+{
+    // 3×3 diagonal PD matrix: LDLT succeeds → ok must be true.
+    Eigen::SparseMatrix<double> A(3, 3);
+    A.insert(0, 0) = 1.0;
+    A.insert(1, 1) = 2.0;
+    A.insert(2, 2) = 3.0;
+    A.makeCompressed();
+
+    Eigen::VectorXd rhs(3);
+    rhs << 1.0, 4.0, 9.0;
+
+    bool fallback = true;
+    bool ok       = false;
+    Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback, &ok);
+
+    EXPECT_TRUE(ok)      << "Full-rank system: ok must be true";
+    EXPECT_FALSE(fallback) << "Full-rank system: LDLT path, no fallback expected";
+    EXPECT_NEAR(x[0], 1.0, 1e-12);
+    EXPECT_NEAR(x[1], 2.0, 1e-12);
+    EXPECT_NEAR(x[2], 3.0, 1e-12);
+}
+
+TEST(SparseQRFallback, OkFlag_TrueOnSingularSystemViaFallback)
+{
+    // Singular matrix (zero row/col 1) — LDLT fails, SparseQR succeeds → ok=true.
+    Eigen::SparseMatrix<double> A(3, 3);
+    A.insert(0, 0) = 2.0;
+    A.insert(2, 2) = 3.0;
+    A.makeCompressed();
+
+    Eigen::VectorXd rhs(3);
+    rhs << 2.0, 0.0, 3.0;
+
+    bool fallback = false;
+    bool ok       = false;
+    Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback, &ok);
+
+    EXPECT_TRUE(ok)      << "Singular-but-consistent system: SparseQR succeeds → ok must be true";
+    EXPECT_TRUE(fallback) << "Singular system: fallback must have been triggered";
+}
+
+TEST(SparseQRFallback, OkFlag_NullPointerIsSafe)
+{
+    // Calling with ok=nullptr must not crash (backward-compatibility).
+    Eigen::SparseMatrix<double> A(2, 2);
+    A.insert(0, 0) = 1.0;
+    A.insert(1, 1) = 1.0;
+    A.makeCompressed();
+
+    Eigen::VectorXd rhs(2);
+    rhs << 3.0, 5.0;
+
+    EXPECT_NO_THROW({
+        Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, nullptr, nullptr);
+        EXPECT_NEAR(x[0], 3.0, 1e-12);
+        EXPECT_NEAR(x[1], 5.0, 1e-12);
+    });
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// H2 / B2 — newton_core reuses the line-search gradient
+//
+// All five solvers now share detail::newton_core.  B2: the gradient computed at
+// the accepted line-search point is carried into the next iteration's
+// convergence check instead of being recomputed at the loop top.  On a
+// unit-Hessian quadratic, exact Newton reaches the minimum in one step, so the
+// total gradient-eval count is exactly 2 (1 initial + 1 accepted trial); the
+// pre-B2 loop would have recomputed G at the top of iteration 1 → 3.
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(NewtonCore, ReusesLineSearchGradient_B2_Convex)
+{
+    const std::vector<double> xstar = {1.0, -2.0, 3.5};
+    int grad_calls = 0;
+
+    auto grad = [&](const std::vector<double>& x) {
+        ++grad_calls;
+        std::vector<double> g(x.size());
+        for (std::size_t i = 0; i < x.size(); ++i) g[i] = x[i] - xstar[i];
+        return g;                                   // G(x) = x − x*, H = +I
+    };
+    auto hess = [&](const std::vector<double>&) {
+        Eigen::SparseMatrix<double> H(3, 3);
+        for (int i = 0; i < 3; ++i) H.insert(i, i) = 1.0;
+        H.makeCompressed();
+        return H;
+    };
+
+    auto res = conformallab::detail::newton_core(
+        std::vector<double>{0.0, 0.0, 0.0}, grad, hess,
+        /*concave=*/false, /*tol=*/1e-9, /*max_iter=*/50);
+
+    ASSERT_TRUE(res.converged);
+    EXPECT_EQ(res.iterations, 1);
+    EXPECT_NEAR(res.x[0],  1.0, 1e-9);
+    EXPECT_NEAR(res.x[1], -2.0, 1e-9);
+    EXPECT_NEAR(res.x[2],  3.5, 1e-9);
+    EXPECT_EQ(grad_calls, 2)
+        << "B2: the loop-top gradient must be reused from the line search";
+}
+
+// Concave path (spherical-style): H is NSD, newton_core factors −H and solves
+// (−H)·Δx = G.  G(x) = x* − x with H = −I makes x* a maximiser; the core must
+// still converge in one step.
+TEST(NewtonCore, ConcavePathConverges)
+{
+    const std::vector<double> xstar = {0.5, -1.5, 2.0};
+    int grad_calls = 0;
+
+    auto grad = [&](const std::vector<double>& x) {
+        ++grad_calls;
+        std::vector<double> g(x.size());
+        for (std::size_t i = 0; i < x.size(); ++i) g[i] = xstar[i] - x[i];
+        return g;                                   // G(x) = x* − x, H = −I
+    };
+    auto hess = [&](const std::vector<double>&) {
+        Eigen::SparseMatrix<double> H(3, 3);
+        for (int i = 0; i < 3; ++i) H.insert(i, i) = -1.0;
+        H.makeCompressed();
+        return H;
+    };
+
+    auto res = conformallab::detail::newton_core(
+        std::vector<double>{0.0, 0.0, 0.0}, grad, hess,
+        /*concave=*/true, /*tol=*/1e-9, /*max_iter=*/50);
+
+    ASSERT_TRUE(res.converged);
+    EXPECT_EQ(res.iterations, 1);
+    EXPECT_NEAR(res.x[0],  0.5, 1e-9);
+    EXPECT_NEAR(res.x[1], -1.5, 1e-9);
+    EXPECT_NEAR(res.x[2],  2.0, 1e-9);
+    EXPECT_EQ(grad_calls, 2);
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// I1 / H1 / N7 — newton_core termination status, iteration count, diagnostics
+// ════════════════════════════════════════════════════════════════════════════
+
+namespace {
+// 1×1 sparse identity / scalar helpers for the synthetic newton_core tests.
+inline Eigen::SparseMatrix<double> diag_sparse(std::initializer_list<double> d)
+{
+    const int n = static_cast<int>(d.size());
+    Eigen::SparseMatrix<double> H(n, n);
+    int i = 0;
+    for (double v : d) { H.insert(i, i) = v; ++i; }
+    H.makeCompressed();
+    return H;
+}
+} // namespace
+
+// I1: a converging solve reports status == Converged and the N7 pivot is the
+// (well-conditioned) unit diagonal; no SparseQR fallback.
+TEST(NewtonCore, Status_Converged_AndDiagnostics_N7)
+{
+    const std::vector<double> xstar = {2.0, -1.0};
+    auto grad = [&](const std::vector<double>& x) {
+        return std::vector<double>{ x[0] - xstar[0], x[1] - xstar[1] };
+    };
+    auto hess = [&](const std::vector<double>&) { return diag_sparse({1.0, 1.0}); };
+
+    auto res = conformallab::detail::newton_core(
+        std::vector<double>{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50);
+
+    EXPECT_TRUE(res.converged);
+    EXPECT_EQ(res.status, conformallab::NewtonStatus::Converged);
+    EXPECT_EQ(res.iterations, 1);
+    EXPECT_FALSE(res.sparse_qr_fallback_used);          // N7
+    EXPECT_NEAR(res.min_ldlt_pivot, 1.0, 1e-12);        // N7: D = I
+}
+
+// I1: a slow (linearly-convergent cubic) problem that does not reach tol within
+// max_iter reports MaxIterations, with iterations == max_iter (H1).
+TEST(NewtonCore, Status_MaxIterations)
+{
+    // G(x) = x³, H = 3x²  →  Newton map x ← (2/3)x (linear convergence).
+    auto grad = [&](const std::vector<double>& x) {
+        return std::vector<double>{ x[0]*x[0]*x[0] };
+    };
+    auto hess = [&](const std::vector<double>& x) {
+        return diag_sparse({ 3.0 * x[0] * x[0] });
+    };
+
+    auto res = conformallab::detail::newton_core(
+        std::vector<double>{1.0}, grad, hess, /*concave=*/false, /*tol=*/1e-8, /*max_iter=*/3);
+
+    EXPECT_FALSE(res.converged);
+    EXPECT_EQ(res.status, conformallab::NewtonStatus::MaxIterations);
+    EXPECT_EQ(res.iterations, 3);
+}
+
+// I1/H1: a constant non-zero gradient cannot be reduced by any step, so the
+// line search stalls on the first iteration → LineSearchStalled, iterations == 0.
+TEST(NewtonCore, Status_LineSearchStalled)
+{
+    auto grad = [&](const std::vector<double>&) {
+        return std::vector<double>{ 1.0, 1.0 };       // constant, independent of x
+    };
+    auto hess = [&](const std::vector<double>&) { return diag_sparse({1.0, 1.0}); };
+
+    auto res = conformallab::detail::newton_core(
+        std::vector<double>{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50);
+
+    EXPECT_FALSE(res.converged);
+    EXPECT_EQ(res.status, conformallab::NewtonStatus::LineSearchStalled);
+    EXPECT_EQ(res.iterations, 0);                       // H1: no step completed
 }

code/tests/cgal/test_spherical_hessian.cpp

@@ -206,3 +206,31 @@ TEST(SphericalHessian, FDCheck_MixedPinnedVertices)
     EXPECT_TRUE(hessian_check_spherical(mesh, x, maps))
         << "FD Hessian check failed for mixed pinned/variable vertices";
 }
+
+// ════════════════════════════════════════════════════════════════════════════
+// I4: spherical_hessian throws on edge DOFs
+//
+// The spherical Hessian only implements the vertex-block cotangent Laplacian;
+// if any edge has a free DOF index (e_idx >= 0) it must fail loudly rather
+// than returning a silently rank-deficient matrix.
+// ════════════════════════════════════════════════════════════════════════════
+
+TEST(SphericalHessian, ThrowsOnEdgeDOF)
+{
+    // Build a tetrahedron and assign one edge as a free DOF.
+    auto mesh = make_tetrahedron();
+    auto maps = setup_spherical_maps(mesh);
+    compute_lambda0_from_mesh(mesh, maps);  // spherical variant (A2 rename pending)
+
+    int idx = 0;
+    for (auto v : mesh.vertices())
+        maps.v_idx[v] = idx++;
+    const int n_v = idx;
+
+    // Give one edge a free DOF index to trigger the guard.
+    auto e0 = *mesh.edges().begin();
+    maps.e_idx[e0] = n_v;  // first free edge DOF
+
+    std::vector<double> x(static_cast<std::size_t>(n_v + 1), 0.0);
+    EXPECT_THROW(spherical_hessian(mesh, x, maps), std::logic_error);
+}

doc/math/geometry-modes.md

@@ -80,8 +80,27 @@ lᵢⱼ(ζ)                             — edge length from ζ value
 βᵢ(l₁, l₂, l₃)                    — vertex angle sum contribution
 ```
 
-**Hessian:** currently a symmetric finite-difference approximation (see `hyper_ideal_hessian.hpp`).
+**Hessian:** a symmetric finite-difference approximation (see `hyper_ideal_hessian.hpp`).
-The analytic Hessian through the `(bᵢ, aₑ) → lᵢⱼ → ζ → αᵢⱼ/βᵢ` chain is Phase 9b.
+The solver uses the **per-face block-FD** variant by default (`O(F)` face evaluations,
+~33–1166× faster than full-FD); the full-FD baseline is retained for cross-validation.
+The analytic Hessian through the `(bᵢ, aₑ) → lᵢⱼ → ζ → αᵢⱼ/βᵢ` chain is Phase 9b-analytic.
+
+**Scale-floor clamp mode (`HyperIdealScaleClamp`).** The vertex scale `bᵢ` must stay
+positive for the geometry to be valid. When a Newton iterate drives `bᵢ` toward or below
+zero, one of two floors is applied — selectable per solver call (last argument of
+`newton_hyper_ideal`, default `HardJava`):
+
+| Mode | Behaviour | When to use |
+|---|---|---|
+| `HardJava` *(default)* | `bᵢ < 0 → HYPER_IDEAL_SCALE_FLOOR` (= 0.01). Bit-for-bit faithful to `HyperIdealFunctional.java`, so the Java golden-oracle parity tests hold **exactly**. It is only C⁰ (in fact value-discontinuous at `bᵢ = 0`), so a step crossing the feasibility boundary sees a kink that can stall Newton. | Default — preserves parity; correct whenever scales stay positive (the normal case). |
+| `SmoothBarrier` | `bᵢ ↦ floor + softplusβ(bᵢ − floor)` with `β = HYPER_IDEAL_SCALE_SHARPNESS` (= 100). C¹ everywhere, ≈ `bᵢ` to machine precision once `β·(bᵢ−floor) ≳ 35`, and decays smoothly to `floor` as `bᵢ → −∞`. Removes the N3 kink, but perturbs values near the boundary and therefore **deviates from the Java oracle** there. | Opt in when a solve is expected to cross the feasibility boundary and convergence robustness matters more than strict Java parity. |
+
+Both modes share `HYPER_IDEAL_SCALE_FLOOR`; the edge clamp (`aₑ < 0 → 0`) is unchanged.
+Away from the floor (the normal `bᵢ ≫ 0.01` regime) the two modes are numerically
+identical — the `LawsonHyperIdeal.SmoothBarrierConvergesToGoldenVector` test confirms the
+smooth mode reaches the same golden solution as `HardJava`. Implemented in
+`clamp_hyper_ideal_scale` (`hyper_ideal_functional.hpp`); origin: numerical-stability
+audit finding **N3**.
 
 **Holonomy:** Möbius isometries Tᵢ ∈ SU(1,1) (orientation-preserving isometries of the Poincaré disk).
 `MobiusMap` in `layout.hpp`: T(z) = (az+b)/(cz+d).

doc/math/references.md

@@ -29,12 +29,14 @@ Java reference implementation: [github.com/varylab/conformallab](https://github.
 | Reference | Used in |
 |---|---|
 | ✅ **Springborn** — *Ideal Hyperbolic Polyhedra and Discrete Uniformization*, Discrete & Computational Geometry **64** (2020), pp. 63–108. DOI: [10.1007/s00454-019-00132-8](https://doi.org/10.1007/s00454-019-00132-8) | `hyper_ideal_geometry.hpp` — ζ₁₃/ζ₁₄/ζ₁₅ functions; `hyper_ideal_functional.hpp` |
+| ✅ **Kolpakov, Mednykh** — *A Formula for the Volume of a Hyperbolic Tetrahedron*, arXiv: [math/0603097](https://arxiv.org/abs/math/0603097) (2006) | Tetrahedron volume with one ideal vertex: `calculateTetrahedronVolumeWithIdealVertexAtGamma` in `hyper_ideal_utility.hpp` (Phase 9b analytic Hessian) |
+| ✅ **Meyerhoff, Ushijima** — *A Note on the Dirichlet Domain*, in: The Epstein Birthday Schrift (2006) | Tetrahedron volume with three ideal vertices: `calculateTetrahedronVolumeFullyIdeal` in `hyper_ideal_utility.hpp` |
 | **Pinkall, Polthier** — *Computing Discrete Minimal Surfaces and Their Conjugates*, Experimental Mathematics (1993) | `euclidean_hessian.hpp` — cotangent Laplacian |
 | **Bobenko, Springborn** — *Variational Principles for Circle Patterns and Koebe's Theorem*, Transactions AMS (2004) | Variational angle-sum framework underlying all three functionals |
 | **Luo** — *Combinatorial Yamabe Flow on Surfaces*, Communications in Contemporary Mathematics (2004) | Inversive-distance functional — **new research** in Phase 9a.2 (no Java original; implemented from this paper + Glickenstein 2011 + Bowers-Stephenson 2004) |
 | **Bowers, Stephenson** — *Uniformizing dessins and Belyĭ maps via circle packing*, Memoirs of the AMS 170(805) (2004) | Introduces **inversive-distance circle packings** (used in Phase 9a.2). *Hinweis:* die zur Initialisierung benutzte Formel I_ij = (ℓ²−r_i²−r_j²)/(2 r_i r_j) ist die **klassische** inversive Distanz (vgl. Glickenstein §5.2: ℓ²=r_i²+r_j²+2r_ir_jη), nicht eine eigene „Bowers-Stephenson-Identität" — B–S liefern die Packungstheorie, nicht diese Formel. |
 | **Glickenstein** — *Discrete conformal variations and scalar curvature on piecewise flat two- and three-dimensional manifolds*, J. Differential Geometry **87**(2) (2011), pp. 201–238 | Analytic Hessian of the inversive-distance functional. ⚠️ *Korrektur:* die Arbeit nummeriert Gleichungen **nicht** im Format „(4.6)" — der Verweis ist durch die **§5.2**-Parametrisierung ℓ²_ij = r²_i + r²_j + 2 r_i r_j η_ij zu ersetzen. Cross-correspondence: η_ij ist die inversive Distanz und entspricht dem Kosinus des **Supplements** des Schnittwinkels (Schnitt bei arccos(−η_ij)) — also I_ij = cos θ_e **nur bis aufs Vorzeichen/Supplement**, nicht wörtlich. |
-| **Bobenko, Pinkall, Springborn** — *Discrete conformal maps and ideal hyperbolic polyhedra*, Geometry & Topology **19**(4) (2015), pp. 2155–2215. arXiv: [1005.2698](https://arxiv.org/abs/1005.2698) | Face-based circle-packing functional (`CPEuclideanFunctional.java` → `cp_euclidean_functional.hpp`, Phase 9a.1) |
+| ✅ **Bobenko, Pinkall, Springborn** — *Discrete conformal maps and ideal hyperbolic polyhedra*, Geometry & Topology **19**(4) (2015), pp. 2155–2215. arXiv: [1005.2698](https://arxiv.org/abs/1005.2698) (first posted 2010) | Face-based circle-packing functional (`CPEuclideanFunctional.java` → `cp_euclidean_functional.hpp`, Phase 9a.1) |
 | **Schläfli** — *On the multiple integral ∫dx dy …*, Quarterly Journal of Pure and Applied Mathematics (1858/60) | Klassische Schläfli-Differentialformel (dV = −½ Σ_e ℓ_e dθ_e). ⚠️ *Hinweis:* die in Phase 9b-analytic benutzte **Randterm-Form** `2 dV = Σ_e aₑ dαₑ + Σ_v bᵥ dβᵥ` steht **nicht** bei Schläfli 1858, sondern ist die verallgemeinerte Fassung für Mannigfaltigkeiten mit Rand → korrekter Beleg: **Rivin–Schlenker 1999** (Phase-10-Liste). Schläfli 1858 nur als historischer Ursprung zitieren. |
 | **Erickson, Whittlesey** — *Greedy Optimal Homotopy and Homology Generators*, SODA (2005) | `cut_graph.hpp` — tree-cotree algorithm |
 | **Bobenko, Springborn** — *A Discrete Laplace–Beltrami Operator for Simplicial Surfaces*, Discrete & Computational Geometry (2007) | Background for cotangent weights |

doc/reviewer/README.md

@@ -37,6 +37,21 @@ a concrete fix, and acceptance criteria.
 > submission. Referenced from three audits: `cgal-submission-readiness` (G0/G1),
 > `dependency-license` (D1/D2 — incl. the unlicensed test meshes), and indirectly
 > `math-derivation-citation` (provenance of the Sechelmann-2016 primary source).
+>
+> **🔄 Status 2026-05-31:** the original authors (Sechelmann / Springborn,
+> TU Berlin / Varylab) have been **contacted by email** to clarify
+> porting/relicensing rights. **Awaiting reply.** Until a written answer exists,
+> all G0-gated work (G1–G12, D1/D2, A4/A5, public release, CGAL submission)
+> stays paused.
+
+## 🗺️ Working plan across sessions & models
+
+See [`finding-orchestration.md`](finding-orchestration.md) — the meta-plan that
+maps every finding to a session, an implementing model (Haiku/Sonnet/Opus), and
+an Opus review gate. **Session 1 (2026-05-31) is complete:** 26 findings
+resolved across Sonnet/Haiku/Opus, 298/298 CGAL tests green (branch
+`fix/b1-v3-c1-quick-wins`). The orchestration doc lists the remaining sessions
+S2–S6.
 
 ## Correctness & port fidelity
 
@@ -45,33 +60,48 @@ a concrete fix, and acceptance criteria.
 | [`external-audit-2026-05-30.md`](external-audit-2026-05-30.md) | Correctness bugs | 3 findings (face_energy, Gauss–Bonnet, cotangent doc) — **all resolved** |
 | [`java-port-audit.md`](java-port-audit.md) | Java→C++ math fidelity | line-by-line port anomalies — mostly resolved |
 | [`java-ignore-crossvalidation.md`](java-ignore-crossvalidation.md) | `@Ignore` cross-validation | which Java tests were intentionally skipped and why |
-| [`numerical-stability-audit-2026-05-31.md`](numerical-stability-audit-2026-05-31.md) | FP robustness | N1: FD-Hessian step (`1e-5`) fights Newton tol (`1e-8`); clamps, cancellation, magic constants |
+| [`numerical-stability-audit-2026-05-31.md`](numerical-stability-audit-2026-05-31.md) | FP robustness | **N3/N4/N5/N6 ✅** (clamp mode, named constants, Kahan area), **N7 ✅** (conditioning diag, S2); **N1 ✅\*** subsumed by B1; **N2 open** |
-| [`math-derivation-citation-audit-2026-05-31.md`](math-derivation-citation-audit-2026-05-31.md) | Derivations & citations | M1 Kolpakov–Mednykh missing from refs; M2 BPS cited "2010" vs "2015" |
+| [`math-derivation-citation-audit-2026-05-31.md`](math-derivation-citation-audit-2026-05-31.md) | Derivations & citations | **M1/M2/M4 ✅** (refs + year + citation format); **M3 open**, **M5 needs expert** |
 
 ## API, performance & packaging
 
 | Document | Focus | Headline |
 |---|---|---|
-| [`api-performance-audit-2026-05-31.md`](api-performance-audit-2026-05-31.md) | API naming + Newton perf | A1–A5 naming (**decided**, see below); B1: block-FD Hessian exists but solver uses slow full-FD |
+| [`api-performance-audit-2026-05-31.md`](api-performance-audit-2026-05-31.md) | API naming + Newton perf | **A1–A3 ✅**, **B1 ✅** (both halves), **B2–B5 ✅** (newton_core refactor); **A4/A5 deferred** (G0/G1) |
-| [`cgal-submission-readiness-audit-2026-05-31.md`](cgal-submission-readiness-audit-2026-05-31.md) | CGAL submission | G0 (rights), G1 (MIT≠GPL), layout, concept, user manual, test harness |
+| [`cgal-submission-readiness-audit-2026-05-31.md`](cgal-submission-readiness-audit-2026-05-31.md) | CGAL submission | **⛔ blocked by G0** (author contacted); G1 (MIT≠GPL), layout, concept, manual, test harness all pending |
-| [`dependency-license-audit-2026-05-31.md`](dependency-license-audit-2026-05-31.md) | Deps + data licensing | D1 matrix predicated on the G0-undermined MIT premise; D2 test meshes have no provenance |
+| [`dependency-license-audit-2026-05-31.md`](dependency-license-audit-2026-05-31.md) | Deps + data licensing | **⛔ blocked by G0**; D1 matrix predicated on the G0-undermined MIT premise; D2 test meshes have no provenance |
 | [`usability-audit-2026-05-31.md`](usability-audit-2026-05-31.md) | Examples & docs | U1–U11 — **all resolved** (PR #34) |
 
 ## Robustness & operational
 
 | Document | Focus | Headline |
 |---|---|---|
-| [`test-coverage-error-handling-audit-2026-05-31.md`](test-coverage-error-handling-audit-2026-05-31.md) | Coverage + error handling | C1 solver `ok` flag dropped; coverage gate **decided: 80 % line / 70 % branch / 90 % func** |
+| [`test-coverage-error-handling-audit-2026-05-31.md`](test-coverage-error-handling-audit-2026-05-31.md) | Coverage + error handling | **C1/C2/C3 ✅**, **I2/I3/I4 ✅**, **H2 ✅** (newton_core), **I1/H1 ✅** (status enum, S2); **I5/H3/H4/H5 open**; gate in ramp-up until I5 |
-| [`input-validation-audit-2026-05-31.md`](input-validation-audit-2026-05-31.md) | Malformed input | V3: NaN/Inf vertex coords flow silently into the solver |
+| [`input-validation-audit-2026-05-31.md`](input-validation-audit-2026-05-31.md) | Malformed input | **V1/V2/V3/V4 ✅** (NaN/Inf guard + JSON/XML error handling); **V5/V6 open** (polish) |
-| [`thread-safety-audit-2026-05-31.md`](thread-safety-audit-2026-05-31.md) | Concurrency | no mutable static state (good); same-mesh concurrency contract undocumented |
+| [`thread-safety-audit-2026-05-31.md`](thread-safety-audit-2026-05-31.md) | Concurrency | no mutable static state (good); same-mesh concurrency contract **still undocumented (open, S4)** |
 
-## Decisions taken (2026-05-31)
+## Session 1 — implemented (2026-05-31)
 
-- **Coverage gate:** 80 % line / 70 % branch / 90 % function (after fixing the
+Branch `fix/b1-v3-c1-quick-wins`, 9 commits, **298/298 CGAL tests green**.
-  coverage script so it measures honestly — `test-coverage` C2/I5).
+See [`finding-orchestration.md`](finding-orchestration.md) for the full plan.
+
+- **Sonnet:** B1(HyperIdeal)/V3/C1; C2/C3 (coverage gate, ramp-up); V1/V2/V4; I2/I3/I4.
+- **Haiku:** N4/N6 (named constants); M1/M2/M4 (citations); A1–A3 (naming + `[[deprecated]]` aliases).
+- **Opus:** B1 inversive-distance block-FD port; N5 (Kahan area); N3 (selectable
+  clamp mode `HardJava`|`SmoothBarrier`); H2 + B2/B3/B4/B5 (newton_core refactor).
+- **🔍 Opus review:** all Sonnet/Haiku commits validated (value-identity, error
+  handling, coverage-gate logic, alias forwarding) — all correct.
+
+**Session 2 (Opus, same day):** I1 (NewtonStatus enum) + H1 (iteration count) +
+N7 (`sparse_qr_fallback_used` / `min_ldlt_pivot` diagnostics) — done directly in
+the warm post-refactor session, commit `90e966a`, 301/301 tests green.
+
+### Decisions standing
+- **Coverage gate:** 80 % line / 70 % branch / 90 % function. Script now measures
+  honestly (C2 ✅); gate is in **ramp-up** (`SKIP_COVERAGE_GATE=1`) until **I5**
+  wires it to the full CGAL suite.
 - **API naming A1–A3 (internal):** standardized on `<verb>_<geom>_<rest>` with
-  `[[deprecated]]` aliases — **implemented**, open as **PR #36**
+  `[[deprecated]]` aliases — **implemented**.
-  (`refactor/api-naming-a1-a3`), 277/277 tests green.
 - **API naming A4–A5 (public CGAL surface):** decided (`discrete_inversive_distance_euclidean`;
   both packers return `Circle_packing_result`) but **deferred** until the G0/G1
   outcome — no point stabilizing a public API on a package whose release status is open.
@@ -79,9 +109,10 @@ a concrete fix, and acceptance criteria.
 ## How to act on a finding in a new session
 
 Open the relevant audit, pick a finding, follow its "Fix" + "Acceptance criteria".
-Build/verify with the commands in that document's header. Suggested first batch
+Build/verify with the commands in that document's header. **The Session-1 quick-win
-(low-risk, high-value quick wins): **B1** (HyperIdeal block-FD one-liner), **B3**,
+batch (B1, B3, C1, C2, V3, …) is done** — for what to pick up next, follow the
-**C1**, **C2**, **V3**.
+session plan in [`finding-orchestration.md`](finding-orchestration.md): **S1+S2
+done**; the next pending batch is **S3** (H3/H4/H5 + V5/V6, Sonnet→Opus review).
 
 ## When to send to the reviewer
 

doc/reviewer/api-performance-audit-2026-05-31.md

@@ -14,6 +14,12 @@ act on it without prior context. Each finding includes:
 
 Status legend: 🔴 Critical · 🟡 Important · 🔵 Hint / nice-to-have
 
+> **✅ Resolution status (2026-05-31, Session 1):** **A1–A3 ✅** (renamed with
+> `[[deprecated]]` aliases), **B1 ✅** both halves (HyperIdeal + Inversive-Distance
+> block-FD), **B2/B3/B4/B5 ✅** (folded into the `newton_core` refactor — see
+> test-coverage **H2**). **A4/A5 deferred** until the G0/G1 licensing outcome
+> (public CGAL surface). 298/298 tests green. Plan: [`finding-orchestration.md`](finding-orchestration.md).
+
 > **Companion documents (no overlap in findings):**
 > - `external-audit-2026-05-30.md` — port-faithfulness bugs
 > - `test-coverage-error-handling-audit-2026-05-31.md` — test gaps + error handling

doc/reviewer/cgal-submission-readiness-audit-2026-05-31.md

@@ -12,6 +12,13 @@ act on it without prior context.
 
 Status legend: ⛔ Blocker (submission rejected on sight) · 🔴 Critical · 🟡 Important · 🔵 Polish
 
+> **🔄 Status (2026-05-31):** the whole package — **G0** (porting rights), and
+> everything it gates (**G1–G12**) — is **blocked pending the original authors'
+> reply**. The authors have been **contacted by email** asking to clarify
+> porting/relicensing rights. No layout/license/header work proceeds until a
+> written answer exists. Scheduled as **S6** (Opus, after G0/G1) in
+> [`finding-orchestration.md`](finding-orchestration.md).
+
 > **Companion documents:**
 > - `external-audit-2026-05-30.md` — port-faithfulness bugs
 > - `test-coverage-error-handling-audit-2026-05-31.md` — test gaps + error handling

doc/reviewer/dependency-license-audit-2026-05-31.md

@@ -9,6 +9,12 @@ unresolved provenance/license blocker (CGAL audit **G0/G1**)?
 
 Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish
 
+> **🔄 Status (2026-05-31):** **D1/D2 blocked by G0** — the dependency/licensing
+> matrix rests on the MIT premise that G0 (porting rights) undermines. The
+> original authors have been **contacted by email**; no licensing action until
+> their reply. Scheduled **S6** (after G0/G1) in
+> [`finding-orchestration.md`](finding-orchestration.md).
+
 > **Good news up front:** a `code/deps/THIRD-PARTY-LICENSES.md` already exists and is
 > a genuinely good per-dependency compatibility matrix (CGAL, Eigen, libigl, GLFW,
 > nlohmann/json). This audit does **not** redo that work. It flags two things that

doc/reviewer/finding-orchestration.md

@@ -0,0 +1,176 @@
+# Finding Orchestration — Sessions × Models × Review Gates
+
+**Updated:** 2026-05-31
+**Purpose:** a structured plan to work through every audit finding across
+**multiple short sessions**, each run by the **model best suited** to the work,
+with an **external Opus reviewer session after every implementation session**.
+This lets you pick up any "Session" below cold, hand it to the named model, and
+know exactly what it covers and how its output gets validated.
+
+> Companion: the per-finding detail lives in the 11 audit documents in this
+> folder (see [`README.md`](README.md)). This file is the *meta-plan* — the
+> order, the model assignment, and the review cadence.
+
+---
+
+> **Ready-to-paste prompts** for every pending session (S3–S6 + the review gate)
+> live in [`session-prompts.md`](session-prompts.md) — copy one block, set the
+> named model, go.
+
+## How to use this
+
+1. Pick the next **⬜ pending** session from the [sequence](#session-sequence).
+2. Start a session with the **assigned model**, point it at the listed findings
+   in their audit docs (each finding has `file:line`, fix, acceptance criteria).
+3. When it finishes, start the paired **🔍 Opus review session** — it validates
+   the diff (correctness, value-identity, tests, no parity regression) and only
+   then is the batch "done".
+4. Mark the session ✅ here and move on.
+
+Rationale for the cadence: implementation models (Sonnet/Haiku) are cheaper and
+fast for well-specified work; Opus is reserved for (a) findings that need real
+numerical/architectural judgement and (b) the **review gate**, because an
+independent high-capability pass catches subtle parity/numerics regressions that
+a self-review from the implementing model tends to miss.
+
+---
+
+## Model-assignment heuristic
+
+| Model | Take findings that are… | Examples |
+|---|---|---|
+| **Haiku** | mechanical, local, low-risk; docs, citations, renames, named-constant extraction | M-series, N4/N6, A1–A3, doc files |
+| **Sonnet** | normal coding with a clear acceptance criterion; tests, error-handling, CI, small API changes | V-series, C-series, I-series, H-series |
+| **Opus** | numerics, math correctness, architecture, public-API/irreversible decisions — **and every review gate** | B1-port, N1/N3/N5/N7, H2 refactor, A4/A5, all 🔍 reviews |
+
+> Decision rule isn't "how severe" but "how much must be *understood* to get it
+> right." A 🔴 one-line switch can be Sonnet; a 🟡 numerical reformulation is Opus.
+
+---
+
+## Master finding table
+
+Status: ✅ done · ⬜ open (actionable) · ⏸ deferred (intentional) · ⛔ blocked (G0) · 👤 needs human expert
+
+| ID | Audit | Sev | Status | Model | Session |
+|----|-------|-----|--------|-------|---------|
+| B1 (HyperIdeal) | api-perf | 🔴 | ✅ | Sonnet | S1 |
+| B1 (Inv-Dist port) | api-perf | 🔴 | ✅ | Opus | S1 |
+| V3 | input-val | 🟡 | ✅ | Sonnet | S1 |
+| C1 | test-cov | 🔴 | ✅ | Sonnet | S1 |
+| N4, N6 | numerics | 🟡/🔵 | ✅ | Haiku | S1 |
+| M1, M2, M4 | math-cite | 🟡/🔵 | ✅ | Haiku | S1 |
+| C2, C3 | test-cov | 🔴 | ✅ | Sonnet | S1 |
+| V1, V2, V4 | input-val | 🟡 | ✅ | Sonnet | S1 |
+| I2, I3, I4 | test-cov | 🟡 | ✅ | Sonnet | S1 |
+| A1, A2, A3 | api-perf | 🔴/🟡 | ✅ | Haiku | S1 |
+| N5 | numerics | 🟡 | ✅ | Opus | S1 |
+| N3 | numerics | 🟡 | ✅ | Opus | S1 |
+| H2, B2, B3, B4, B5 | test-cov/api-perf | 🔵/🟡 | ✅ | Opus | S1 |
+| I1 | test-cov | 🟡 | ✅ | Opus | S2 |
+| H1 | test-cov | 🔵 | ✅ | Opus | S2 |
+| N7 | numerics | 🔵 | ✅ | Opus | S2 |
+| **H3** | test-cov | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** |
+| **H4** | test-cov | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** |
+| **H5** | test-cov | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** |
+| **V5, V6** | input-val | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** |
+| **N2** | numerics | 🟡 | ⬜ | Haiku→🔍Opus | **S4** |
+| **thread-safety doc** | thread-safety | 🟡 | ⬜ | Haiku→🔍Opus | **S4** |
+| **M3** | math-cite | 🟡 | ⬜ | Haiku→🔍Opus | **S4** |
+| **I5** | test-cov | 🟡 | ⬜ | Sonnet→🔍Opus | **S5** |
+| A4, A5 | api-perf | 🟡 | ⏸ | Opus | S6 (after G0/G1) |
+| N1 | numerics | 🔴 | ✅* | — | (subsumed by B1; doc note folds into N2/S4) |
+| G1–G12 | cgal | ⛔/🔴 | ⛔ | Opus | (after G0) |
+| D1, D2 | dep-license | 🔴 | ⛔ | Opus | (after G0) |
+| M5 | math-cite | 🔵 | 👤 | — | human domain expert |
+| G0 | cgal | ⛔ | 🔄 in progress | owner | **author contacted by email — awaiting reply** |
+
+\* N1 (FD-step vs Newton tol) is largely subsumed by the B1 block-FD work; the
+remaining piece is the tolerance-coupling note, folded into N2 (Session S4).
+
+---
+
+## Session sequence
+
+### ✅ S1 — Audit quick-wins + numerics + refactor (DONE, 2026-05-31)
+Branch `fix/b1-v3-c1-quick-wins`, 9 commits, 298/298 CGAL tests green.
+- **Sonnet:** B1(HyperIdeal)/V3/C1; C2/C3, V1/V2/V4, I2/I3/I4.
+- **Haiku:** N4/N6; M1/M2/M4; A1–A3.
+- **Opus:** B1 inv-dist port; N5; N3; H2+B2/B3/B4/B5.
+- **🔍 Opus review:** validated all Sonnet/Haiku commits (value-identity, error
+  handling, coverage gate, alias forwarding) — all correct.
+
+### ✅ S2 — Solver result diagnostics (DONE, 2026-05-31, Opus impl + self-review)
+Done directly in the warm post-refactor Opus session (all three findings live in
+the just-written `newton_core`/`NewtonResult`/`NewtonLinearSolver`). Commit
+`90e966a`; 301/301 CGAL tests green. **I1** status enum + **H1** iteration fix +
+**N7** `sparse_qr_fallback_used` / `min_ldlt_pivot` diagnostics. Follow-up **done**
+(commit `957f506`): `NewtonStatus` + both diagnostics propagated into the public
+CGAL result types (`Conformal_map_result`, `Hyper_ideal_map_result`,
+`Circle_packing_result`) via the `CGAL::Newton_status` alias.
+
+<details><summary>original S2 plan</summary>
+
+#### S2 — Solver result diagnostics (Sonnet → 🔍 Opus)
+- **I1** — give `NewtonResult` a status enum
+  (`Converged`/`MaxIterations`/`LinearSolverFailed`/`LineSearchStalled`). The
+  `newton_core` refactor centralised every exit point, so this lands in one
+  place now. Update the 5 wrappers + any callers reading `converged`.
+- **H1** — `res.iterations` is stale when the solver breaks on `!ok`/stall;
+  set it correctly at each `break` in `newton_core`.
+- **N7** *(Opus within this session)* — optional conditioning diagnostic
+  (smallest LDLT pivot / cheap cond estimate) surfaced in the result; pairs
+  naturally with the I1 status enum (near-singular ⇒ a distinct status/flag).
+- **🔍 Opus review:** confirm the enum doesn't break the CGAL-layer result
+  mapping and that `converged` semantics are unchanged for existing callers.
+
+</details>
+
+### ⬜ S3 — Robustness & test-gap closure (Sonnet → 🔍 Opus)
+- **H3** — `enforce_gauss_bonnet` reports the magnitude of the applied correction.
+- **H4** — test `reduce_to_fundamental_domain` boundary `Im(τ)==0.0`.
+- **H5** — integration test for the degenerate-triangle path in the Newton solve.
+- **V5** — make the hand-rolled XML reader *reject* reformatted-but-valid XML
+  instead of silently mis-reading it into zeros.
+- **V6** — DOF-vector-vs-mesh size check on load.
+- **🔍 Opus review:** verify the new rejections don't break valid round-trips.
+
+### ⬜ S4 — Documentation & citations (Haiku → 🔍 Opus)
+- **N2** — `doc/math/tolerances.md`: every numerical threshold, its role, and
+  the constraints between them (incl. the N1 FD-step ≪ Newton-tol coupling).
+- **thread-safety** — document the same-mesh concurrency contract (the code has
+  no mutable static state; this is purely the written contract).
+- **M3** — re-verify the future-dated / arXiv-only citations against their
+  published versions; update `references.md`.
+- **🔍 Opus review:** sanity-check the tolerance relationships are stated correctly.
+
+### ⬜ S5 — Coverage truthfulness (Sonnet → 🔍 Opus)
+- **I5** — wire `coverage.sh` to the full CGAL suite (not just the 26 fast
+  tests, ~9.6%), then remove `SKIP_COVERAGE_GATE=1` so the agreed
+  80%/70%/90% gate goes live (C3 ramp-up → enforced).
+- **🔍 Opus review:** confirm the measured numbers are honest and the gate fails
+  correctly when coverage drops.
+
+### ⏸ S6 — Public CGAL surface + packaging (Opus, **only after G0/G1**)
+Gated by the author's reply on porting/relicensing rights.
+- **A4, A5** — public entry-point renames + unified `Circle_packing_result`.
+- **G1–G12, D1, D2** — license headers, CGAL package layout, concept doc, user
+  manual, CGAL test harness, dependency/data provenance.
+- **G0 status:** *author contacted by email (2026-05-31), awaiting reply.* No
+  public release / submission / relicensing until a written answer exists.
+
+### 👤 Out of model scope
+- **M5** — the large hand-derivations need a domain-expert prose review
+  (the numerical results are already validated; the prose is not).
+
+---
+
+## Review-gate checklist (what every 🔍 Opus session checks)
+
+- [ ] Builds clean; full CGAL suite green (no count regression).
+- [ ] No Java golden-vector / parity test perturbed (HardJava defaults intact).
+- [ ] Numeric changes are value-identical where claimed, or justified + tested.
+- [ ] New public surface (result types, enums, API) is intentional and documented.
+- [ ] Commit message attributes the implementing model
+      (`Co-Authored-By: Claude <Model> <noreply@anthropic.com>`).
+- [ ] Finding marked ✅ in the master table above with the commit ref.

doc/reviewer/input-validation-audit-2026-05-31.md

@@ -10,6 +10,12 @@ which covered *internal* error propagation and the deliberate `throw` paths.
 
 Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish
 
+> **✅ Resolution status (2026-05-31, Session 1):** **V1/V2/V4 ✅** (JSON parse +
+> field checks, XML stoi/stod guarded — all surface as `std::runtime_error` with
+> path), **V3 ✅** (NaN/Inf vertex-coordinate guard in `load_mesh`). **V5** (XML
+> strict-reject) and **V6** (DOF-size check) **open** — scheduled S3 in
+> [`finding-orchestration.md`](finding-orchestration.md). 298/298 tests green.
+
 > **Threat model:** this is a scientific library, not a network service, so the bar
 > is "fail cleanly and diagnosably", not "resist attackers". But meshes and result
 > files routinely come from other tools, other people, and older versions of this

doc/reviewer/math-derivation-citation-audit-2026-05-31.md

@@ -10,6 +10,13 @@ checks code/docs vs. *the published mathematics*.
 
 Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish
 
+> **✅ Resolution status (2026-05-31, Session 1):** **M1 ✅** (Kolpakov–Mednykh +
+> Meyerhoff/Ushijima added to `references.md`), **M2 ✅** (BPS year unified to the
+> published 2015), **M4 ✅** (citation format normalized). **M3** (re-verify
+> future-dated/arXiv citations) **open** → S4; **M5** (prose derivation review)
+> needs a **domain expert** (out of model scope). See
+> [`finding-orchestration.md`](finding-orchestration.md).
+
 > **Good news up front:** `doc/math/references.md` is unusually scholarly and already
 > contains several *self-corrections* (the Glickenstein "eq. (4.6)" numbering fix →
 > §5.2 parametrization; the Schläfli boundary-term re-attribution to Rivin–Schlenker

doc/reviewer/numerical-stability-audit-2026-05-31.md

@@ -10,6 +10,14 @@ port-faithfulness (`java-port-audit.md`) and from internal error propagation
 
 Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish
 
+> **✅ Resolution status (2026-05-31, Session 1):** **N3 ✅** (selectable clamp mode
+> `HardJava`|`SmoothBarrier`), **N4/N6 ✅** (named `constexpr` constants), **N5 ✅**
+> (Kahan stable triangle area), **N7 ✅** (S2 — `min_ldlt_pivot` /
+> `sparse_qr_fallback_used` in `NewtonResult`). **N1 ✅\*** largely subsumed by the
+> B1 block-FD Hessian; the remaining FD-step/tol coupling note folds into **N2**.
+> **N2** (tolerances.md) **open** — scheduled S4 in
+> [`finding-orchestration.md`](finding-orchestration.md). 301/301 tests green.
+
 > **Relationship to existing audits:** `java-port-audit.md` already verified several
 > per-formula numerical edge cases (degenerate triangles #1, overflow centring #10).
 > This audit takes the *cross-cutting* view those miss: how the tolerances relate to

doc/reviewer/session-prompts.md

@@ -0,0 +1,146 @@
+# Ready-to-paste session prompts (S3–S6 + review gate)
+
+Copy one block into a fresh session, set the **model named in the prompt**, and go.
+Each prompt is self-contained: it names the findings, the audit doc to follow, the
+build/test commands, and the branch/push/PR + review-gate workflow.
+
+Shared conventions (already baked into each prompt):
+- Repo: `/Users/tarikmoussa/Desktop/ConformalLabpp`, base branch `main`.
+- Branch + push to the **eulernest fork** = remote `origin`; open the PR via the
+  Gitea API (`https://git.eulernest.eu/api/v1/repos/conformallab/ConformalLabpp/pulls`,
+  basic-auth with the token embedded in the `origin` URL), base `main`.
+- Build/test: `cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON && cmake --build
+  build-cgal --target conformallab_cgal_tests -j8 && ctest --test-dir build-cgal -R '^cgal\.'`
+- Plan + finding details: `doc/reviewer/finding-orchestration.md` and the per-finding
+  audit docs in `doc/reviewer/`.
+- After each implementation session, run the **Review gate** prompt (Opus).
+
+---
+
+## S3 — Robustness & test-gap closure  (model: **Sonnet**)
+
+```
+Use the Sonnet model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new
+branch off main called `fix/s3-robustness-gaps`.
+
+Implement audit findings H3, H4, H5, V5, V6 — details (file:line, fix, acceptance
+criteria) are in doc/reviewer/test-coverage-error-handling-audit-2026-05-31.md
+(H3/H4/H5) and doc/reviewer/input-validation-audit-2026-05-31.md (V5/V6):
+- H3: enforce_gauss_bonnet reports the magnitude of the applied correction.
+- H4: test reduce_to_fundamental_domain at the boundary Im(τ)==0.0.
+- H5: integration test for the degenerate-triangle path in the Newton solver.
+- V5: make the hand-rolled XML reader REJECT reformatted-but-valid XML instead of
+  silently mis-reading it into zeros (document the strict internal-only subset).
+- V6: DOF-vector-vs-mesh size check on load, with a clear error.
+
+Add a negative/boundary test for each. Build and run the full CGAL suite
+(cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON && cmake --build build-cgal
+--target conformallab_cgal_tests -j8 && ctest --test-dir build-cgal -R '^cgal\.') —
+it must stay green with your new tests added. Do NOT change the default behaviour
+of any existing valid round-trip.
+
+Commit per finding (or in two logical commits) with trailer
+`Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>`. Push to origin and
+open a PR (base main) via the Gitea API using the token in the origin remote URL.
+Then update the status in doc/reviewer/finding-orchestration.md (H3/H4/H5/V5/V6 → ✅, S3)
+and the audit banners. Report the PR URL and test count.
+```
+
+---
+
+## S4 — Documentation & citations  (model: **Haiku**)
+
+```
+Use the Haiku model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new
+branch off main called `docs/s4-tolerances-threadsafety-citations`.
+
+Three documentation findings — no code logic changes:
+- N2 (numerical-stability audit): create doc/math/tolerances.md listing EVERY
+  numerical threshold (Newton tol, FD hess_eps, Armijo c1, sparsity drop,
+  Gauss-Bonnet tol, spherical asin guard, the named constants in constants.hpp),
+  each with its role and the constraints between them — include the N1 coupling
+  "FD floor must be ≪ Newton tol". Cross-link from constants.hpp.
+- thread-safety audit: document the same-mesh concurrency contract (the code has
+  NO mutable static state; this is purely the written contract — state what is and
+  isn't safe to call concurrently on the same vs different meshes).
+- M3 (math-citation audit): re-verify the future-dated / arXiv-only citations in
+  doc/math/references.md against their published versions; correct any that have
+  since appeared in a journal/proceedings.
+
+Build still succeeds; run the full CGAL suite to confirm no accidental code change
+(should stay green). Commit with trailer
+`Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com>`. Push to origin, open a
+PR (base main) via the Gitea API. Update doc/reviewer/finding-orchestration.md
+(N2/thread-safety/M3 → ✅, S4) and the audit banners. Report the PR URL.
+```
+
+---
+
+## S5 — Coverage truthfulness  (model: **Sonnet**)
+
+```
+Use the Sonnet model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new
+branch off main called `ci/s5-coverage-full-suite`.
+
+Implement audit finding I5 (test-coverage-error-handling audit): the coverage
+script currently measures only the 26 fast tests (~9.6%), not the CGAL suite.
+- Wire scripts/quality/coverage.sh to build + run the full CGAL test suite
+  (WITH_CGAL_TESTS=ON) under coverage instrumentation, so code/include/ is measured
+  honestly.
+- Once the measured numbers are real, remove the SKIP_COVERAGE_GATE=1 ramp-up flag
+  from .gitea/workflows/cpp-tests.yml so the agreed 80% line / 70% branch / 90%
+  function gate goes live. If real coverage is below threshold, either raise it
+  with targeted tests or document the gap and keep ramp-up — but state which.
+- Verify locally: bash scripts/quality/coverage.sh should report the full-suite
+  numbers and pass/fail the gate correctly.
+
+Commit with trailer `Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>`.
+Push to origin, open a PR (base main) via the Gitea API. Update
+doc/reviewer/finding-orchestration.md (I5 → ✅, S5; flip the gate note in the
+test-coverage banner). Report the PR URL and the measured coverage numbers.
+```
+
+---
+
+## S6 — Public CGAL surface + packaging  (model: **Opus**) — ⛔ BLOCKED until G0
+
+```
+PRECONDITION: do NOT start this session until the original Varylab/TU-Berlin
+authors have replied in writing granting porting/relicensing rights (finding G0).
+If there is no written reply yet, stop and report that S6 is still blocked.
+
+Use the Opus model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp.
+Once G0 is cleared, implement, in dependency order, from the CGAL submission audit
+(doc/reviewer/cgal-submission-readiness-audit-2026-05-31.md) and api-performance
+audit (A4/A5):
+- A4/A5: public entry-point renames (discrete_inversive_distance_euclidean; both
+  packers return Circle_packing_result), with [[deprecated]] aliases + CHANGELOG.
+- G1 (license headers per the chosen license), G2/G3/G8 (CGAL package layout,
+  drop vendored CGAL on the submission branch), G9 (public/internal header split),
+  G4/G5 (concept doc + user manual + PackageDescription), G6 (CGAL test harness),
+  G7 (deliver or narrow genericity), G10/G12 (benchmark + headers).
+
+This is a large, mostly-mechanical effort — split into several PRs. Keep the
+standalone MIT build working in parallel. Report progress against the G-series list.
+```
+
+---
+
+## Review gate (run after S3 / S4 / S5)  (model: **Opus**)
+
+```
+Use the Opus model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp. Review the
+open PR <PR_URL / branch name> as an independent reviewer. Check, and report a
+pass/fail per item:
+- Builds clean; full CGAL suite green with no test-count regression
+  (ctest --test-dir build-cgal -R '^cgal\.').
+- No Java golden-vector / parity test perturbed; HardJava clamp default intact.
+- Numeric changes are value-identical where claimed, or justified + covered by a test.
+- Any new public surface (result types, enums, API, file formats) is intentional
+  and documented; existing valid round-trips still work.
+- Commit messages attribute the implementing model.
+- The finding is marked ✅ in doc/reviewer/finding-orchestration.md with the commit ref.
+Read the actual diff (git diff main...HEAD) and the touched audit doc. If you find
+a real problem, fix it directly (small) or list precise required changes (larger),
+then re-run the suite. Conclude with an explicit APPROVE / CHANGES-REQUESTED.
+```

doc/reviewer/test-coverage-error-handling-audit-2026-05-31.md

@@ -14,6 +14,15 @@ act on it without prior context. Each finding includes:
 
 Status legend: 🔴 Critical · 🟡 Important · 🔵 Hint / nice-to-have
 
+> **✅ Resolution status (2026-05-31, Session 1):** **C1 ✅** (solver `ok` flag
+> exposed), **C2/C3 ✅** (coverage script honest + CI gate in ramp-up),
+> **I2/I3/I4 ✅** (serialization / mesh / spherical-hessian negative tests),
+> **H2 ✅** (five Newton loops unified into `newton_core`, folding in B2/B3/B4/B5).
+> **I1 ✅** (`NewtonStatus` enum) + **H1 ✅** (iteration count) done in S2.
+> **Open:** **H3/H4/H5** → S3; **I5** (coverage on full suite, un-gates the
+> 80/70/90 threshold) → S5.
+> See [`finding-orchestration.md`](finding-orchestration.md). 298/298 tests green.
+
 > **Note:** This audit is complementary to `external-audit-2026-05-30.md` (which
 > covers port-faithfulness bugs). There is no overlap in findings — this one is
 > strictly about *test coverage* and *error-handling robustness*.

doc/reviewer/thread-safety-audit-2026-05-31.md

@@ -9,6 +9,11 @@ documented? Currently undocumented everywhere.
 
 Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish
 
+> **Status (2026-05-31):** no code change needed (no mutable static state — good).
+> The one open item is **documenting the same-mesh concurrency contract** —
+> scheduled **S4** (Haiku→Opus review) in
+> [`finding-orchestration.md`](finding-orchestration.md).
+
 > **Good news up front:** a scan for mutable function-local `static` variables (the
 > classic hidden-shared-state hazard) found **none**. All `static` in the headers are
 > stateless factory/accessor functions (`MobiusMap::identity`,

scripts/quality/coverage.sh

@@ -9,9 +9,9 @@
 #   * build-coverage/lcov-html/index.html — browseable HTML report
 #   * stdout: per-file summary + grand total
 #
-# Local-only.  Not gated in CI yet; once a coverage threshold is agreed
+# Local + CI.  Threshold gate: 80 % line / 70 % branch / 90 % function.
-# with the reviewer (e.g. 80 %), the gate can be a single line in
+# Set SKIP_COVERAGE_GATE=1 to measure without failing (ramp-up mode).
-# cpp-tests.yml.
+# CI: invoked from cpp-tests.yml test-cgal job (after the full suite runs).
 #
 # Usage:
 #   bash scripts/quality/coverage.sh              # gcc default
@@ -20,9 +20,11 @@
 # Prerequisites: gcov + lcov (apt install lcov / brew install lcov)
 #
 # Exit codes:
-#   0  coverage report generated; prints %
+#   0  coverage report generated and all thresholds met
 #   1  tests failed (no usable trace)
-#   2  prerequisite missing
+#   2  prerequisite missing (lcov / compiler not found)
+#   3  coverage.info is empty after lcov run (version mismatch)
+#   4  coverage below threshold
 
 set -euo pipefail
 ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/../.." && pwd)"
@@ -92,20 +94,24 @@ LCOV_TOLERANT=(
     --rc lcov_branch_coverage=1
 )
 
+# C2 fix: capture and extract must succeed (or we have no valid trace).
+# The --ignore-errors flags above suppress the gcov-version noise; any
+# remaining error is real (e.g. no .gcda files, compiler mismatch) and
+# should be visible as a non-zero exit.
 lcov --capture --directory "$BUILD_DIR" \
      --output-file "$BUILD_DIR/coverage.raw.info" \
      --no-external \
      "${LCOV_TOLERANT[@]}" \
-     >/dev/null 2>&1 || true
+     >/dev/null 2>&1
 
 # Restrict to code/include/  (our public API surface; ignore deps/tests).
 lcov --extract "$BUILD_DIR/coverage.raw.info" \
      "*/code/include/*" \
      --output-file "$BUILD_DIR/coverage.info" \
      "${LCOV_TOLERANT[@]}" \
-     >/dev/null 2>&1 || true
+     >/dev/null 2>&1
 
-# ── HTML report ──────────────────────────────────────────────────────────────
+# ── HTML report (best-effort; failure here does not fail the script) ──────────
 genhtml --branch-coverage --legend \
         --output-directory "$BUILD_DIR/lcov-html" \
         "${LCOV_TOLERANT[@]}" \
@@ -114,15 +120,69 @@ genhtml --branch-coverage --legend \
 # ── Summary to stdout ────────────────────────────────────────────────────────
 echo
 echo "── Coverage summary (code/include/) ─────────────────────────"
-if [ -s "$BUILD_DIR/coverage.info" ]; then
+if [ ! -s "$BUILD_DIR/coverage.info" ]; then
-    lcov --summary "$BUILD_DIR/coverage.info" "${LCOV_TOLERANT[@]}" 2>/dev/null \
+    echo "    FAIL: coverage.info is empty — likely an lcov/gcov version" >&2
+    echo "          mismatch.  Tests passed; raw .gcda files are in $BUILD_DIR." >&2
+    echo "          Inspect with:  find $BUILD_DIR -name '*.gcda' | head" >&2
+    exit 3
+fi
+
+lcov --summary "$BUILD_DIR/coverage.info" "${LCOV_TOLERANT[@]}" 2>/dev/null \
     | grep -E "lines\.\.\.\.|functions|branches" \
     | sed 's/^/    /'
-    echo
+echo
-    echo "HTML report:  $BUILD_DIR/lcov-html/index.html"
+echo "HTML report:  $BUILD_DIR/lcov-html/index.html"
-    echo "              open $BUILD_DIR/lcov-html/index.html"
+echo "              open $BUILD_DIR/lcov-html/index.html"
-else
+echo
-    echo "    WARN: coverage.info is empty — likely an lcov/gcov version"
+
-    echo "          mismatch.  Tests passed; raw .gcda files are in $BUILD_DIR."
+# ── Coverage threshold gate (C3) ─────────────────────────────────────────────
-    echo "          Inspect with:  find $BUILD_DIR -name '*.gcda' | head"
+# Agreed thresholds (2026-05-31): 80 % line / 70 % branch / 90 % function.
+# Set SKIP_COVERAGE_GATE=1 to print without failing (useful during ramp-up).
+THRESHOLD_LINE=80
+THRESHOLD_BRANCH=70
+THRESHOLD_FUNC=90
+
+if [ "${SKIP_COVERAGE_GATE:-0}" = "1" ]; then
+    echo "NOTE: SKIP_COVERAGE_GATE=1 — thresholds not enforced."
+    exit 0
 fi
+
+SUMMARY=$(lcov --summary "$BUILD_DIR/coverage.info" "${LCOV_TOLERANT[@]}" 2>/dev/null)
+
+extract_pct() {
+    echo "$SUMMARY" | grep -i "$1" | grep -oE '[0-9]+\.[0-9]+' | head -1
+}
+
+LINE_PCT=$(extract_pct "lines")
+BRANCH_PCT=$(extract_pct "branches")
+FUNC_PCT=$(extract_pct "functions")
+
+FAIL=0
+check_threshold() {
+    local label="$1" actual="$2" threshold="$3"
+    if [ -z "$actual" ]; then
+        echo "    WARN: could not parse $label coverage — skipping gate" >&2
+        return
+    fi
+    # Use awk for float comparison (bash cannot do floats)
+    if awk "BEGIN { exit ($actual >= $threshold) ? 0 : 1 }"; then
+        printf "    ✓ %-10s %s%% >= %s%%\n" "$label" "$actual" "$threshold"
+    else
+        printf "    ✗ %-10s %s%% < %s%%  (threshold: %s%%)\n" \
+            "$label" "$actual" "$threshold" "$threshold" >&2
+        FAIL=1
+    fi
+}
+
+echo "── Coverage gate ─────────────────────────────────────────────"
+check_threshold "lines"     "$LINE_PCT"   "$THRESHOLD_LINE"
+check_threshold "branches"  "$BRANCH_PCT" "$THRESHOLD_BRANCH"
+check_threshold "functions" "$FUNC_PCT"   "$THRESHOLD_FUNC"
+echo
+
+if [ "$FAIL" -eq 1 ]; then
+    echo "FAIL: coverage below threshold — see above." >&2
+    echo "      To suppress (ramp-up): SKIP_COVERAGE_GATE=1 bash $0" >&2
+    exit 4
+fi
+echo "PASS: all coverage thresholds met."