ConformalLabpp/doc/reviewer/external-audit-2026-05-30.md

External Code Audit — ConformalLabpp v0.10.0

Date: 2026-05-30
Auditor: External reviewer (Claude Sonnet 4.6)
Branch: review/external-audit-2026-05-30
Base: docs/fix-test-count-post-merge (HEAD at audit time)

This document is self-contained. A new session can pick up any finding below and act on it without prior context. Each finding includes:

• exact file path + line numbers (verified by direct file read)
• a minimal reproduction of the problematic code
• the correct fix or recommended action
• acceptance criteria for "done"

Status legend: 🔴 Bug · 🟡 API/Doc error · 🟠 Test gap · 🔵 Architectural risk

---

How to read this document in a new session

# 1. Check out the audit branch
git checkout review/external-audit-2026-05-30

# 2. Build the CGAL test suite (needed for verification)
cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON
cmake --build build-cgal --target conformallab_cgal_tests -j$(nproc)

# 3. Run the full suite before making any change
ctest --test-dir build-cgal -R '^cgal\.' --output-on-failure
# Expected: 246 passed, 0 failed

# 4. Pick a finding below, apply the fix, re-run ctest, then commit.

The Java reference implementation lives at:

/Users/tarikmoussa/Desktop/conformallab/src/de/varylab/discreteconformal/

Consult it for any port-faithfulness question.

---

FINDING-A — 🔴 CRITICAL BUG: face_energy() silently wrong for mixed ideal/hyper-ideal configurations

Location

code/include/hyper_ideal_functional.hpp lines 319–344

Problem

The face_energy() function handles only two cases: "all three vertices hyper-ideal" and "exactly one vertex ideal (pinned)". When two or all three vertices are ideal (v?b = m.v_idx[v?] < 0), the cascade falls through to a branch that calls the one-ideal-vertex volume formula — which is mathematically wrong for two or three ideal vertices.

// CURRENT (broken for >= 2 ideal vertices):
static double face_energy(const FaceAngles& fa)
{
    ...
    if (fa.v1b && fa.v2b && fa.v3b) {          // all hyper-ideal  ✓
        V = calculateTetrahedronVolume(
                fa.beta1, fa.beta2, fa.beta3,
                fa.alpha23, fa.alpha31, fa.alpha12);
    } else if (!fa.v1b) {                        // BUG: enters even when !v1b && !v2b
        V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
                fa.beta1, fa.alpha31, fa.alpha12,
                fa.alpha23, fa.beta2, fa.beta3);
    } else if (!fa.v2b) {                        // BUG: enters even when !v2b && !v3b
        V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
                fa.beta2, fa.alpha12, fa.alpha23,
                fa.alpha31, fa.beta3, fa.beta1);
    } else {   // !v3b only                      // only correct for exactly 1 ideal
        V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
                fa.beta3, fa.alpha23, fa.alpha31,
                fa.alpha12, fa.beta1, fa.beta2);
    }
    return aa + bb + 2.0 * V;
}

The same structural error exists in face_angles_from_local_dofs() at lines 192–215:

if (l12 > l23 + l31) {
    o.beta1 = 0.0; o.beta2 = 0.0; o.beta3 = PI;
    o.alpha12 = PI; o.alpha23 = 0.0; o.alpha31 = 0.0;
} else if (l23 > l12 + l31) { ... }
else if (l31 > l12 + l23) { ... }
else { /* normal */ }

(this part is actually correct — no bug here; the degenerate case structure is standard. Keeping the note for completeness.)

Trigger condition

Only triggered when the HyperIdealMaps has some vertices pinned (v_idx[v] = -1) and some variable. The default workflow assign_all_dof_indices(mesh, maps) makes ALL vertices variable (v?b = true everywhere), which always takes the first branch — safe. The bug only fires for mixed ideal/hyper-ideal configurations.

Fix

The missing cases are the fully-ideal (all three ideal) and two-ideal-vertex cases. The Java reference HyperIdealFunctional.java must be consulted to find the correct volume formula for two ideal vertices. The typical fix structure is:

