Dirac delta function $\delta(f(x))$ of function $f$ with a higher-order zero Dirac delta function have this property:
\begin{equation}
\delta(f(x))=\textstyle \sum_i\frac{\delta(x-a_i)}{\lvert f^\prime(a_i)\rvert}.
\end{equation}
And its derivation is:
\begin{eqnarray}
\int_{-\infty}^{\infty}g(x)\delta(f(x))&=&\sum_i\int_{a_i-\epsilon}^{a_i+\epsilon}g(x)\delta(f(x)),\qquad f(a_i)=0\cr
&=&\sum_i\int_{f(a_i-\epsilon)}^{f(a_i+\epsilon)}g(f^{-1}(y))\delta(y)\frac{1}{\lvert f^\prime(f^{-1}(y))\rvert}dy,\qquad x=f^{-1}(y)\cr
&=&\sum_i\frac{a_i}{\lvert f^\prime(a_i)\rvert}.\cr
&&\therefore\ \delta(f(x))=\sum_i\frac{\delta(x-a_i)}{\lvert f^\prime(a_i)\rvert}.
\end{eqnarray}
Here are questions:


*

*How about if $f^\prime(a_i)$ is 0?

*I think $\int_{-\infty}^{\infty}g(x)\delta(f(x))=\sum_i g(a_i)$ because Dirac delta function $\delta(f(x))$ in LHS is zero at $a_i$, and so $\int\delta(x)f(x)\ dx$ is just $f(x)$. What's wrong with this?

 A: *

*Firstly, OP is essentially asking the following question. 

What is $$I~=~\int_{\mathbb{R}}\! dx~g(x)~\delta(f(x)) \tag{A}$$ if the function $f$ only has isolated zeros $a_i$ with multiplicity $n_i$, where the index $i$ labels the zeros?

That is a good question. Below we will take a heuristic approach, and leave it to the reader to insert adequate assumptions to justify the various steps in the calculation. The answer is then

$$I~=~\sum_i \frac{(n_i-1)!}{|f^{(n_i)}(a_i)|} \lim_{y\to 0^+}\frac{g(a_i+\sqrt[n_i]{y})+g(a_i-\sqrt[n_i]{y})}{2y^{1-\frac{1}{n_i}}}. \tag{B}$$

Obviously, the limits in eq. (B) might not exist.
Sketched proof of formula (B). Since each zero is isolated, it is enough to just consider one zero $a$. By translational symmetry of the Lebesgue measure, we can assume that the zero is at the origin $a=0$. It follows that it is enough to consider $$f(x)~=~Ax^n, \tag{C}$$  where $$A~=~\frac{f^{(n)}(0)}{n!}~\neq~ 0.\tag{D}$$ We can throw the constant $A$ outside of the Dirac delta distribution as $\frac{1}{|A|}$. 
Then
$$I~=~\frac{1}{|A|}\int_{\mathbb{R}}\! dx~g(x)\delta(x^n)
~=~\frac{1}{|A|}\int_{\mathbb{R}_+}\! dx~[g(x)+g(-x)]\delta(x^n)$$ 
$$~\stackrel{y=x^n}{=}~\frac{1}{|A|}\int_{\mathbb{R}_+}\! \frac{dy}{ny^{1-\frac{1}{n}}}~[g(\sqrt[n]{y})+g(-\sqrt[n]{y})]\delta(y)
~=~\lim_{y\to 0^+}\frac{g(\sqrt[n]{y})+g(-\sqrt[n]{y})}{2n|A|y^{1-\frac{1}{n}}}, \tag{E}$$
which reproduces the familiar result 
$$I~=~\frac{g(0)}{|f^{\prime}(0)|}\quad\text{for}\quad n~=~1.\tag{F}$$ 
$\Box$

*Secondly, OP is asking what is wrong with his second calculation? For starters, it seems that OP is forgetting the Jacobian factor under integration by substitution.
A: If $f'(a_i) = 0$ then the substitution $x = f^{-1}(y)$ is invalid because $f^{-1}$ is not differentiable in $a_i$.
