conformallab/ConformalLabpp
3
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

doc(reviewer): add external audit 2026-05-30

Browse Source 
Full external code review of v0.10.0 against all math-critical headers
(direct source read, not automated scan). Documents 5 open bugs/API errors,
3 test gaps inherited from java-port-audit.md, 1 architectural CI risk,
and 5 minor findings. Each finding is self-contained with file+line
references, a concrete fix proposal, and acceptance criteria so a new
session can pick up any single finding without prior context.

Key findings:
- FINDING-A (critical): face_energy() wrong for mixed ideal/hyper-ideal configs
- FINDING-B (medium):   Gauss–Bonnet API conceptually wrong for HyperIdeal
- FINDING-C (medium):   cotangent formula in euclidean_hessian.hpp header is wrong
- FINDING-D (medium):   DOF-assignment doc falsely claims "pin before" works
- FINDING-E (medium):   cp_euclidean gradient_check uses absolute not relative error
- FINDING-F/G/H:        three open test gaps from java-port-audit.md
- FINDING-I (arch):     246/272 CGAL tests not gated in CI

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
2026-05-30 23:44:58 +02:00
committed by Tarik Moussa
parent 6b3497f451
commit 7f1a077269
1 changed files with 694 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

694
doc/reviewer/external-audit-2026-05-30.md Normal file
View File

@@ -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 <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.
```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<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.
```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<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:
```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).