// PROPOSED FIX — verify against Java reference before applying:
if (fa.v1b && fa.v2b && fa.v3b) {
    // all hyper-ideal
    V = calculateTetrahedronVolume(
            fa.beta1, fa.beta2, fa.beta3,
            fa.alpha23, fa.alpha31, fa.alpha12);
} else if (fa.v1b && fa.v2b && !fa.v3b) {
    // ideal at v3 only
    V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
            fa.beta3, fa.alpha23, fa.alpha31,
            fa.alpha12, fa.beta1, fa.beta2);
} else if (fa.v1b && !fa.v2b && fa.v3b) {
    // ideal at v2 only
    V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
            fa.beta2, fa.alpha12, fa.alpha23,
            fa.alpha31, fa.beta3, fa.beta1);
} else if (!fa.v1b && fa.v2b && fa.v3b) {
    // ideal at v1 only
    V = calculateTetrahedronVolumeWithIdealVertexAtGamma(
            fa.beta1, fa.alpha31, fa.alpha12,
            fa.alpha23, fa.beta2, fa.beta3);
} else if (!fa.v1b && !fa.v2b && fa.v3b) {
    // two ideal: v1, v2  — FORMULA NEEDED (see Java reference)
    V = 0.0; // TODO: implement two-ideal-vertex formula
} else if (!fa.v1b && fa.v2b && !fa.v3b) {
    // two ideal: v1, v3  — FORMULA NEEDED
    V = 0.0; // TODO
} else if (fa.v1b && !fa.v2b && !fa.v3b) {
    // two ideal: v2, v3  — FORMULA NEEDED
    V = 0.0; // TODO
} else {
    // all three ideal: volume = sum of Lobachevsky only
    V = 0.0; // TODO: verify correct formula
}

Check Java: HyperIdealFunctional.java → triangleEnergyAndAlphas() for all cases.

Resolution (2026-05-30)

Root cause clarified: The Java reference (HyperIdealFunctional.java lines 219–233) has the identical cascade structure — it silently applies the one-ideal-vertex formula to the first ideal vertex found, ignoring additional ideal vertices in the same face. The C++ was a faithful port; it did not introduce a new divergence.

The fix therefore does NOT change the calculation (which would diverge from Java). Instead it:

1. Adds <stdexcept> include to hyper_ideal_functional.hpp.
2. Adds an ideal_count guard at the top of face_energy() that counts ideal vertices (!v?b) and throws std::logic_error for ideal_count >= 2, replacing the silent wrong result with a clear diagnostic message that names the limitation and points to the audit document.
3. Adds three new GTest cases to test_hyper_ideal_functional.cpp:
   • MultiIdealGuard_TwoIdealVertices_Throws — two ideal vertices → throw ✅
   • MultiIdealGuard_AllThreeIdealVertices_Throws — all ideal → throw ✅
   • MultiIdealGuard_ExactlyOneIdeal_DoesNotThrow — one ideal → no throw ✅
4. Adds a doc comment above face_energy() explaining the supported configurations (0-ideal / 1-ideal), the Java reference limitation, and why multi-ideal faces are not implemented (requires new research beyond the Java reference).

Test result: 262/262 CGAL tests pass (was 246 before + 3 new guard tests + 13 from prior suite growth). No regressions.

Further work: Implementing the correct volume formulas for 2-ideal and 3-ideal faces is a separate research item (not a port — Java does not have them either). It is now tracked in doc/roadmap/research-track.md under "Hyper-ideal volume formulas for 2- and 3-ideal-vertex faces (Phase 9b+)". The throw remains the correct safe behaviour until that research item is resolved.

Acceptance criteria

• All 8 cases (2³ combinations) are either correct or throw explicitly
• Three new tests cover the three guard cases
• 262 CGAL tests pass, 0 failed

---

FINDING-B — 🟡 API CONCEPTUAL ERROR: Gauss–Bonnet check silently wrong for HyperIdeal

Location

code/include/gauss_bonnet.hpp lines 87–88, 128–134

Problem

The function gauss_bonnet_sum() is overloaded for all five map types including HyperIdealMaps, and check_gauss_bonnet() is templated so it accepts any Maps:

// gauss_bonnet.hpp:87
inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps& mp)
    { return gauss_bonnet_sum(m, mp.theta_v); }

// gauss_bonnet.hpp:128-134
template <typename Maps>
inline void check_gauss_bonnet(const ConformalMesh& mesh, const Maps& maps,
                                double tol = 1e-8)
{
    check_gauss_bonnet(mesh, gauss_bonnet_sum(mesh, maps), tol);  // always throws for HyperIdeal!
}

