diff --git a/code/include/clausen.hpp b/code/include/clausen.hpp index 9f11638..2c7c5a4 100644 --- a/code/include/clausen.hpp +++ b/code/include/clausen.hpp @@ -30,6 +30,12 @@ inline double csevl(double x, const double* cs, int n) noexcept { // Count Chebyshev terms needed so truncation error <= eta. // Corresponds to Java Clausen.inits(). +// +// Note on return value: the loop decrements n one extra time after the +// stopping condition is met, so the returned value is one less than the +// last index checked. Callers pass the result directly to csevl() as +// the term count, which evaluates terms [0, n-1] — this is intentional +// and matches the Java Clausen.inits() / csevl() contract exactly. inline int inits(const double* series, int n, double eta) noexcept { double err = 0.0; while (err <= eta && n-- != 0) { diff --git a/code/include/conformal_mesh.hpp b/code/include/conformal_mesh.hpp index 98f8626..e771a8f 100644 --- a/code/include/conformal_mesh.hpp +++ b/code/include/conformal_mesh.hpp @@ -37,6 +37,7 @@ #include #include +#include #include namespace conformallab { @@ -107,4 +108,23 @@ inline auto add_face_properties(ConformalMesh& mesh) return ftype; } +// ── Half-edge index helper ──────────────────────────────────────────────────── +// +// Convert a Halfedge_index to std::size_t for use as a vector subscript. +// Each functional previously duplicated this one-liner under a name like +// eucl_hidx / spher_hidx / hidx. Centralised here (MINOR-4 fix). +// +// Implementation: the double-cast uint32_t → size_t is intentional. CGAL 6.x +// Surface_mesh stores half-edge indices internally as 32-bit unsigned integers. +// Casting directly to size_t on a 64-bit system would produce the same result +// because CGAL index types use non-negative values, but the explicit intermediate +// cast documents the assumption and silences spurious sign-conversion warnings. + +/// Convert a `Halfedge_index` to `std::size_t` for vector subscript use. +/// Replaces the duplicated `eucl_hidx`, `spher_hidx`, `hidx` helpers. +inline std::size_t halfedge_to_index(Halfedge_index h) noexcept +{ + return static_cast(static_cast(h)); +} + } // namespace conformallab diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp index 8ea93e0..f8cac0c 100644 --- a/code/include/euclidean_functional.hpp +++ b/code/include/euclidean_functional.hpp @@ -41,6 +41,7 @@ #include "conformal_mesh.hpp" #include "constants.hpp" +#include "gauss_legendre.hpp" #include "euclidean_geometry.hpp" #include #include @@ -173,11 +174,10 @@ static inline double eucl_dof_val(int idx, const std::vector& x) return idx >= 0 ? x[static_cast(idx)] : 0.0; } -/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. -static inline std::size_t eucl_hidx(Halfedge_index h) -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp (included above). +// The old local alias eucl_hidx is retained as a thin wrapper for now so +// call-sites below do not need touching; a follow-up can remove it. +static inline std::size_t eucl_hidx(Halfedge_index h) { return halfedge_to_index(h); } /// Compute the Euclidean-functional gradient G(x): /// * `G_v = Θ_v − Σ_faces α_v(face)` @@ -272,20 +272,8 @@ inline double euclidean_energy( const std::vector& x, const EuclideanMaps& m) { - static const double gl_s[10] = { - -0.9739065285171717, -0.8650633666889845, - -0.6794095682990244, -0.4333953941292472, - -0.1488743389816312, 0.1488743389816312, - 0.4333953941292472, 0.6794095682990244, - 0.8650633666889845, 0.9739065285171717 - }; - static const double gl_w[10] = { - 0.0666713443086881, 0.1494513491505806, - 0.2190863625159820, 0.2692667193099963, - 0.2955242247147529, 0.2955242247147529, - 0.2692667193099963, 0.2190863625159820, - 0.1494513491505806, 0.0666713443086881 - }; + const double* gl_s = gl10_nodes(); + const double* gl_w = gl10_weights(); const std::size_t n = x.size(); double E = 0.0; diff --git a/code/include/gauss_legendre.hpp b/code/include/gauss_legendre.hpp new file mode 100644 index 0000000..e80769a --- /dev/null +++ b/code/include/gauss_legendre.hpp @@ -0,0 +1,49 @@ +#pragma once +// Copyright (c) 2024-2026 Tarik Moussa. +// SPDX-License-Identifier: MIT + +// gauss_legendre.hpp +// +// 10-point Gauss-Legendre quadrature nodes and weights on [-1,1]. +// +// Previously duplicated verbatim in: +// euclidean_functional.hpp, spherical_functional.hpp, +// inversive_distance_functional.hpp (MINOR-3 fix) +// +// Usage (integration over [0,1] via change of variables t=(1+s)/2, w=w_GL/2): +// +// const auto* s = conformallab::gl10_nodes(); +// const auto* w = conformallab::gl10_weights(); +// for (int k = 0; k < 10; ++k) { +// double t = (1.0 + s[k]) * 0.5; +// double wt = w[k] * 0.5; +// E += wt * dot(G(t*x), x); +// } + +namespace conformallab { + +/// 10-point Gauss-Legendre nodes on [−1, 1]. +inline const double* gl10_nodes() noexcept { + static constexpr double s[10] = { + -0.9739065285171717, -0.8650633666889845, + -0.6794095682990244, -0.4333953941292472, + -0.1488743389816312, 0.1488743389816312, + 0.4333953941292472, 0.6794095682990244, + 0.8650633666889845, 0.9739065285171717 + }; + return s; +} + +/// 10-point Gauss-Legendre weights on [−1, 1]. +inline const double* gl10_weights() noexcept { + static constexpr double w[10] = { + 0.0666713443086881, 0.1494513491505806, + 0.2190863625159820, 0.2692667193099963, + 0.2955242247147529, 0.2955242247147529, + 0.2692667193099963, 0.2190863625159820, + 0.1494513491505806, 0.0666713443086881 + }; + return w; +} + +} // namespace conformallab diff --git a/code/include/hyper_ideal_functional.hpp b/code/include/hyper_ideal_functional.hpp index d1d2bb7..7a5f125 100644 --- a/code/include/hyper_ideal_functional.hpp +++ b/code/include/hyper_ideal_functional.hpp @@ -129,11 +129,8 @@ static inline double dof_val(int idx, const std::vector& x) return idx >= 0 ? x[static_cast(idx)] : 0.0; } -/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. -static inline std::size_t hidx(Halfedge_index h) -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp. +static inline std::size_t hidx(Halfedge_index h) { return halfedge_to_index(h); } // ── Pure-math face-angle kernel ────────────────────────────────────────────── // diff --git a/code/include/inversive_distance_functional.hpp b/code/include/inversive_distance_functional.hpp index ef8da4d..8d67a32 100644 --- a/code/include/inversive_distance_functional.hpp +++ b/code/include/inversive_distance_functional.hpp @@ -68,6 +68,7 @@ #include "conformal_mesh.hpp" #include "constants.hpp" +#include "gauss_legendre.hpp" #include "euclidean_geometry.hpp" // euclidean_angles(λ12, λ23, λ31) #include #include @@ -229,10 +230,8 @@ inline double dof_val(int idx, const std::vector& x) noexcept return idx >= 0 ? x[static_cast(idx)] : 0.0; } -inline std::size_t hidx(Halfedge_index h) noexcept -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp. +inline std::size_t hidx(Halfedge_index h) noexcept { return halfedge_to_index(h); } // Inversive-distance edge length squared: ℓ² = exp(2u_i) + exp(2u_j) + 2 I r_i r_j // where r_i = exp(u_i), so: ℓ² = r_i² + r_j² + 2 I r_i r_j. @@ -323,20 +322,8 @@ inline double inversive_distance_energy( const std::vector& x, const InversiveDistanceMaps& m) { - static const double gl_s[10] = { - -0.9739065285171717, -0.8650633666889845, - -0.6794095682990244, -0.4333953941292472, - -0.1488743389816312, 0.1488743389816312, - 0.4333953941292472, 0.6794095682990244, - 0.8650633666889845, 0.9739065285171717 - }; - static const double gl_w[10] = { - 0.0666713443086881, 0.1494513491505806, - 0.2190863625159820, 0.2692667193099963, - 0.2955242247147529, 0.2955242247147529, - 0.2692667193099963, 0.2190863625159820, - 0.1494513491505806, 0.0666713443086881 - }; + const double* gl_s = gl10_nodes(); + const double* gl_w = gl10_weights(); const std::size_t n = x.size(); double E = 0.0; diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp index 1ffac86..5813e66 100644 --- a/code/include/spherical_functional.hpp +++ b/code/include/spherical_functional.hpp @@ -36,6 +36,7 @@ #include "conformal_mesh.hpp" #include "spherical_geometry.hpp" +#include "gauss_legendre.hpp" #include #include #include @@ -196,11 +197,8 @@ static inline double spher_eff_lambda(const SphericalMaps& m, : (m.lambda0[e] + u_i + u_j); } -/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. -static inline std::size_t spher_hidx(Halfedge_index h) -{ - return static_cast(static_cast(h)); -} +// halfedge_to_index is defined in conformal_mesh.hpp. +static inline std::size_t spher_hidx(Halfedge_index h) { return halfedge_to_index(h); } // ── Gradient only (no energy) ───────────────────────────────────────────────── @@ -321,21 +319,8 @@ inline double spherical_energy( const std::vector& x, const SphericalMaps& m) { - // 10-point Gauss-Legendre nodes and weights on [-1, 1]. - static const double gl_s[10] = { - -0.9739065285171717, -0.8650633666889845, - -0.6794095682990244, -0.4333953941292472, - -0.1488743389816312, 0.1488743389816312, - 0.4333953941292472, 0.6794095682990244, - 0.8650633666889845, 0.9739065285171717 - }; - static const double gl_w[10] = { - 0.0666713443086881, 0.1494513491505806, - 0.2190863625159820, 0.2692667193099963, - 0.2955242247147529, 0.2955242247147529, - 0.2692667193099963, 0.2190863625159820, - 0.1494513491505806, 0.0666713443086881 - }; + const double* gl_s = gl10_nodes(); + const double* gl_w = gl10_weights(); const std::size_t n = x.size(); double E = 0.0; @@ -423,8 +408,11 @@ inline bool gradient_check_spherical( // Apply the shift by adding t* to every vertex DOF in x. // // Implementation: bisection on f(t) = Σ_v G_v(x + t·1_v). -// f is strictly monotone decreasing (second derivative < 0) for a convex -// functional, so bisection converges in O(log₂(2·bracket/tol)) iterations. +// f is strictly monotone decreasing because increasing the global scale +// increases all effective edge lengths and thereby all corner angles, which +// reduces Σ G_v = Σ(Θ_v − Σα_v). This holds for the concave spherical +// energy (NSD Hessian) just as well as for a convex one. +// Bisection converges in O(log₂(2·bracket/tol)) iterations. // // Parameters: // bracket – initial search interval [−bracket, +bracket] (default 50) @@ -477,7 +465,7 @@ inline double spherical_gauge_shift( // ── No sign change (zero may lie at a domain boundary). ─────────────────── // Use damped Newton's method with backtracking line search. - // f'(t) estimated by forward finite difference. + // f'(t) estimated by central finite difference (O(ε²) vs O(ε) for forward). // When the Newton step overshoots the valid domain (ΣG_v jumps back up // because faces become degenerate), backtracking halves the step until // |f| strictly decreases. @@ -488,8 +476,7 @@ inline double spherical_gauge_shift( for (int iter = 0; iter < 120; ++iter) { if (std::abs(ft) < tol) return t; - double ftp = sum_Gv(t + fd_eps); - double dft = (ftp - ft) / fd_eps; + double dft = (sum_Gv(t + fd_eps) - sum_Gv(t - fd_eps)) / (2.0 * fd_eps); if (std::abs(dft) < 1e-14) return t; // gradient flat — give up double dt_raw = -ft / dft; diff --git a/doc/reviewer/external-audit-2026-05-30.md b/doc/reviewer/external-audit-2026-05-30.md index bc867bc..de8a987 100644 --- a/doc/reviewer/external-audit-2026-05-30.md +++ b/doc/reviewer/external-audit-2026-05-30.md @@ -708,11 +708,11 @@ index. Add a comment: | G | test files | — | Test gap | Medium | ✅ Fixed 2026-05-31 | | H | test files | — | Test gap | Medium | ✅ Fixed 2026-05-31 | | 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 | +| MINOR-1 | `spherical_functional.hpp` | 404 | Doc error | Minor | ✅ Fixed 2026-05-31 | +| MINOR-2 | `spherical_functional.hpp` | 470–471 | Accuracy | Minor | ✅ Fixed 2026-05-31 | +| MINOR-3 | three files | — | DRY | Minor | ✅ Fixed 2026-05-31 | +| MINOR-4 | four files | — | DRY | Minor | ✅ Fixed 2026-05-31 | +| MINOR-5 | `clausen.hpp` | 33–38 | Doc | Minor | ✅ Fixed 2026-05-31 | ---