Correlation of random variables with joint PDF proportional to $x^{a-1}y^{b-1}(1-x-y)^{c-1} $ The random variables $X$ and $Y$ have joint PDF $$f(x,y)= \frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)}x^{a-1}y^{b-1}(1-x-y)^{c-1} $$ where $0 \leq x \leq 1 , 0 \leq y \leq 1, x+y < 1 $ where $a , b,c$ are positive real. Note that $\Gamma$ is the gamma function. 

Find the correlation coefficient of $X$ and $Y$. 

The final answer is supposed to be $$- \sqrt{\frac{ab}{(b+c)(a+c)}}$$
I first tried to find $$E(X)=\frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)} \int_0^1\int_0^1 x^{a}y^{b-1}(1-x-y)^{c-1} dxdy.$$ But I have no clue what the integral on the RHS is.
 A: To avoid some lengthy computations, one can make a systematical use of the fact that one is told that $f$ is a PDF hence one knows that:

For every positive $(a,b,c)$, the function $f_{a,b,c}$ defined by $$f_{a,b,c}(x,y)= K(a,b,c)x^{a-1}y^{b-1}(1-x-y)^{c-1}\mathbf 1_{x\gt0,y\gt0,x+y\lt1},$$ with $$K(a,b,c)=\frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)},$$ integrates to $1$.

In particular, for every $(n,m)$, $$E(X^nY^m)=\iint x^ny^mf_{a,b,c}(x,y)\mathrm dx\mathrm dy=\iint \frac{K(a,b,c)}{K(a+n,b+m,c)}f_{a+n,b+m,c}(x,y)\mathrm dx\mathrm dy,$$ that is, $$E(X^nY^m)=\frac{K(a,b,c)}{K(a+n,b+m,c)}=\frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)}\frac{\Gamma(a+n)\Gamma(b+m)}{\Gamma(a+b+c+n+m)}.$$ 
(User @DilipSarwate suggested this approach in a comment, before I wrote this answer.)
Thus, for every nonnegative integers $(n,m)$, introducing the notation $s=a+b+c$, one gets

$$E(X^nY^m)=\frac{(a+n-1)\cdots(a+1)a\times (b+m-1)\cdots(b+1)b}{(s+n+m-1)\cdots(s+1)s}.$$ 

This remark allows to compute without pain the moments $E(X)$, $E(Y)$, $E(XY)$, $E(X^2)$ and $E(Y^2)$, and to deduce the correlation coefficient of $X$ and $Y$.
One gets $$E(X)=\frac{a}s,\ E(Y)=\frac{b}s,\ E(X^2)=\frac{a(a+1)}{s(s+1)},\ E(Y^2)=\frac{b(b+1)}{s(s+1)},\ E(XY)=\frac{ab}{s(s+1)},$$ hence $$\mathrm{var}(X)=\frac{a(s-a)}{s^2(s+1)},\qquad \mathrm{var}(Y)=\frac{b(s-b)}{s^2(s+1)},$$ and finally the formula in your post 

$$\varrho_{X,Y}=-\sqrt{\frac{ab}{(s-a)(s-b)}}.$$

Safety-check: When $c=0$, one gets $\varrho_{X,Y}=-1$. Why is this a safety-check?
A: $\begin{align}
\text{You domain is:}
\\
\left.\begin{array}{c} 0\leq x\leq 1 \\ 0\leq y\leq 1 \\ x+y< 1 \end{array}\right\} & \equiv x\in [0,1], y\in [0, x)
\\[2ex]\text{So the expectation is:}
\\[2ex] \mathsf E[X] & = \frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)}\int_0^1 \int_0^{1-x} x\; f_{X,Y}(x,y)\,\operatorname{d}y \operatorname{d} x 
\\[1ex] & = \frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)}\int_0^1 \int_0^{1-x} x^a y^{b-1}(1-x-y)^{c-1}\operatorname{d}y\operatorname{d}x 
\\[2ex] & = \frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)}\int_0^1 x^a \int_0^{1-x}y^{b-1}((1-x)-y)^{c-1}\operatorname{d}y\operatorname{d}x 
\\[4ex]
\text{Then use the definition of a p.d.f.}
\\[2ex]
1 & = \frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)}\int_0^1 x^{a-1} \int_0^{1-x} y^{b-1}(1-x-y)^{c-1}\operatorname{d}y\operatorname{d}x
\\[2ex]\text{By substitution}
\\[2ex]
\therefore \mathsf E[X] & = \frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)}\frac{\Gamma(1+a)\Gamma(b)\Gamma(c)}{\Gamma(1+a+b+c)}
\\ & = \frac{a}{a+b+c}
\\[2ex]\text{And so on...}
\end{align}$
