docs(roadmap): add phase orchestration system mirroring the reviewer audit workflow

Brings doc/roadmap/ to the same operational level as doc/reviewer/ by adding
the two missing structural files and updating existing docs to close the gap
identified in the system review.

New files:
- doc/roadmap/phase-orchestration.md  — master phase table (model assignments,
  session IDs, status, review-gate checklist); mirrors finding-orchestration.md
- doc/roadmap/session-prompts.md      — copy-paste-ready prompts for P1–P4
  (Haiku/Sonnet/Opus) + universal Opus review gate; mirrors reviewer/session-prompts.md

Updated files:
- CLAUDE.md: roadmap section now lists porting-status.md, phase-orchestration.md,
  session-prompts.md; reviewer section adds finding-orchestration.md and
  session-prompts.md; agentic-workflow section has direct "S3 is next / P1 is next"
  entry points so agents don't need to derive the next action from scratch
- doc/roadmap/phases.md: Current-focus table at the top (7 tracks, next session
  per track, gating); cross-links to geometry-central-comparison.md,
  software-landscape.md, complexity.md added in GC and 9b-analytic sections
- doc/roadmap/porting-status.md: snapshot date updated to 2026-05-31 (post S1+S2);
  test count replaced by reference to doc/api/tests.md; §4 gains four new solver
  rows (newton_core refactor, NewtonStatus enum, diagnostics, selectable clamp mode);
  §7 gains four new research-extension rows from S1/S2
- doc/roadmap/research-track.md: header companion-docs section links to
  novelty-statement.md, software-landscape.md, and phase-orchestration.md

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

.gitea/workflows/cpp-tests.yml

@@ -3,11 +3,13 @@ name: C++ Tests
 # Trigger keywords in commit message (checked via head_commit.message):
 #   /test-cgal      — full CGAL test suite (277 tests, ~5 min build)
 #   /quality-gates  — license, codespell, shellcheck, CGAL conventions
-#   /ci-all         — all of the above + /docs + /links (across all workflows)
+#
+# ⚠  /ci-all was removed: the Pi runner (3-4 GB RAM) cannot sustain
+#    multiple Docker containers simultaneously.  Use one keyword per commit.
 #
 # Examples:
 #   git commit -m "fix: correct angle formula /test-cgal"
-#   git commit -m "release prep /ci-all"
+#   git commit -m "chore: update headers /quality-gates"
 #
 # test-fast always runs on every push — it is fast (< 5 s) and cheap.
 
@@ -15,7 +17,6 @@ on:
   push:
     branches:
       - main
-      - dev
       - "claude/**"
       - "feature/**"
       - "review/**"
@@ -31,6 +32,7 @@ jobs:
     runs-on: eulernest
     container:
       image: git.eulernest.eu/conformallab/ci-cpp:latest
+      options: "--memory=800m --memory-swap=1200m"
 
     steps:
       - uses: actions/checkout@v4
@@ -39,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: >
@@ -76,9 +82,7 @@ jobs:
 # ─────────────────────────────────────────────────────────────────────────────
   test-cgal:
     needs: test-fast
-    if: |
-      contains(github.event.head_commit.message, '/test-cgal') ||
-      contains(github.event.head_commit.message, '/ci-all')
+    if: contains(github.event.head_commit.message, '/test-cgal')
     runs-on: eulernest
     container:
       image: git.eulernest.eu/conformallab/ci-cpp:latest
@@ -133,12 +137,11 @@ jobs:
 # Cheap (~30 s): license headers, CGAL conventions, codespell, shellcheck.
 # ─────────────────────────────────────────────────────────────────────────────
   quality-gates:
-    if: |
-      contains(github.event.head_commit.message, '/quality-gates') ||
-      contains(github.event.head_commit.message, '/ci-all')
+    if: contains(github.event.head_commit.message, '/quality-gates')
     runs-on: eulernest
     container:
       image: git.eulernest.eu/conformallab/ci-cpp:latest
+      options: "--memory=600m --memory-swap=900m"
 
     steps:
       - uses: actions/checkout@v4
@@ -161,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/)"

.gitea/workflows/doc-build.yaml

@@ -10,7 +10,6 @@ on:
   push:
     branches:
       - main
-      - dev
       - "claude/**"
       - "feature/**"
       - "review/**"
@@ -27,11 +26,11 @@ jobs:
   doc-build:
     if: |
       github.event_name == 'workflow_dispatch' ||
-      contains(github.event.head_commit.message, '/docs') ||
-      contains(github.event.head_commit.message, '/ci-all')
+      contains(github.event.head_commit.message, '/docs')
     runs-on: eulernest
     container:
       image: git.eulernest.eu/conformallab/ci-cpp:latest
+      options: "--memory=800m --memory-swap=1200m"
     continue-on-error: true   # never block the merge
     steps:
       - uses: actions/checkout@v4

.gitea/workflows/markdown-links.yml

@@ -14,7 +14,6 @@ on:
   push:
     branches:
       - main
-      - dev
       - "claude/**"
       - "feature/**"
       - "review/**"
@@ -27,11 +26,11 @@ jobs:
     if: |
       github.event_name == 'schedule' ||
       github.event_name == 'workflow_dispatch' ||
-      contains(github.event.head_commit.message, '/links') ||
-      contains(github.event.head_commit.message, '/ci-all')
+      contains(github.event.head_commit.message, '/links')
     runs-on: eulernest
     container:
       image: git.eulernest.eu/conformallab/ci-cpp:latest
+      options: "--memory=400m --memory-swap=600m"
 
     steps:
       - uses: actions/checkout@v4

CLAUDE.md

@@ -270,6 +270,11 @@ Three jobs in `.gitea/workflows/cpp-tests.yml`:
 | `doc-build` | *(none)* | Doxygen | `/docs` in commit message or `workflow_dispatch` | **active** |
 | `markdown-links` | *(none)* | python3 | `/links` in commit message, weekly cron, `workflow_dispatch` | **active** |
 
+> **⚠ Pi runner limit:** use **one keyword per commit**. The Pi (3-4 GB RAM) cannot run
+> multiple Docker containers simultaneously. `/ci-all` was removed for this reason.
+> Typical workflow: one commit with `/test-cgal`, then if green a separate commit with
+> `/quality-gates`.
+
 Runner: `eulernest` — self-hosted Raspberry Pi, ARM64, Ubuntu 22.04. Docker image: `git.eulernest.eu/conformallab/ci-cpp:latest`. `test-cgal` and `quality-gates` both need `test-fast` to pass first (`needs: test-fast`).
 
 `quality-gates` runs four required structural gates: `license-headers.sh`, `cgal-conventions.py`, `codespell.sh`, `shellcheck.sh --strict`.  Seven more gates (clang-format, cmake-format, cppcheck, sanitizers, clang-tidy, multi-compiler, reproducible-build) are local-only — see `scripts/quality/README.md`.
@@ -370,9 +375,12 @@ Full line-by-line audit of all math-critical headers against `de.varylab.discret
 
 | Question | Document |
 |---|---|
-| Phases 1–10 with status and sub-tasks | `doc/roadmap/phases.md` |
+| Phases 1–13 with status, effort, and sub-tasks | `doc/roadmap/phases.md` |
+| Operational truth of which Java maths is in C++ today | `doc/roadmap/porting-status.md` |
 | Which Java classes are ported, which are planned, which are skipped? | `doc/roadmap/java-parity.md` |
 | New research items (beyond Java) — citations, acceptance criteria | `doc/roadmap/research-track.md` |
+| Phase orchestration — model assignments, priority, phase-session mapping | `doc/roadmap/phase-orchestration.md` |
+| Ready-to-paste session prompts for upcoming phases (P1–P4) | `doc/roadmap/session-prompts.md` |
 
 ### Tutorials & onboarding
 
@@ -392,6 +400,8 @@ Full line-by-line audit of all math-critical headers against `de.varylab.discret
 | Internal meeting agenda | `doc/reviewer/agenda.md` |
 | Reviewer landing index | `doc/reviewer/README.md` |
 | Hand-curated reviewer landing page (HTML, source-of-truth for codeberg pages) | `doc/reviewer/hub.html` |
+| Audit orchestration — model assignments, session sequence, status tracker | `doc/reviewer/finding-orchestration.md` |
+| Ready-to-paste session prompts for pending audit sessions (S3–S6) | `doc/reviewer/session-prompts.md` |
 
 The published hub lives at https://tmoussa.codeberg.page/ConformalLabpp/ (Doxygen index at `/doxygen.html`).  See the "Codeberg pages" quirk below for how it is republished.
 
@@ -414,6 +424,8 @@ Recommended loops when working in this repo. Prefer the cheapest gate that catch
 - **Before any commit**: run the four required gates locally — they mirror CI exactly and are seconds-cheap: `bash scripts/quality/license-headers.sh`, `python3 scripts/quality/cgal-conventions.py`, `bash scripts/quality/codespell.sh`, `bash scripts/quality/shellcheck.sh --strict`.
 - **Before tagging a release**: also run the two now-un-gated structural gates (test-cgal is disabled in CI): `BUILD_DIR=build bash scripts/check-test-counts.sh` and `bash scripts/try_it.sh`.  Update `CHANGELOG.md`, `CITATION.cff`, and the `doc/api/tests.md` counts (single source of truth).
 - **Touching public-API headers**: rebuild Doxygen (`cmake --build build --target doc`) and re-check coverage (`bash scripts/doxygen-coverage.sh --threshold 100`); regenerate `doc/api/headers.md` via `python3 scripts/gen-headers-md.py` (or `bash scripts/regen-docs.sh`).
+- **Starting work on an audit finding**: open `doc/reviewer/finding-orchestration.md`, pick the next ⬜ pending session, copy its prompt from `doc/reviewer/session-prompts.md`, set the named model, and go. **S3 is next** (H3/H4/H5/V5/V6, Sonnet → Opus review).
+- **Starting work on a roadmap phase**: open `doc/roadmap/phase-orchestration.md`, pick the next ⬜ pending phase-session, copy its prompt from `doc/roadmap/session-prompts.md`, set the named model, and go. **P1 is next** (9g.1 + 9h.1 + 9h.2 + 9d.3, Haiku → Opus review).
 - **Landing to `main`** (origin is protected): branch → push to `origin` → open PR via `gh`/Gitea API → merge via API → also push `codeberg/main` directly → keep both remotes in sync.
 - **Republishing the reviewer hub**: see the Codeberg `pages` quirk below — manual force-push of an orphan branch; verify the live URL with a cache-bust query.
 - **Delegation**: this repo's heavy builds are slow on the ARM64 runner — when a task is genuinely parallelisable and independent, consider a background agent; otherwise handle inline.  Always verify an agent's actual diff, not just its summary.

README.md

