Critisism of Previous Answers
My personal feeling is that I am not 100 % happy with any of the answers posted as each of them either does something slightly mathematically dubious or in some cases I may not be familiar with a particular niche area of mathematics and therefore was not able to understand and interpret the answers provided.
I do not wish to make this statement without justifying it, and so I have the following comments in regards to previous answers.
KevinG's answer is the closest to the one I shall present here, however the original question was not answered to completion with rigour. The proof I show here explains the final steps summarized in KevinG's statement "Combine these with shifting to each of the desired roots, and you can obtain the equality you're looking for".
Kurkyl's answer is just a description of where the attempt in the original question "goes wrong".
I do not understand the answer by jbc. This appears to be a branch of mathematics with which I am not familiar.
robjohn's solution requires f(x) be invertable and bijective, which is not true in general.
Finally I personally feel J. Heller's answer is unnecessarily complicated - the use of non-standard or obscure definitions of the delta function may well be correct however is not likely to be clear to the average reader as it was not clear to me.
I hope this feedback is taken positively, it is not meant to be overly critical, however I thought I should justify my reasons for posting a new answer considering the number of answers already available.
As an aside, the context for the solution I will provide comes from a problem in Quantum Field Theory, and so without further ado I personally feel this method is the clearest.
#Prerequisite Proof
We shall require proof that
$$\int^{\pm\infty}\delta(ax)\;\mathrm{d}x=\frac{1}{\left|a\right|}\int^{\pm\infty}\delta(x)\;\mathrm{d}x$$
We proceed as follows; firstly define $y=ax$, hence $\mathrm{d}y=a\;\mathrm{d}x$ we obtain
$$\int^{x=\pm\infty}\delta(ax)\;\mathrm{d}x=\int^{y=a\cdot(\pm\infty)}\delta(y)\;\frac{\mathrm{d}y}{a}=
\frac{1}{a}\int^{\pm\infty}\delta(y)\;\mathrm{d}y$$
under the assumption that $a>0$. ($a$ is real and non-negative.) Now consider $a$ negative. We define $b=-a$, hence $b>0$ and the integral above is instead
$$\int^{x=\pm\infty}\delta(-bx)\;\mathrm{d}x=\int^{y=-b\cdot(\pm\infty)}\delta(y)\;\frac{\mathrm{d}y}{-b}=\int^{y=b\cdot(\mp\infty)}\delta(y)\;\frac{\mathrm{d}y}{-b}=
\frac{1}{b}\int^{\pm\infty}\delta(y)\;\mathrm{d}y$$
where we have exchanged the limits of integration in the final step. This integral can be written in terms of $a$, $b=\left|a\right|$
$$\frac{1}{\left|a\right|}\int^{\pm\infty}\delta(y)\;\mathrm{d}y$$
hence we have shown that for $a\in\Re$
$$\int^{\pm\infty}\delta(ax)\;\mathrm{d}x=\frac{1}{\left|a\right|}\int^{\pm\infty}\delta(x)\;\mathrm{d}x$$
Origional Question Proof
We wish to prove
$$\int^{\pm\infty}g(x)\delta(f(x))\;\mathrm{d}x=\sum_i\frac{g(a_i)}{\left|f^\prime(a_i)\right|}$$
where $f^\prime(a_i)=\left[\frac{\mathrm{d}f(x)}{\mathrm{d}x}\right]_{a_i}$, and $a_i$ are the zeros of $f(x)$.
The Taylor expansion of f(x) about the point $a_i$ is
$$f(x)=f(a_i)+f^\prime(a_i)(x-a_i)+\mathcal{O}(x^2)=f^\prime(a_i)(x-a_i)+\mathcal{O}(x^2)$$
as $f(a_i)=0$.
Inserting the expansion into the integral
$$\int^{\pm\infty}g(x)\delta\left(f^\prime(a_i)(x-a_i)\right)\;\mathrm{d}x=\sum_i\int^{a_i\pm\epsilon}g(x)\delta\left(f^\prime(a_i)(x-a_i)\right)\;\mathrm{d}x$$
To be explicit, introduce the new variable $y=x-a_i$, $\mathrm{d}y=\mathrm{d}x$,
$$\rightarrow\quad\sum_i\int^{y=0\pm\epsilon}g(y+a_i)\delta\left(f^\prime(a_i)y\right)\;\mathrm{d}y$$
using the expression from the above proof
$$\rightarrow\quad\sum_i\int^{y=0\pm\epsilon}\frac{g(y+a_i)}{\left|f^\prime(a_i)\right|}\delta\left(y\right)\;\mathrm{d}y=\sum_i\frac{g(a_i)}{\left|f^\prime(a_i)\right|}$$
hence we have obtained the desired result. I personally prefer not to introduce the change of variables from $x$ to $y$, in my original notes I have the expression
$$\sum_i\int^{a_i\pm\epsilon}g(x)\delta\left(f^\prime(a_i)(x-a_i)\right)\;\mathrm{d}x=\sum_i\int^{a_i\pm\epsilon}\frac{g(x)}{\left|f^\prime(a_i)\right|}\delta(x-a_i)\;\mathrm{d}x$$
this can be proved by extending the "Prerequisite Proof" above.
Answer to question in comments...
It has been 5 years since I posted this question and some years before I have done any particularly serious math.
Question was why are the higher order terms in the Taylor expansion ignored.
For the second order term in $x^2$, we have:
$$\int^{a_i\pm\epsilon}\delta(f^{\prime\prime}(a_i)(x-a_i)^2)\;\mathrm{d}x=\int^{a_i\pm\epsilon}\frac{1}{f^{\prime\prime}(a_i)}\delta((x-a_i)^2)\;\mathrm{d}x$$
$$\int^{a_i\pm\epsilon}\delta((x-a_i)^2)\;\mathrm{d}x=0$$
The last line can be proven by substituting variables (twice) $z=y^2$, $y=x-a_i$. Or alternatively you consider that the dirac delta is zero everywhere except at 0, and since we have an integral with a limit which approaches zero from above (positive) and then moves away heading in the positive direction, the integral must cancel itself. (Since there is a square.)
This raises questions about the cubic order term. But the change of variables will always bring in a factor of
$$z^{\frac{1}{n}}$$
where $n$ is the order of the term, so this will always be zero.
(I think)
