#pragma once
// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// newton_solver.hpp
//
// Phase 4a — Newton solver for all three discrete conformal functionals.
//
// Solves G(x) = 0 where G is the gradient of the discrete conformal energy.
//
// ┌──────────────────────────────────────────────────────────────────────────┐
// │  Gradient sign conventions                                               │
// │    Euclidean / Spherical:  G_v = Θ_v − Σ α_v    (target − actual)      │
// │    HyperIdeal:             G_v = Σ β_v − Θ_v    (actual − target)      │
// │                                                                          │
// │  All solvers use the same Newton step  Δx = −H⁻¹·G                     │
// │                                                                          │
// │  Hessian sign at equilibrium                                             │
// │    Euclidean:   H PSD  → SimplicialLDLT on  H                          │
// │    Spherical:   H NSD  → SimplicialLDLT on −H  (solve (−H)Δx = G)      │
// │    HyperIdeal:  H PSD  → SimplicialLDLT on  H  (analytical H: future)  │
// └──────────────────────────────────────────────────────────────────────────┘
//
// SparseQR fallback:
//   When SimplicialLDLT reports a failure (e.g. singular H on a closed mesh
//   without a pinned vertex), the solver automatically retries with
//   Eigen::SparseQR, which finds the minimum-norm Newton step orthogonal to
//   the null space.  This handles the gauge mode on closed surfaces without
//   requiring the caller to pin a vertex explicitly.
//
// Requires:
//   Eigen::SimplicialLDLT, Eigen::SparseQR (Eigen sparse module)

#include "euclidean_hessian.hpp"
#include "spherical_hessian.hpp"
#include "hyper_ideal_hessian.hpp"
#include "cp_euclidean_functional.hpp"
#include "inversive_distance_functional.hpp"
#include "inversive_distance_hessian.hpp"
#include <Eigen/SparseCholesky>
#include <Eigen/SparseQR>
#include <Eigen/OrderingMethods>
#include <Eigen/Dense>
#include <algorithm>
#include <cmath>

