  We show that in a certain, angle-averaged squeezed limit, the $N$-point
function of matter is related to the response of the matter power spectrum to a
long-wavelength density perturbation,
$P^{-1}d^nP(k|\delta_L)/d\delta_L^n|_{\delta_L=0}$, with $n=N-2$. By performing
N-body simulations with a homogeneous overdensity superimposed on a flat
Friedmann-Robertson-Lema\^itre-Walker (FRLW) universe using the \emph{separate
universe} approach, we obtain measurements of the nonlinear matter power
spectrum response up to $n=3$, which is equivalent to measuring the fully
nonlinear matter $3-$ to $5-$point function in this squeezed limit. The
sub-percent to few percent accuracy of those measurements is unprecedented. We
then test the hypothesis that nonlinear $N$-point functions at a given time are
a function of the linear power spectrum at that time, which is predicted by
standard perturbation theory (SPT) and its variants that are based on the ideal
pressureless fluid equations. Specifically, we compare the responses computed
from the separate universe simulations and simulations with a rescaled initial
(linear) power spectrum amplitude. We find discrepancies of 10\% at $k\simeq
0.2 - 0.5 \,h\,{\rm Mpc}^{-1}$ for $5-$ to $3-$point functions at $z=0$. The
discrepancy occurs at higher wavenumbers at $z=2$. Thus, SPT and its variants,
carried out to arbitrarily high order, are guaranteed to fail to describe
matter $N$-point functions ($N>2$) around that scale.
