  We investigate the precision with which the Lyman-$\alpha$ damping wing
signature imprinted on the spectra of high-redshift quasars (QSOs) by the
foreground neutral intergalactic medium (IGM) can measure the history of cosmic
reionization. We leverage a novel inference pipeline based on a generative
probabilistic model for the entire spectrum (both red- and blueward of the
Lyman-$\alpha$ line), accounting for all relevant sources of uncertainty - the
stochasticity caused by patchy reionization, the impact of the quasar's
ionizing radiation on the IGM, it's unknown intrinsic spectrum, and spectral
noise. Performing fast JAX-based Hamiltonian Monte-Carlo (HMC) parameter
inference, we precisely measure the underlying global IGM neutral fraction as
well as the lifetime of the quasar. Running a battery of tests on over a
thousand mocks, we find optimal precision when running the pipeline with a six
parameter PCA continuum model (five coefficients and a normalization) on
$\mathrm{S}/\mathrm{N} \sim 10$ spectra, binned to a $\sim
500\,\mathrm{km}/\mathrm{s}$ velocity pixel scale, and extending at least out
to the C IV $\lambda\,1549\,\text{\r{A}}$ emission line. After marginalizing
out nuisance parameters associated with the quasar continuum, a single spectrum
constrains the IGM neutral fraction to $28.0_{-8.8}^{+8.2}\,\%$ and the quasar
lifetime to $0.80_{-0.55}^{+0.22}\,\mathrm{dex}$, improving notably towards
spectra with a stronger IGM damping wing imprint. Higher precision can be
achieved by averaging over statistical quasar samples. We identify two primary
sources of uncertainty that contribute approximately equally to the total error
budget: the uncertain quasar continuum model and the stochastic distribution of
neutral regions arising from both the reionization topology and the location of
the quasar's ionization front.