namespace conformallab {

// ── Result ────────────────────────────────────────────────────────────────────

/// Result of `newton_solve(...)` — converged DOF vector + diagnostics.
struct NewtonResult {
    std::vector<double> x;            ///< DOF vector at termination.
    int                 iterations;   ///< Newton steps taken.
    double              grad_inf_norm;///< max |Gᵢ| at termination.
    bool                converged;    ///< `true` iff `grad_inf_norm < tol`.
};

// ── Internal helpers ──────────────────────────────────────────────────────────

namespace detail {

// Solve A·Δx = rhs with SimplicialLDLT; on failure fall back to SparseQR.
// Returns Δx.  ok is set to false only if both solvers fail.
// If fallback_used is non-null, it is set to true iff SparseQR was needed.
inline Eigen::VectorXd solve_with_fallback(
    const Eigen::SparseMatrix<double>& A,
    const Eigen::VectorXd&             rhs,
    bool&                              ok,
    bool*                              fallback_used = nullptr)
{
    if (fallback_used) *fallback_used = false;

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt(A);
    if (ldlt.info() == Eigen::Success) {
        Eigen::VectorXd dx = ldlt.solve(rhs);
        if (ldlt.info() == Eigen::Success) { ok = true; return dx; }
    }
    // Fallback: SparseQR — handles singular/rank-deficient H (gauge modes).
    if (fallback_used) *fallback_used = true;
    Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr(A);
    if (qr.info() == Eigen::Success) {
        Eigen::VectorXd dx = qr.solve(rhs);
        if (qr.info() == Eigen::Success) { ok = true; return dx; }
    }
    ok = false;
    return Eigen::VectorXd::Zero(rhs.size());
}

} // namespace detail

// ── Public linear-system solver (SparseQR fallback) ──────────────────────────
//
// Solve A·x = rhs with Eigen::SimplicialLDLT; if that fails (singular or
// rank-deficient A), retry with Eigen::SparseQR which finds the minimum-norm
// solution orthogonal to the null space.
//
// This is the same primitive used internally by all three Newton solvers.
// Exposing it publicly lets callers (tests, downstream code) reuse the logic
// and — via the optional fallback_used pointer — verify which code path ran.
//
//   fallback_used  – if non-null, set to true iff SparseQR was invoked
//   Returns Eigen::VectorXd::Zero(rhs.size()) if both solvers fail.
/// Solve `A·x = rhs` with the same SimplicialLDLT → SparseQR fallback
/// strategy used inside all three Newton solvers.  If `fallback_used`
/// is non-null, it is set to `true` iff the SparseQR fallback ran.
/// If `ok` is non-null, it is set to `true` iff at least one solver succeeded.
/// Returns `Eigen::VectorXd::Zero(rhs.size())` if both solvers fail.
inline Eigen::VectorXd solve_linear_system(
    const Eigen::SparseMatrix<double>& A,
    const Eigen::VectorXd&             rhs,
    bool*                              fallback_used = nullptr,
    bool*                              ok            = nullptr)
{
    bool ok_local = false;
    auto x = detail::solve_with_fallback(A, rhs, ok_local, fallback_used);
    if (ok) *ok = ok_local;
    return x;
}

namespace detail {   // re-open for the remaining helpers

// Globalised line search for the Newton system G(x) = 0.
//
// Merit function  f(x) = ½‖G(x)‖²₂.  Driving f to its minimum drives the
// residual G to zero; because the merit only depends on ‖G‖ it is sign-agnostic
// and works identically for the convex Euclidean / HyperIdeal / CP / InvDist
// energies and the concave Spherical energy.
//
//   Phase 1 — Newton direction `dx` (satisfies H·dx = −G, so the merit slope
//             ∇f·dx = GᵀH·dx = −‖G‖² < 0 — always a descent direction).
//             Backtrack α ∈ {1, ½, ¼, …} until the Armijo sufficient-decrease
//             condition holds:
//                   ‖G(x + α·dx)‖²  ≤  (1 − 2·c1·α)·‖G(x)‖²
//
//   Phase 2 — if Phase 1 exhausts its halvings, fall back to the steepest-
//             descent direction of the merit, `d_sd = −H·G` (slope
//             ∇f·d_sd = −‖H·G‖² ≤ 0 regardless of the definiteness of H), with
//             the analogous Armijo test:
//                   ‖G(x + α·d_sd)‖²  ≤  ‖G(x)‖² − 2·c1·α·‖d_sd‖²
//
// If neither phase satisfies Armijo, return the best (smallest-residual) point
// visited.  If nothing beat ‖G(x)‖, return x unchanged and set *improved=false,
// so the caller can stop cleanly instead of taking the old divergent full step.
//
// B2 (api-performance audit): the gradient computed at the accepted trial point
// is normally discarded (only its norm is kept), forcing the next Newton
// iteration to recompute G(x) from scratch.  When `accepted_grad` is non-null it
// receives the gradient vector at the returned point, so the caller can reuse it
// as the next iteration's convergence gradient.  It is only written when the
// returned point differs from `x` (i.e. `*improved == true`); on a pure stall
// (`*improved == false`) the caller must recompute G(x) itself.
template <typename GradFn>
inline std::vector<double> line_search(
    const std::vector<double>& x,
    const Eigen::VectorXd&     dx,
    const Eigen::VectorXd&     d_sd,
    double                     norm0,
    GradFn&&                   grad_fn,
    bool*                      improved      = nullptr,
    std::vector<double>*       accepted_grad = nullptr,
    int                        max_halvings  = 20,
    double                     c1            = 1e-4)
{
    const int    n        = static_cast<int>(x.size());
    const double norm0_sq = norm0 * norm0;

    std::vector<double> xnew(static_cast<std::size_t>(n));
    std::vector<double> best_x    = x;
    double              best_norm = norm0;

    std::vector<double> trial_grad;   // gradient at the most recent trial point
    std::vector<double> best_grad;    // gradient at best_x (B2)

    // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew` and the
    // gradient there in `trial_grad`.
    auto eval = [&](const Eigen::VectorXd& dir, double alpha) -> double {
        for (int i = 0; i < n; ++i)
            xnew[static_cast<std::size_t>(i)] =
                x[static_cast<std::size_t>(i)] + alpha * dir[i];
        trial_grad = grad_fn(xnew);
        double s = 0.0;
        for (double v : trial_grad) s += v * v;
        return std::sqrt(s);
    };

    // ── Phase 1: Newton direction, Armijo backtracking ────────────────────────
    double alpha = 1.0;
    for (int ls = 0; ls < max_halvings; ++ls) {
        double norm_new = eval(dx, alpha);
        if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; }
        if (norm_new * norm_new <= (1.0 - 2.0 * c1 * alpha) * norm0_sq) {
            if (improved) *improved = true;
            if (accepted_grad) *accepted_grad = trial_grad;
            return xnew;
        }
        alpha *= 0.5;
    }