@@ -14,7 +14,7 @@ Algorithmic foundation:
 > DOI: [10.14279/depositonce-5415](https://depositonce.tu-berlin.de/items/8e2988b2-d991-45b5-aad5-9fb7988f3b2f) · CC BY-SA 4.0 ·
 > [Java original](https://github.com/varylab/conformallab) · [sechel.de](https://sechel.de/)
 
-**Status:** v0.9.0 — Phases 1–9a complete, Phase 8b-Lite CGAL API surface. Newton solvers for **five** DCE models (Euclidean / Spherical / HyperIdeal / CP-Euclidean / Inversive-Distance), priority-BFS layout in ℝ²/S²/Poincaré disk, Gauss–Bonnet, tree-cotree cut graph, Möbius holonomy, period matrix (genus 1), fundamental domain, halfedge_uv texture atlas, JSON/XML serialisation, CLI app.  Full test suite passing, 0 skipped — see [`doc/api/tests.md`](doc/api/tests.md) for the per-suite breakdown.
+**Status:** v0.10.0 — Phases 1–9b complete, Phase 8b-Lite CGAL API surface. Newton solvers for **five** DCE models (Euclidean / Spherical / HyperIdeal / CP-Euclidean / Inversive-Distance), priority-BFS layout in ℝ²/S²/Poincaré disk, Gauss–Bonnet, tree-cotree cut graph, Möbius holonomy, period matrix (genus 1), fundamental domain, halfedge_uv texture atlas, JSON/XML serialisation, CLI app.  277 tests passing, 0 skipped — see [`doc/api/tests.md`](doc/api/tests.md) for the per-suite breakdown.
 
 ---
 
@@ -75,6 +75,12 @@ cmake -S code -B build -DWITH_CGAL_TESTS=ON -DCONFORMALLAB_DEV_BUILD=ON
 # build at the cost of 5–15× slower test RUN.  Neutral on macOS.
 cmake -S code -B build -DWITH_CGAL_TESTS=ON -DCONFORMALLAB_FAST_TEST_BUILD=ON
 
+# Low-memory build: -O0, no PCH, unity batch 1, --no-keep-memory linker.
+# Drops cc1plus peak from ~700 MB to ~150 MB per TU.  Use on Raspberry Pi
+# or any runner with ≤ 4 GB RAM.  Always build with -j1.
+cmake -S code -B build -DWITH_CGAL_TESTS=ON -DCONFORMALLAB_LOW_MEMORY_BUILD=ON
+cmake --build build --target conformallab_cgal_tests -j1
+
 # Pristine measurement: disable both performance levers, e.g. for
 # scripts/quality/coverage.sh that needs every TU compiled fresh.
 cmake -S code -B build -DWITH_CGAL_TESTS=ON \
@@ -103,14 +109,12 @@ ConformalMesh mesh = load_mesh("input.off");
 EuclideanMaps maps = setup_euclidean_maps(mesh);
 compute_euclidean_lambda0_from_mesh(mesh, maps);
 
-// Assign DOFs — pin first vertex (gauge fix)
-auto vit = mesh.vertices().begin();
-maps.v_idx[*vit++] = -1;
-int idx = 0;
-for (; vit != mesh.vertices().end(); ++vit) maps.v_idx[*vit] = idx++;
+// Assign DOFs — pin first vertex as gauge fix, index the rest 0..n-1
+auto gauge = *mesh.vertices().begin();
+int n = assign_euclidean_vertex_dof_indices(mesh, maps, gauge);
 
 // Natural equilibrium target: x* = 0 by construction
-std::vector<double> x0(idx, 0.0);
+std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
 auto G0 = euclidean_gradient(mesh, x0, maps);
 for (auto v : mesh.vertices())
     if (maps.v_idx[v] >= 0) maps.theta_v[v] -= G0[maps.v_idx[v]];

code/examples/CMakeLists.txt

@@ -42,6 +42,22 @@ target_include_directories(example_layout SYSTEM PRIVATE ${EXAMPLE_INCLUDES})
 target_include_directories(example_layout PRIVATE        ${EXAMPLE_PRIVATE_INCLUDES})
 target_compile_definitions(example_layout PRIVATE        ${EXAMPLE_DEFS})
 
+# ── example_flatten (primary use case: real conformal flattening) ─────────────
+# Demonstrates Θ_v = 2π → non-trivial conformal map; contrast with
+# example_euclidean which uses "natural theta" (u_v ≈ 0, testing trick).
+add_executable(example_flatten example_flatten.cpp)
+target_include_directories(example_flatten SYSTEM PRIVATE ${EXAMPLE_INCLUDES})
+target_include_directories(example_flatten PRIVATE        ${EXAMPLE_PRIVATE_INCLUDES})
+target_compile_definitions(example_flatten PRIVATE        ${EXAMPLE_DEFS})
+
+# ── example_cgal_api (CGAL public API: one-call interface) ────────────────────
+# Demonstrates CGAL::discrete_conformal_map_euclidean from
+# <CGAL/Discrete_conformal_map.h>; contrast with example_flatten (internal API).
+add_executable(example_cgal_api example_cgal_api.cpp)
+target_include_directories(example_cgal_api SYSTEM PRIVATE ${EXAMPLE_INCLUDES})
+target_include_directories(example_cgal_api PRIVATE        ${EXAMPLE_PRIVATE_INCLUDES})
+target_compile_definitions(example_cgal_api PRIVATE        ${EXAMPLE_DEFS})
+
 # ── example_viewer (requires WITH_VIEWER) ─────────────────────────────────────
 if(WITH_VIEWER)
     add_executable(example_viewer example_viewer.cpp)

code/examples/example_cgal_api.cpp

@@ -0,0 +1,171 @@
+// example_cgal_api.cpp
+//
+// conformallab++ — CGAL public API example
+//
+// This example demonstrates the HIGH-LEVEL CGAL API defined in
+// <CGAL/Discrete_conformal_map.h>.  It is the recommended entry point for
+// users who want a simple one-call interface without managing Maps bundles,
+// DOF assignment, or Newton solver details.
+//
+// Contrast with example_euclidean.cpp / example_flatten.cpp which use the
+// INTERNAL API (setup_euclidean_maps + newton_euclidean + euclidean_layout).
+//
+// ┌─────────────────────────────────────────────────────────────────────┐
+// │  When to use the CGAL API vs the internal API                       │
+// │                                                                     │
+// │  CGAL API (this file):                                              │
+// │    • One call, sensible defaults                                    │
+// │    • Named parameters for tuning (tolerance, cone angles, …)       │
+// │    • Result type: Conformal_map_result<FT>  (u_per_vertex indexed  │
+// │      by raw vertex index)                                           │
+// │    • Layout via CGAL::euclidean_layout (Conformal_layout.h)        │
+// │                                                                     │
+// │  Internal API (example_flatten.cpp, example_layout.cpp):           │
+// │    • Full control over every pipeline step                         │
+// │    • Access to holonomy, period matrix, cut graph                  │
+// │    • Result type: NewtonResult  (x indexed by DOF index)           │
+// │    • Required for closed surfaces (genus ≥ 1)                     │
+// └─────────────────────────────────────────────────────────────────────┘
+//
+// Build (requires -DWITH_CGAL=ON):
+//   cmake -S code -B build -DWITH_CGAL=ON
+//   cmake --build build --target example_cgal_api
+//
+// Run:
+//   ./build/examples/example_cgal_api                        # built-in mesh
+//   ./build/examples/example_cgal_api code/data/obj/cathead.obj flat.off
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Discrete_conformal_map.h>   // CGAL public API
+#include <CGAL/Conformal_layout.h>         // CGAL layout wrapper
+
+// For loading meshes and saving results we still use the internal helpers.
+#include "mesh_io.hpp"
+#include "mesh_builder.hpp"
+#include "layout.hpp"
+
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <algorithm>
+
+int main(int argc, char* argv[])
+{
+    using Kernel = CGAL::Simple_cartesian<double>;
+    using Mesh   = CGAL::Surface_mesh<Kernel::Point_3>;
+
+    // ── Step 1: obtain mesh ───────────────────────────────────────────────────
+    Mesh mesh;
+    std::string input_path  = (argc > 1) ? argv[1] : "";
+    std::string output_path = (argc > 2) ? argv[2] : "/tmp/conformallab_cgal_api_out.off";
+
+    if (input_path.empty()) {
+        std::cout << "[example_cgal_api] No input given.  Using make_quad_strip().\n"
+                  << "  For a non-trivial result, supply a 3-D mesh:\n"
+                  << "    ./example_cgal_api code/data/obj/cathead.obj\n\n";
+        mesh = conformallab::make_quad_strip();
+    } else {
+        std::cout << "[example_cgal_api] Loading mesh: " << input_path << "\n";
+        try { mesh = conformallab::load_mesh(input_path); }
+        catch (const std::exception& e) {
+            std::cerr << "Error loading mesh: " << e.what() << "\n";
+            return 1;
+        }
+    }
+
+    std::cout << "[example_cgal_api] Mesh: "
+              << mesh.number_of_vertices() << " vertices, "
+              << mesh.number_of_faces()    << " faces.\n";
+
+    // ── Step 2: CGAL API call ─────────────────────────────────────────────────
+    //
+    // Default invocation — one function call, no explicit target angles:
+    //   • No vertex_curvature_map supplied → "natural theta" default:
+    //     Θ_v is set to the ACTUAL angle sum at x = 0.  This makes x* = 0
+    //     the equilibrium, so Newton converges in 0 iterations with u_v = 0.
+    //
+    // ⚠  Natural theta is a TESTING CONVENTION, not a conformal flattening.
+    //    The output u_v = 0 means "no deformation" — the map is identity.
+    //    For real conformal flattening use:
+    //      • example_flatten.cpp (internal API, recommended for open meshes)
+    //      • Supply vertex_curvature_map(theta_map) with Θ_v = 2π for a
+    //        closed mesh (must satisfy Gauss-Bonnet; see below)
+    //
+    // The function sets up maps, computes λ°, assigns DOFs, runs Newton,
+    // and returns the converged result.
+    auto result = CGAL::discrete_conformal_map_euclidean(mesh);
+
+    std::cout << "[example_cgal_api] CGAL API result (natural theta — identity map):\n"
+              << "  converged  = " << std::boolalpha << result.converged << "\n"
+              << "  iterations = " << result.iterations
+              << "  (0 = natural theta, trivially at equilibrium)\n"
+              << "  |grad|_inf = " << std::scientific << std::setprecision(2)
+              << result.gradient_norm << "\n";
+
+    // ── Step 3: inspect the result ────────────────────────────────────────────
+    //
+    // result.u_per_vertex is indexed by raw vertex index (v.idx()), NOT by
+    // DOF index.  Length = num_vertices(mesh).  Pinned vertices have u = 0.
+    //
+    // This differs from the internal API (NewtonResult.x) which is indexed
+    // by DOF index (v_idx[v]).  The CGAL API handles the mapping internally.
+    double u_min = 0.0, u_max = 0.0;
+    for (auto v : mesh.vertices()) {
+        double u = result.u_per_vertex[static_cast<std::size_t>(v.idx())];
+        u_min = std::min(u_min, u);
+        u_max = std::max(u_max, u);
+    }
+    std::cout << "  u_v range  = [" << std::fixed << std::setprecision(4)
+              << u_min << ", " << u_max << "]\n\n";
+
+    // Print a few per-vertex values
+    std::cout << "[example_cgal_api] First 6 per-vertex scale factors u_v:\n";
+    int shown = 0;
+    for (auto v : mesh.vertices()) {
+        if (shown++ >= 6) break;
+        double u = result.u_per_vertex[static_cast<std::size_t>(v.idx())];
+        std::cout << "  v" << v.idx() << "  u = " << std::fixed
+                  << std::setprecision(6) << u << "\n";
+    }
+
+    // ── Step 4: tuning via named parameters ───────────────────────────────────
+    //
+    // The CGAL API accepts named parameters for fine-grained control.
+    // Example: tighter tolerance and more iterations.
+    //
+    //   auto result2 = CGAL::discrete_conformal_map_euclidean(
+    //       mesh,
+    //       CGAL::parameters::gradient_tolerance(1e-12)
+    //                        .max_iterations(500));
+    //
+    // Example: supply explicit cone angles (requires Gauss-Bonnet to hold):
+    //
+    //   auto angle_map = mesh.add_property_map<Mesh::Vertex_index, double>(
+    //       "my:angles", 2.0 * M_PI).first;
+    //   // … set angle_map[v] for cone vertices …
+    //   auto result3 = CGAL::discrete_conformal_map_euclidean(
+    //       mesh,
+    //       CGAL::parameters::vertex_curvature_map(angle_map));
+    //
+    // See <CGAL/Discrete_conformal_map.h> for all named parameters.
+
+    // ── Step 5: compute layout via the CGAL layout wrapper ────────────────────
+    //
+    // The CGAL API does not return a layout directly; call CGAL::euclidean_layout
+    // (Conformal_layout.h) with the converged DOF vector and the internal maps.
+    // Note: to access the DOF vector and maps from a CGAL-API call, use the
+    // result.x field (if exposed) or call the internal API directly.
+    //
+    // For Phase 8b-Lite, the cleanest way to get a layout after the CGAL call
+    // is to use the internal API (example_flatten.cpp) — the CGAL result carries
+    // u_per_vertex for inspection but the full pipeline (layout, holonomy, period
+    // matrix) still requires the internal Maps bundle.
+
+    if (result.converged)
+        std::cout << "\n[example_cgal_api] ✓ Conformal map computed successfully.\n"
+                  << "  For UV layout output, see example_flatten.cpp (internal API)\n"
+                  << "  which provides direct access to euclidean_layout().\n";
+
+    return result.converged ? 0 : 1;
+}

code/examples/example_euclidean.cpp

@@ -59,20 +59,25 @@ int main(int argc, char* argv[])
     compute_euclidean_lambda0_from_mesh(mesh, maps);
 
     // ── Step 3: pin the first vertex (gauge fix) ──────────────────────────
-    auto vit = mesh.vertices().begin();
-    Vertex_index v_pinned = *vit++;
-    maps.v_idx[v_pinned] = -1;   // pinned: u[v_pinned] = 0 (fixed)
-    int idx = 0;
-    for (; vit != mesh.vertices().end(); ++vit)
-        maps.v_idx[*vit] = idx++;
-    const int n = idx;
+    // The gauge-vertex overload pins the chosen vertex (v_idx = -1) and
+    // assigns sequential indices 0..n-1 to the rest in a single call.
+    Vertex_index v_pinned = *mesh.vertices().begin();
+    const int n = assign_euclidean_vertex_dof_indices(mesh, maps, v_pinned);
 
     std::cout << "[example_euclidean] DOFs: " << n << " (1 vertex pinned).\n";
 
     // ── Step 4: natural equilibrium — set theta_v = actual angle sum at x=0 ─
-    //    After this step x* = 0 is the equilibrium (no deformation).
-    //    In a real application you would set theta_v = desired angle (e.g. 2π
-    //    for flat disks, or the cone angles for a cone metric).
+    //
+    // ⚠  TESTING CONVENTION — NOT A REAL CONFORMAL MAP
+    //
+    // "Natural theta" sets Θ_v = actual angle sum at x = 0, so x* = 0 is the
+    // equilibrium by construction.  The solver converges in 0–1 iterations and
+    // u_v ≈ 0 everywhere.  This is useful for testing the solver pipeline but
+    // produces NO conformal deformation.
+    //
+    // For a REAL conformal flattening (the primary use case):
+    //   → see example_flatten.cpp
+    //   which sets Θ_v = 2π (flat interior target) and produces non-trivial u_v.
     {
         std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
         auto G0 = euclidean_gradient(mesh, x0, maps);

code/examples/example_flatten.cpp

@@ -0,0 +1,206 @@
+// example_flatten.cpp
+//
+// conformallab++ — Real conformal flattening example
+//
+// This example demonstrates the PRIMARY USE CASE of the library:
+// conformally flatten a 3-D surface mesh to the plane with minimal
+// angle distortion.  Every interior vertex is assigned a target cone
+// angle of 2π (a regular flat vertex); the solver finds the unique
+// conformal factor u_v that realises this target.
+//
+// Contrast with the other examples (example_euclidean, example_layout)
+// which use the "natural theta" testing trick that produces x* = 0 —
+// a valid solver test but NOT a conformal flattening.
+//
+// ┌─────────────────────────────────────────────────────────────────────┐
+// │  Pipeline for conformal flattening of an OPEN mesh (disk topology) │
+// │                                                                     │
+// │  1. Load mesh                                                       │
+// │  2. Setup maps + compute λ° from geometry                          │
+// │  3. Pin all boundary vertices (they define the boundary of the UV) │
+// │     Set Θ_v = 2π for all interior vertices (flat target)           │
+// │     ← NO Gauss–Bonnet check needed for open meshes                 │
+// │  4. Solve Newton                                                    │
+// │  5. Compute planar layout — the conformal UV parameterisation       │
+// │  6. Save UV-mapped OFF + report distortion                         │
+// └─────────────────────────────────────────────────────────────────────┘
+//
+// Build (requires -DWITH_CGAL=ON):
+//   cmake -S code -B build -DWITH_CGAL=ON
+//   cmake --build build --target example_flatten
+//
+// Run:
+//   # With a real 3-D surface mesh (open, disk topology):
+//   ./build/examples/example_flatten code/data/obj/cathead.obj flat.off
+//
+//   # Without arguments: uses a built-in synthetic open mesh (6 vertices)
+//   ./build/examples/example_flatten
+
+#include "conformal_mesh.hpp"
+#include "mesh_builder.hpp"
+#include "mesh_io.hpp"
+#include "euclidean_functional.hpp"
+#include "gauss_bonnet.hpp"
+#include "newton_solver.hpp"
+#include "layout.hpp"
+#include "constants.hpp"
+#include <CGAL/boost/graph/iterator.h>
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <algorithm>
+
+using namespace conformallab;
+
+int main(int argc, char* argv[])
+{
+    // ── Step 1: load or synthesise a mesh ────────────────────────────────────
+    ConformalMesh mesh;
+    std::string input_path  = (argc > 1) ? argv[1] : "";
+    std::string output_path = (argc > 2) ? argv[2] : "/tmp/conformallab_flatten_out.off";
+
+    if (input_path.empty()) {
+        // Fallback: load cathead.obj from the standard data location if
+        // it exists alongside the executable; otherwise use a synthetic mesh.
+        // For a meaningful non-trivial flattening, supply a real 3-D mesh:
+        //   ./example_flatten code/data/obj/cathead.obj
+        std::cout << "[example_flatten] No input given.  Using make_quad_strip() "
+                     "(flat synthetic mesh — u_v will be near-zero).\n"
+                  << "  For a non-trivial flattening, provide a 3-D mesh:\n"
+                  << "    ./example_flatten code/data/obj/cathead.obj\n\n";
+        mesh = make_quad_strip();
+    } else {
+        std::cout << "[example_flatten] Loading mesh: " << input_path << "\n";
+        try { mesh = load_mesh(input_path); }
+        catch (const std::exception& e) {
+            std::cerr << "Error loading mesh: " << e.what() << "\n";
+            return 1;
+        }
+    }
+
+    const int V = static_cast<int>(mesh.number_of_vertices());
+    const int F = static_cast<int>(mesh.number_of_faces());
+    std::cout << "[example_flatten] Mesh: " << V << " vertices, " << F << " faces\n";
+
+    // ── Step 2: setup maps + compute initial edge lengths ─────────────────────
+    auto maps = setup_euclidean_maps(mesh);
+    compute_euclidean_lambda0_from_mesh(mesh, maps);
+
+    // ── Step 3: DOF assignment — pin boundary, free interior ──────────────────
+    //
+    // For an OPEN mesh (disk topology):
+    //   • Boundary vertices are pinned (v_idx = -1): they define the
+    //     boundary of the UV domain and are not optimised.
+    //   • Interior vertices get a DOF (v_idx ≥ 0) and target Θ_v = 2π.
+    //     This means "make every interior point look like a flat plane vertex".
+    //
+    // For a CLOSED mesh (e.g. a sphere or torus), use the Euclidean pipeline
+    // on a cut-open mesh (see example_layout.cpp + cut_graph.hpp), or use the
+    // spherical / hyper-ideal functional instead.
+    //
+    // ⚠  This is the KEY difference from example_euclidean.cpp:
+    //    There, "natural theta" sets Θ_v = actual angle sum → x* = 0 (trivial).
+    //    Here, Θ_v = 2π → the solver finds the REAL conformal map.
+    int idx = 0;
+    int n_boundary = 0, n_interior = 0;
+
+    for (auto v : mesh.vertices()) {
+        // CGAL: a vertex is on the boundary iff it has an incident border halfedge.
+        bool is_bnd = false;
+        for (auto h : CGAL::halfedges_around_target(v, mesh))
+            if (mesh.is_border(h)) { is_bnd = true; break; }
+
+        if (is_bnd) {
+            maps.v_idx[v] = -1;      // pinned — boundary defines the UV border
+            ++n_boundary;
+        } else {
+            maps.theta_v[v] = TWO_PI;  // flat interior target — the actual goal
+            maps.v_idx[v]   = idx++;
+            ++n_interior;
+        }
+    }
+
+    std::cout << "[example_flatten] Boundary vertices (pinned): " << n_boundary
+              << "  Interior (free DOFs): " << n_interior << "\n";
+
+    if (n_interior == 0) {
+        std::cerr << "[example_flatten] No interior vertices — mesh has no free DOFs.\n"
+                  << "  Use a mesh with interior vertices (e.g. cathead.obj).\n";
+        return 1;
+    }
+
+    // For open meshes the Gauss–Bonnet identity holds in a different form and
+    // does NOT need to be checked before calling newton_euclidean.  The solver
+    // converges as long as at least one interior vertex exists.
+    // (For CLOSED meshes: call enforce_gauss_bonnet(mesh, maps) here.)
+
+    // ── Step 4: Newton ────────────────────────────────────────────────────────
+    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
+
+    std::cout << "[example_flatten] Running Newton (tol = 1e-9)…\n";
+    auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-9);
+
+    if (res.converged)
+        std::cout << "[example_flatten] Converged in " << res.iterations
+                  << " iterations  ||G||_inf = " << std::scientific
+                  << std::setprecision(2) << res.grad_inf_norm << "\n";
+    else
+        std::cout << "[example_flatten] WARNING: did not converge after "
+                  << res.iterations << " iterations  ||G||_inf = "
+                  << res.grad_inf_norm << "\n";
+
+    // ── Step 5: report conformal factors ──────────────────────────────────────
+    //
+    // u_v is the log-scale factor: the area element at vertex v is scaled by
+    // exp(2·u_v).  For a non-trivial mesh the values are non-zero.
+    double u_min = 0.0, u_max = 0.0;
+    for (auto v : mesh.vertices()) {
+        int iv = maps.v_idx[v];
+        if (iv >= 0) {
+            double u = res.x[static_cast<std::size_t>(iv)];
+            u_min = std::min(u_min, u);
+            u_max = std::max(u_max, u);
+        }
+    }
+    std::cout << "[example_flatten] Conformal factors u_v: "
+              << "min = " << std::fixed << std::setprecision(4) << u_min
+              << "  max = " << u_max << "\n";
+
+    if (std::abs(u_max - u_min) < 1e-6)
+        std::cout << "[example_flatten] Note: u_v ≈ 0 everywhere — the mesh is "
+                     "already flat in the\n"
+                     "  Euclidean sense (e.g. a planar mesh).  Use a 3-D surface "
+                     "for non-trivial output.\n";
+    else
+        std::cout << "[example_flatten] Non-trivial u_v range = "
+                  << (u_max - u_min) << " — real conformal deformation computed.\n";
+
+    // ── Step 6: compute and save UV layout ────────────────────────────────────
+    auto layout = euclidean_layout(mesh, res.x, maps);
+
+    if (layout.success) {
+        // Find the bounding box of the UV coords
+        double xmin = 1e30, xmax = -1e30, ymin = 1e30, ymax = -1e30;
+        for (auto v : mesh.vertices()) {
+            if (static_cast<std::size_t>(v.idx()) < layout.uv.size()) {
+                const auto& p = layout.uv[static_cast<std::size_t>(v.idx())];
+                xmin = std::min(xmin, p.x()); xmax = std::max(xmax, p.x());
+                ymin = std::min(ymin, p.y()); ymax = std::max(ymax, p.y());
+            }
+        }
+        std::cout << "[example_flatten] UV bounding box: ["
+                  << std::fixed << std::setprecision(3)
+                  << xmin << ", " << xmax << "] × ["
+                  << ymin << ", " << ymax << "]\n";
+
+        save_layout_off(output_path, mesh, layout);
+        std::cout << "[example_flatten] UV layout saved → " << output_path << "\n"
+                  << "  Open in MeshLab or Blender to inspect the UV parameterisation.\n";
+    } else {
+        std::cerr << "[example_flatten] Layout failed.\n";
+    }
+
+    return (res.converged && layout.success) ? 0 : 1;
+}

code/examples/example_layout.cpp

@@ -37,6 +37,11 @@
 using namespace conformallab;
 
 // ── Helper: natural target angles so x* = 0 is the equilibrium ───────────────
+//
+// ⚠  TESTING CONVENTION — NOT A REAL CONFORMAL MAP
+// Natural theta sets Θ_v = actual angle sum at x=0, making x* = 0 trivially
+// the equilibrium (u_v ≈ 0, no deformation).  For real conformal flattening,
+// see example_flatten.cpp which uses Θ_v = 2π (flat interior target).
 static void set_natural_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n)
 {
     std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
@@ -49,14 +54,13 @@ static void set_natural_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n)
 }
 
 // ── Helper: pin vertex 0, assign DOF indices 0..n-1 to the rest ──────────────
