  We consider Ising mixed $p$-spin glasses at high-temperature and without
external field, and study the problem of sampling from the Gibbs distribution
$\mu$ in polynomial time. We develop a new sampling algorithm with complexity
of the same order as evaluating the gradient of the Hamiltonian and, in
particular, at most linear in the input size. We prove that, at sufficiently
high-temperature, it produces samples from a distribution $\mu^{alg}$ which is
close in normalized Wasserstein distance to $\mu$. Namely, there exists a
coupling of $\mu$ and $\mu^{alg}$ such that if $({\boldsymbol x},{\boldsymbol
x}^{alg})\in\{-1,+1\}^n\times \{-1,+1\}^n$ is a pair drawn from this coupling,
then $n^{-1}{\mathbb E}\{\|{\boldsymbol x}-{\boldsymbol
x}^{alg}\|_2^2\}=o_n(1)$. For the case of the Sherrington-Kirkpatrick model,
our algorithm succeeds in the full replica-symmetric phase.
  We complement this result with a negative one for sampling algorithms
satisfying a certain `stability' property, which is verified by many standard
techniques.
  No stable algorithm can approximately sample at temperatures below the onset
of shattering, even under the normalized Wasserstein metric. Further, no
algorithm can sample at temperatures below the onset of replica symmetry
breaking.
  Our sampling method implements a discretized version of a diffusion process
that has become recently popular in machine learning under the name of
`denoising diffusion.' We derive the same process from the general construction
of stochastic localization. Implementing the diffusion process requires to
efficiently approximate the mean of the tilted measure. To this end, we use an
approximate message passing algorithm that, as we prove, achieves sufficiently
accurate mean estimation.