The Gauss–Bonnet identity checked here is:

Σ_v (2π − Θ_v) = 2π · χ(M)     ← Euclidean / flat case only

For a hyperbolic metric (which is what HyperIdeal computes), the correct identity is:

Σ_v (2π − Θ_v) − Area(M) = 2π · χ(M)

For a regular (no cone singularities) genus-2 surface with Θ_v = 2π for all v:

• LHS: Σ(2π − 2π) = 0
• RHS expected by check: 2π·χ = 2π·(2−2·2) = −4π

So |0 − (−4π)| = 4π ≫ tol → check_gauss_bonnet always throws for valid hyperbolic configurations. The overload exists but calling it produces a wrong result.

Fix options

Option A (recommended): Delete the HyperIdealMaps overload of gauss_bonnet_sum and add a compile-time or doc-level guard that prevents check_gauss_bonnet from being instantiated with HyperIdealMaps. Add a comment explaining why.

// DELETE this overload:
// inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps& mp)

// ADD to check_gauss_bonnet doc:
// NOTE: Do NOT call with HyperIdealMaps — the hyperbolic Gauss–Bonnet identity
// includes an Area term absent from this check. For HyperIdeal configurations
// there is no simple global angle constraint analogous to the Euclidean case.

Option B: Keep the overload but make it compute the correct hyperbolic quantity Σ(2π−Θ_v) − Area(M) and add a separate check_hyperbolic_gauss_bonnet that compares against 2π·χ. This requires computing the area from the HyperIdeal metric, which is non-trivial.

Resolution (2026-05-31)

1. gauss_bonnet_sum(mesh, HyperIdealMaps) deleted — replaced with = delete overload and a multi-line comment explaining the Area-term discrepancy. Attempting to call this is now a compile error.
2. enforce_gauss_bonnet(mesh, HyperIdealMaps&) deleted — explicit = delete overload prevents the generic template from being silently instantiated with HyperIdealMaps.
3. Header comment block rewritten — now has a clear box explaining which geometries are supported (Euclidean/Spherical) and which are not (HyperIdeal), with the correct hyperbolic Gauss–Bonnet identity shown.
4. One new GTest in test_phase6.cpp:
   • HyperIdeal_EuclideanSumDiscrepancy_DocumentsWhyCheckIsDeleted — verifies numerically that for a regular (Θ_v=2π) tetrahedron the Euclidean sum = 0 while 2π·χ = 4π, i.e. the Euclidean check would produce deficit = −4π for a valid HyperIdeal target.
5. Three compile-time static_asserts in the test file (SFINAE-based) confirm that gauss_bonnet_sum is NOT invocable with HyperIdealMaps but IS invocable with EuclideanMaps and SphericalMaps.

Test result: 263/263 CGAL tests pass. No regressions.

Acceptance criteria

• Calling check_gauss_bonnet(mesh, hyper_ideal_maps) → compile error
• Calling enforce_gauss_bonnet(mesh, hyper_ideal_maps) → compile error
• Comments explain why Euclidean G-B does not apply to HyperIdeal
• 263 CGAL tests pass, 0 failed

---

FINDING-C — 🟡 DOCUMENTATION BUG: Cotangent formula in header comment is wrong

Location

code/include/euclidean_hessian.hpp line 26–27 (box comment) and line 57–58 (comment above euclidean_cot_weights)

Problem

The box comment at the top of the file states:

│    cot_k = (t_adj1·l123 − t_adj2·t_opp) / denom2                        │
│          = cotangent of the angle αk at vertex k

And the comment above euclidean_cot_weights:

// cot_k = (t_adj·l123 − t_opp·t_other) / (8·Area)

Both are mathematically wrong. The correct formula (verified numerically) is:

cot_k = (t_opp · l123 − t_adj1 · t_adj2) / (8·Area)

Numerical proof (3-4-5 right triangle, l23=3, l31=4, l12=5):