    // ── Phase 2: steepest-descent fallback (−H·G), Armijo backtracking ────────
    const double dsd_sq = d_sd.squaredNorm();
    if (dsd_sq > 0.0) {
        alpha = 1.0;
        for (int ls = 0; ls < max_halvings; ++ls) {
            double thresh_sq = norm0_sq - 2.0 * c1 * alpha * dsd_sq;
            double norm_new  = eval(d_sd, alpha);
            if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; }
            if (thresh_sq >= 0.0 && norm_new * norm_new <= thresh_sq) {
                if (improved) *improved = true;
                if (accepted_grad) *accepted_grad = trial_grad;
                return xnew;
            }
            alpha *= 0.5;
        }
    }

    // ── Both phases failed Armijo — never take the divergent full step. ───────
    //    Return the best point seen; if none improved, stay put and signal stall.
    const bool any_improvement = (best_norm < norm0);
    if (improved) *improved = any_improvement;
    if (accepted_grad && any_improvement) *accepted_grad = best_grad;
    return best_x;
}

// ── Cached linear solver (B4) ─────────────────────────────────────────────────
//
// solve_with_fallback rebuilds a fresh SimplicialLDLT every call, re-running the
// symbolic analysis (fill-reducing reordering / COLAMD) each Newton iteration —
// even though the Hessian sparsity pattern is iteration-invariant for a fixed
// mesh + DOF layout.  This object keeps one persistent LDLT and runs
// `analyzePattern` once, then only `factorize` per iteration.
//
// FD-built Hessians (hyper-ideal, inversive-distance) prune near-zero triplets,
// so their pattern can shift slightly between iterations; we detect that via a
// change in `nonZeros()` and re-analyze, which keeps the cache correct and
// self-healing.  The SparseQR fallback for singular/rank-deficient systems is
// preserved verbatim (built fresh on the rare path).  The fallback semantics are
// identical to `solve_with_fallback`, so behaviour is unchanged.
struct NewtonLinearSolver {
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
    bool          analyzed = false;
    Eigen::Index  last_nnz = -1;
    bool          fallback_used = false;   ///< true iff the last solve used SparseQR

    Eigen::VectorXd solve(const Eigen::SparseMatrix<double>& A,
                          const Eigen::VectorXd&             rhs,
                          bool&                              ok)
    {
        fallback_used = false;
        if (!analyzed || A.nonZeros() != last_nnz) {
            ldlt.analyzePattern(A);
            analyzed = true;
            last_nnz = A.nonZeros();
        }
        ldlt.factorize(A);
        if (ldlt.info() == Eigen::Success) {
            Eigen::VectorXd x = ldlt.solve(rhs);
            if (ldlt.info() == Eigen::Success) { ok = true; return x; }
        }
        // Fallback: SparseQR — handles singular/rank-deficient A (gauge modes).
        fallback_used = true;
        Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr(A);
        if (qr.info() == Eigen::Success) {
            Eigen::VectorXd x = qr.solve(rhs);
            if (qr.info() == Eigen::Success) { ok = true; return x; }
        }
        ok = false;
        return Eigen::VectorXd::Zero(rhs.size());
    }
};