+// Uses the gauge-vertex overload introduced in external-audit Finding-D:
+// assign_euclidean_vertex_dof_indices(mesh, maps, gauge) pins the chosen
+// vertex and assigns sequential indices in a single pass.
 static int pin_first(ConformalMesh& mesh, EuclideanMaps& maps)
 {
-    auto vit = mesh.vertices().begin();
-    maps.v_idx[*vit++] = -1;   // pinned
-    int idx = 0;
-    for (; vit != mesh.vertices().end(); ++vit)
-        maps.v_idx[*vit] = idx++;
-    return idx;
+    return assign_euclidean_vertex_dof_indices(
+        mesh, maps, *mesh.vertices().begin());
 }
 
 // ── Main ─────────────────────────────────────────────────────────────────────
@@ -124,7 +128,9 @@ int main(int argc, char* argv[])
     if (layout.success) {
         std::cout << "UV coordinates:\n";
         for (auto v : mesh.vertices()) {
-            auto& p = layout.uv[v.idx()];
+            // layout.uv is indexed by v.idx() (raw integer vertex index).
+            // Valid as long as no vertices were removed/compacted after loading.
+            auto& p = layout.uv[static_cast<std::size_t>(v.idx())];
             std::cout << "  v" << v.idx()
                       << ": (" << std::fixed << std::setprecision(6)
                       << p.x() << ", " << p.y() << ")\n";

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

@@ -7,13 +7,17 @@
 \file CGAL/Discrete_conformal_map.h
 \ingroup PkgConformalMapRef
 
-User-facing entry point for the Discrete_conformal_map package.
+User-facing entry points for the Discrete_conformal_map package (Phase 8b-Lite).
 
-This header provides a single function — `discrete_conformal_map_euclidean`
-— that computes a Euclidean discrete-conformal flattening of an open or
-closed triangle mesh.  Spherical and hyperbolic variants are scheduled
-for Phase 8b.2 once the Euclidean pattern is validated by Phase 9a
-(Inversive-Distance functional).
+This header provides three functions covering all three DCE geometries:
+
+- `CGAL::discrete_conformal_map_euclidean`   — flat conformal map (ℝ²), open or closed mesh
+- `CGAL::discrete_conformal_map_spherical`   — spherical uniformisation (S²), genus-0 mesh
+- `CGAL::discrete_conformal_map_hyper_ideal` — hyperbolic conformal map (H²), genus ≥ 1
+
+For circle-packing models see the companion headers:
+- `<CGAL/Discrete_circle_packing.h>`       — `discrete_circle_packing_euclidean` (CP-Euclidean)
+- `<CGAL/Discrete_inversive_distance.h>`   — `discrete_inversive_distance_map` (Luo 2004)
 
 \section Example Simplest usage
 
@@ -74,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
 
@@ -96,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);
 };
 
 
@@ -268,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) ──────────────────
     //
@@ -406,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(
@@ -455,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);
 };
 
 /*!
@@ -534,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 {
@@ -85,12 +87,26 @@ inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31)
     if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0)
         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};
+    const double l123 = l12 + l23 + l31;
 
-    // denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area
-    const double denom2 = 2.0 * std::sqrt(denom2_sq);
+    // Triangle area via Kahan's numerically stable formula.  Sort the side
+    // 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;
-    if (v2b && b2  < 0.0) b2 = 0.01;
-    if (v3b && b3  < 0.0) b3 = 0.01;
+    b1 = clamp_hyper_ideal_scale(b1, v1b, clamp);
+    b2 = clamp_hyper_ideal_scale(b2, v2b, clamp);
+    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;
-    if (fa.v2b && fa.b2  < 0.0)           fa.b2  = 0.01;
-    if (fa.v3b && fa.b3  < 0.0)           fa.b3  = 0.01;
+    fa.b1 = clamp_hyper_ideal_scale(fa.b1, fa.v1b, clamp);
+    fa.b2 = clamp_hyper_ideal_scale(fa.b2, fa.v2b, clamp);
+    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 Gm = evaluate_hyper_ideal(mesh, xm, 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, /*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/layout.hpp

@@ -138,18 +138,38 @@ struct MobiusMap {
 /// per-vertex UV coordinates plus a per-half-edge UV atlas for seamed
 /// textures.
 struct Layout2D {
-    /// uv[v.idx()] — primary 2-D position (first / shallowest-BFS-depth visit).
+    /// Primary 2-D position per vertex (first / shallowest-BFS-depth visit).
+    ///
+    /// **Indexing:** `uv[v.idx()]` — indexed by the raw integer vertex index.
+    /// **Length:** `mesh.number_of_vertices()`.
+    ///
+    /// \pre  No vertices have been removed from `mesh` after loading (i.e.
+    ///       `mesh.is_valid()` and no compaction was performed).  On a fresh
+    ///       `Surface_mesh` loaded from file, `v.idx()` is always in
+    ///       `[0, number_of_vertices())` and contiguous.  If vertices were
+    ///       deleted and `mesh.collect_garbage()` was called, re-run the
+    ///       layout — indices will have shifted.
+    ///
+    /// Access pattern:
+    /// ```cpp
+    /// for (auto v : mesh.vertices())
+    ///     Eigen::Vector2d uv_v = layout.uv[v.idx()];
+    /// ```
     std::vector<Eigen::Vector2d> uv;
 
-    /// halfedge_uv[h.idx()] — UV of source(h) as seen from face(h).
+    /// UV of `source(h)` as seen from `face(h)`, indexed by `h.idx()`.
     ///
-    /// For interior (non-seam) halfedges: equals uv[source(h).idx()].
-    /// For seam halfedges:  carries the UV from the virtual unfolding across
-    ///   the cut (i.e. the trilaterated position that was NOT used as the
-    ///   primary uv).  This gives each face its own copy of a seam vertex,
-    ///   enabling a proper GPU texture atlas without vertex duplication.
+    /// **Indexing:** `halfedge_uv[h.idx()]` — raw integer halfedge index.
+    /// **Length:** `mesh.number_of_halfedges()`.  Same no-compaction precondition
+    /// as `uv` (see above).
     ///
-    /// Size = mesh.number_of_halfedges().  Border halfedges = (0,0).
+    /// For interior (non-seam) halfedges: equals `uv[source(h).idx()]`.
+    /// For seam halfedges: carries the UV from the virtual unfolding across the
+    ///   cut — the trilaterated position that was NOT used as the primary `uv`.
+    ///   This gives each face its own copy of a seam vertex, enabling a proper
+    ///   per-halfedge GPU texture atlas without vertex duplication.
+    ///
+    /// Border halfedges (outer face) hold `(0, 0)`.
     std::vector<Eigen::Vector2d> halfedge_uv;
 
     bool success  = false; ///< `true` iff the BFS placed every vertex.

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;
-    return detail::solve_with_fallback(A, rhs, ok, fallback_used);
+    bool ok_local = false;
+    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,
@@ -144,9 +187,10 @@ inline std::vector<double> line_search(
     const Eigen::VectorXd&     d_sd,
     double                     norm0,
     GradFn&&                   grad_fn,
-    bool*                      improved     = nullptr,
-    int                        max_halvings = 20,
-    double                     c1           = 1e-4)
+    bool*                      improved      = nullptr,
+    std::vector<double>*       accepted_grad = nullptr,
+    int                        max_halvings  = 20,
+    double                     c1            = 1e-4)
 {
     const int    n        = static_cast<int>(x.size());
     const double norm0_sq = norm0 * norm0;
@@ -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;
-    const int n = static_cast<int>(x.size());
+    // Layout is loop-invariant — decide the Hessian variant once (B3): cyclic
+    // layout (edge DOFs present) → analytic cyclic Hessian, which covers the
+    // vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian omits.
+    bool has_edge_dof = false;
+    for (auto e : mesh.edges())
+        if (m.e_idx[e] >= 0) { has_edge_dof = true; break; }
 
-    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) ───────
-        double norm0 = G.norm();
-        Eigen::VectorXd d_sd = -(H * G);   // merit-function steepest descent
-        bool improved = true;
-        x = detail::line_search(x, dx, d_sd, norm0,
-            [&](const std::vector<double>& xnew) {
-                return euclidean_gradient(mesh, xnew, m);
-            }, &improved);
-        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;
+    return detail::newton_core(
+        std::move(x0),
+        [&](const std::vector<double>& xc) { return euclidean_gradient(mesh, xc, m); },
+        [&](const std::vector<double>& xc) {
+            return has_edge_dof ? euclidean_hessian_analytic(mesh, xc, m)
+                                : euclidean_hessian(mesh, xc, m);
+        },
+        /*concave=*/false, tol, max_iter);
 }
 
 // ── Spherical Newton solver ───────────────────────────────────────────────────
