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 $β<β_{\min}:=\min\{β_c,β_G\}$, 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 $β_{\min}$; 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 $β_c$, which is the static phase transition point; while for $A$ non-concave, if $β_G<β_c$, the algorithms work up to the known algorithmic threshold $β_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.