ConformalLabpp/doc/math/validation.md

Mathematical Validation

This document lists analytically known results and explains how to verify them against conformallab++ output. It is the primary tool for an independent mathematician to check the correctness of the implementation.

---

How to run the examples

cmake -S code -B build -DWITH_CGAL=ON -DCMAKE_BUILD_TYPE=Release
cmake --build build --target conformallab_cgal_tests
ctest --test-dir build -R cgal --output-on-failure

All tests pass, 0 skipped (per-suite breakdown: doc/api/tests.md).

---

1 — Gauss–Bonnet (topology)

Theorem. For any closed triangulated surface M,

Σᵥ (2π − Θᵥ) = 2π · χ(M)

where χ(M) = 2 − 2g is the Euler characteristic.

Surface	g	χ	Σ(2π − Θᵥ)
Sphere (tetrahedron, cube, …)	0	2	4π
Torus	1	0	0
Double torus	2	−2	−4π

How to check:

#include "gauss_bonnet.hpp"
auto defect = gauss_bonnet_sum(mesh, maps);   // Σ(2π − Θᵥ)
auto chi    = euler_characteristic(mesh);     // free function, not a member
EXPECT_NEAR(defect, 2.0 * M_PI * chi, 1e-10);

Covered by: cgal.GaussBonnet.* tests in test_phase6.cpp.

---

2 — Period matrix: fundamental domain invariants

Theorem (SL(2,ℤ)-reduction). Every lattice τ ∈ ℍ has a unique representative in the standard fundamental domain

F = { τ ∈ ℍ :  |τ| ≥ 1,  |Re(τ)| ≤ 1/2,  Im(τ) > 0 }

After calling compute_period_matrix(hol), the returned τ must satisfy:

Condition	Invariant
pd.tau_reduced.imag() > 0	τ lies in the upper half-plane
std::abs(pd.tau_reduced) >= 1.0 - 1e-10	τ outside unit disk
std::abs(pd.tau_reduced.real()) <= 0.5 + 1e-10	τ in vertical strip

These three conditions hold for any closed genus-1 triangulated surface processed through Euclidean uniformization — they are topology, not geometry.

Covered by: cgal.PeriodMatrix.TauInFundamentalDomain_* tests in test_phase7.cpp.

---

3 — Torus of revolution: conformal modulus

Setup. The bundled torus meshes are surfaces of revolution (a tube of minor radius r swept around a major circle of radius R > r), not abstract square/hexagonal flat tori:

mesh	major R	minor r	cross section
torus_4x4.off	2	1	square (4-gon)
torus_hex_6x6.off	3	1	hexagon (6-gon)
torus_8x8.off	3	1	octagon (8-gon)

Expected. The induced metric ds² = (R + r cos φ)² dθ² + r² dφ² is made flat by the conformal change of variable dψ = r/(R + r cos φ) dφ. The ψ-period is ∮ r/(R + r cos φ) dφ = 2πr/√(R²−r²), so the flat torus is the rectangular lattice 2π·ℤ × (2πr/√(R²−r²))·ℤ. Its period ratio, reduced so that |τ| ≥ 1, is therefore purely imaginary:

Re(τ) = 0                  (meridian ⟂ longitude reflection symmetry)
Im(τ) = √(R² − r²) / r     (reduced conformal modulus)

giving the analytic targets

mesh	analytic τ = i·√(R²−r²)/r
torus_4x4.off	i·√3 ≈ 1.732 i
torus_hex_6x6.off	i·√8 ≈ 2.828 i
torus_8x8.off	i·√8 ≈ 2.828 i

A common pitfall is to expect τ = i (or the order-6 fixed point e^{iπ/3}) from the 4-fold (6-fold) symmetry. That reasoning is wrong here: a torus of revolution's rotational symmetry is a rotation about the axis, which acts on the surface as a fixed-point-free translation along the longitude — it is not an order-4 conformal automorphism with a fixed point, so it does not pin τ to a modular fixed point. The correct invariant is the rectangular modulus above.

Reproduced end-to-end (solve → compute_cut_graph → euclidean_layout(…, &cg, &hol) → compute_period_matrix). The coarse polygonal cross sections approximate the circular modulus from above; the gap shrinks as the cross section gains sides:

mesh	analytic	computed τ	rel. error
torus_4x4.off	1.732 i	≈ 1.79 i	~3 % (4-gon)
torus_hex_6x6.off	2.828 i	≈ 2.85 i	~0.8 % (6-gon)
torus_8x8.off	2.828 i	≈ 2.84 i	~0.4 % (8-gon)

Covered by: cgal.HolonomyEndToEnd.Torus*_TauMatchesRevolutionModulus in test_phase7.cpp.

Conformal flattening of the torus is wired end-to-end and converges in a handful of Newton steps (3–4 on the bundled meshes):

./bin/conformallab_core -i code/data/off/torus_4x4.off -g euclidean -v -o lay.off
# → topology: closed, free DOFs=15, genus=1
# → Euclidean: converged=yes  iter=3  |grad|_inf≈1e-12

Equivalent C++ (current API — note newton_euclidean takes an x0 vector and the DOF indices must be assigned first):

ConformalMesh mesh = load_mesh("code/data/off/torus_4x4.off");
EuclideanMaps  maps = setup_euclidean_maps(mesh);   // Θ_v = 2π (flat target)
compute_euclidean_lambda0_from_mesh(mesh, maps);

// Pin one vertex (scale gauge), free the rest; make the flat target
// Gauss-Bonnet-consistent.
int idx = 0; bool pinned = false;
for (auto v : mesh.vertices())
    maps.v_idx[v] = (!pinned ? (pinned = true, -1) : idx++);
enforce_gauss_bonnet(mesh, maps);

std::vector<double> x0(idx, 0.0);
auto res = newton_euclidean(mesh, x0, maps);        // converged after ~3 iters

The CLI reports the period ratio for genus-1 inputs, e.g.

./bin/conformallab_core -i code/data/off/torus_4x4.off -g euclidean -v
# → period ratio τ = 0.000000 + 1.793... i  (genus 1, reduced to fundamental domain)

---

4 — Higher-resolution cross sections converge to the circular modulus

torus_hex_6x6.off (R=3, r=1, hexagon) and torus_8x8.off (R=3, r=1, octagon) share the same analytic modulus i·√8 ≈ 2.828 i (see §3). Because both approximate the circular tube cross section, the computed τ approaches the analytic value as the polygon gains sides: the 8-gon (≈ 2.84 i, ~0.4 %) is closer than the 6-gon (≈ 2.85 i, ~0.8 %), which is closer than the 4-gon of torus_4x4.off (~3 %). All three are asserted in cgal.HolonomyEndToEnd.Torus*_TauMatchesRevolutionModulus.

---

5 — Newton convergence rate

Theorem. Because the Euclidean and hyper-ideal energies are strictly convex (after gauge-fixing), Newton's method converges quadratically near the optimum.

Expected: for any mesh with up to a few hundred faces, Newton converges in fewer than 30 iterations starting from u = 0.

std::vector<double> x0(n_dofs, 0.0);
auto res = newton_euclidean(mesh, x0, maps);
EXPECT_LT(res.iterations, 30);
EXPECT_LT(res.grad_inf_norm, 1e-8);   // field is grad_inf_norm, default tol 1e-8

Covered by: cgal.NewtonSolver.* (Euclidean ×3, Spherical ×4, HyperIdeal ×4) and cgal.NewtonPhase9a.* (the two circle-packing solvers).

---

6 — Gradient check (finite differences)

For each functional F(u), the gradient G = ∂F/∂u is verified by:

|G(u)ᵢ  −  (F(u + εeᵢ) − F(u − εeᵢ)) / (2ε)| < 1e-6

with ε = 1e-5. This check is run inside the test suite for all three geometries (Euclidean, Spherical, HyperIdeal) at u = 0 and at random u.

Relevant test suites:

cgal.EuclideanFunctional.GradientCheck_*
cgal.SphericalFunctional.GradientCheck_*
cgal.HyperIdealFunctional.GradientCheck_*

A failing gradient check means the energy and its derivative are inconsistent — the Newton solver will converge to the wrong point.

---

