Algorithm Overview
==================

**Root solver**: The inverter processes each output grid point independently (parallelized with OpenMP).
For a target point :math:`\mathbf{x}_\text{target} = (x_0, x_1, \ldots)`, the inversion solves a nonlinear system of the form:

.. math::

   \mathbf{F}(T, \mu_B, \mu_S, \mu_Q) =
   \begin{bmatrix}
   x_0(T,\mu) - x_{0,\text{target}} \\
   x_1(T,\mu) - x_{1,\text{target}} \\
   \vdots
   \end{bmatrix}
   = 0,

where :math:`(x_0, x_1, \ldots)` are the user-configured output grid axes (e.g. :math:`s, n_B, n_S, n_Q` or :math:`e, n_B, y_S, y_Q`).

The solver uses GSL's derivative modified-Powell hybrid method (``gsl_multiroot_fsolver_hybrids``) to iteratively refine :math:`(T, \mu_B, \mu_S, \mu_Q)` until the residual norm :math:`|\mathbf{F}|` falls below tolerance.
This method combines a Newton-type step with a trust-region (dogleg) globalization strategy. The Jacobian

.. math::

   J = \frac{\partial \mathbf{F}}{\partial (T, \mu_B, \mu_S, \mu_Q)},

is approximated internally via finite differences of :math:`\mathbf{F}`.

**Seed generation**: The seeds (or initial guesses) are generated for each EoS branch. The algorithm retrieves the :math:`k` nearest
table rows (in output space, measured by scaled Euclidean distance) as initial guesses.
The last accepted solution for that branch is also tried as a warm-start seed.

**Solution validation**: A converged solution is accepted only if
(a) all coordinates lie within the EoS domain,
(b) optional physicality filters pass (:math:`0 \leq c_s^2 \leq 1`, bounded :math:`|\mu_i|/T`), and
(c) the output-space envelope check passes (when ``enforce_output_envelope`` is enabled).

**Phase selection**: When multiple phase branches of the table each return a valid
solution, the inverter selects the thermodynamically preferred branch: the mixed
(unstable or meta-stable) phase is chosen when exactly one such solution exists;
if the table has no phase data, the single stable solution is used directly.

**Output**: The accepted :math:`(T, \mu_B, \mu_S, \mu_Q)` and all derived quantities
are written to the output table. Failed points are filled with NaN or excluded,
depending on ``output.on_inversion_failure``.