This problem arises from the following property of Dirac $\delta-$function: $$\delta(f(x))=\sum_{a_i\in Z(f)}\frac{\delta(x-a_i)}{|\frac{df}{dx}(a_i)|} $$ where $Z(f):=\{x\in dom(f)|\,f(x)=0\}$, the zero-set of $f$.
Here, many people use the Taylor expansion of $f(x)$ to prove the property. Specifically, they expand $f(x)$ around those zero point of $f$, for example, $f(x)=f(a_i)+f'(a_i)(x-a_i)+\mathcal{O}((x-a_i)^2)$, where $f(a_i)=0$, and say that $\mathcal{O}((x-a_i)^2)$ is higher order infinitesimal to $(x-a_i)$ therefore they are neglegible, which sounds reasonable.
However, when $f(x)$ is the power of $x$, for example, $f(x)=x^2$. Then the above form tell us $\delta(x^2)=\frac{\delta(x-0)}{2\cdot0}$ which is ill-defined. Similarly to $x^n$.
One may define the extended version such that $\frac{1}{0}:=\infty$, then the whole $\delta(x^2)$ still make sense. However, if one tries to evaluate the following integral:
$$\int_{-\infty}^{\infty}x\delta(x^2)dx $$ by the above formula, it will give us $$\int_{-\infty}^{\infty}x\delta(x^2)dx=\int_{-\infty}^{\infty}x\frac{\delta(x)}{2\cdot 0}dx=\frac{0}{2\cdot 0} $$ which is an ill-defined result.
On the other hand, notice $xdx=\frac{1}{2}d(x^2)$
$$\int_{-\infty}^{\infty}x\delta(x^2)dx=\frac{1}{2}\int_{-\infty}^{\infty}\delta(x^2)d(x^2)=\frac{1}{2}$$ which I think this is just a coincidence. (the reason is one can have the integrand to be $x^3\delta(x^2)$, again use $xdx=\frac{1}{2}d(x^2)$, this time the first way generate the same wield formula, while the second way gives you $0$).
Can anyone explain it to me? (Notice, whenever the denominator is not zero, everything looks fine)
The formula \begin{equation} \delta(f(x))=\sum_{a_i\in Z(f)}\frac{\delta(x-a_i)}{\left|\frac{df}{dx}(a_i)\right|} \tag{1} \label{1} \end{equation} is only a formula for expressing the composition between a Dirac distribution and a function belonging to a well defined class of ordinary function. Its "ill-definiteness" simply due to the fact that the formula is valid not valid for the function $f(x)=x^2$. Precisely, equation \eqref{1} can be proved only under the following two hypotheses ([1], §1.9 pp. 22-23), both implicitly used in all the proposed answers to question "Dirac Delta Function of a Function":
The general definition of $\delta(f(x))$ in $\mathscr{D}^\prime$
The problem of defining of $\delta(f(x))$ for a general $f$ is a particular instance of the problem of defining the composition of an ordinary function with a distribution: this in turn is one of the motivating problems which gave rise to several "nonlinear theories" of generalized functions. The approach followed by Vladimirov ([1], §1.9 p. 22) for defining this composition in In $\mathscr{D}^\prime$ is the following one: $$ \delta(f(x))=\lim_{\varepsilon\to0+}\omega_\varepsilon(f(x))\,\text{ in }\mathscr{D}^\prime\iff \langle\delta(f),\varphi\rangle=\lim_{\varepsilon\to0+}\!\int\limits_{Z(f)\cap[a,b]}\!\!\!\!\omega_\varepsilon(f(x))\varphi(x)\mathrm{d}x\tag{2}\label{2} \quad\forall\varphi\in\mathscr{D}([a,b])$$ where $\mathscr{D}([a,b])$ is the space of $C^\infty$ functions whose support is contained in any interval $[a,b]$ of interest, and \begin{equation} \omega_\varepsilon(x)= \begin{cases} C_\varepsilon e^{-\frac{\varepsilon^2}{\varepsilon^2-x^2}} & |x|\leq\varepsilon\\ &\\ 0 & |x|\geq\varepsilon \end{cases} \qquad C_\varepsilon=\frac{1}{\varepsilon \int\limits_{|\xi|<1}e^{-\frac{1}{1-\xi^2}}\mathrm{d}\xi}\qquad\forall\varepsilon>0 \end{equation}
is the standard $\delta$-sequence converging to the Dirac distribution (I made explicit the domain of integration in this formula respect to the one found in reference [1] in order to show clearly what is the rôle of $Z(f)$).
If $f(x)$ does not satisfies the requirement 1 and 2, the the limit \eqref{2} is not equal to formula \eqref{1} and may or may not exists, i.e $\delta(f(x))$ may or may not be a distribution.
Finally, let's consider the case where $Z(f)$ is not made of only isolated points: if, for example $f(x)=0$ for all $x$ in a given interval $[a,b]$, the limit \eqref{2} is $\infty$ since $\omega_\varepsilon(f(x))$ goes to $\infty$ for $\varepsilon\to0$ on a set of finite measure $b-a$, thus again does not defines a distribution.
[1] Vladimirov, V. S. (2002), Methods of the theory of generalized functions, Analytical Methods and Special Functions, 6, London–New York: Taylor & Francis, pp. XII+353, ISBN 0-415-27356-0, MR 2012831, Zbl 1078.46029.