  The continuous random energy model (CREM) is a toy model of spin glasses on
$\{0,1\}^N$ that, in the limit, exhibits an infinitely hierarchical correlation
structure. We give two polynomial-time algorithms to approximately sample from
the Gibbs distribution of the CREM in the high-temperature regime, based on a
Markov chain and a sequential sampler. The running time depends algebraically
on the desired TV distance and failure probability and exponentially in
$(1/g')^{O(1)}$, where $g'$ is the gap to a certain inverse temperature
threshold; this contrasts with previous results which only attain $o(N)$
accuracy in KL divergence. If the covariance function $A$ of the CREM is
concave, the algorithms work up to the critical threshold $\beta_c$, which is
the static phase transition point; moreover, for certain $A$, the algorithms
work up to the known algorithmic threshold $\beta_G$ proposed in Addario-Berry
and Maillard (2020) for non-trivial sampling guarantees. Our result depends on
quantitative bounds for the fluctuation of the partition function and a new
contiguity result of the ``tilted" CREM obtained from sampling, which is of
independent interest. We also show that the spectral gap is exponentially small
with high probability, suggesting that the algebraic dependence is unavoidable
with a Markov chain approach.
