Dirac delta function $\delta(f(x))$ of function $f$ with a higher-order zero
Solution 1:
-
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.
Solution 2:
If $f'(a_i) = 0$ then the substitution $x = f^{-1}(y)$ is invalid because $f^{-1}$ is not differentiable in $a_i$.