I'm trying to follow this answer, and got stuck in the use of the gamma function in the third step:
\begin{align} \operatorname{E}[XY] &= \sum_{y=0}^n \int_{x=0}^1 xy \binom{n}{y} x^y (1-x)^{n-y} \, dx \\ &= \sum_{y=0}^n y \binom{n}{y} \int_{x=0}^1 x^{y+1} (1-x)^{n-y} \, dx \\ &= \sum_{y=0}^n y \binom{n}{y} \frac{\Gamma(y+2)\Gamma(n-y+1)}{\Gamma(n+3)} \\ \end{align}
My guess is that this is a well-known method to eliminate the integral, and then greatly simplify things by considering that $y$ is a natural number.
The equation, though, is rather intimidating. I see for instance that $\Gamma(y+2) =\int_0^\infty x^{(y+1)}e^{-x}dx,$ which is promising, but there is that pesky exponential, plus the limits of integration are different, let alone the fact that there are now three integrals at play.