• t12=2, t23=6, t31=4, l123=12, 8·Area=48
• cot(α₁) = 4/3 (angle at v₁ opposite l23=3)
• Code result: (t23·l123 − t31·t12)/48 = (6·12 − 4·2)/48 = 64/48 = 4/3 ✓
• Header formula: (t_adj1·l123 − t_adj2·t_opp)/48 = (t12·l123 − t31·t23)/48 = (2·12 − 4·6)/48 = 0 ✗

The implementation in euclidean_cot_weights (lines 91–96) is correct. Only the documentation is wrong.

Fix

// REPLACE both comment instances with the correct formula:

// cot_k = (t_opp · l123 − t_adj1 · t_adj2) / (8·Area)
//
// where for vertex k:
//   t_opp  = t-value of the edge OPPOSITE to k   (= 2(s − l_opp))
//   t_adj1, t_adj2 = t-values of the two edges ADJACENT to k
//   l123   = l12 + l23 + l31  (perimeter)
//   8·Area = denom2

Acceptance criteria

• Both comment blocks corrected (box comment + function-level comment)
• A one-line comment in the function body confirms which t-value is t_opp for each return value, e.g.: // cot1: t_opp=t23 (opposite v1)
• No code changes — only documentation

---

FINDING-D — 🟡 DOCUMENTATION BUG: DOF-assignment functions claim "pin before" works, but it doesn't

Location

• code/include/euclidean_functional.hpp lines 97–107 (assign_euclidean_vertex_dof_indices)
• code/include/spherical_functional.hpp lines 90–96 (assign_vertex_dof_indices)
• code/include/inversive_distance_functional.hpp lines 123–139 (assign_inversive_distance_vertex_dof_indices)

Problem

All three functions iterate unconditionally over every vertex and assign a sequential index, overwriting any previously set -1 pin:

// euclidean_functional.hpp:104-107
inline int assign_euclidean_vertex_dof_indices(ConformalMesh& mesh, EuclideanMaps& m)
{
    int idx = 0;
    for (auto v : mesh.vertices()) m.v_idx[v] = idx++;   // overwrites ALL
    return idx;
}

The doc comment above this function says:

/// **Note:** does NOT pin a gauge vertex.  For closed meshes the caller
/// must set one `m.v_idx[v] = -1` either before or after this call

"Before" is wrong — the loop overwrites any pre-set pin.

The inversive-distance version is worse:

// inversive_distance_functional.hpp:126-130
/// 1. set one `m.v_idx[v] = -1` *before* calling this function (then
///    the call is a no-op for that vertex) — OR —

Explicitly claims the call is "a no-op" for a pre-pinned vertex, which is false.

Consequence

A caller who pins m.v_idx[first_vertex] = -1 and then calls assign_euclidean_vertex_dof_indices() will get a fully-free system with no gauge fix. On a closed mesh the Hessian is singular (gauge mode). SimplicialLDLT silently falls back to SparseQR, which finds a minimum-norm solution — no error is reported, but the result is not the intended pinned solution.

Fix

Option A (minimal): Fix the doc comments only. Remove "before" as an option; only "after" works.

/// NOTE: this function assigns indices to ALL vertices unconditionally.
/// To pin a gauge vertex, set `m.v_idx[v] = -1` AFTER calling this function.

Option B (API improvement): Add an overload that accepts a gauge vertex:

inline int assign_euclidean_vertex_dof_indices(
    ConformalMesh& mesh, EuclideanMaps& m, Vertex_index gauge)
{
    int idx = 0;
    for (auto v : mesh.vertices())
        m.v_idx[v] = (v == gauge) ? -1 : idx++;
    return idx;
}

Acceptance criteria

• Doc comments corrected in all three files to say "AFTER" only
• Optionally: overload with explicit gauge vertex added
• No code changes required for correctness (existing callers set pin after, which already works)

---

FINDING-E — 🟡 INCONSISTENCY: gradient_check_cp_euclidean uses absolute error, all others use relative

Location

code/include/cp_euclidean_functional.hpp lines 338–349 (gradient_check_cp_euclidean) and lines 373–387 (hessian_check_cp_euclidean)

Problem

Every gradient check in the library normalises the error by the gradient magnitude:

// euclidean_functional.hpp:343 — RELATIVE error
double scale = std::max(1.0, std::abs(G[si]));
if (err / scale > tol) ok = false;

But gradient_check_cp_euclidean uses absolute error:

