Let $E = \{(x,y) \in \mathbb{R}^2 \mid x^2+y^2 \leq 1 \}$. Compute $$ \int_E x^n y^m \ dx \ dy$$ for all $n,m \in \mathbb{N}$.
What is the best way to do this? One can, of course, transform it using polar coordinates which leaves some trigonometric integral that can be solved with integration by parts and trigonometric formulas.
That approach seems rather nasty, though. Is there a nicer way?
You can do it in many ways, as the resulting beta-function has many equivalent representations.
As already was noted in the comments we can assume both $n\mapsto2n$ and $m\mapsto2m$ to be even and restrict the integration domain to the first quadrant.
For example: $$\begin{array}{} \frac I4&=\int\limits_0^1x^{2n}dx\int\limits_0^\sqrt{1-x^2}y^{2m}dy =\frac{1}{2m+1}\int\limits_0^1 x^{2n}(1-x^2)^{m+\frac12}dx\\ &\stackrel{x^2\mapsto t}=\frac{1}{2(2m+1)}\int\limits_0^1 t^{n-\frac12}(1-t)^{m+\frac12}dt=\frac{B(n+\frac12,m+\frac32)}{2(2m+1)}, \end{array} $$ or: $$\begin{array}{} \frac I4&=\int\limits_0^1 r^{2n+2m+1}dr\int\limits_0^{\pi/2}\cos^{2n}\phi\sin^{2m}\phi\;d\phi =\frac{B(n+\frac12,m+\frac12)}{4(n+m+1)}. \end{array} $$
The results seemingly differ but in fact they are the same due to: $$ B(x,y+1)=B(x,y)\frac{y}{x+y}, $$ and can be expressed in terms of more "elementary" functions as: $$\begin{array}{} I(n,m)=\frac{B(n+\frac12,m+\frac12)}{n+m+1} =\frac1{n+m+1}\frac{\Gamma(n+\frac12)\Gamma(m+\frac12)}{\Gamma(n+m+1)} =\frac{(2n-1)!!(2m-1)!!}{2^{n+m}(n+m+1)!}\pi.\end{array} $$