@@ -322,55 +489,14 @@ inline NewtonResult newton_spherical(
     double                     tol      = 1e-8,
     int                        max_iter = 200)
 {
-    std::vector<double> x = x0;
-    const int n = static_cast<int>(x.size());
-
-    NewtonResult res;
-    res.converged      = false;
-    res.iterations     = 0;
-    res.grad_inf_norm  = 0.0;
-
-    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;
+    // Concave energy: the Hessian H is NSD, so newton_core factors −H (PSD) and
+    // 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.
+    return detail::newton_core(
+        std::move(x0),
+        [&](const std::vector<double>& xc) { return spherical_gradient(mesh, xc, m); },
+        [&](const std::vector<double>& xc) { return spherical_hessian(mesh, xc, m); },
+        /*concave=*/true, tol, max_iter);
 }
 
 // ── 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;
-    const int n = static_cast<int>(x.size());
-
-    NewtonResult res;
-    res.converged      = false;
-    res.iterations     = 0;
-    res.grad_inf_norm  = 0.0;
-
-    for (int iter = 0; iter < max_iter; ++iter) {
-        // ── Gradient ──────────────────────────────────────────────────────────
-        auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false).gradient;
-        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;
+    // block-FD Hessian (~33–1166× faster than full-FD); `clamp` selects the
+    // vertex-scale floor mode (N3).  Convex energy → factor H directly.
+    return detail::newton_core(
+        std::move(x0),
+        [&](const std::vector<double>& xc) {
+            return evaluate_hyper_ideal(mesh, xc, m, /*energy=*/false, /*grad=*/true, clamp).gradient;
+        },
+        [&](const std::vector<double>& xc) {
+            return hyper_ideal_hessian_block_fd_sym(mesh, xc, m, hess_eps, clamp);
+        },
+        /*concave=*/false, tol, max_iter);
 }
 
 // ── 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;
-    const int n = static_cast<int>(x.size());
-
-    NewtonResult res;
-    res.converged      = false;
-    res.iterations     = 0;
-    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;
+    // Convex energy with an exact analytic Hessian (BPS-2015 2×2-per-edge).
+    return detail::newton_core(
+        std::move(x0),
+        [&](const std::vector<double>& xc) { return cp_euclidean_gradient(mesh, xc, m); },
+        [&](const std::vector<double>& xc) { return cp_euclidean_hessian(mesh, xc, m); },
+        /*concave=*/false, tol, max_iter);
 }
 
 // ── 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 analytic gradient (same pattern as the Phase 4a HyperIdeal solver).
-/// An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in
-/// `doc/roadmap/research-track.md` as Phase 9a.2-analytic.
+/// The Hessian is computed by **per-face block finite differences**
+/// (`inversive_distance_hessian_block_fd_sym`, same pattern as the HyperIdeal
+/// solver after Phase 9b) — O(F) face evaluations instead of the O(n·F) of the
+/// 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;
-    const int n = static_cast<int>(x.size());
-
-    NewtonResult res;
-    res.converged      = false;
-    res.iterations     = 0;
-    res.grad_inf_norm  = 0.0;
-
-    // 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;
+    // Convex energy; block-FD Hessian (~n/6× faster than full-FD).
+    return detail::newton_core(
+        std::move(x0),
+        [&](const std::vector<double>& xc) { return inversive_distance_gradient(mesh, xc, m); },
+        [&](const std::vector<double>& xc) {
+            return inversive_distance_hessian_block_fd_sym(mesh, xc, m, hess_eps);
+        },
+        /*concave=*/false, tol, max_iter);
 }
 
 } // namespace conformallab

code/include/serialization.hpp

@@ -93,31 +93,60 @@ inline std::vector<double> load_result_json(
 {
     using json = nlohmann::json;
     std::ifstream ifs(path);
-    if (!ifs) throw std::runtime_error("Cannot open: " + path);
-    json j; ifs >> j;
+    if (!ifs) throw std::runtime_error("conformallab: cannot open: " + path);
 
-    if (geom && j.contains("geometry"))
-        *geom = j["geometry"].get<std::string>();
-
-    std::vector<double> x = j.at("dof_vector").get<std::vector<double>>();
-
-    if (res) {
-        res->x             = x;
-        res->converged     = j["solver"]["converged"].get<bool>();
-        res->iterations    = j["solver"]["iterations"].get<int>();
-        res->grad_inf_norm = j["solver"]["grad_inf_norm"].get<double>();
+    // 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());
     }
 
-    if (layout2d && j.contains("layout") && j["layout"]["dim"] == 2) {
-        auto uv = j["layout"]["uv"];
-        layout2d->uv.resize(uv.size());
-        for (std::size_t i = 0; i < uv.size(); ++i)
-            layout2d->uv[i] = { uv[i][0].get<double>(), uv[i][1].get<double>() };
-        layout2d->has_seam = j["layout"].value("has_seam", false);
-        layout2d->success  = true;
-    }
+    try {
+        if (geom && j.contains("geometry"))
+            *geom = j["geometry"].get<std::string>();
 
-    return x;
+        // 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     = s.at("converged").get<bool>();
+            res->iterations    = s.at("iterations").get<int>();
+            res->grad_inf_norm = s.at("grad_inf_norm").get<double>();
+        }
+
+        if (layout2d && j.contains("layout") && j["layout"]["dim"] == 2) {
+            auto uv = j["layout"]["uv"];
+            layout2d->uv.resize(uv.size());
+            for (std::size_t i = 0; i < uv.size(); ++i)
+                layout2d->uv[i] = { uv[i][0].get<double>(), uv[i][1].get<double>() };
+            layout2d->has_seam = j["layout"].value("has_seam", false);
+            layout2d->success  = true;
+        }
+
+        return x;
+    } catch (const json::exception& e) {
+        throw std::runtime_error(
+            "conformallab: malformed JSON in " + path + ": " + e.what());
+    }
 }
 
 // ════════════════════════════════════════════════════════════════════════════
