Note: this is fixed in version 9.
When I perform the double integral in Mathematica,
Integrate[(x (1 - x))^z (y (1 - y))^z, {x, 0, 1}, {y, 0, 1}]
which should give
$$B(z+1,z+1)^2 = \frac{\Gamma(z+1)^4}{\Gamma\left(2(z+1)\right)^2}$$
where $B(x,y)$ is the Beta function and $\Gamma(z)$ is the Gamma function because the integral is a product of two Beta functions, I instead get this ratio of Gamma functions times the extra factor $(-1)^{2z}$. What is going on here? To make matters stranger, if I do the integral instead using two nested calls to Integrate (one to integrate out $x$ and one to integrate out $y$), I get the ratio of Gamma functions without the incorrect extra factor.
Assuming[Re[z]>0 && Im[z]==0, Integrate[(x (1 - x))^z (y (1 - y))^z, {x, 0, 1}, {y, 0, 1}]. $(-1)^{2z}$ should definitely not be there. $\endgroup$