Integrating $x(xy+(1-x)(1-y))^{n-1}$ over Unit Square I'm trying to integrate the following function ( which comes out of calculations regarding reliability polynomials of graphs ):
$$
x\,\left[\, xy + \left(\, 1 - x\,\right)\left(\, 1 - y\,\right)\,\right]
^{\, n - 1}
$$
over the unit square
$\left(\vphantom{\large A}\,\mbox{i.e.,}\ 0 \le x \le 1,\,\,\, 0 \le y \le 1\,\right)$. I haven't found anything through searching, and evaluating the expression in WolframAlpha does not finish within the standard time limit. 
However, when I ask Mathematica, there is an answer, which is quite complicated:

Is there analytical reasoning as to why this is the result? 
Edit: since we can assume $n \ge 0$ and is an integer, we can use Assuming in Mathematica:
f[n_] := Integrate[x (x*y + (1 - x) (1 - y))^(n - 1), {x, 0, 1}, {y, 0, 1}]
Assuming[n >= 0 && Element[n, Integers], f[n]]

to achieve:
$$-\,2^{-2 - n}\,\,\,
\frac{\left(\, -1\,\right)^{n}\,\beta\left(\,-1,-n,1 + n\,\right) + \beta\left(\, 2,1 + n,0\,\right) +
\pi\left[\,\mathrm{i} + \cot\left(\, n\pi\,\right)\right]}{n}
$$
where $\beta\left(\, a,b,c\,\right)$ is the incomplete beta function. I would still like to know how this is expression is derived.
 A: Put yourself in baricentric coordinates (which is always good when you are integrating over a symmetric domain).
$$
x = u + 1/2,\;\;y = v + 1/2
$$
then 
$$
\begin{gathered}
  f(u,v,n) = \left( {1/2 + u} \right)\left( {\left( {1/2 + u} \right)\left( {1/2 + v} \right) + \left( {1/2 - u} \right)\left( {1/2 - v} \right)} \right)^{\,n - 1}  =  \hfill \\
   = \left( {1/2 + u} \right)\left( {1/2 + 2\,u\,v} \right)^{\,n - 1}  =  \hfill \\
   = 1/2\left( {1 + \frac{\partial }
{{\partial \,v}}\left( {2\,u\,v} \right)} \right)\left( {1/2 + 2\,u\,v} \right)^{\,n - 1}  \hfill \\
   = 1/2\left( {\left( {1/2 + 2\,u\,v} \right)^{\,n - 1}  + \frac{1}
{n}\frac{\partial }
{{\partial \,v}}\left( {1/2 + 2\,u\,v} \right)^{\,n} } \right) \hfill \\ 
\end{gathered} 
$$
and from here integrate by parts and then - or straightly - apply binomial expansion
considering that $uv$ is odd in both variables, and so are its odd powers.
A: Approximation: You may use that either $xy\leq 1/4$ or $(1-x)(1-y)\leq 1/4$. For large $n$ you may then use the approximation $$(xy+(1-x)(1-y))^{n-1} \approx (xy)^{n-1} + ((1-x)(1-y))^{n-1}$$ so that the integral becomes:
$$ \int_0^1\int_0^1 x  (xy)^{n-1} +  (1 - (1-x))(1-x)^{n-1}(1-y))^{n-1} \; dx \;dy=  \int_0^1\int_0^1 x^{n-1}y^{n-1} \; dx\;dy = \frac{1}{n^2} $$ You may try to compare with numerics (or the Wolfram result) for different $n$ values.