// ── Unified Newton core (H2) ──────────────────────────────────────────────────
//
// All five DCE Newton solvers (Euclidean, Spherical, HyperIdeal, CP-Euclidean,
// Inversive-Distance) ran a near-identical loop that differed only in (a) the
// gradient function, (b) the Hessian function, and (c) the spherical sign quirk
// (the concave spherical energy factorises −H).  This template captures that one
// loop so a fix lands once instead of five times.
//
// `grad(x) -> std::vector<double>`        : the energy gradient G(x).
// `hess(x) -> Eigen::SparseMatrix<double>`: the TRUE Hessian H(x) (un-negated).
// `concave`                               : if true, solve (−H)·Δx = G (spherical);
//                                           otherwise H·Δx = −G.  Both are the same
//                                           Newton system; `concave` only selects
//                                           the PSD matrix actually factorised.
//
// Folds in B2 (reuse the line-search gradient as the next convergence gradient),
// B4 (cached symbolic factorization), and B5 (no dead solver variable).
template <typename GradFn, typename HessFn>
inline NewtonResult newton_core(
    std::vector<double> x,
    GradFn&&            grad,
    HessFn&&            hess,
    bool                concave,
    double              tol,
    int                 max_iter)
{
    const int n = static_cast<int>(x.size());

    NewtonResult res;
    res.converged     = false;
    res.iterations    = 0;
    res.grad_inf_norm = 0.0;

    NewtonLinearSolver solver;

    // B2: gradient carried across the loop; the first one is the only "extra"
    // evaluation — every later iteration reuses the line-search gradient.
    std::vector<double> G_std = grad(x);

    for (int iter = 0; iter < max_iter; ++iter) {
        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);

        const 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;
        }

        Eigen::SparseMatrix<double> H = hess(x);

        // Newton system: factor the PSD matrix (H for convex, −H for concave).
        Eigen::SparseMatrix<double> A = concave ? Eigen::SparseMatrix<double>(-H) : H;
        Eigen::VectorXd rhs = concave ? Eigen::VectorXd(G) : Eigen::VectorXd(-G);

        bool ok = false;
        Eigen::VectorXd dx = solver.solve(A, rhs, ok);
        if (!ok) break;

        // Globalised line search (Armijo + steepest-descent fallback).
        // d_sd = −(H·G) is the merit-function steepest descent for f = ½‖G‖²,
        // using the TRUE (un-negated) Hessian for both signs.
        const double norm0 = G.norm();
        Eigen::VectorXd d_sd = -(H * G);
        bool improved = true;
        std::vector<double> G_next;
        x = detail::line_search(x, dx, d_sd, norm0, grad, &improved, &G_next);
        if (!improved) break;            // line search stalled — stop cleanly

        G_std = std::move(G_next);       // B2: reuse for the next convergence check
        res.iterations = iter + 1;
    }

    // Final gradient norm (reached only on max-iter / break — recomputed to be
    // correct after a stalled or failed step).
    auto G_final = grad(x);
    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 detail

// ── Euclidean Newton solver ────────────────────────────────────────────────────