// cp_euclidean_functional.hpp:345 — ABSOLUTE error (inconsistent)
if (std::abs(G[i] - fd) > tol) {
    std::cerr << "[cp-euclidean] FD gradient mismatch ...";
    return false;
}

Same issue in hessian_check_cp_euclidean line 378:

if (std::abs(an - fd) > tol) {   // absolute, not relative

The default tol = 1e-6 is acceptable for unit-scale problems but will produce false failures if the gradient values grow large (e.g., many faces, large ρ).

Fix

// gradient_check_cp_euclidean — replace the comparison:
double err   = std::abs(G[i] - fd);
double scale = std::max(1.0, std::abs(G[i]));
if (err / scale > tol) {
    std::cerr << "[cp-euclidean] FD gradient mismatch at DOF " << i
              << ": analytic=" << G[i] << "  FD=" << fd
              << "  rel-err=" << (err / scale) << "\n";
    return false;
}

// hessian_check_cp_euclidean — same pattern:
double err   = std::abs(an - fd);
double scale = std::max(1.0, std::abs(an));
if (err / scale > tol) { ... }

Acceptance criteria

• Both gradient_check_cp_euclidean and hessian_check_cp_euclidean use relative error (normalised by max(1.0, |analytic|))
• The existing CP-Euclidean gradient check tests still pass

---

FINDING-F — 🟠 TEST GAP: Degenerate-triangle gradient (Finding 1 from java-port-audit.md)

Location

code/tests/cgal/test_euclidean_functional.cpp and code/tests/cgal/test_spherical_functional.cpp

Problem (from java-port-audit.md item 1, still open)

The degenerate-triangle fix (Finding 1 in java-port-audit.md) made euclidean_angles() and spherical_angles() return the limiting angles (π opposite the over-long edge, 0/0 for the others) instead of {0,0,0}. This change is untested: no test builds a triangle with one edge longer than the sum of the others and asserts the π corner.

Required test

// Add to test_euclidean_functional.cpp:
TEST(EuclideanGeometry, DegenerateTriangle_LimitingAngles) {
    // l12 = 10, l23 = 1, l31 = 1 → l23+l31=2 < l12=10 → degenerate
    // Expected: alpha3 = π (vertex v3 opposite l12), alpha1=alpha2=0
    auto fa = euclidean_angles_from_lengths(10.0, 1.0, 1.0);
    EXPECT_FALSE(fa.valid);
    EXPECT_NEAR(fa.alpha3, conformallab::PI, 1e-12);
    EXPECT_NEAR(fa.alpha1, 0.0, 1e-12);
    EXPECT_NEAR(fa.alpha2, 0.0, 1e-12);
}

// Similarly for spherical_angles in test_spherical_functional.cpp

Acceptance criteria

• Test added for Euclidean degenerate triangle → alpha3 = π
• Test added for Spherical degenerate triangle → same
• Both tests pass

---

FINDING-G — 🟠 TEST GAP: euclidean_hessian edge-DOF guard must throw (Finding 2 from java-port-audit.md)

Location

code/tests/cgal/test_euclidean_hessian.cpp (or new file)

Problem (from java-port-audit.md item 2, still open)

Finding 2 added a throw std::logic_error guard to euclidean_hessian() when edge DOFs are present. No test asserts this throw fires.

Required test

TEST(EuclideanHessian, EdgeDOFGuard_Throws) {
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    assign_euclidean_all_dof_indices(mesh, maps);  // assigns edge DOFs
    std::vector<double> x(euclidean_dimension(mesh, maps), 0.0);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    EXPECT_THROW(
        euclidean_hessian(mesh, x, maps),
        std::logic_error
    );
}

Acceptance criteria

• Test added and passes (throw confirmed)
• Test is in the cgal suite (build with -DWITH_CGAL_TESTS=ON)

---

FINDING-H — 🟠 TEST GAP: Missing end-to-end torus with Re(τ) < 0 before reduction (java-port-audit.md item 7)

Location

code/tests/cgal/test_pipeline.cpp or code/tests/cgal/test_phase7.cpp

Problem (from java-port-audit.md item 7, partly open)

compute_period_matrix now calls normalizeModulus which maps τ into the half-strip 0 ≤ Re(τ) ≤ ½. This is tested at the function level by the Java oracle (NormalizeModulus_GoldenJava). But there is no end-to-end test where the initial τ has Re(τ) < 0 and the pipeline is verified to fold it into Re(τ) ≥ 0.

If someone reverts compute_period_matrix to call reduce_to_fundamental_domain instead of normalizeModulus, the function-level test still passes but the pipeline output would silently diverge from the Java oracle.

Required test

A torus whose geometry produces Re(τ) < 0 before SL(2,ℤ) reduction. One approach: construct a torus with an asymmetric lattice (e.g., parallelogram with obtuse angle on the left side) where the natural τ has negative real part.

TEST(PeriodMatrix, EndToEnd_NegativeReTau_FoldedToPositive) {
    // Build or load a torus mesh whose natural τ has Re < 0.
    // Run the full pipeline: newton_euclidean → layout → compute_period_matrix.
    // Assert: Re(result.tau) >= 0  (normalizeModulus was applied)
    // Assert: Im(result.tau) >= 0
    // Assert: |result.tau| >= 1
}

Acceptance criteria

• End-to-end test added that exercises a mesh with Re(τ) < 0 pre-reduction
• Test asserts Re(τ) ≥ 0 after compute_period_matrix
• 246 CGAL tests still pass with new test included

---

FINDING-I — 🔵 ARCHITECTURAL RISK: 246/272 CGAL tests are not gated in CI

Location

.gitea/workflows/cpp-tests.yml line (see if: false block), CLAUDE.md lines 267–274

Problem

The CI pipeline has three jobs:

1. test-fast — 26 pure-math tests (no CGAL), active
2. test-cgal — 246 CGAL tests, disabled (if: false since 2026-05-26)
3. quality-gates — structural linting, active

This means the Newton solvers, layout, holonomy, period matrix, cut graph, CGAL public API, and all five DCE models are not regression-tested on any push.

Root cause: the CGAL build OOMs on the 1.6 GB ARM64 Raspberry Pi runner at -j parallel compilation.

Recommended actions (in priority order)

1. Immediate (low risk): Run test-cgal with -j1 (serial build) to avoid OOM. The wall time increases but correctness is not compromised.

   cmake --build build-cgal --target conformallab_cgal_tests -j1
   

2. Short term: Split the CGAL test binary into subsets so a failing compilation is localised. Add a minimal subset (e.g. Newton + gradient checks only) as a new CI job that runs on every push.

3. Medium term: Add a GitHub Actions job (arm64 runner, 7 GB RAM) to mirror CI. Self-hosted Raspberry Pi is not suitable for a library targeting CGAL submission.

4. Document the gap explicitly in CHANGELOG.md and doc/release-policy.md: "v0.10.0 ships with CGAL CI disabled — run ctest -R '^cgal\.' locally before tagging any release."

Acceptance criteria

• At least one of the above options implemented
• The CGAL test suite runs in CI on every PR (not just locally)
• CHANGELOG.md documents the current CI limitation

---

MINOR FINDINGS (quick fixes, no architectural impact)

MINOR-1 — spherical_gauge_shift misleading comment [spherical_functional.hpp:404]

// CURRENT (wrong):
// f is strictly monotone decreasing (second derivative < 0) for a convex functional

// FIX:
// f is strictly monotone decreasing (sum of vertex angles increases with scale)
// for any conservative gradient, including the concave spherical energy.

MINOR-2 — spherical_gauge_shift uses forward FD, should use central [spherical_functional.hpp:470–471]

// CURRENT (forward difference, O(ε)):
double ftp = sum_Gv(t + fd_eps);
double dft = (ftp - ft) / fd_eps;

// FIX (central difference, O(ε²), same cost in this context):
double dft = (sum_Gv(t + fd_eps) - sum_Gv(t - fd_eps)) / (2.0 * fd_eps);

MINOR-3 — Gauss-Legendre constants duplicated in three files

gl_s[10] and gl_w[10] are identically copy-pasted in:

• code/include/euclidean_functional.hpp lines 257–270
• code/include/spherical_functional.hpp lines 304–317
• code/include/inversive_distance_functional.hpp lines 312–325

Fix: Extract to a new header gauss_legendre.hpp:

// code/include/gauss_legendre.hpp (new file)
namespace conformallab::detail {
    inline const double* gl10_nodes() { static const double s[10] = {...}; return s; }
    inline const double* gl10_weights() { static const double w[10] = {...}; return w; }
}

MINOR-4 — hidx() function duplicated in four files

Same one-liner in euclidean_functional.hpp, spherical_functional.hpp, hyper_ideal_functional.hpp, inversive_distance_functional.hpp. Move to conformal_mesh.hpp as:

// code/include/conformal_mesh.hpp (add):
inline std::size_t halfedge_idx(Halfedge_index h) noexcept {
    return static_cast<std::size_t>(static_cast<std::uint32_t>(h));
}

MINOR-5 — inits() in clausen.hpp has unintuitive off-by-one semantics

The function returns n after decrement, which is one less than the last-checked index. Add a comment:

// Returns the index one below the last term needed to reach the requested
// accuracy — callers use this as the `n` argument to csevl(), which then
// evaluates terms [0..n-1].  Matches Java Clausen.inits() semantics exactly.

---

Summary table

ID	File	Lines	Type	Severity	Status
A	hyper_ideal_functional.hpp	319–344	Bug	Critical (mixed config)	✅ Fixed 2026-05-30
B	gauss_bonnet.hpp	87–88, 128–134	API error	Medium	✅ Fixed 2026-05-31
C	euclidean_hessian.hpp	26–27, 57–58	Doc error	Medium	✅ Fixed 2026-05-31
D	euclidean_functional.hpp + 2 others	97–107	Doc error	Medium	🟡 Open
E	cp_euclidean_functional.hpp	338–349, 373–387	Inconsistency	Medium	🟡 Open
F	test files	—	Test gap	Medium	🟠 Open
G	test files	—	Test gap	Medium	🟠 Open
H	test files	—	Test gap	Medium	🟠 Open
I	CI workflow	—	Arch risk	High	🔵 Open
MINOR-1	spherical_functional.hpp	404	Doc error	Minor	🟡 Open
MINOR-2	spherical_functional.hpp	470–471	Accuracy	Minor	🟡 Open
MINOR-3	three files	—	DRY	Minor	🟡 Open
MINOR-4	four files	—	DRY	Minor	🟡 Open
MINOR-5	clausen.hpp	33–38	Doc	Minor	🟡 Open

---

What was verified as correct (no action needed)

The following were carefully audited and found faithful:

• clausen.hpp — Chebyshev expansion, csevl, inits, clausen2, Lobachevsky, ImLi2. All match Java oracle to 1e-12.
• euclidean_geometry.hpp — euclidean_angles() and euclidean_angles_from_lengths(). Degenerate handling correct (π/0/0), centering trick correct.
• spherical_geometry.hpp — spherical_l(), spherical_angles(). Degenerate handling correct. Half-angle formula numerically stable.
• euclidean_functional.hpp — gradient G_v = Θ_v − Σα_v, edge gradient G_e = α_opp⁺ + α_opp⁻ − φ_e, energy via Gauss-Legendre path integral: all correct.
• spherical_functional.hpp — vertex path: correct. Edge-DOF replacement parameterization (Finding 3 in java-port-audit.md): correct.
• euclidean_hessian.hpp — euclidean_cot_weights() formula and sign are correct (only the comment is wrong, see Finding C). Analytic cyclic Hessian ∂α_i/∂s_j formula is correct and internally consistent.
• hyper_ideal_geometry.hpp — lij, sigma_i, sigma_ij, alpha_ij, zeta, zeta13/14/15: all match Java oracle.
• newton_solver.hpp — merit f = ½‖G‖², Armijo conditions Phase 1 and 2, steepest-descent fallback d_sd = −H·G (correct descent direction even for NSD H), spherical sign flip (−H)·Δx = G: all mathematically correct.
• gauss_bonnet.hpp — χ = V−E+F, g = (2−χ)/2, enforce shift δ = (lhs−rhs)/V: correct for Euclidean and Spherical. Wrong for HyperIdeal (Finding B).
• cp_euclidean_functional.hpp — p_function, energy, gradient: match Java BPS-2010 formula. Hessian h_jk = sin θ / (cosh Δρ − cos θ): correct.
• inversive_distance_functional.hpp — edge_length_squared, gradient (Θ − Σα pattern), degenerate-face limiting angles (Finding 9 fix correct).