Explicit computation of singular part of two-loop sunrise diagram For $\phi^4$, there is two-loop self-energy contribution from sunrise (sunset) diagram.

The integration is 
$$
   I(p)=\int\frac{d^D p_1}{(2\pi)^D}\frac{d^Dp_2}{(2\pi)^D}\frac{1}{(p_1^2+m^2)(p_2+m^2)[(p-p_1-p_2)^2+m^2]}.
$$ 
I try to use trick like Feynman Parameter, but still cannot get the explicit result. I saw this question Two-loop regularization
and I found that the exact solution of above integration is very complex and need a lot of fancy technique. However in order to do the renormalization I need only singular part of $I(p)$ ($\epsilon$-expansion with $D=4-\epsilon$).
My question:


*

*Is there some easy and direct method just to compute the singular part of $I(p)$ i.e. the coefficient of $1/\epsilon^2$ and $1/\epsilon$? 


By the way, I saw this answer https://physics.stackexchange.com/a/79236/169288 But it only give an explicit computation of singular of $I^\prime(p^2=m^2)$ which is not my requirement. 
You just need to provide the literature or textbooks which show the explicit computation details. I can rarely find textbooks covering explicit two-loop computation.
PS: The following is my method by ordinary Feynman parameter trick.
$$  I(p)=\int\frac{d^D p_1}{(2\pi)^D}\frac{d^Dp_2}{(2\pi)^D}\frac{1}{(p_1^2+m^2)(p_2+m^2)[(p-p_1-p_2)^2+m^2]}$$
$$I(p)=\int\frac{d^D p_1}{(2\pi)^D}\frac{d^Dp_2}{(2\pi)^D} \int dx dy dz \delta(x+y+z-1) \frac{2}{\mathcal{D}^3}$$
with
$$\mathcal{D}=x p_1^2 +y p_2^2 +z(p-p_1-p_2)^2 +m^2 = \alpha k_1^2 + \beta k_2^2 + \gamma p^2 +m^2 $$
with 
$$\alpha = x+z $$
$$\beta = \frac{xy + yz +zx}{x+z}$$
$$\gamma = \frac{xyz}{xy + yz +zx}$$
$$k_1= p_1 + \frac{z}{x+z}(p_2-p)$$
$$k_2 =p_2 - \frac{xz}{xy + yz +zx} p.$$
And the Jacobian $\frac{\partial(k_1,k_2)}{\partial(p_1,p_2)}=1.$
$$I(p)=\int_0^1 dx dy dz \delta(x+y+z-1) \int\frac{d^D k_1}{(2\pi)^D}\frac{d^Dk_2}{(2\pi)^D} \frac{2}{(\alpha k_1^2 + \beta k_2^2 + \gamma p^2 +m^2 )^3} $$
$$I(p)=\int dx dy dz \delta(x+y+z-1) \int\frac{d^D k_1}{(2\pi)^D}\frac{d^Dk_2}{(2\pi)^D} \int_0^{+\infty} dt t^2 e^{-t(\alpha k_1^2 + \beta k_2^2 + \gamma p^2 +m^2) }.$$
Gaussian integral of $k_1$ and $k_2$
$$I(p)= \int dx dy dz \delta(x+y+z-1)\int_0^{+\infty} dt \frac{t^{2-D}}{ (4\pi)^D (\alpha \beta)^{D/2}}e^{-t(\gamma p^2 +m^2)}$$
$$I(p)=\int dx dy dz \delta(x+y+z-1) \frac{\Gamma(3-D)}{(4\pi)^D(\alpha \beta)^{D/2}(\gamma p^2 +m^2)^{3-D}}.$$
Use $D= 4-\epsilon$
$$I(p)=\frac{1}{(4 \pi)^4}\int dx dy dz \delta(x+y+z-1) \frac{\gamma p^2 +m^2}{(\alpha \beta)^2}\Gamma(-1+\epsilon) \left(\frac{\sqrt{\alpha \beta}}{\gamma p^2 +m^2}\right)^\epsilon$$
$$\Gamma(-1+\epsilon)= -\frac{1}{\epsilon} +\gamma_E-1 +\mathcal{O}(\epsilon)$$
$$\left(\frac{\sqrt{\alpha \beta}}{\gamma p^2 +m^2}\right)^\epsilon = 1 + \epsilon \ln\left(\frac{\sqrt{\alpha \beta}}{\gamma p^2 +m^2}\right)+\mathcal{O}(\epsilon^2).$$
Up to $0$th order of $\epsilon$
$$I(p)=\frac{1}{(4 \pi)^4}\int dx dy dz \delta(x+y+z-1) \frac{\gamma p^2 +m^2}{(\alpha \beta)^2}\left(-\frac{1}{\epsilon}- \ln\left(\frac{\sqrt{\alpha \beta}}{\gamma p^2 +m^2}\right) +\gamma_E-1\right) .$$
There is two integral I need to compute
$$I_1 = \int_0^1 dx dy dz \delta(x+y+z-1) \frac{\gamma p^2 +m^2}{(\alpha \beta)^2} $$
$$I_2 =  \int_0^1 dx dy dz \delta(x+y+z-1) \frac{\gamma p^2 +m^2}{(\alpha \beta)^2}\ln\left(\frac{\sqrt{\alpha \beta}}{\gamma p^2 +m^2}\right).$$
I can't find the explicit result of these two.
 A: If you can be satsfied with the $m=0$ case, the integral is easy provided you work in configuaration space rather that momentum space. The $x$-space propagator in $n$ dimensions  is 
$$
g(x,x') = \frac{1}{ (n-2)S_{n-1}} \left(\frac {1}{|x-x'|}\right)^{n-2}
$$
where $S_{n-1} = 2\pi^{n/2}/\Gamma(n/2)$ is the surface  area of the $n$-ball.
Your  Feynman diagram in configuration space  is
$$
I(p)= \int {d^nx} e^{ipx}[g(x,0)]^3 
$$ 
which you can evaluate (after re-naming  the following integrals integration variable $k$ as  $x$ and its  $x$ to the external momentum $p$) using the standard Fourier  integral
$$
\int \frac{d^n k}{(2\pi)^n} e^{ikx} |k^2|^s = \frac{4^s}{\pi^{n/2}}\frac{\Gamma(s+n/2)}{\Gamma(-s)} \frac{1}{|x|^{2s+n}}.
$$ 
The $m=0$ case is good enough  to get the $\beta$ function and the wavefunction renormalization $Z$
as these do not depend on $m$ anyway.  If the $m$ dependence is important to you you will need to evaluate the Fourier trasnform of the cube of a Bessel "K" function 
