1. Since $x\gg p$, we see that $\sin(px)$ is highly oscillatory. In fact, the integral becomes
$$\int_0 ^\infty \mathrm{d}p\ p \sin px \ e^{-it\sqrt{p^2 +m^2}}\sim \int_{-\infty} ^\infty \mathrm{d}p\ p\ e^{ipx-it\sqrt{p^2 +m^2}}$$
modulo some factor of $\pm2/i$. Observe now this integral resembles $\int f(p)\exp(g(p))\,\mathrm{d}p$. We find the point $\tilde{p}$ such that
$$g'(\tilde{p})=0.$$
Then just replace $g(p)$ with $g(\tilde{p})+\frac{1}{2}g''(\tilde{p})(p-\tilde{p})^{2}$ and carry out the integral as a moment of a Gaussian. For more on this approximation, see e.g. the relevant chapter of Hunter and Nachtergaele's book (freely and legally available).
2. This is just a Fourier transform of $p\,(p^{2}+m^{2})^{-1/2}$.
3. I assume you are referring to Eq (2.51) on page 27. We write the integral as
$$I(t) = \int^{\infty}_{m}\sqrt{E^{2}-m^{2}}e^{-iEt}\,\mathrm{d}E.$$
Peskin and Schroeder consider this integral as $t\to\infty$. If we consider a change of variables to
$$E^{2}-m^{2}=\mu^{2}\quad\Longrightarrow\quad \mathrm{d}E = \frac{\mu}{\sqrt{m^{2}+\mu^{2}}}\mathrm{d}\mu$$
We have
\begin{align}
I(t) &= \int^{\infty}_{0} \frac{\mu^{2}}{\sqrt{m^{2}+\mu^{2}}}e^{-it\sqrt{m^{2}+\mu^{2}}}\mathrm{d}\mu \\
&=\frac{1}{2}\int^{\infty}_{-\infty} \frac{\mu^{2}}{\sqrt{m^{2}+\mu^{2}}}e^{-it\sqrt{m^{2}+\mu^{2}}}\mathrm{d}\mu
\end{align}
Let
$$f(\mu) = \frac{\mu^{2}}{\sqrt{m^{2}+\mu^{2}}},\quad\mbox{and}\quad
\phi(\mu) = \sqrt{m^{2}+\mu^{2}}$$
so
$$I(t)=\frac{1}{2}\int^{\infty}_{-\infty}f(\mu)e^{-it\phi(\mu)}\mathrm{d}\mu.$$
Observe
$$f(\mu)=\mu\phi'(\mu).$$
As $t\to\infty$, the integral becomes highly oscillatory.
There are two ways to approach the problem from here. The first, unforgivably handwavy but faster: take the stationary phase approximation, and pretend that $f(\mu_{\text{crit}})$ is some arbitrary constant.
The critical points for $\phi$ are $\mu_{0}=0$ and $\mu_{\pm}=\pm im$. We only care about the real $\mu$, so we Taylor expand about $\mu_0$ to second order:
\begin{align}
\phi(\mu)&=\phi(0)+\frac{1}{2!}\phi''(0)\mu^{2}\\
&=m + \frac{1}{2m}\mu^{2}
\end{align}
We now approximate the integral as
$$I(t) \sim \int^{\infty}_{-\infty}f(c)e^{-itm}e^{-it\mu^{2}/2m}\mathrm{d}\mu \approx f(c)e^{-itm}\sqrt{\frac{4\pi m}{t}}.\tag{1}$$
The other approximation doesn't fix $f$. Observe $f(\mu)\sim|\mu|$, so we have
$$I(t) \sim e^{-itm} \int^{\infty}_{0}\mu e^{-it\mu^{2}/2m}\mathrm{d}\mu.$$
We have (using Fresnel integrals)
$$\int^{\infty}_{0}\mu e^{-it\mu^{2}/2m}\mathrm{d}\mu\sim \frac{im}{t}.$$
Hence
$$I(t)\sim\frac{im}{t}e^{-imt}.\tag{2}$$
This can be seen by partial integration
$$\frac{\partial}{\partial \rho}\sqrt{\rho^2-m^2}=\frac{\rho}{\sqrt{\rho^2-m^2}}$$
OP edit: More explicitly, we use this to write $(3)$ as
\begin{align}
D(x-y) &= \frac{1}{(2\pi)^2r}\int^\infty_m d\rho \frac{\partial}{\partial \rho}\sqrt{\rho^2-m^2} e^{-\rho r} \\
&= \frac{1}{(2\pi)^2r}\left[\sqrt{\rho^2-m^2} e^{-\rho r}\right]^\infty_m-\frac{1}{(2\pi)^2r}\int^\infty_m d\rho \sqrt{\rho^2-m^2} \frac{\partial}{\partial \rho}e^{-\rho r}\\
&= \frac{1}{(2\pi)^2}\int^\infty_m d\rho \sqrt{\rho^2-m^2} e^{-\rho r}\\
&= \frac{m}{(2\pi)^2r}K_1(mr)
\end{align}
Best Answer
In those paragraphs, P&S want to obtain the behaviour of the amplitude $D(x-y)$ in the limit $r \to \infty$. Starting from the integral $$ \tag{1} \label{int} \frac{-i}{2(2\pi)^2 r} \int_{-\infty}^\infty \text{d}p \frac{p e^{ipr}}{\sqrt{p^2 +m^2}}, $$ you can see that a limit for $r \to \infty$ is not well defined, since the exponential has an oscillating behaviour. Also, the final form of the integral can be recognized as a modified Bessel function of the second kind. See this Wikipedia page for the asymptotic expansion of this function. The modification of the integration contour is a way to obtain a more straightforward integral. I hope this answers your question 1.
About point 2., in order to evaluate this integral, P&S apply the Cauchy's theorem defining the integrand as a function of the complex variable $p$. The square root in the denominator causes the branch cut in Figure 2.3. When rewriting the contour integral as $$ \tag{2} \label{1} \int_{i \infty}^{i m} \dots + \int_{i m}^{i \infty} \dots $$ we have to be careful in writing the integrands, since the branch cut generates a discontinuity between the two sides of the cut. In fact the square root in the complex plane is a multi-valued function, and the argument of $p^2+m^2$ shifts of $2\pi$ when passing from the right of the branch cut to the left. Thus $$ \sqrt{p^2 + m^2} = \cases{|p^2+m^2|^{1/2} \exp\left(i\frac{\arg(p^2+m^2)}{2}\right), \\ |p^2+m^2|^{1/2} \exp\left(i\frac{\arg(p^2+m^2)+ 2\pi}{2}\right) = - |p^2+m^2|^{1/2} \exp\left(i\frac{\arg(p^2+m^2)}{2}\right),} $$ where the first value is at the right and the second at the left of the branch cut. You see that the difference is a minus sign.
Therefore, making the substitution $p = i \rho$, \eqref{1} becomes $$ i\int_\infty^m \text{d} \rho \frac{\rho e^{-r \rho}}{- \sqrt{\rho^2 - m^2}} + i\int_m^\infty \text{d} \rho \frac{\rho e^{-r \rho}}{\sqrt{\rho^2 - m^2}}, $$ which finally gives $$ 2i \int_m^\infty \text{d} \rho \frac{\rho e^{-r \rho}}{\sqrt{\rho^2 - m^2}}. $$ Substituting this in \eqref{int} you find your result. Note that if you had chosen a different contour, for example the analogous of the "pushed" contour of Figure 2.3 but with the U-turn at $p = 0 + i 0$, you would not have the phase difference in the integrands on the lines $(im, 0)$ and $(0, im)$, so these two pieces would have cancelled each other (see \eqref{1}) and you would have obtained the same result with $m$ as lower integration limit.