7 — Holonomy composition (Möbius maps)

For a closed surface, the composition of holonomies around any contractible cycle must be the identity. In genus 1 with a single handle:

T₁ · T₂ · T₁⁻¹ · T₂⁻¹ = Id    (commutator = Id for a torus)

because π₁(T²) = ℤ × ℤ is abelian.

For genus g ≥ 2, the fundamental group is non-abelian and this check does not hold, but the representation ρ: π₁(Σ_g) → SU(1,1) must still satisfy the relation

[T₁, T₂] · [T₃, T₄] · … = Id    (product of g commutators = Id)

The Möbius arithmetic these checks rely on (identity, inverse, composition, from_three) is verified in test_phase7.cpp under cgal.MobiusMap.*. The Euclidean end-to-end holonomy extraction is now validated against the analytic torus-of-revolution modulus (§3, cgal.HolonomyEndToEnd.*). A standalone cgal.HolonomyData.* commutator-closes-up suite for the hyperbolic (Möbius) holonomy is still future work.

---

9 — Cross-validation with geometry-central (optional / hypothetical)

Note: This section describes a possible external cross-validation that is not a prerequisite for the correctness of the implementation. It is of interest because geometry-central implements the same mathematical core (Gillespie, Springborn, Crane — SIGGRAPH 2021, building on Springborn 2020), but with a different algorithmic strategy (Ptolemaic flips + intrinsic triangulations instead of Newton on the original triangulation).

Which outputs are comparable?

Output	conformallab++	geometry-central	Comparable?
u-vector (scale parameters)	res.x	u after Yamabe flow	✓ after normalisation
UV coordinates	layout.uv[v]	conformal parameterisation	✓ up to Möbius transformation
Gauss-Bonnet deficit	gauss_bonnet_sum()	implicit via curvature flow	✓ (analytically identical)
Number of Newton iterations	res.iterations	Yamabe steps	~ (different algorithm)
Period matrix τ	pd.tau_reduced	not available	✗
Möbius holonomy	hol.T_a, T_b	not available	✗

Normalisation alignment

The u-vector in conformallab++ has one degree of freedom (global additive constant — gauge freedom after pin-fixing). geometry-central may use a different convention. Normalise before comparing:

// conformallab++: centre u
double mean_u = std::accumulate(x.begin(), x.end(), 0.0) / x.size();
std::vector<double> x_norm(x.size());
for (int i = 0; i < x.size(); ++i) x_norm[i] = x[i] - mean_u;

// Then compare with the geometry-central u-vector (also centred):
// max|x_norm[i] - gc_u[i]| < 1e-8  →  identical convergence point

When is the comparison useful?

Point in time	What is possible
Now (Phase 7)	Manual comparison using the same .off/.obj test meshes
After Phase 8	Automated comparison script (Python or separate C++ binary)
Phase 10 (research)	Algorithm comparison: Newton vs. Ptolemaic flips on difficult meshes

Connection to the literature

The Springborn 2020 paper ("Ideal Hyperbolic Polyhedra and Discrete Uniformization") is already implemented in conformallab++ — it is the mathematical foundation for the HyperIdeal geometry mode (Phase 2/3). The geometry-central implementation is based on the extension by Gillespie, Springborn & Crane (2021), which uses the same variational principle of Bobenko–Springborn 2004 but additionally applies Ptolemaic flips to improve the triangulation during optimisation — an idea not yet implemented in conformallab++ (→ GC-2 in the phase roadmap).

---

8 — Checklist for an independent reviewer

Run these in order to validate the implementation:

• ctest --test-dir build -R cgal --output-on-failure → all pass, 0 skipped (count: doc/api/tests.md)
• cgal.GaussBonnet.* all pass → topology is correctly read from mesh
• cgal.EuclideanFunctional.GradientCheck_* pass → energy = integral of gradient
• cgal.PeriodMatrix.* pass → SL(2,ℤ) reduction correct (on prescribed holonomy)
• cgal.MobiusMap.* pass → Möbius arithmetic (identity, inverse, compose) correct

All of the above are deterministic, analytic tests — no mesh loading, no file I/O, no floating-point non-determinism beyond standard IEEE-754.