Integration of the splitting function I have a problem performing the following integration provided in the paper by Catani and Seymour (arXiv: hep-ph/9605323) page 27. Given is the integral 
$$
\mathcal{V}=\int_0^1 (z(1-z))^{-\epsilon} \int_0^1 (1-y)^{1-2\epsilon}y^{-1-\epsilon}V(z;y) dydz
$$
with 
$$
V(z;y)=\frac{2}{1-z(1-y)}-(1+z)-\epsilon(1-z).
$$
In the paper the exact result and an approximation is stated as follows: 
$$
\mathcal{V}=\frac{\Gamma(1-\epsilon)^3}{\Gamma(1-3\epsilon)}\left[\frac{1}{\epsilon^2}+\frac{1}{\epsilon}\frac{3+\epsilon}{2(1-3\epsilon)}\right]=\frac{1}{\epsilon^2}+\frac{3}{2\epsilon}+5-\frac{\pi^2}{2}+\mathcal{O}(\epsilon).
$$
To me it is obvious how to deal with the second and the third summand of $V$ under the integral. I simply expand and employ the definition of the Euler Beta function 
$$
B(a,b)=\int_0^1 t^{a-1}(1-t)^{b-1} dt=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}.
$$
Hence 
\begin{align}
\int_0^1 (z(1-z))^{-\epsilon} \int_0^1 (1-y)^{1-2\epsilon}y^{-1-\epsilon}\left[-(1+z)-\epsilon(1-z)\right] dydz
& = \frac{\Gamma(1-\epsilon)^3}{\Gamma(1-3\epsilon)}\left[\frac{1}{\epsilon}\frac{3+\epsilon}{2(1-3\epsilon)}\right]
\\ & = \frac{3}{2\epsilon}+5+\mathcal{O}(\epsilon)
\end{align}
for the simple part. 
Do you have any suggestions how to treat the singular term of $V$ in order to get the exact result? In particular how to solve 
$$
\int_0^1 (z(1-z))^{-\epsilon} \int_0^1 (1-y)^{1-2\epsilon}y^{-1-\epsilon}\left[\frac{2}{1-z(1-y)}\right] dydz=\frac{\Gamma(1-\epsilon)^3}{\Gamma(1-3\epsilon)}\left[\frac{1}{\epsilon^2}\right]=\frac{1}{\epsilon^2}-\frac{\pi^2}{2}+\mathcal{O}(\epsilon).
$$
Thank you in advance. 
Edit: After some effort in reverse engineering (basically rewriting the product of Gamma-functions in the exact result in terms of two Euler Beta functions) the initial problem boils down to showing the equality of
\begin{align}
& \int_0^1 (z(1-z))^{-\epsilon} \int_0^1 y^{1-2\epsilon}(1-y)^{-1-\epsilon}\left[\frac{2}{1-zy}\right] dydz \\
&=\int_0^1 (t(1-t))^{-\epsilon} \int_0^1 s^{1-2\epsilon}(1-s)^{-1-\epsilon}\left[\frac{(1-s)}{t(1-t)s^2}\right] dsdt.
\end{align}
I guess the solution is an integral transformation. Does anyone know how to transform in order to show the equality of the upper expressions?
 A: I'm not sure the extent to which you would like a complete analytical solution to the integral but since the integrations over $y$ and $z$ decouple you can proceed by identifying, say the integration over $z$, as a Hypergeometric2F1 which has the following integral representation $$_2F_1(a,b,c;t) = \frac{\Gamma(c)}{\Gamma(b) \Gamma(c-b)} \int_0^1 \mathrm{d}z \,z^{b-1} (1-z)^{c-b-1} (1-zt)^{-a}\,\,\,\text{for}\,\,\,\text{Re}(c)> \text{Re}(b) > 0.$$ Let $I$ be the integral of interest. In your case, with $t = 1-y$, you have $$I =  \frac{2\Gamma(1-\epsilon)^2}{\Gamma(2-2\epsilon)} \int_0^1 \mathrm{d}y\, (1-y)^{1-2\epsilon}y^{-1-\epsilon}\, _2F_1(1,1-\epsilon,2-2\epsilon;1-y).$$ Now, as $|t|<1$ for all $y \in (0,1)$, one is within the radius of convergence to permit use of the power series representation for the Hypergeometric2F1 in terms of so called Pochhammer symbols $$_2F_1(a,b,c;t) = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n}{(c)_n} \frac{t^n}{n!},$$ where $(a)_n = \Gamma(a+n)/\Gamma(a)$ with $(1)_n := n!$. It follows $$I = 2\Gamma(1-\epsilon) \sum_{n=0}^{\infty} \frac{\Gamma(1-\epsilon+n)}{\Gamma(2-2\epsilon+n)} \int_0^1 \mathrm{d}y\, (1-y)^{1-2\epsilon+n}\,y^{-1-\epsilon}.$$ The integral over $y$ is an Euler Beta type so you can readily write $$I = 2\Gamma(1-\epsilon)\Gamma(-\epsilon) \sum_{n=0}^{\infty}\frac{\Gamma(1-\epsilon+n)}{\Gamma(2-2\epsilon+n)}\frac{\Gamma(2-2\epsilon+n)}{\Gamma(2+n-3\epsilon)} = 2\Gamma(1-\epsilon)\Gamma(-\epsilon) \sum_{n=0}^{\infty}\frac{\Gamma(1-\epsilon+n)}{\Gamma(2+n-3\epsilon)}$$
At this point the explicit summation over $n$ can be input into say mma, returning a combination of Gamma functions containing another $\Gamma(-\epsilon)$ which overlaps with the $\Gamma(-\epsilon)$ already present to provide the double pole in $\epsilon$. Indeed you find $$I = \frac{1}{\epsilon^2} -\frac{\pi^2}{2} + O(\epsilon).$$
