diff --git a/doc/reviewer/external-audit-2026-05-30.md b/doc/reviewer/external-audit-2026-05-30.md new file mode 100644 index 0000000..c5599d0 --- /dev/null +++ b/doc/reviewer/external-audit-2026-05-30.md @@ -0,0 +1,694 @@ +# 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 + +```bash +# 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. + +```cpp +// 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: +```cpp +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: + +```cpp +// 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. + +### Acceptance criteria +- [ ] All 8 cases (2Β³ combinations of v1b/v2b/v3b) are handled with explicit branches +- [ ] A test exercises a mixed configuration (e.g., one vertex pinned) and checks the + gradient via `gradient_check()` β€” the result should match the FD check +- [ ] 246 CGAL tests still pass + +--- + +## 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: + +```cpp +// 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 +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. + +```cpp +// 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. + +### Acceptance criteria +- [ ] Calling `check_gauss_bonnet(mesh, hyper_ideal_maps)` on a valid genus-2 mesh + with Θ_v = 2Ο€ either (a) does not compile, (b) throws with a meaningful message + distinguishing "wrong formula" from "actual violation", or (c) is removed +- [ ] A comment explains why Euclidean G-B does not apply to HyperIdeal +- [ ] 246 CGAL tests still pass + +--- + +## 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`: + +```cpp +// 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 + +```cpp +// 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: + +```cpp +// 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: + +```cpp +// 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. + +```cpp +/// 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: + +```cpp +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: + +```cpp +// 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: + +```cpp +// 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: + +```cpp +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 + +```cpp +// 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 + +```cpp +// 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 + +```cpp +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 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. + +```cpp +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. + ```yaml + 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`] + +```cpp +// 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`] + +```cpp +// 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`: +```cpp +// 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: +```cpp +// code/include/conformal_mesh.hpp (add): +inline std::size_t halfedge_idx(Halfedge_index h) noexcept { + return static_cast(static_cast(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: + +```cpp +// 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) | πŸ”΄ Open | +| B | `gauss_bonnet.hpp` | 87–88, 128–134 | API error | Medium | 🟑 Open | +| C | `euclidean_hessian.hpp` | 26–27, 57–58 | Doc error | Medium | 🟑 Open | +| 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).