A four-dimensional integral in Peskin & Schroeder

That's equivalent simply to $c\int dx/x$. Switch to the Euclidean spacetime, $k_0=ik_4$ where $(k_1,\dots k_4)$ is $k_E$; i.e. analytically continue in $k_0$ (Wick rotation). The integral is $$\int \frac{i\cdot d^4 k_E}{(2\pi)^4} \frac{1}{(k_E^2)^2} \exp(ik\cdot \epsilon)$$ So it's proportional to the Fourier transform of $1/k_E^4$. The original function is $SO(4)$ symmetric, so the Fourier transform must be symmetric as well and depend on $\epsilon^2$ only. Dimensional analysis implies that the result is dimensionless i.e. it must be a combination of a constant and $\ln(\epsilon^2)$. The logarithm is there with a nonzero coefficient so the constant only determines how to take the logarithm: it should properly be written as $\ln(\epsilon^2/\epsilon_0^2)$ for some constant $\epsilon_0$ with the same dimension.

The only remaining unknown is the coefficient and one gets $4\pi^2$ from the remaining integral. It's a sort of waste of resources to compute this special integral; it's better to compute the more general integrals in appendix A.4, see especially formulae (A.44)-(A.49) on page 807, which I won't copy here because that's why Peskin and Schroeder wrote the textbook.

I will give another approach to this identity. First, we notice that

$$\int\frac{d^4k}{(2\pi)^4}\,\frac{1}{(k^2-m^2)^2}e^{ik\cdot\epsilon}=-i\frac{\partial}{\partial m^2}D_F(x)\big|_{x=\epsilon}$$

For space-like vector $\epsilon^2=-r^2<0$, we have

$$D_F(x)=\frac{m}{4\pi^2r}K_1(mr)$$

whose derivation refers to Weinberg's book vol. 1, page 202. For $r\rightarrow 0$, the following expansion holds

$$ K_1(mr)=\frac{1}{mr}+\frac{mr}{2}\log\frac{mr}{2}$$

With this conditions, we finally obtain

$$\int\frac{d^4k}{(2\pi)^4}\,\frac{1}{(k^2-m^2)^2}e^{ik\cdot\epsilon}=\frac{i}{16\pi^2}\log\frac{1}{\epsilon^2}$$

Further another approach.
After Wick rotation($k^0=ik_E^0,\,\epsilon^0=i\epsilon_E^0$) the integral is
$$ I_1 \equiv \int \frac{i\cdot d^4 k_E}{(2\pi)^4} \frac{1}{(k_E^2)^2} \exp(-ik_E\cdot \epsilon_E) . $$ In (19.43) we want the quntity $$ \frac{\partial}{\partial \epsilon^{\gamma}}I_1. $$ Note that $$ \frac{1}{k_E^2}=\int_0^{\infty}du e^{-k_E^2 u}, $$ $$ \frac{1}{(k_E^2)^2}=\int_0^{\infty}du\int_0^{\infty}dv e^{-k_E^2 (u+v)}, $$ $$ I_1=\frac{i}{(2\pi)^4} \int d^4k_E \int_0^{\infty}du\int_0^{\infty}dv e^{-(u+v)k_E^2-i\epsilon_E \cdot k_E}. $$ $ \displaystyle \int_{-\infty}^{\infty} dk_E^i \exp[-(u+v)(k_E^i)^2-i\epsilon_E^i k_E^i] = \sqrt{\frac{\pi}{u+v}} \exp\left[ -\frac{(\epsilon_E^i)^2}{4(u+v)} \right] $ \begin{alignat}{2} \therefore I_1&=&& \frac{i}{(2\pi)^4} \int_0^{\infty}du\int_0^{\infty}dv \frac{\pi^2}{(u+v)^2} \exp\left[-\frac{\epsilon_E^2}{4(u+v)} \right] \\ &=&& \frac{i}{16\pi^2}I_2 \left(\frac{\epsilon_E^2}{4}\right) \end{alignat} where $\displaystyle I_2(x) \equiv \int_0^{\infty}du\int_0^{\infty}dv\, \frac{1}{(u+v)^2}\exp\left(-\frac{x}{u+v}\right) $. The calculation of $I_2$ is here. \begin{alignat}{2} I_1 &=&& \frac{i}{16\pi^2}\left( -\log\left(\frac{\epsilon_E^2}{4}\right) +\gamma-1+\lim_{M\to \infty}\log M +\cal{O}(\epsilon_E^2) \right) \\ &=&& \frac{i}{16\pi^2}\left( -\log\left(-\frac{\epsilon^2}{4}\right) +\gamma-1+\lim_{M\to \infty}\log M +\cal{O}(\epsilon^2) \right) \end{alignat} After $\epsilon \to 0$, we have
$$ \frac{\partial}{\partial \epsilon^{\gamma}}I_1 = \frac{i}{16\pi^2}\frac{\partial}{\partial \epsilon^{\gamma}}\log \frac{1}{\epsilon^2}. $$