# Evaluating $\int_{0}^{1}\int_{0}^{1}\frac{r^{i+j}(1-r)^{k+l}s^{2m-i-j}(1-s)^{2m-k-l}}{(r+s)^{m}(2-r-s)^{m}}drds$

I'm trying to compute a closed form expression for the integral $$\int_{0}^{1}\int_{0}^{1}\frac{r^{i+j}(1-r)^{k+l}s^{2m-i-j}(1-s)^{2m-k-l}}{(r+s)^{m}(2-r-s)^{m}}drds \quad i,j,k,l \in\{0,1,\ldots,m\},$$ which occurs in calculating $$\mathrm{Ex}\left[\frac{R^{i}(1-R)^{k}S^{m-i}(1-S)^{m-k}}{(R+S)^{m}(2-R-S)^{m}}\right]$$ where $R$ and $S$ are independent random variables with $R\sim Beta(j+1,l+1)$, $S\sim Beta(m-j+1,m-l+1)$ and $i,j,k,l \in\{0,1,\ldots,m\}$.

Right now the only idea I have is to try to find a partial fraction decomposition that will allow me to compute the integral of the resulting terms using integration by parts. However, this would be a lengthy and tedious calculation. I would really appreciate any ideas or suggestions for computing this integral in a more elegant way.

-

Denote $a=i+j$ and $b=k+l$: $$I_{m}(a,b) = \int_0^1 \int_0^1 \frac{r^a (1-r)^b s^{2m-a} (1-s)^{2m-b}}{(r+s)^m (2-r-s)^m} \mathrm{d}r \,\mathrm{d}s$$ By interchanging $(r,s) \to (1-r,1-s)$ we get $I_m(a,b)=I_m(b,a)$. Similarly by $(r,s) \to (s,r)$, we get $I_m(a,b) = I_m(2m-a,2m-b)$. It therefore suffices to consider $0 \leqslant a \leqslant b \leqslant m$.
Quick experimentation with Mathematica suggests that the answer: $$I_m(a,b) = R_m(a,b) + Q_m(a,b) \ln(2)$$ where $R_m$ and $Q_m$ are rational numbers:
This seems to be do to a simple structure of the anti-derivative: $$\begin{eqnarray} \int \int \frac{r^a (1-r)^b s^{2m-a} (1-s)^{2m-b}}{(r+s)^m (2-r-s)^m} \mathrm{d}r \,\mathrm{d}s &=& \frac{P_1(r,s)}{(r+s)^{m-2}(2-r-s)^{m-2}} \\ && + \left(P_2(r) + P_3(s)\right) \log(r+s) \\ && + \left(P_4(r) + P_5(s) \right) \log(2-r-s) + C \end{eqnarray}$$ where $P_i$ are polynomials. It is very likely that the creative telescoping method can be used to establish recurrence equations satisfied by $I_m(a,b)$. Quite possibly the software packages developed at RISC-Linz may be useful to do this.