/// Solve the Euclidean discrete conformal problem: find u ∈ ℝ^V such that
/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v.
///
/// Starting from x0, Newton's method minimises E(u) (the Euclidean DCE energy,
/// which is convex) by iterating  u ← u − H⁻¹·G  with backtracking line search.
/// The Hessian H is the cotangent Laplacian — PSD with one zero eigenvalue on
/// closed surfaces (gauge mode).  A SparseQR fallback handles this automatically.
///
/// \param mesh      Input triangulated surface (edges must carry lambda0 + theta_v).
/// \param x0        Initial DOF vector (length = number of free vertices).
///                  Pass all-zeros for a flat start (typical).
/// \param m         EuclideanMaps: lambda0[e], theta_v[v], v_idx[v] must be set.
///                  Call setup_euclidean_maps() + compute_euclidean_lambda0_from_mesh()
///                  + enforce_gauss_bonnet() before passing here.
/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
/// \param max_iter  Maximum Newton iterations.  Default: 200.
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \note On closed meshes without a pinned vertex, SimplicialLDLT detects the
///       gauge singularity and falls back to SparseQR automatically.
///
/// \see doc/math/discrete-conformal-theory.md §3 for the mathematical background.
inline NewtonResult newton_euclidean(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const EuclideanMaps&       m,
    double                     tol      = 1e-8,
    int                        max_iter = 200)
{
    // Layout is loop-invariant — decide the Hessian variant once (B3): cyclic
    // layout (edge DOFs present) → analytic cyclic Hessian, which covers the
    // vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian omits.
    bool has_edge_dof = false;
    for (auto e : mesh.edges())
        if (m.e_idx[e] >= 0) { has_edge_dof = true; break; }

    return detail::newton_core(
        std::move(x0),
        [&](const std::vector<double>& xc) { return euclidean_gradient(mesh, xc, m); },
        [&](const std::vector<double>& xc) {
            return has_edge_dof ? euclidean_hessian_analytic(mesh, xc, m)
                                : euclidean_hessian(mesh, xc, m);
        },
        /*concave=*/false, tol, max_iter);
}

// ── Spherical Newton solver ───────────────────────────────────────────────────

/// Solve the spherical discrete conformal problem: find u ∈ ℝ^V such that
/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v (genus-0 / sphere-like surfaces).
///
/// The spherical DCE energy is *concave*, so the Hessian H is NSD at the solution.
/// The solver factorises −H (which is PSD) and solves (−H)·Δx = G.
/// A gauge vertex must be pinned (set v_idx = -1) to remove the rotational mode.
///
/// \param mesh      Input triangulated surface, genus 0.
/// \param x0        Initial DOF vector (length = free vertices, excluding gauge_vertex).
///                  All-zeros is a good start.
/// \param m         SphericalMaps: lambda0[e], theta_v[v], v_idx[v], gauge_vertex set.
///                  Call setup_spherical_maps() + compute_spherical_lambda0_from_mesh()
///                  + enforce_gauss_bonnet() (checks Σ(2π-Θ) > 0) before passing here.
/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
/// \param max_iter  Maximum Newton iterations.  Default: 200.
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \note Unlike the Euclidean solver, the spherical solver does NOT need a SparseQR
///       fallback — the gauge vertex pins the null mode directly.
///
/// \see doc/math/geometry-modes.md §Spherical for sign-convention details.
inline NewtonResult newton_spherical(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const SphericalMaps&       m,
    double                     tol      = 1e-8,
    int                        max_iter = 200)
{
    // Concave energy: the Hessian H is NSD, so newton_core factors −H (PSD) and
    // solves (−H)·Δx = G — the same Newton system, just the sign that keeps LDLT
    // well-defined.  The merit steepest-descent inside still uses the true H.
    return detail::newton_core(
        std::move(x0),
        [&](const std::vector<double>& xc) { return spherical_gradient(mesh, xc, m); },
        [&](const std::vector<double>& xc) { return spherical_hessian(mesh, xc, m); },
        /*concave=*/true, tol, max_iter);
}

// ── HyperIdeal Newton solver ──────────────────────────────────────────────────

