fix(minor): all five MINOR findings — doc, accuracy, and DRY

MINOR-1 (spherical_functional.hpp:426)
  Wrong comment said "second derivative < 0 for a convex functional".
  The spherical energy is *concave* (NSD Hessian); the monotone-f argument
  applies to both convex and concave functionals equally.  Comment rewritten
  to explain the actual physics: increasing scale increases all angles and
  thus reduces Σ G_v.

MINOR-2 (spherical_functional.hpp:494-495)
  Forward finite difference O(ε) → central finite difference O(ε²):
    old: dft = (sum_Gv(t + fd_eps) - ft) / fd_eps
    new: dft = (sum_Gv(t + fd_eps) - sum_Gv(t - fd_eps)) / (2*fd_eps)
  Same cost when the extra sum_Gv(t - fd_eps) replaces the cached ft.

MINOR-3 (euclidean_functional.hpp, spherical_functional.hpp,
         inversive_distance_functional.hpp)
  New header gauss_legendre.hpp centralises the 10-point Gauss-Legendre
  nodes and weights (gl10_nodes() / gl10_weights()).  The three energy
  functions now use the shared accessors instead of duplicated local
  static arrays.

MINOR-4 (euclidean_functional.hpp, spherical_functional.hpp,
         hyper_ideal_functional.hpp, inversive_distance_functional.hpp)
  halfedge_to_index() centralised in conformal_mesh.hpp.  All four local
  aliases (eucl_hidx, spher_hidx, hidx, id_detail::hidx) now delegate to
  it as one-line wrappers; the aliases are kept for now to avoid a larger
  call-site churn, clearly documented as thin wrappers.

MINOR-5 (clausen.hpp:33-38)
  Added a comment above inits() explaining the intentional off-by-one
  return value and how it interacts with csevl() — matching the Java
  Clausen.inits() / csevl() contract.

277/277 CGAL + 26/26 pure-math tests pass, 0 failed.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

code/include/spherical_functional.hpp

@@ -36,6 +36,7 @@
 
 #include "conformal_mesh.hpp"
 #include "spherical_geometry.hpp"
+#include "gauss_legendre.hpp"
 #include <CGAL/boost/graph/iterator.h>
 #include <vector>
 #include <cmath>
@@ -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<std::size_t>(static_cast<std::uint32_t>(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<double>& 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;