I am having some issues understanding Dirac's delta function/distribution behaviour under change of coordinates. There is a statement, if $(x_1,\ldots,x_n)$ are cartesian coordinates and $y_1,\ldots,y_n$ are some general coordinates in Euclidean $\mathbb{E}^n$ space, then: $$ x_1=x_1(y_1,\ldots,y_n),\\ \vdots \\ x_n=x_n(y_1,\ldots,y_n)$$ and the Jacobian: $$ J(y_1,\ldots,y_n) = \begin{vmatrix} \frac{\partial x_1}{\partial y_1} & \cdots & \frac{\partial x_1}{\partial y_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial x_n}{\partial y_1} & \cdots & \frac{\partial x_n}{\partial y_n} \end{vmatrix} $$
define the transformation. Provided that $J(y_1^{(P)},\ldots,y_n^{(P)})\ne0$ for $(y_1^{(P)},\ldots,y_n^{(P)})$ coordinates of some point $P$ coresponding to cartesian coordinates $(x_1^{(P)},\ldots,x_n^{(P)})$, then:
$$\delta(x_1-x_1^{(P)})\cdots\delta(x_n-x_n^{(P)}) = \frac{1}{|J(y_1^{(P)},\ldots,y_n^{(P)})|} \delta(y_1-x_1^{(P)})\cdots\delta(x_n-y_n^{(P)})$$.
I wonder how to show this fact? If I consider $\delta$ in the original idea by Dirac, then:
$$f(x_1^{(P)},\ldots,x_n^{(P)})=\int_{\Omega}\delta(x_1-x_1^{(P)})\cdots\delta(x_n-x_n^{(P)}) f(x_1,\ldots,x_n)\,{\rm d}x_1\ldots{\rm d}x_n=\\= \int\ldots\int\frac{\delta(y_1-y_1^{(P)})\cdots\delta(y_n-y_n^{(P)})}{|J(y_1^{(P)},\ldots,y_n^{(P)})|} f(y_1,\ldots,y_n) |J(y_1^{(P)},\ldots,y_n^{(P)})|\,{\rm d}y_1\ldots{\rm d}y_n=\\ = \begin{cases} f(y_1^{(P)},\ldots,y_n^{(P)}) & \quad \text{all integration bounds contain the point $P$}\\ 0 & \quad \text{Any one of the integration bounds does not contain the point $P$} \end{cases} $$
And almost all seems fine but I am not sure if such proof is sufficient for this interpretation.
When I consider $\delta$ as a distribution, I do not really understand how to show that this holds not even whether it is true. I remember that for a distribution, one can define a multiplication with a $C^\infty$ function, but I still cannot see how to show that the above holds, since in terms of distributions $\delta$ is defined only as a continuous linear functional that does $<\delta,\phi>=\phi(x)$ to any $C^\infty$ test function with compact support in some set (there is really no integration involved).