/// Solve the hyper-ideal discrete conformal problem: find (b, a) ∈ ℝ^{V+E} such that
/// Σ β_v(b,a) = Θ_v  and  Σ α_e(b,a) = θ_e  for all vertices v and edges e.
///
/// Used for genus-g surfaces (g ≥ 1) under hyperbolic cone metrics.
/// The energy is *strictly convex* (Springborn 2020, Theorem 1.3), so Newton
/// converges globally from any starting point.
///
/// DOF layout: first V_free entries are vertex variables b_v (hyper-ideal radii),
/// followed by E entries for edge variables a_e (intersection angles).
/// Use assign_hyper_ideal_all_dof_indices(mesh, maps) to set v_idx and e_idx automatically —
/// no vertex needs to be pinned.
///
/// \param mesh      Input triangulated surface, genus g ≥ 1.
/// \param x0        Initial DOF vector (length = V + E).  All-zeros typical.
/// \param m         HyperIdealMaps: lambda0[e], theta_v[v], v_idx[v], e_idx[e] set.
///                  Call setup_hyper_ideal_maps() + compute_hyper_ideal_lambda0_from_mesh().
/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
/// \param max_iter  Maximum Newton iterations.  Default: 200.
/// \param hess_eps  Finite-difference step for Hessian approximation.  Default: 1e-5.
///                  (Phase 9b will replace this with an analytic Hessian.)
/// \param clamp     Vertex-scale floor mode (N3 audit).  Default `HardJava`
///                  reproduces the Java oracle's hard `b<0 → floor` snap
///                  (C⁰, parity-faithful).  `SmoothBarrier` uses the C¹
///                  softplus floor — smoother near the feasibility boundary
///                  but deviates from the Java golden values.  See
///                  `HyperIdealScaleClamp` in hyper_ideal_functional.hpp.
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \see Springborn (2020), Theorem 1.3 for the strict convexity proof.
/// \see doc/math/geometry-modes.md §Hyper-ideal for DOF layout details.
inline NewtonResult newton_hyper_ideal(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const HyperIdealMaps&      m,
    double                     tol      = 1e-8,
    int                        max_iter = 200,
    double                     hess_eps = 1e-5,
    HyperIdealScaleClamp       clamp    = HyperIdealScaleClamp::HardJava)
{
    // block-FD Hessian (~33–1166× faster than full-FD); `clamp` selects the
    // vertex-scale floor mode (N3).  Convex energy → factor H directly.
    return detail::newton_core(
        std::move(x0),
        [&](const std::vector<double>& xc) {
            return evaluate_hyper_ideal(mesh, xc, m, /*energy=*/false, /*grad=*/true, clamp).gradient;
        },
        [&](const std::vector<double>& xc) {
            return hyper_ideal_hessian_block_fd_sym(mesh, xc, m, hess_eps, clamp);
        },
        /*concave=*/false, tol, max_iter);
}

// ── 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 2015 §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-2015 mapping.
inline NewtonResult newton_cp_euclidean(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const CPEuclideanMaps&     m,
    double                     tol      = 1e-8,
    int                        max_iter = 200)
{
    // Convex energy with an exact analytic Hessian (BPS-2015 2×2-per-edge).
    return detail::newton_core(
        std::move(x0),
        [&](const std::vector<double>& xc) { return cp_euclidean_gradient(mesh, xc, m); },
        [&](const std::vector<double>& xc) { return cp_euclidean_hessian(mesh, xc, m); },
        /*concave=*/false, tol, max_iter);
}

// ── 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.
///
/// The Hessian is computed by **per-face block finite differences**
/// (`inversive_distance_hessian_block_fd_sym`, same pattern as the HyperIdeal
/// solver after Phase 9b) — O(F) face evaluations instead of the O(n·F) of the
/// full-FD baseline.  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)
{
    // Convex energy; block-FD Hessian (~n/6× faster than full-FD).
    return detail::newton_core(
        std::move(x0),
        [&](const std::vector<double>& xc) { return inversive_distance_gradient(mesh, xc, m); },
        [&](const std::vector<double>& xc) {
            return inversive_distance_hessian_block_fd_sym(mesh, xc, m, hess_eps);
        },
        /*concave=*/false, tol, max_iter);
}

} // namespace conformallab