Integrating over Dirac delta TL;DR: How does one show that $(\clubsuit)$ holds.

Some context and how I arrived at my problem. In a Thermodynamics problem set I was asked to calculate the partition function of the ideal gas in the microcanonic ensemble using the Hamiltionian
$$\mathcal{H} = \frac{1}{2m}\sum_{i=1}^N p_i^2.$$
The problem seems quite straight forward, but I fail at a certain point when trying to solve the integral:
$$\begin{align}Z &= \int_{\mathbb{R}^{6N}} \frac{d^{3N}pd^{3N}q}{h^{3N}N!}\delta\left(E-\mathcal{H}\right)\\
&= \frac{1}{h^{3N}N!}\int_{\mathbb{R}^{3N}}d^{3N}q \int_{\mathbb{R}^{3N}} d^{3N}p\, \delta\left(E-\sum_{i=1}^N \frac{p_i^2}{2m}\right)\\ 
&= \frac{V^N(2m)^{3N/2}}{h^{3N}N!}\int_{\mathbb{R}^{3N}}d^{3N}u\, \delta \left(E-\sum_{i=1}^Nu_i^2\right)\\ 
&= \frac{V^N(2m)^{3N/2}}{h^{3N}N!} \int_0^\infty \delta (E-r^2)r^{3N-1}dr \int d^{3N-1}\Omega\\ 
&{=} \frac{V^N(2m)^{3N/2}}{h^{3N}N!}\frac { 2 \pi ^ { 3 N / 2 } } { \Gamma \left( \frac { 3 N } { 2 } \right) }\int_0^\infty \delta (E-r^2)r^{3N-1}dr.\end{align}$$
And here is where the trouble starts. I'm not really sure on how to proceed from here on. It seems the right thing to do would be a coordinate transformation to get something of the form $\delta(x-r)$, but I just can't figure out how to do that. Comparing my current result with the solution I expect something like
$$\int_0^\infty \delta (E-r^2)r^{3N-1}dr=\frac { 1 } { 2 } E ^ { 3 N / 2 - 1 }.\tag{$\clubsuit$}$$
 A: You have the integral$$\int_0^\infty \delta (E-r^2)r^{3N-1}dr\tag1$$
Then let $g(r)=-r^2+E$. Of course, $g(r)$ is only zero when $r=\sqrt E$ ($E>0$ since it is an energy, and we don't care about $r=-\sqrt E$ since this is not within out integration bounds), so we use the definition for the composition of the Dirac delta function with a different function given on Wikipedia to get the delta function into a more managable form.
$$\delta(g(x))=\sum_i\frac{\delta(x-x_i)}{|g'(x_i)|}$$ where the sum is over all the roots of $g(x)$.
For our function $g(r)$, this gives $$g(-r^2+E)=\frac{\delta(r-\sqrt E)}{2\sqrt E}$$
Plugging this into $(1)$, $$\int_0^\infty \delta (E-r^2)r^{3N-1}dr=\frac1{2\sqrt E}\int_0^\infty\delta(r-\sqrt E)r^{3N-1}dr=\frac12 E^{-1/2}E^{3N/2} E^{-1/2}=\frac12E^{{3N}/2-1}$$ as required.
A: User John Doe has already given a good answer. Here is a slightly different calculation:
$$ \int_{\mathbb{R}_+}\! \mathrm{d}r~r^{3N-1}\delta(r^2-E)
~\stackrel{u=r^2}{=}~
\frac{1}{2}\int_{\mathbb{R}_+}\! \mathrm{d}u~u^{3N/2-1}\delta(u-E)
~=~\frac{1}{2}\theta(E)E^{3N/2-1},\tag{$\clubsuit$}$$
where $\theta$ denotes the Heaviside step function. 