@@ -255,9 +284,23 @@ inline std::vector<double> load_result_xml(
         // Solver metadata
         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"));
-                res->grad_inf_norm = std::stod(detail_xml::xml_get_attr(line, "grad_inf_norm"));
+                res->converged = (detail_xml::xml_get_attr(line, "converged") == "true");
+                // V2: stoi/stod throw std::invalid_argument on empty or non-numeric
+                // 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/api/contracts.md

@@ -13,11 +13,11 @@ of all preceding units.
 | `load_mesh()` | Valid file path, supported format (OFF/OBJ/PLY) | Manifold, oriented, triangulated `ConformalMesh` |
 | `setup_*_maps()` | Triangulated mesh | Initialised property maps; `lambda0` zeroed |
 | `compute_*_lambda0_from_mesh()` | `setup_*_maps()` called | `lambda0[e]` set from 3-D edge lengths |
-| `check_gauss_bonnet()` | `theta_v[v]` set for all vertices | Throws `std::runtime_error` if Σ(2π−Θᵥ) ≠ 2π·χ(M) |
-| `enforce_gauss_bonnet()` | `theta_v[v]` set | Σ(2π−Θᵥ) = 2π·χ(M) guaranteed; `theta_v` modified |
+| `check_gauss_bonnet()` | `theta_v[v]` set · **EuclideanMaps or SphericalMaps only** | Throws `std::runtime_error` if Σ(2π−Θᵥ) ≠ 2π·χ(M). Deleted overload for `HyperIdealMaps` — see note below. |
+| `enforce_gauss_bonnet()` | `theta_v[v]` set · **EuclideanMaps or SphericalMaps only** | Σ(2π−Θᵥ) = 2π·χ(M) guaranteed; `theta_v` modified. Deleted overload for `HyperIdealMaps` — see note below. |
 | `newton_euclidean()` | GB satisfied · DOFs assigned · `lambda0` initialised | `NewtonResult.x` — converged scale factors; `.converged`, `.iterations`, `.grad_inf_norm` |
 | `newton_spherical()` | GB satisfied · DOFs assigned · gauge vertex pinned | Same as above |
-| `newton_hyper_ideal()` | `assign_all_dof_indices()` called · `lambda0` initialised | Same as above |
+| `newton_hyper_ideal()` | `assign_all_dof_indices()` called · **no `lambda0` needed** · **no GB pre-check** | Same as above. HyperIdeal computes lengths internally from b_v, a_e (via ζ₁₃/ζ₁₄/ζ₁₅); `lambda0` is irrelevant. Energy is strictly convex → Newton converges for any targets without a pre-check. |
 | `compute_cut_graph()` | Closed, orientable, triangulated mesh | `CutGraph.cut_edge_flags` — 2g seam edges; `.genus` |
 | `euclidean_layout()` | `newton_euclidean()` converged | `Layout2D.uv[v]` · `.halfedge_uv[h]` · `HolonomyData.translations` |
 | `spherical_layout()` | `newton_spherical()` converged | `Layout3D.xyz[v]` |
@@ -53,17 +53,28 @@ The solver never writes to pinned DOFs. All indexing is 0-based and contiguous.
 
 ---
 
-## Gauss–Bonnet — the most common source of failure
+## Gauss–Bonnet — the most common source of failure for Euclidean and Spherical
 
 Prescribing angles that violate Gauss–Bonnet means no conformal factor can realise the
 target metric — the Newton solver will iterate indefinitely without converging.
 
 ```cpp
-// Option A: verify before solving (throws on violation)
-check_gauss_bonnet(mesh, maps);
+// Option A: verify before solving — throws std::runtime_error on violation.
+// Works with EuclideanMaps and SphericalMaps only.
+check_gauss_bonnet(mesh, euclidean_maps);  // or spherical_maps
+check_gauss_bonnet(mesh, spherical_maps);
 
-// Option B: auto-correct (redistributes defect uniformly across all vertices)
-enforce_gauss_bonnet(mesh, maps);
+// Option B: auto-correct — redistributes the defect uniformly across all vertices.
+enforce_gauss_bonnet(mesh, euclidean_maps);
+enforce_gauss_bonnet(mesh, spherical_maps);
+
+// ⚠  HyperIdealMaps: both functions have deleted overloads — calling them is a
+// compile error.  Reason: the correct hyperbolic Gauss–Bonnet identity includes
+// a surface area term: Σ(2π−Θᵥ) − Area(M) = 2π·χ(M), not Σ(2π−Θᵥ) = 2π·χ(M).
+// No pre-check is needed for HyperIdeal: the energy is strictly convex
+// (Springborn 2020 Theorem 1.3), so newton_hyper_ideal converges for any targets.
 ```
 
 The target violation is `|Σ(2π−Θᵥ) − 2π·χ(M)| > tol`. Default tolerance: `1e-10`.
+Applicable to closed meshes only — for open meshes, pin the boundary directly
+and skip the Gauss–Bonnet check (the identity does not hold for meshes with boundary).

doc/getting-started.md

@@ -75,6 +75,21 @@ cmake --build build -j$(nproc)
 > **Note:** `-DWITH_CGAL=ON` implies `-DWITH_VIEWER=ON`. Do not use this in headless
 > environments — it will fail with `Failed to find wayland-scanner`.
 
+### Mode 4 — Low-memory build (RAM-constrained CI / Raspberry Pi ≤ 4 GB)
+
+For machines where the CGAL build OOMs (peak ~700 MB per compilation unit at `-O3`):
+
+```bash
+cmake -S code -B build -DWITH_CGAL_TESTS=ON \
+      -DCONFORMALLAB_LOW_MEMORY_BUILD=ON
+cmake --build build --target conformallab_cgal_tests -j1
+```
+
+`LOW_MEMORY_BUILD` applies four measures: `-O0` (no debug info), PCH off,
+unity batch size 1, `--no-keep-memory` linker flag.  Drops cc1plus peak from
+~700 MB to ~150-200 MB per TU.  Tests run ~15× slower at `-O0` but all pass.
+**Always use `-j1`** — parallel compilation would defeat the memory savings.
+
 ---
 
 ## Running a single test
@@ -114,16 +129,53 @@ After a full build (`-DWITH_CGAL=ON`):
 ./bin/conformallab_core --help
 ```
 
+### CLI parameter reference
+
+| Flag | Default | Description |
+|------|---------|-------------|
+| `-i / --input` | *required* | Input mesh file (OFF / OBJ / PLY) |
+| `-o / --output` | *(none)* | Save layout as OFF file |
+| `-j / --json` | *(none)* | Serialise solver result + UV to JSON |
+| `-x / --xml` | *(none)* | Serialise solver result + UV to XML |
+| `-g / --geometry` | `euclidean` | Target geometry: `euclidean` · `spherical` · `hyper_ideal` |
+| `-s / --show` | `false` | Open the input mesh in the interactive viewer |
+| `-v / --verbose` | `false` | Print mesh topology, DOF counts, convergence details |
+
+> **Tip:** `./bin/conformallab_core --help` always shows the canonical up-to-date
+> list generated by CLI11.  The table above matches `conformallab_cli.cpp`
+> as of v0.10.0.
+
 ## Example programs
 
 ```bash
+# PRIMARY USE CASE: conformally flatten a mesh to the plane
+./build/examples/example_flatten   code/data/obj/cathead.obj flat.off
+#  → non-trivial u_v (e.g. range ≈ 2.96), real conformal deformation
+
+# CGAL public API: one-call interface (natural theta by default)
+./build/examples/example_cgal_api  [input.off]
+
+# Full pipeline with JSON/XML serialisation and round-trip test
 ./build/examples/example_layout    [input.off] [layout.off] [result.json]
+
+# Solver test (natural theta — u_v ≈ 0, used for pipeline validation)
 ./build/examples/example_euclidean [input.off] [output.off]
+
+# Hyper-ideal (hyperbolic) functional
 ./build/examples/example_hyper_ideal [input.off] [output.off]
-./build/examples/example_viewer    [input.off]   # interactive, requires WITH_VIEWER
+
+# Interactive viewer (requires WITH_VIEWER)
+./build/examples/example_viewer    [input.off]
 ```
 
-`example_layout.cpp` is the best starting point — it shows the complete pipeline in ~120 lines.
+**Start here:** `example_flatten.cpp` shows the primary use case — real conformal
+flattening with `Θ_v = 2π`.  `example_layout.cpp` adds JSON/XML serialisation.
+
+> **Note on "natural theta":** `example_euclidean` and `example_layout` use the
+> "natural theta" testing trick (`Θ_v = actual angle sum at x=0`), which makes
+> `x* = 0` trivially the equilibrium.  The output `u_v ≈ 0` is expected and
+> correct for a pipeline test, but means **no conformal deformation was applied**.
+> For real UV parameterisation, use `example_flatten.cpp`.
 
 **Expected output of `example_euclidean` on the built-in quad-strip mesh:**
 ```

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`).
-The analytic Hessian through the `(bᵢ, aₑ) → lᵢⱼ → ζ → αᵢⱼ/βᵢ` chain is Phase 9b.
+**Hessian:** a symmetric finite-difference approximation (see `hyper_ideal_hessian.hpp`).
+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

@@ -18,6 +18,102 @@ After the meeting, a fourth file `2026-XX-XX-meeting-notes.md` should
 land here too — the memo template at the bottom of `agenda.md` is the
 suggested structure.
 
+---
+
+# External audit documents (2026-05-30 / 2026-05-31)
+
+Twelve self-contained audit documents covering the project from every angle.
+Each is written so a **fresh session can pick up any finding and act on it**
+without prior context: every finding has an exact `file:line`, a code snippet,
+a concrete fix, and acceptance criteria.
+
+## ⛔ Cross-cutting blocker — resolve before any public release or CGAL submission
+
+> **G0 — Porting rights.** The original Varylab/TU-Berlin project
+> (`gitlab.discretization.de:varylab/conformallab`) carries **no license**
+> (all rights reserved), and conformallab++ is a documented **port**
+> (derivative work). The right to port + (re)license must be clarified with the
+> original authors **before** any license decision, public release, or CGAL
+> 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
+
+| Document | Focus | Headline |
+|---|---|---|
+| [`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 | **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/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–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 | **⛔ 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 | **⛔ 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/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 | **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 **still undocumented (open, S4)** |
+
+## Session 1 — implemented (2026-05-31)
+
+Branch `fix/b1-v3-c1-quick-wins`, 9 commits, **298/298 CGAL tests green**.
+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**.
+- **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.
+
+## 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. **The Session-1 quick-win
+batch (B1, B3, C1, C2, V3, …) is done** — for what to pick up next, follow the
+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
 
 - **briefing.md** + **questions.md**: as part of the meeting-confirmation

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`,

doc/reviewer/usability-audit-2026-05-31.md

@@ -0,0 +1,556 @@
+# Usability & Documentation Audit — ConformalLabpp v0.10.0
+
+**Date:** 2026-05-31
+**Auditor:** External reviewer (Claude Sonnet 4.6)
+**Branch:** `review/usability-audit-2026-05-31`
+**Base:** `main` (post external-audit-2026-05-30 merge)
+**Focus:** Documentation quality, API usability, new-user experience
+
+This document is self-contained. A new session can pick up any finding
+and act on it without prior context. Each finding includes exact file paths,
+the concrete problem, and a precise fix with acceptance criteria.
+
+Status legend: 🔴 Usability bug · 🟡 Doc error/stale · 🟠 Missing content
+
+---
+
+## How to read this document in a new session
+
+```bash
+git checkout review/usability-audit-2026-05-31
+
+# No build needed for most fixes — documentation and example changes only.
+# For changes that touch compilable code (examples), verify with:
+cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON -DCONFORMALLAB_LOW_MEMORY_BUILD=ON
+cmake --build build-cgal --target conformallab_cgal_tests -j1
+ctest --test-dir build-cgal -R '^cgal\.' --output-on-failure
+```
+
+---
+
+## FINDING-U1 — 🔴 CRITICAL: All examples show the trivial identity map, not real conformal flattening
+
+### Location
+- `README.md` lines 100–121 ("Minimal usage")
+- `code/examples/example_euclidean.cpp` lines 72–83 (Step 4)
+- `code/examples/example_layout.cpp` lines 40–49 (`set_natural_theta`)
+- `code/examples/example_hyper_ideal.cpp` lines 89–98 (Step 3)
+
+### Problem
+
+Every example uses the "natural theta" trick:
+```cpp
+auto G0 = euclidean_gradient(mesh, x0, maps);
+for (auto v : mesh.vertices())
+    if (maps.v_idx[v] >= 0)
+        maps.theta_v[v] -= G0[maps.v_idx[v]];
+```
+
+This sets Θ_v to the *actual angle sums at x=0*, making x*=0 the equilibrium.
+The result: `u_v ≈ 0` everywhere — **no deformation**. The map is identity.
+
+Natural theta is a valid test pattern (it proves Newton converges in 0 or 1
+steps and guards against solver regressions). But as the *only* pattern shown to
+a new user, it gives the false impression that the library produces trivial
+output. The real use case — **conformally flatten a surface to the plane** —
+is never demonstrated.
+
+### What is missing
+
+A "real-world" example that shows the actual use case:
+
+```cpp
+// Conformally flatten a mesh to the plane.
+// Set Θ_v = 2π (regular interior vertex, zero cone angle defect) for all
+// free vertices; enforce_gauss_bonnet makes the assignment topologically
+// consistent.  The result is a flat conformal map with minimal angle distortion.
+for (auto v : mesh.vertices())
+    maps.theta_v[v] = conformallab::TWO_PI;
+conformallab::enforce_gauss_bonnet(mesh, maps);
+auto res = conformallab::newton_euclidean(mesh, x0, maps);
+auto layout = conformallab::euclidean_layout(mesh, res.x, maps);
+// layout.uv[v.idx()] now holds 2D coordinates of the flattening
+```
+
+### Fix
+
+1. **Add a new example** `code/examples/example_flatten.cpp`:
+   - Takes a real mesh from `code/data/` (e.g. `cathead.obj` or `torus_8x8.off`)
+   - Sets Θ_v = 2π, calls `enforce_gauss_bonnet`, solves Newton, computes layout
+   - Saves a UV-mapped OFF file
+   - Shows non-trivial u_v values at output
+
+2. **Update README "Minimal usage"** to show the flatten use case instead of
+   (or in addition to) natural theta.  Add a comment explaining what natural
+   theta is and when to use it.
+
+3. **Add a comment** to all existing examples explaining that natural theta
+   is a *testing convenience*, not a typical use case.
+
+### Acceptance criteria
+- [ ] `example_flatten.cpp` compiles and runs with a real mesh input
+- [ ] Output shows non-trivial `u_v ≠ 0` and a visible UV parameterisation
+- [ ] README "Minimal usage" demonstrates real flattening (not identity)
+- [ ] Existing examples have a comment: "Note: natural theta → x* = 0 (testing
+      pattern). For real flattening, set theta_v = TWO_PI and call enforce_gauss_bonnet."
+
+---
+
+## FINDING-U2 — 🔴 CRITICAL: No example for the CGAL public API
+
+### Location
+- `code/include/CGAL/Discrete_conformal_map.h` (the "user-facing entry")
+- `code/examples/` (no example uses CGAL::discrete_conformal_map_euclidean)
+
+### Problem
+
+`Discrete_conformal_map.h` is documented as "User-facing entry point for the
+Discrete_conformal_map package."  The Doxygen doc shows a minimal usage snippet:
+
+```cpp
+auto result = CGAL::discrete_conformal_map_euclidean(mesh);
+```
+
+But this snippet only exists in a Doxygen comment — no compilable example
+demonstrates it. All four `code/examples/*.cpp` use the *internal* API
+(`setup_euclidean_maps + newton_euclidean`). A CGAL user who reaches this
+library via the CGAL ecosystem has no runnable starting point.
+
+The CGAL API and the internal API also produce *different result types*:
+- Internal: `NewtonResult.x` — DOF-index-indexed `std::vector<double>`
+- CGAL: `Conformal_map_result.u_per_vertex` — vertex-index-indexed `std::vector<FT>`
+
+This discrepancy is not explained anywhere for users.
+
+### Fix
+
+Add `code/examples/example_cgal_api.cpp`:
+```cpp
+// Demonstrates the CGAL public API (Discrete_conformal_map.h).
+// Contrast with example_euclidean.cpp which uses the internal API directly.
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Discrete_conformal_map.h>
+// ...
+auto result = CGAL::discrete_conformal_map_euclidean(mesh);
+// result.u_per_vertex[v.idx()] — per-vertex scale factor
+// result.converged, result.iterations, result.gradient_norm
+```
+
+Also add to README: a short paragraph explaining when to use CGAL API
+vs internal API (CGAL API: simpler, less control; internal API: full pipeline,
+holonomy, period matrix).
+
+### Acceptance criteria
+- [ ] `example_cgal_api.cpp` compiles with `-DWITH_CGAL=ON`
+- [ ] Example uses `CGAL::discrete_conformal_map_euclidean` (not internal API)
+- [ ] README explains the two API levels and when to use each
+
+---
+
+## FINDING-U3 — 🟡 DOC ERROR: `contracts.md` incorrect after Finding-B
+
+### Location
+`doc/api/contracts.md` line 16
+
+### Problem
+
+```markdown
+| `check_gauss_bonnet()` | `theta_v[v]` set for all vertices |
+  Throws `std::runtime_error` if Σ(2π−Θᵥ) ≠ 2π·χ(M) |
+```
+
+After the external-audit-2026-05-30 Finding-B fix, calling
+`check_gauss_bonnet(mesh, hyper_ideal_maps)` is a **compile error**
+(the overload is `= delete`).  The contract table implies it works for
+all map types — which is now false.
+
+Also missing from the table: the HyperIdeal functional needs no
+Gauss-Bonnet pre-check because its energy is strictly convex
+(Springborn 2020 Theorem 1.3) — Newton converges for any target angles.
+
+### Fix
+
+```markdown
+| `check_gauss_bonnet()` | `theta_v[v]` set · **Euclidean or Spherical maps only** |
+  Throws if Σ(2π−Θᵥ) ≠ 2π·χ(M). **Not applicable to HyperIdealMaps** —
+  deleted overload; HyperIdeal needs no pre-check (strictly convex energy). |
+```
+
+Also add a row for `enforce_gauss_bonnet`:
+```markdown
+| `enforce_gauss_bonnet()` | `theta_v[v]` set · **Euclidean or Spherical maps only** |
+  Shifts all Θᵥ by δ = (lhs−rhs)/V so that GB holds exactly.
+  **Not applicable to HyperIdealMaps** — deleted overload. |
+```
+
+### Acceptance criteria
+- [ ] `check_gauss_bonnet` row notes "Euclidean/Spherical only, not HyperIdeal"
+- [ ] `enforce_gauss_bonnet` row added with same note
+- [ ] A brief explanation why HyperIdeal doesn't need GB: "strictly convex energy"
+
+---
+
+## FINDING-U4 — 🟡 DOC ERROR: Stale version in README (v0.9.0)
+
+### Location
+`README.md` line 17
+
+### Problem
+
+```markdown
+**Status:** v0.9.0 — Phases 1–9a complete, Phase 8b-Lite CGAL API surface.
+```
+
+The current release is **v0.10.0** (tagged on `main`; CLAUDE.md confirms this).
+CLAUDE.md was updated to v0.10.0 but the README was not.
+
+### Fix
+
+```markdown
+**Status:** v0.10.0 — Phases 1–9b complete. Newton solvers for five DCE
+models (Euclidean / Spherical / HyperIdeal / CP-Euclidean / Inversive-Distance),
+priority-BFS layout in ℝ²/S²/Poincaré disk, Gauss–Bonnet, tree-cotree cut
+graph, Möbius holonomy, period matrix (genus 1), fundamental domain,
+halfedge_uv texture atlas, JSON/XML serialisation, CLI app.
+277 tests passing, 0 skipped — see [`doc/api/tests.md`](doc/api/tests.md).
+```
+
+### Acceptance criteria
+- [ ] README status line shows v0.10.0 and the correct feature list
+- [ ] Test count updated to 277
+
+---
+
+## FINDING-U5 — 🟡 DOC ERROR: CGAL header comment describes superseded state
+
+### Location
+`code/include/CGAL/Discrete_conformal_map.h` lines 14–16
+
+### Problem
+
+```cpp
+/// This header provides a single function — `discrete_conformal_map_euclidean`
+/// — that computes a Euclidean discrete-conformal flattening … Spherical and
+/// hyperbolic variants are scheduled for Phase 8b.2 once the Euclidean
+/// pattern is validated
+```
+
+But the file already contains `discrete_conformal_map_spherical()` and
+`discrete_conformal_map_hyper_ideal()` (and the other models).  The comment
+describes Phase 8a state; the file is at Phase 8b-Lite.
+
+### Fix
+
+Update the `\file` Doxygen block at the top of `Discrete_conformal_map.h`:
+
+```cpp
+/*!
+\file CGAL/Discrete_conformal_map.h
+\ingroup PkgConformalMapRef
+
+User-facing entry points for the Discrete_conformal_map package.
+
+This header provides five functions covering all three DCE geometries:
+  - `CGAL::discrete_conformal_map_euclidean`   (flat / ℝ²)
+  - `CGAL::discrete_conformal_map_spherical`   (S², genus 0)
+  - `CGAL::discrete_conformal_map_hyper_ideal` (H², genus ≥ 1)
+  - `CGAL::discrete_circle_packing_euclidean`  → Discrete_circle_packing.h
+  - `CGAL::discrete_inversive_distance_map`    → Discrete_inversive_distance.h
+...
+*/
+```
+
+### Acceptance criteria
+- [ ] `\file` doc block lists all five entry functions
+- [ ] No forward-reference to "planned Phase 8b.2" — it is already implemented
+
+---
+
+## FINDING-U6 — 🟡 STALE EXAMPLES: New gauge-vertex overload not reflected
+
+### Location
+- `README.md` lines 107–110
+- `code/examples/example_euclidean.cpp` lines 62–68
+- `code/examples/example_layout.cpp` lines 52–59 (`pin_first` helper)
+
+### Problem
+
+Finding-D (external-audit-2026-05-30) added a clean gauge-vertex overload:
+```cpp
+assign_euclidean_vertex_dof_indices(mesh, maps, gauge_vertex);
+```
+
+But all examples still use the old verbose loop:
+```cpp
+auto vit = mesh.vertices().begin();
+maps.v_idx[*vit++] = -1;   // pinned
+int idx = 0;
+for (; vit != mesh.vertices().end(); ++vit)
+    maps.v_idx[*vit] = idx++;
+```
+
+A new user copying from the examples will write the old error-prone pattern
+instead of the new clean overload — defeating the purpose of the fix.
+
+### Fix
+
+Replace the old loop in all three locations with the new overload:
+```cpp
+// Pin the first vertex, assign sequential DOF indices to the rest.
+int n = assign_euclidean_vertex_dof_indices(mesh, maps,
+                                            *mesh.vertices().begin());
+```
+
+Same for `assign_vertex_dof_indices` (spherical) and
+`assign_inversive_distance_vertex_dof_indices` wherever they appear in examples.
+
+### Acceptance criteria
+- [ ] README "Minimal usage" uses the overload
+- [ ] `example_euclidean.cpp` uses the overload
+- [ ] `example_layout.cpp` `pin_first` helper replaced with the overload
+- [ ] Old manual loop removed from all user-facing code
+
+---
+
+## FINDING-U7 — 🟠 MISSING: `CONFORMALLAB_LOW_MEMORY_BUILD` absent from `getting-started.md`
+
+### Location
+`doc/getting-started.md` — "Build modes" section
+
+### Problem
+
+`getting-started.md` documents five compile-time modes but does not list
+`CONFORMALLAB_LOW_MEMORY_BUILD=ON`.  A Raspberry Pi or low-RAM user who
+reads `getting-started.md` (the natural first stop) cannot find the flag.
+It only appears in CLAUDE.md and the README table.
+
+### Fix
+
+Add to `getting-started.md` after the `CONFORMALLAB_FAST_TEST_BUILD` entry:
+
+```markdown
+### Low-memory build (RAM-constrained CI / Raspberry Pi)
+
+For machines with ≤ 4 GB RAM, the default CGAL build (PCH + -O3) OOMs
+the compiler.  `LOW_MEMORY_BUILD` uses `-O0`, disables PCH, and sets
+unity batch size to 1, keeping `cc1plus` peak at ~150-200 MB per TU:
+
+```bash
+cmake -S code -B build -DWITH_CGAL_TESTS=ON \
+      -DCONFORMALLAB_LOW_MEMORY_BUILD=ON
+cmake --build build --target conformallab_cgal_tests -j1
+```
+
+Tests run ~15× slower at -O0 but all pass.  Useful for CI runners with
+limited memory.
+```
+
+### Acceptance criteria
+- [ ] `getting-started.md` documents `CONFORMALLAB_LOW_MEMORY_BUILD`
+- [ ] Example command shows `-j1` (required for memory safety)
+
+---
+
+## FINDING-U8 — 🟠 MISSING: `LOW_MEMORY_BUILD` absent from README compile-time table
+
+### Location
+`README.md` lines 53–82 ("Compile-time workflow modes")
+
+### Problem
+
+The README code block shows five compile-time patterns but omits
+`CONFORMALLAB_LOW_MEMORY_BUILD`.  This is the flag that re-enabled CI
+and is directly relevant to anyone building on resource-constrained hardware.
+
+### Fix
+
+Add to the README compile-time examples:
+```bash
+# Low-memory build: -O0, no PCH, unity batch 1 (~150 MB peak per cc1plus).
+# Use on Raspberry Pi / CI runners with ≤ 4 GB RAM.  Tests run ~15× slower.
+cmake -S code -B build -DWITH_CGAL_TESTS=ON -DCONFORMALLAB_LOW_MEMORY_BUILD=ON
+cmake --build build --target conformallab_cgal_tests -j1
+```
+
+### Acceptance criteria
+- [ ] README shows `LOW_MEMORY_BUILD` usage with a one-line explanation
+
+---
+
+## FINDING-U9 — 🟠 MISSING: `layout.uv` indexing semantics not documented
+
+### Location
+`code/examples/example_layout.cpp` line 127
+`code/include/layout.hpp` — `Layout2D` struct definition
+
+### Problem
+
+```cpp
+auto& p = layout.uv[v.idx()];  // example_layout.cpp:127
+```
+
+`Layout2D.uv` is a `std::vector` indexed by the raw integer `.idx()` of
+`Vertex_index`.  This is undocumented: nowhere is it stated that the vector
+length equals `num_vertices(mesh)` and that the index is contiguous.  In CGAL,
+`vertex.idx()` is contiguous on a fresh `Surface_mesh` but may have gaps after
+vertex removals.
+
+The `halfedge_uv` field (also in `Layout2D`) has the same issue for halfedge
+indices.
+
+### Fix
+
+In `layout.hpp`, add explicit documentation to the `Layout2D` struct:
+
+```cpp
+struct Layout2D {
+    /// UV coordinate of vertex v in the layout plane.
+    /// Indexed by `v.idx()` — length equals `mesh.number_of_vertices()`.
+    /// \pre  No vertices have been removed from `mesh` since construction
+    ///       (i.e., `mesh.is_valid(false)` and no compaction was performed).
+    std::vector<Eigen::Vector2d> uv;
+
+    /// UV of source(h) as seen from face(h), indexed by `h.idx()`.
+    /// At seam halfedges the two opposite halfedges carry different UV values,
+    /// enabling proper per-halfedge UV for GPU texture atlasing.
+    std::vector<Eigen::Vector2d> halfedge_uv;
+    ...
+};
+```
+
+Also update `example_layout.cpp` to add a comment:
+```cpp
+// layout.uv is indexed by v.idx() — valid as long as no vertices were
+// removed from the mesh after loading.
+auto& p = layout.uv[v.idx()];
+```
+
+### Acceptance criteria
+- [ ] `Layout2D.uv` and `Layout2D.halfedge_uv` have explicit Doxygen docs
+      stating the index semantics and the no-compaction precondition
+- [ ] `example_layout.cpp` has an inline comment at the access site
+
+---
+
+## FINDING-U10 — 🟠 MISSING: CLI parameter reference not in README or getting-started
+
+### Location
+- `README.md` (has CLI examples but no parameter table)
+- `doc/getting-started.md` (same)
+- `code/src/apps/v0/conformallab_cli.cpp` lines 333–340 (source of truth)
+
+### Problem
+
+The CLI app (`conformallab_core`) is the primary non-programmatic entry point
+for the library.  The README shows a few usage examples but no parameter
+reference.  A user running `--help` gets CLI11's auto-generated output, but
+that output is not reproduced anywhere in the documentation.
+
+Current parameters (from `conformallab_cli.cpp`):
+
+| Flag | Default | Description |
+|------|---------|-------------|
+| `-i / --input` | *required* | Input mesh (OFF / OBJ / PLY) |
+| `-o / --output` | *(none)* | Output layout OFF file |
+| `-j / --json` | *(none)* | Serialise result to JSON |
+| `-x / --xml` | *(none)* | Serialise result to XML |
+| `-g / --geometry` | `euclidean` | Target geometry: `euclidean` \| `spherical` \| `hyper_ideal` |
+| `-s / --show` | `false` | Open input mesh in interactive viewer |
+| `-v / --verbose` | `false` | Print topology / DOF stats |
+
+### Fix
+
+Add a concise parameter table to `doc/getting-started.md` in the
+"First run — CLI app" section.  The README already shows the most common
+invocations; the table can live in getting-started to avoid README bloat.
+
+### Acceptance criteria
+- [ ] `doc/getting-started.md` has a parameter table matching the source
+- [ ] Each flag has a one-line description and notes the default
+- [ ] Cross-reference to `--help` for the canonical up-to-date list
+
+---
+
+## FINDING-U11 — 🟠 ROADMAP: CLI missing three useful parameters + two models
+
+### Location
+`code/src/apps/v0/conformallab_cli.cpp`
+
+### Problem (not a bug — a gap for a future session to fill)
+
+**Missing solver-tuning parameters** (low effort, high value):
+
+The Newton tolerance and iteration limit are hardcoded at their defaults
+(`tol=1e-8`, `max_iter=200`).  Any user who wants to solve a stiff mesh
+more carefully or stop early has to recompile.
+
+```bash
+# Not yet possible:
+./conformallab_core -i hard_mesh.off -g euclidean --tol 1e-12 --max-iter 1000
+./conformallab_core -i quick_check.off --tol 1e-4 --max-iter 20
+```
+
+**Missing geometry modes** (medium effort):
+
+The two Phase-9a functionals are not reachable from the CLI at all:
+
+```bash
+# Not yet possible — needs run_cp_euclidean() + run_inversive_distance() helpers:
+./conformallab_core -i mesh.off -g cp_euclidean
+./conformallab_core -i mesh.off -g inversive_distance
+```
+
+CP-Euclidean (`Discrete_circle_packing.h`) and Inversive-Distance
+(`Discrete_inversive_distance.h`) are fully implemented in the library and
+exposed via the CGAL public API, but a CLI user cannot reach them.
+
+### Suggested implementation (for a future session)
+
+1. **`--tol <double>` and `--max-iter <int>`:** one `app.add_option` each,
+   thread through all three `run_*` functions.  ~15 lines.
+
+2. **`-g cp_euclidean` and `-g inversive_distance`:** add two new
+   `run_cp_euclidean()` / `run_inversive_distance()` helpers following
+   the existing `run_euclidean()` pattern.  ~60 lines each.  Add to the
+   `-g` `CLI::IsMember` list.  Register the two new geometry strings in
+   the dispatch block.
+
+3. **Update README + getting-started** to list all five geometry modes.
+
+### Effort estimate
+- `--tol` + `--max-iter`: ~30 min
+- Phase-9a CLI exposure: ~2 hours (mostly copy-paste + adaptation)
+- Tests: add two CLI smoke tests (geometry converges on tetrahedron)
+
+### Acceptance criteria (when implemented)
+- [ ] `--tol` and `--max-iter` are accepted and forwarded to all three/five solvers
+- [ ] `./conformallab_core -i mesh.off -g cp_euclidean -o out.off` works
+- [ ] `./conformallab_core -i mesh.off -g inversive_distance -o out.off` works
+- [ ] `--help` output lists all five geometry modes
+
+---
+
+## Summary table
+
+| ID | File | Type | Severity | Status |
+|----|------|------|----------|--------|
+| U1 | README + 4 examples | Usability | 🔴 Critical | ✅ Fixed 2026-05-31 |
+| U2 | examples/ (missing) | Usability | 🔴 Critical | ✅ Fixed 2026-05-31 |
+| U3 | `contracts.md:16` | Doc error | 🟡 Medium | ✅ Fixed 2026-05-31 |
+| U4 | `README.md:17` | Stale | 🟡 Medium | ✅ Fixed 2026-05-31 |
+| U5 | `Discrete_conformal_map.h:14` | Stale | 🟡 Medium | ✅ Fixed 2026-05-31 |
+| U6 | README + 2 examples | Stale | 🟡 Medium | ✅ Fixed 2026-05-31 |
+| U7 | `getting-started.md` | Missing | 🟠 Minor | ✅ Fixed 2026-05-31 |
+| U8 | `README.md` | Missing | 🟠 Minor | ✅ Fixed 2026-05-31 |
+| U9 | `layout.hpp` + example | Missing | 🟠 Minor | ✅ Fixed 2026-05-31 |
+| U10 | `getting-started.md` | Missing | 🟠 Minor | ✅ Fixed 2026-05-31 |
+| U11 | `conformallab_cli.cpp` | Roadmap | 🟠 Minor | ✅ Added to phases.md (9h) |
+
+**Priority order for fixing:**
+1. ~~U1 + U2~~ ✅ done
+2. U3 (correctness — contracts.md wrong after Finding-B)
+3. U4 + U5 + U6 (stale content — quick fixes)
+4. U7 + U8 + U9 + U10 (missing docs — minor)
+5. U11 (roadmap — CLI extensions, not urgent)

doc/roadmap/phase-orchestration.md

@@ -0,0 +1,163 @@
+# Phase Orchestration — Sessions × Models × Review Gates
+
+**Updated:** 2026-06-01
+**Purpose:** structured plan to work through roadmap phases across multiple
+short sessions, each run by the **model best suited** to the work, with an
+**external Opus reviewer session after every implementation session**.
+Mirrors the audit workflow in
+[`doc/reviewer/finding-orchestration.md`](../reviewer/finding-orchestration.md).
+
+> **Ready-to-paste prompts** for every pending session (P1–P4 + 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 phase details in
+   `phases.md` and (for research items) `research-track.md`. Each phase entry
+   has a Java reference, a mathematical source, and acceptance criteria.
+3. When it finishes, start the paired **🔍 Opus review session** — it validates
+   the diff (correctness, gradient checks, no parity regression) and only then
+   is the batch "done".
+4. Mark the session ✅ here and move on.
+
+**Rationale:** same as the audit system — Haiku/Sonnet for well-specified work,
+Opus reserved for math/architecture work and the review gate.
+
+---
+
+## Model-assignment heuristic
+
+| Model | Take phases that are… | Examples |
+|---|---|---|
+| **Haiku** | mechanical Java ports, local CLI additions, doc-only, no new math | 9g.1 (quality measures), 9h.1/9h.2 (CLI), 9d.3 (Stereographic port) |
+| **Sonnet** | medium Java ports or research with explicit acceptance criteria | 9e (CirclePatternLayout), Phase 12 (Decorated DCE) |
+| **Opus** | math derivations, architecture decisions, all 🔍 review gates | 9b-analytic (Schläfli Hessian), 9c (genus > 1), all reviews |
+
+> Decision rule isn't "how large" but "how much must be *understood* to get it
+> right." A mechanical 600-line Java port can be Haiku; a 200-line numerical
+> reformulation needs Opus.
+
+---
+
+## Master phase table
+
+Status: ✅ done · ⬜ open (actionable) · ⏸ deferred/gated · ⛔ blocked · 👤 needs reviewer input
+
+| Phase | Focus | Effort | Sev | Status | Model | Session |
+|-------|-------|--------|-----|--------|-------|---------|
+| 9h.1 | Newton CLI params (`--tol`, `--max-iter`) | 30 min | 🔵 | ⬜ | Haiku→🔍Opus | **P1** |
+| 9h.2 | Phase 9a models in CLI (`-g cp_euclidean`, `-g inversive_distance`) | 2–4 h | 🔵 | ⬜ | Haiku→🔍Opus | **P1** |
+| 9g.1 | Conformal quality measures (`conformal_quality.hpp`) | ~3 days | 🟡 | ⬜ | Haiku→🔍Opus | **P1** |
+| 9d.3 | StereographicUnwrapper port | ~3 days | 🟡 | ⬜ | Haiku→🔍Opus | **P1** |
+| 12 | Decorated DCE & geometric transition | medium | 🔴 | ⬜ | Sonnet→🔍Opus | **P2** |
+| 9e | CirclePatternLayout + CirclePatternUtility | medium | 🟡 | ⬜ | Sonnet→🔍Opus | **P3** |
+| 9d.4 | MobiusCenteringFunctional (variational Möbius centring) | medium | 🟡 | ⬜ | Sonnet→🔍Opus | **P3** |
+| 9g.2 | Period-matrix convergence study (experiment test) | ~3 days | 🔵 | ⬜ | Sonnet→🔍Opus | **P3** |
+| GC-1 | geometry-central cross-validation | small | 🔵 | ⏸ | Sonnet→🔍Opus | (after reviewer Q5) |
+| 9d.1 | ConesUtility port (Euclidean cone singularities) | medium | 🟡 | ⏸ | Sonnet→🔍Opus | (after reviewer Q1) |
+| 9d.2 | Non-Euclidean cone extensions (RESEARCH) | large | 🟡 | ⏸ | Opus | (after reviewer Q1) |
+| 9f | Polygon Laplacian (RESEARCH, no Java parent) | medium | 🟡 | ⏸ | Sonnet/Opus | (after reviewer Q1) |
+| 9b-analytic | Analytic HyperIdeal Hessian via Schläfli | 10–14 days | 🔴 | ⏸ | Opus | **P4** (after reviewer Q3) |
+| 9c | 4g-polygon fundamental domain (genus g > 1) | ~5 weeks | 🔴 | ⏸ | Opus | **P5** (after reviewer Q4 + G0) |
+| 10a | Holomorphic/harmonic 1-forms + DEC layer | ~6 weeks | 🔴 | ⏸ | Opus | **P6** (after 9c) |
+| 10b | Siegel period matrix Ω ∈ H_g | ~1 week | 🟡 | ⏸ | Opus | **P6** (after 10a) |
+| 10c | Full genus-g≥2 uniformization | large | 🔴 | ⏸ | Opus | **P7** (after 10b) |
+| 13 | Canonical tessellations + polyhedral realisation (Chain B capstone) | very large | 🔴 | ⏸ | Opus | **P8** (after 10c) |
+| A4, A5 (public naming) | CGAL entry-point renames + `Circle_packing_result` unification | 🟡 | ⏸ | Opus | S6 in audit system (after G0) |
+
+---
+
+## Session sequence
+
+### ⬜ P1 — Quick wins (Haiku → 🔍 Opus)
+
+Four independent items, all with well-specified acceptance criteria and no new
+math. Can run in one batch or as two sub-sessions (9h.1+9h.2 first as a
+~half-day warmup, then 9g.1+9d.3 as a ~1-week batch).
+
+- **9h.1** — `--tol` and `--max-iter` CLI options in `conformallab_cli.cpp`;
+  thread through the three `run_*()` helpers; update `doc/getting-started.md`.
+- **9h.2** — `-g cp_euclidean` and `-g inversive_distance` routes in the CLI;
+  add both to the `CLI::IsMember` validator; update README + getting-started.
+- **9g.1** — `code/include/conformal_quality.hpp` with `IsothermicityMeasure`,
+  `DiscreteConformalEquivalenceMeasure`, `FlippedTriangles`, `LengthCrossRatio`,
+  `ConvergenceUtility` measures; one test each in `code/tests/cgal/`.
+- **9d.3** — `code/include/stereographic_layout.hpp` porting
+  `StereographicUnwrapper.java` (266 LoC); round-trip test required.
+- **🔍 Opus review:** gradient-check test for quality measures, CLI smoke test
+  on a real mesh, no count regression.
+
+### ⬜ P2 — Decorated DCE transition (Sonnet → 🔍 Opus)
+
+- **Phase 12** — Penner-coordinate decoration layer + geometric transition driver
+  + validation harness. Mathematical reference: Bobenko–Lutz 2025 §3.
+  Acceptance criteria (all four must pass — from `phases.md`):
+  - Decoration round-trip `I_ij ↔ (r_i, r_j, ℓ)` at machine precision.
+  - At κ=0: bit-for-bit match with existing euclidean/inversive path.
+  - Gauss–Bonnet holds per geometry; invariant constant across κ-transition.
+- **🔍 Opus review:** validate the three acceptance criteria numerically;
+  verify the cross-geometry invariant test; confirm no parity regression.
+
+### ⬜ P3 — Circle pattern embedding + Möbius centring + convergence study (Sonnet → 🔍 Opus)
+
+- **9e** — `code/include/circle_pattern_layout.hpp` porting
+  `CirclePatternLayout` + `CirclePatternUtility` (Phase 9a.1 prerequisite ✅).
+  Test: given ρ values from a solved CP-Euclidean system, embedded vertex
+  positions are self-consistent (intersection angles match the input).
+- **9d.4** — upgrade `normalise_hyperbolic()` in `layout.hpp` via
+  `MobiusCenteringFunctional` (Lorentz energy; Java: 289 LoC). Retain
+  Fréchet mean as fallback. Test: compare old vs new centering on brezel.obj.
+- **9g.2** — `code/tests/cgal/test_period_matrix_convergence.cpp` (experiment):
+  genus-1 mesh with known analytic τ, igl::loop subdivision, noise, assert
+  |τ_computed − τ_expected| decreases with refinement.
+- **🔍 Opus review:** cross-validate circle pattern positions against solver
+  output; verify convergence-study direction; confirm no count regression.
+
+### ⏸ P4 — Analytic HyperIdeal Hessian (Opus) — awaiting reviewer Q3
+
+**Precondition:** reviewer Q3 answer must confirm the ~6× speedup over
+block-FD is worth ~2 weeks at their typical mesh sizes. If Q3 is "no" or
+"later", phase stays ⏸. If "yes" or "above V > X", proceed.
+
+- **Phase 9b-analytic** — full analytic Hessian via Schläfli identity.
+  Derivation: `doc/math/hyperideal-hessian-derivation.md`.
+  Chain: `(bᵢ, aₑ) → ℓᵢⱼ → ζ₁₃/ζ₁₄/ζ₁₅ → αᵢⱼ/βᵢ → ∂²E/∂u²`.
+  Block-FD path must be retained as a compile-time cross-check.
+  Acceptance: analytic and block-FD agree to 1e-6; golden-value tests pass
+  bit-for-bit; benchmark: analytic faster on V > 500.
+- **No separate review gate** (Opus implements + self-reviews; block-FD is the
+  independent cross-check).
+
+### ⏸ P5 — Genus > 1 fundamental domain (Opus) — awaiting reviewer Q4 + G0
+
+**Preconditions:** (a) reviewer Q4 answer (port-literal vs. re-derive from
+Springborn 2020 §5); (b) G0 porting rights confirmed in writing.
+
+- **Phase 9c** — 4g-polygon boundary walk; high-precision substrate
+  (`boost::multiprecision::cpp_dec_float_50`, local to uniformization module).
+  See `phases.md` §9c for the full Java source list and the precision
+  prerequisite note.
+
+### ⏸ P6–P8 — Higher-genus research chain (Opus, sequential)
+
+Phases 10a → 10b → 10c → 13. Each gates the next. Start only after P5 (9c) is
+complete. Full scope in `phases.md` and `research-track.md`.
+
+---
+
+## Review-gate checklist (every 🔍 Opus session checks)
+
+- [ ] Builds clean; full CGAL suite green; count matches `doc/api/tests.md`.
+- [ ] Any new functional has a gradient-check test (see `CLAUDE.md` §Test design patterns).
+- [ ] No Java golden-vector / parity test perturbed (HardJava clamp defaults intact).
+- [ ] Numeric changes are value-identical where claimed, or justified + covered by a test.
+- [ ] New public surface (result types, enums, CGAL headers) is intentional;
+      documented in `doc/api/headers.md` and `doc/api/contracts.md`.
+- [ ] Commit message attributes the implementing model
+      (`Co-Authored-By: Claude <Model> <noreply@anthropic.com>`).
+- [ ] Phase marked ✅ in the master table above with the commit ref.
+- [ ] `porting-status.md` §7 updated if the phase adds a C++-only capability.

doc/roadmap/phases.md

@@ -10,6 +10,26 @@
 
 ---
 
+## ◼ Current focus (as of 2026-06-01)
+
+> **Agent entry point.** Read this box first. For model assignments, priority
+> order, and ready-to-paste prompts, see
+> [`phase-orchestration.md`](phase-orchestration.md) and
+> [`session-prompts.md`](session-prompts.md).
+> Audit-finding sessions run in parallel from `doc/reviewer/`.
+
+| Track | Next session | Gating |
+|---|---|---|
+| **Audit findings** | S3: H3/H4/H5/V5/V6 — Sonnet → Opus review | None — start now |
+| **Phase quick-wins** | P1: 9g.1 + 9h.1 + 9h.2 + 9d.3 — Haiku → Opus review | None — start now |
+| **Research: Decorated DCE** | P2: Phase 12 — Sonnet → Opus review | None — start now |
+| **Circle pattern + convergence** | P3: 9e + 9d.4 + 9g.2 — Sonnet → Opus review | None — start now |
+| **Research: Analytic Hessian** | P4: Phase 9b-analytic — Opus | ⏸ awaiting reviewer Q3 |
+| **Java port: genus > 1** | P5: Phase 9c — Opus | ⏸ awaiting reviewer Q4 + G0 |
+| **CGAL packaging** | S6 (audit system) — Opus | ⛔ blocked on G0 (author reply pending) |
+
+---
+
 ## ◼ Porting complete — Phases 1–7
 
 ```
@@ -107,14 +127,16 @@ mesh type.
          Status: 🟡 PR #9 open, 7 tests passing, ~96× speed-up measured.
 
 9b-analytic   Full analytic HyperIdeal Hessian via Schläfli identity
-         → planned, see research-track.md
+         → planned, see research-track.md + session-prompts.md §P4
          Mathematical source: Springborn 2020 §4 + Schläfli 1858/60
                               + Rivin, Schlenker 1999 "The Schläfli formula in
                                 Einstein manifolds with boundary" (ERA-AMS 5, 18–23)
                               + Cho-Kim 1999 + Glickenstein 2011 §4
          Algorithm: explicit chain rule through (bᵢ,aₑ) → ℓᵢⱼ → ζ₁₃/ζ₁₄/ζ₁₅ → αᵢⱼ/βᵢ
-         Includes: short LaTeX correctness note in doc/math/.
+         Includes: 805-line LaTeX derivation in doc/math/hyperideal-hessian-derivation.md.
          Effort: 10–14 days net.  Trigger: profiling on V > 5000.
+         Complexity / scaling context: doc/math/complexity.md.
+         ⏸ GATED: awaiting reviewer Q3 answer ("worth ~2 weeks at your mesh sizes?").
 
 9c — Genus g > 1 fundamental domain  (Java port + research extensions)
 ──────────────────────────────────────────────────────────────────────
@@ -287,6 +309,37 @@ mesh type.
            no new library surface).  Effort: small (~2–3 days).
 ```
 
+9h — CLI usability extensions  (infrastructure — no Java equivalent)
+──────────────────────────────────────────────────────────────────
+
+Identified by usability-audit-2026-05-31 (Finding-U11).  Two independent
+sub-tasks; either can land independently.
+
+```
+9h.1   Newton solver tuning parameters  (~30 min)
+         Currently hardcoded in all three run_*() helpers:
+           tol = 1e-8, max_iter = 200
+         Add to conformallab_cli.cpp:
+           app.add_option("--tol",      tol,      "Newton gradient tolerance [1e-8]");
+           app.add_option("--max-iter", max_iter, "Newton iteration limit    [200]");
+         Thread both through run_euclidean / run_spherical / run_hyper_ideal.
+         Update getting-started.md CLI parameter table.
+         Status: 🔲 planned.  Effort: ~30 min.  No dependencies.
+
+9h.2   Phase-9a models in CLI  (~2–4 hours)
+         The CP-Euclidean and Inversive-Distance functionals are fully
+         implemented in the library and exposed via the CGAL public API
+         (Discrete_circle_packing.h, Discrete_inversive_distance.h), but
+         a CLI user cannot reach them.  Add:
+           -g cp_euclidean       → run_cp_euclidean()
+           -g inversive_distance → run_inversive_distance()
+         Following the existing run_euclidean() pattern (~60 lines each).
+         Add both geometry strings to the CLI::IsMember validator.
+         Update README + getting-started.md parameter table.
+         Status: 🔲 planned.  Effort: ~2–4 h.  Requires: 9a complete ✅.
+         Java reference: none (both functionals are research / literature ports).
+```
+
 ---
 
 ## ◼ New research directions — Phases 10d–10g  (2026 library scan)
@@ -383,6 +436,11 @@ poor input triangulations. Both share the same
 mathematical core (discrete conformal equivalence, Gauss–Bonnet,
 variational principle of Bobenko–Springborn 2004).
 
+Full architecture-level comparison (overlap analysis, adoption rationale,
+scientific value): [`doc/architecture/geometry-central-comparison.md`](../architecture/geometry-central-comparison.md).
+Software landscape (libigl, CGAL, pmp-library, geometry-central):
+[`doc/math/software-landscape.md`](../math/software-landscape.md).
+
 ---
 
 ## ◼ Phase 10 — Genus g ≥ 2 (research with partial Java support)

doc/roadmap/porting-status.md

@@ -1,6 +1,8 @@
 # Porting status overview
 
-> **Snapshot date:** 2026-05-22 (commit graph at v0.9.0 + 2 open PRs).
+> **Snapshot date:** 2026-05-31 (post audit-sessions S1 + S2 on branch
+> `fix/b1-v3-c1-quick-wins`; 9 commits; 301/301 CGAL tests green).
+> For the authoritative live count see [`doc/api/tests.md`](../api/tests.md).
 >
 > **Audience.** External collaborators evaluating whether to use, extend,
 > or contribute to conformallab++.  This document is the **operational
@@ -55,8 +57,8 @@ legacy API (`code/include/*.hpp`) and the CGAL public API
 | **CP-Euclidean** (BPS) | face circles  | **face**     | CPEuclidean 260 LoC | analytic 2×2-per-edge  | ✅        | `discrete_circle_packing_euclidean`         | ⛔ N/A     |
 | **Inversive Distance** | vertex circles| vertex       | ❌ no Java (Luo 2004 + Glickenstein 2011 from literature) | FD (analytic planned) | ✅ | `discrete_inversive_distance_map`           | ⛔ pending |
 
-Total: 250+ tests covering all five models, 0 skipped.  Per-suite
-breakdown: [`doc/api/tests.md`](../api/tests.md).
+Total test count: see [`doc/api/tests.md`](../api/tests.md) — single source
+of truth (counts change as sessions land; do not hardcode them here).
 
 ### What "UV out" means
 
@@ -99,6 +101,10 @@ into `pmap` — no separate user code needed.  See
 | Newton with line search                         | ✅     | `newton_solver.hpp`, all five solvers |
 | SimplicialLDLT + SparseQR fallback              | ✅     | gauge-singular meshes handled automatically |
 | Block-FD Hessian framework                      | ✅     | shipped for HyperIdeal (Phase 9b); 96× speed-up |
+| `newton_core` refactor — single exit path       | ✅     | shipped S1 (2026-05-31); prerequisite for clean diagnostic propagation |
+| `NewtonStatus` enum (`Converged` / `MaxIterations` / `LinearSolverFailed` / `LineSearchStalled`) | ✅ | shipped S2 (2026-05-31); propagated into all three CGAL result types via `CGAL::Newton_status` alias |
+| Solver diagnostics: `sparse_qr_fallback_used`, `min_ldlt_pivot` | ✅ | shipped S2 (2026-05-31); available on all public CGAL result types |
+| Selectable Newton clamp mode (`HardJava` \| `SmoothBarrier`) | ✅ | shipped S1 (2026-05-31); `HardJava` is the default (Java-parity) |
 | Analytic Hessian via Schläfli                   | 🔲    | Phase 9b-analytic; derivation in [`hyperideal-hessian-derivation.md`](../math/hyperideal-hessian-derivation.md) |
 
 ---
@@ -216,8 +222,12 @@ These are conformallab++ contributions beyond porting — the
 | Genus-2 test mesh + brezel2.obj scalability  | ✅ shipped |
 | Memory-safe layout via `halfedge_uv` storage | ✅ shipped |
 | `doc/release-policy.md` formal release policy | ✅ shipped (PR #13) |
-| Analytic HyperIdeal Hessian via Schläfli     | 🔲 derivation written (`hyperideal-hessian-derivation.md`); implementation Phase 9b-analytic |
 | Output UV map integrated into wrappers       | ✅ shipped (PR #14) for 3 of 5 models |
+| `[[deprecated]]` aliases for pre-S1 API names (A1–A3 rename → `<verb>_<geom>_<rest>`) | ✅ shipped S1 (2026-05-31) |
+| Named numeric constants in `constants.hpp` (all thresholds / FD steps / guard values) | ✅ shipped S1 (2026-05-31) |
+| Kahan-compensated triangle area in `enforce_gauss_bonnet` | ✅ shipped S1 (2026-05-31) |
+| `NewtonStatus` enum + solver diagnostics on all public CGAL result types | ✅ shipped S2 (2026-05-31) |
+| Analytic HyperIdeal Hessian via Schläfli     | 🔲 derivation written (`hyperideal-hessian-derivation.md`); implementation Phase 9b-analytic (P4) |
 | Full uniformisation for genus g ≥ 2           | 🔲 Phase 10c — fully new research |
 
 ---

doc/roadmap/research-track.md

@@ -7,6 +7,15 @@
 > port-tracking sheet `doc/roadmap/java-parity.md` so that the porting
 > work and the research work can be planned independently.
 >
+> **Companion documents:**
+> - [`doc/math/novelty-statement.md`](../math/novelty-statement.md) — why these
+>   contributions are novel and who the target audience is.
+> - [`doc/math/software-landscape.md`](../math/software-landscape.md) — how
+>   conformallab++ relates to libigl, geometry-central, and CGAL (relevant for
+>   deciding research vs. duplication at the boundary cases).
+> - [`phase-orchestration.md`](phase-orchestration.md) — model assignments and
+>   session prompts for implementing items in this document.
+>
 > **Created:** 2026-05-21, after a full doc audit that identified four
 > items previously mislabelled as "ports".  This document corrects the
 > record and extends it with the explicit research plan for Phase

doc/roadmap/session-prompts.md

@@ -0,0 +1,208 @@
+# Ready-to-paste session prompts (P1–P4 + 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 phase(s), the detail doc to follow,
+the build/test commands, and the branch/push/PR + review-gate workflow.
+
+Shared conventions (baked into each prompt):
+- Repo root: `/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 command (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\.'`
+- Phase details: `doc/roadmap/phases.md` (per-phase plan);
+  `doc/roadmap/research-track.md` (research items with acceptance criteria).
+- Build flags reference: `CLAUDE.md` §Build commands (canonical source; use
+  `-DCONFORMALLAB_LOW_MEMORY_BUILD=ON -j1` if on the Raspberry Pi runner).
+- After each implementation session, run the **Review gate** prompt (Opus).
+
+---
+
+## P1 — Quick wins  (model: **Haiku**)
+
+```
+Use the Haiku model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new
+branch off main called `feat/p1-quick-wins`.
+
+Implement four independent additions — full details (Java references, math
+references, acceptance criteria) in doc/roadmap/phases.md at the sections
+labelled 9h.1, 9h.2, 9g.1, 9d.3:
+
+- 9h.1  Add --tol and --max-iter CLI options in code/src/conformallab_cli.cpp.
+  Thread both through run_euclidean / run_spherical / run_hyper_ideal.
+  Update doc/getting-started.md CLI parameter table.
+
+- 9h.2  Add -g cp_euclidean and -g inversive_distance routes in the CLI,
+  following the existing run_euclidean() pattern (~60 lines each).
+  Add both geometry strings to the CLI::IsMember validator.
+  Update README + doc/getting-started.md.
+
+- 9g.1  Create code/include/conformal_quality.hpp implementing:
+    IsothermicityMeasure, DiscreteConformalEquivalenceMeasure,
+    FlippedTriangles, LengthCrossRatio, ConvergenceUtility measures.
+  Java source classes are listed in phases.md §9g.1.
+  Each function must have at least one sanity test in code/tests/cgal/
+  (e.g. FlippedTriangles returns 0 on a valid layout; LengthCrossRatio
+  is 1.0 on an equilateral triangle). Register in code/tests/cgal/CMakeLists.txt.
+
+- 9d.3  Create code/include/stereographic_layout.hpp porting
+  StereographicUnwrapper.java (266 LoC). See phases.md §9d.3 for the math
+  (stereographic projection S²→ℂ∪{∞} + Möbius centring for genus-0 surfaces).
+  Add at least one round-trip test (north pole → ∞; a unit-sphere point →
+  expected complex value).
+
+Run the full CGAL suite after all four are implemented:
+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.
+
+Commit per phase (or in two logical commits) with trailer
+`Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com>`.
+Push to origin and open a PR (base main) via the Gitea API using the token
+in the origin remote URL. Update doc/roadmap/phase-orchestration.md (mark each
+completed phase ✅ with the commit ref). Report the PR URL and test count.
+```
+
+---
+
+## P2 — Decorated DCE transition  (model: **Sonnet**)
+
+```
+Use the Sonnet model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new
+branch off main called `feat/p2-decorated-dce`.
+
+Implement Phase 12 — Decorated DCE & geometric transition.
+Full details and acceptance criteria in:
+  doc/roadmap/phases.md §Phase 12
+  doc/roadmap/research-track.md §Phase 12
+
+Mathematical reference: Bobenko, Lutz 2025 "Decorated Discrete Conformal
+Equivalence in Non-Euclidean Geometries" (Discrete & Comput. Geom.;
+arXiv:2310.17529) §3 — Penner-coordinate decoration unifying the three
+background geometries.
+
+Scope (from phases.md):
+1. Decoration layer — per-vertex circle/horocycle radius as Penner coordinate;
+   implement the map ↔ existing inversive distance I_ij via
+   ℓ² = r_i² + r_j² + 2 r_i r_j η.
+   → Create code/include/decorated_dce.hpp.
+2. Transition driver — deform background curvature κ ∈ {+,0,−} while holding
+   the discrete conformal invariant fixed; call the three existing solvers
+   (euclidean, spherical, hyper_ideal).
+3. Validation harness — code/tests/cgal/test_decorated_dce.cpp.
+
+All four acceptance criteria from phases.md must be met:
+- Decoration round-trip I_ij ↔ (r_i, r_j, ℓ) at machine precision.
+- At κ=0: bit-for-bit match with existing euclidean / inversive path.
+- Gauss-Bonnet holds per geometry across the κ-transition.
+- Invariant (I_ij) is constant across the three-geometry transition to tol.
+
+Run the full CGAL suite (must stay green with new tests).
+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/roadmap/phase-orchestration.md (Phase 12 → ✅, session P2 + commit ref).
+Report the PR URL.
+```
+
+---
+
+## P3 — Circle pattern embedding + Möbius centring + convergence study  (model: **Sonnet**)
+
+```
+Use the Sonnet model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new
+branch off main called `feat/p3-circle-pattern-convergence`.
+
+Three items — full details in doc/roadmap/phases.md:
+
+- 9e  Create code/include/circle_pattern_layout.hpp porting
+  CirclePatternLayout + CirclePatternUtility. Phase 9a.1 (the CP-Euclidean
+  energy + solver) is a prerequisite and is already landed on main.
+  Java references: unwrapper/circlepattern/{CirclePatternLayout,
+  CirclePatternUtility,CPEuclideanRotation}.java (phases.md §9e).
+  Required test: given ρ values from a solved CP-Euclidean system, verify
+  that the embedded vertex positions are self-consistent (each face's three
+  circle-intersection points form the correct intersection angles to tol).
+
+- 9d.4  Upgrade normalise_hyperbolic() in code/include/layout.hpp to use the
+  variational MobiusCenteringFunctional (Lorentz energy
+  E = Σ log(-⟨x,p⟩/√(-⟨x,x⟩)), gradient + Hessian).
+  Java reference: functional/MobiusCenteringFunctional.java (289 LoC).
+  Retain the existing Fréchet mean as a fallback if Newton fails.
+  Required test: compare old vs new centring output on brezel.obj; both
+  must place the centroid within tol of the origin.
+
+- 9g.2  Add code/tests/cgal/test_period_matrix_convergence.cpp (experiment,
+  not a library feature — see phases.md §9g.2):
+  Generate a genus-1 elliptic mesh with a known analytic τ; subdivide via
+  igl::loop; add per-vertex Gaussian noise; compute |τ_discrete − τ_analytic|
+  after each step. Assert that the residual decreases monotonically with
+  refinement (the discrete period matrix converges).
+
+Run the full CGAL suite (must stay green).
+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/roadmap/phase-orchestration.md (9e/9d.4/9g.2 → ✅, P3 + commit ref).
+Report the PR URL and the convergence-study output.
+```
+
+---
+
+## P4 — Analytic HyperIdeal Hessian  (model: **Opus**) — ⏸ GATED on reviewer Q3
+
+```
+PRECONDITION: do NOT start this session until the reviewer has answered Q3
+("Is the ~6× speedup over block-FD worth ~2 weeks at your typical mesh sizes?").
+  If the answer is "no" or "not a priority", stop — Phase 9b-analytic stays ⏸.
+  If the answer is "yes" or "above V > X", proceed.
+
+Use the Opus model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new
+branch off main called `feat/p4-analytic-hessian`.
+
+Implement Phase 9b-analytic — the full analytic HyperIdeal Hessian via the
+Schläfli identity. The complete chain-rule derivation (805-line LaTeX note) and
+the implementation plan are in:
+  doc/math/hyperideal-hessian-derivation.md
+  doc/roadmap/phases.md §9b-analytic
+  doc/roadmap/research-track.md §Phase 9b-analytic
+
+Chain: (bᵢ, aₑ) → ℓᵢⱼ → ζ₁₃/ζ₁₄/ζ₁₅ → αᵢⱼ/βᵢ → ∂²E/∂u².
+Replace the block-FD path in code/include/hyper_ideal_hessian.hpp with the
+analytic Hessian. Retain the block-FD path available as a compile-time flag
+(-DCONFORMALLAB_HYPER_IDEAL_FD_CHECK or runtime enum) for cross-validation.
+
+Acceptance criteria:
+- All existing HyperIdeal golden-value tests pass bit-for-bit (HardJava clamp).
+- New test: analytic and block-FD Hessians agree to 1e-6 on the tetrahedron.
+- Benchmark: measure the analytic vs block-FD wall-time on the largest test
+  mesh; report the ratio. Analytic must be faster for V > 500.
+
+Run the full CGAL suite (must stay green). Commit; push; open PR via Gitea API.
+Update doc/roadmap/phase-orchestration.md (9b-analytic → ✅, P4 + commit ref).
+Report PR URL and the measured speed ratio.
+```
+
+---
+
+## Review gate (run after P1 / P2 / P3)  (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 count regression
+  (ctest --test-dir build-cgal -R '^cgal\.'); count matches doc/api/tests.md.
+- Any new functional has a gradient-check test (pattern in CLAUDE.md §Test design patterns).
+- No Java golden-vector / parity test perturbed; HardJava clamp default intact.
+- Numeric changes are value-identical where claimed, or justified + covered by a test.
+- New public surface (result types, enums, CGAL headers) is intentional and
+  documented in doc/api/headers.md and doc/api/contracts.md.
+- Commit messages attribute the implementing model (Co-Authored-By trailer).
+- Phase(s) marked ✅ in doc/roadmap/phase-orchestration.md with the commit ref.
+
+Read the actual diff (git diff main...HEAD -- code/include/ code/tests/ doc/)
+and the phases.md entry for the phase(s) involved. 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.
+```

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
-# with the reviewer (e.g. 80 %), the gate can be a single line in
-# cpp-tests.yml.
+# Local + CI.  Threshold gate: 80 % line / 70 % branch / 90 % function.
+# Set SKIP_COVERAGE_GATE=1 to measure without failing (ramp-up mode).
+# 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
-    lcov --summary "$BUILD_DIR/coverage.info" "${LCOV_TOLERANT[@]}" 2>/dev/null \
-        | grep -E "lines\.\.\.\.|functions|branches" \
-        | sed 's/^/    /'
-    echo
-    echo "HTML report:  $BUILD_DIR/lcov-html/index.html"
-    echo "              open $BUILD_DIR/lcov-html/index.html"
-else
-    echo "    WARN: coverage.info is empty — likely an lcov/gcov version"
-    echo "          mismatch.  Tests passed; raw .gcda files are in $BUILD_DIR."
-    echo "          Inspect with:  find $BUILD_DIR -name '*.gcda' | head"
+if [ ! -s "$BUILD_DIR/coverage.info" ]; then
+    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 "HTML report:  $BUILD_DIR/lcov-html/index.html"
+echo "              open $BUILD_DIR/lcov-html/index.html"
+echo
+
+# ── Coverage threshold gate (C3) ─────────────────────────────────────────────
+# 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."