Phase 9a-Newton: newton_cp_euclidean + newton_inversive_distance

Wires the two Phase-9a functionals into the Newton-solver layer so
they are operational end-to-end.  CGAL test count: 212 → 219 (+7).

Solvers
───────
* newton_cp_euclidean(mesh, x0, m, tol, max_iter)
    - Uses cp_euclidean_hessian — analytic 2×2-per-edge BPS-2010
      formula h_jk = sin θ / (cosh Δρ − cos θ).
    - SparseQR fallback handles the gauge-singular case when no face
      is pinned (caller error, but we recover gracefully).
    - Strictly-convex energy ⇒ quadratic convergence near optimum.

* newton_inversive_distance(mesh, x0, m, tol, max_iter, hess_eps)
    - Uses an inline FD Hessian (n × gradient evaluations per step) —
      mirrors the Phase 4a HyperIdeal solver in spirit.
    - Analytic alternative via Glickenstein 2011 eq. (4.6) is tracked
      in doc/roadmap/research-track.md as Phase 9a.2-analytic.
    - Sensitive to initial point; the test suite always starts from
      a natural-theta setup (u = 0 is the equilibrium when
      compute_inversive_distance_init_from_mesh was called).

Tests (test_newton_phase9a.cpp, 7 cases)
────────────────────────────────────────
* CPEuclidean_NaturalPhi_ClosedTetrahedron_ConvergesInZeroIterations
* CPEuclidean_PerturbedStart_ConvergesBackToEquilibrium
* CPEuclidean_OpenTetrahedron_NaturalPhi_Converges
* InversiveDistance_NaturalTheta_Triangle_ConvergesInZero
* InversiveDistance_PerturbedQuadStrip_Converges
* InversiveDistance_PerturbedTetrahedron_Converges
* CPEuclidean_UsesAnalyticHessian
    Regression guard: 3-DOF problem converges in ≤ 10 iterations even
    with strong perturbation, confirming the analytic Hessian path is
    actually used.

All seven tests pass.  Full CGAL suite: 219/219 PASSED, 0 SKIPPED.

Roadmap additions (`doc/roadmap/phases.md`)
───────────────────────────────────────────
New Phase 11+ section flags two Java sub-packages as optional/deferred
ports, recorded for project memory but not roadmap commitments:

* 11a — Schottky uniformisation (Java plugin/schottky/*, ~3000 LoC)
        Hyperbolic loxodromic group acting on S²; complement of the
        Phase 10c Fuchsian-group representation in H².  Requires
        Phase 10b period matrix + Möbius-group machinery from Phase 7.
        Effort: very large (4-6 weeks).

* 11b — Riemann maps (Java plugin/riemannmap/*, ~1500 LoC)
        Discrete Riemann mapping theorem; texture mapping of bounded
        planar regions, classical conformal mapping for engineering.
        Requires Phase 10b' quasi-isothermic or Phase 9a.1 CP-Euclidean.
        Effort: large (3-4 weeks).

Both are explicitly NOT roadmap commitments — they live in the doc so
they aren't re-discovered.

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

code/include/newton_solver.hpp

@@ -31,6 +31,8 @@
 #include "euclidean_hessian.hpp"
 #include "spherical_hessian.hpp"
 #include "hyper_ideal_hessian.hpp"
+#include "cp_euclidean_functional.hpp"
+#include "inversive_distance_functional.hpp"
 #include <Eigen/SparseCholesky>
 #include <Eigen/SparseQR>
 #include <Eigen/OrderingMethods>
@@ -375,4 +377,187 @@ inline NewtonResult newton_hyper_ideal(
     return res;
 }
 
+// ── CP-Euclidean Newton solver (Phase 9a.1) ───────────────────────────────────
+
+/// Solve the CP-Euclidean circle-packing problem: find ρ ∈ ℝ^F such that the
+/// per-face angle sums match φ_f at every free face.
+///
+/// The CP-Euclidean energy (Bobenko-Pinkall-Springborn 2010 §6) is strictly
+/// convex on its open domain of validity, so the Hessian H is PSD and the
+/// solution is unique up to the gauge mode pinned by `f_idx == −1`.
+/// `cp_euclidean_hessian` provides the analytic 2×2-per-edge formula
+/// `h_jk = sin θ / (cosh Δρ − cos θ)`; no FD machinery is required.
+///
+/// \param mesh      Input triangle mesh (closed or with boundary).
+/// \param x0        Initial DOF vector (length = number of free faces).
+///                  All-zeros is a valid start.
+/// \param m         CPEuclideanMaps: f_idx must have one pinned face
+///                  (`f_idx[f0] == −1`); theta_e and phi_f set by the caller.
+/// \param tol       Convergence threshold on `‖G‖∞`.  Default: 1e-8.
+/// \param max_iter  Newton iteration limit.  Default: 200.
+/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
+///
+/// \note  Unlike the Euclidean solver, the CP-Euclidean Hessian is exact
+///        (analytic), so the SparseQR fallback only triggers in genuine
+///        gauge-singular situations (no pinned face).
+/// \see   doc/architecture/phase-9a-validation.md §1 for the BPS-2010 mapping.
+inline NewtonResult newton_cp_euclidean(
+    ConformalMesh&             mesh,
+    std::vector<double>        x0,
+    const CPEuclideanMaps&     m,
+    double                     tol      = 1e-8,
+    int                        max_iter = 200)
+{
+    std::vector<double> x = x0;
+    const int n = static_cast<int>(x.size());
+
+    NewtonResult res;
+    res.converged      = false;
+    res.iterations     = 0;
+    res.grad_inf_norm  = 0.0;
+
+    for (int iter = 0; iter < max_iter; ++iter) {
+        auto G_std = cp_euclidean_gradient(mesh, x, m);
+        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
+
+        double inf_norm = G.cwiseAbs().maxCoeff();
+        if (inf_norm < tol) {
+            res.converged     = true;
+            res.grad_inf_norm = inf_norm;
+            res.iterations    = iter;
+            res.x             = x;
+            return res;
+        }
+
+        auto H = cp_euclidean_hessian(mesh, x, m);
+        bool ok = false;
+        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
+        if (!ok) break;
+
+        double norm0 = G.norm();
+        x = detail::line_search(x, dx, norm0,
+            [&](const std::vector<double>& xnew) {
+                return cp_euclidean_gradient(mesh, xnew, m);
+            });
+
+        res.iterations = iter + 1;
+    }
+
+    auto G_final = cp_euclidean_gradient(mesh, x, m);
+    double inf_final = 0.0;
+    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
+    res.grad_inf_norm = inf_final;
+    res.x             = x;
+    return res;
+}
+
+// ── Inversive-Distance Newton solver (Phase 9a.2) ─────────────────────────────
+
+/// Solve the inversive-distance circle-packing problem: find u ∈ ℝ^V such that
+/// Σ_{faces adj v} α_v(u) = Θ_v at every free vertex (Luo 2004 Lemma 3.1).
+///
+/// The inversive-distance energy is (locally) strictly convex on the open
+/// domain where every triangle satisfies the inequalities.  Luo's 1-form is
+/// closed there, so the path-integral energy is well-defined.
+///
+/// MVP implementation: the Hessian is computed by **finite differences** of
+/// the analytic gradient (same pattern as the Phase 4a HyperIdeal solver).
+/// An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in
+/// `doc/roadmap/research-track.md` as Phase 9a.2-analytic.
+///
+/// \param mesh      Input triangle mesh.
+/// \param x0        Initial DOF vector (length = number of free vertices).
+/// \param m         InversiveDistanceMaps: v_idx has at least one pinned
+///                  vertex; I_e and r0 set by compute_inversive_distance_init.
+/// \param tol       Convergence threshold on `‖G‖∞`.  Default: 1e-8.
+/// \param max_iter  Newton iteration limit.  Default: 200.
+/// \param hess_eps  FD step size for the Hessian.  Default: 1e-5.
+/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
+///
+/// \note  Convergence is sensitive to the initial point: u = 0 is the
+///        natural choice when `compute_inversive_distance_init_from_mesh`
+///        has been called, since the Bowers-Stephenson identity reconstructs
+///        the input edge lengths at u = 0.
+inline NewtonResult newton_inversive_distance(
+    ConformalMesh&                mesh,
+    std::vector<double>           x0,
+    const InversiveDistanceMaps&  m,
+    double                        tol      = 1e-8,
+    int                           max_iter = 200,
+    double                        hess_eps = 1e-5)
+{
+    std::vector<double> x = x0;
+    const int n = static_cast<int>(x.size());
+
+    NewtonResult res;
+    res.converged      = false;
+    res.iterations     = 0;
+    res.grad_inf_norm  = 0.0;
+
+    // Local FD Hessian builder — n × (cost of gradient eval).
+    auto build_hessian = [&](const std::vector<double>& xc) -> Eigen::SparseMatrix<double> {
+        std::vector<Eigen::Triplet<double>> trips;
+        trips.reserve(static_cast<std::size_t>(n) * 16);   // sparse heuristic
+
+        std::vector<double> xp = xc, xm = xc;
+        for (int j = 0; j < n; ++j) {
+            const std::size_t sj = static_cast<std::size_t>(j);
+            xp[sj] = xc[sj] + hess_eps;
+            xm[sj] = xc[sj] - hess_eps;
+
+            auto Gp = inversive_distance_gradient(mesh, xp, m);
+            auto Gm = inversive_distance_gradient(mesh, xm, m);
+
+            xp[sj] = xm[sj] = xc[sj];   // restore
+
+            for (int i = 0; i < n; ++i) {
+                double val = (Gp[static_cast<std::size_t>(i)]
+                            - Gm[static_cast<std::size_t>(i)])
+                             / (2.0 * hess_eps);
+                if (std::abs(val) > 1e-15)
+                    trips.emplace_back(i, j, val);
+            }
+        }
+        Eigen::SparseMatrix<double> H(n, n);
+        H.setFromTriplets(trips.begin(), trips.end());
+        // Symmetrise — FD rounding may introduce tiny asymmetries.
+        Eigen::SparseMatrix<double> Ht = H.transpose();
+        return (H + Ht) * 0.5;
+    };
+
+    for (int iter = 0; iter < max_iter; ++iter) {
+        auto G_std = inversive_distance_gradient(mesh, x, m);
+        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
+
+        double inf_norm = G.cwiseAbs().maxCoeff();
+        if (inf_norm < tol) {
+            res.converged     = true;
+            res.grad_inf_norm = inf_norm;
+            res.iterations    = iter;
+            res.x             = x;
+            return res;
+        }
+
+        auto H = build_hessian(x);
+        bool ok = false;
+        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
+        if (!ok) break;
+
+        double norm0 = G.norm();
+        x = detail::line_search(x, dx, norm0,
+            [&](const std::vector<double>& xnew) {
+                return inversive_distance_gradient(mesh, xnew, m);
+            });
+
+        res.iterations = iter + 1;
+    }
+
+    auto G_final = inversive_distance_gradient(mesh, x, m);
+    double inf_final = 0.0;
+    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
+    res.grad_inf_norm = inf_final;
+    res.x             = x;
+    return res;
+}
+
 } // namespace conformallab