What is the pdf of $\frac{|x-y|}{(x+y)(2-x-y)}$ when $x,y$ are i.i.d uniform on $[0,1]$? If $x,y$ are i.i.d uniform random variables on $[0,1]$.
I know that the PDF of $|x-y|$ is:
$$f(z) = \begin{cases}
2(1-z) & \text{for $0 < z < 1$} \\
0 & \text{otherwise.}
\end{cases}$$
I know that the PDF of $x+y$ is
$$f(t)=f(z) = \begin{cases}
z & \text{for $0 < z < 1$} \\
2-z & \text{for $1 \le z < 2$} \\
0 & \text{otherwise.}
\end{cases}$$
I am trying to calculate the PDF of $\frac{|x-y|}{(x+y)(2-x-y)}$
What I'm trying to do is to calculate first the PDF of $(x+y)(2-x-y)$ by doing the following: since $(x+y)$ and $(2-x-y)$ have the same PDFs, the PDF of the product should be:
$$f_{(x+y)^2}(z)=\int_{0}^1 f_{x+y}(t) f_{x+y}(z/t) \frac{1}{|t|}\,dt=\int_{0}^1z/t \,dt=z-z\,\log z$$
if $0<z<1$ and
$$f_{(x+y)^2}(z)=\int_{1}^2 f_{x+y}(t) f_{x+y}(z/t) \frac{1}{|t|}\,dt=$$
$$=\int_{1}^2 \frac{(2-t)(2-z/t)}{t} \,dt=-2 - z + z \log 2 + \log 16$$
if $1<z<2$ and $0$ elsewhere.
However, this is not a density function because its area is not 1. I cannot find the mistake.

Let $$Z = \frac{|X-Y|}{(X+Y)(2-X-Y)}.$$  We first determine the support of $Z$ for $X, Y \sim \operatorname{Uniform}(0,1)$.  A quick glance shows that the minimum value is $0$; the maximum is attained at $(X,Y) \in \{ (1,0), (0,1) \}$; therefore, $Z \in [0,1]$.  Then we want $$\Pr[Z \le z] = \int_{x=0}^1 \int_{y=0}^1 \mathbb 1\left(\frac{|x-y|}{(x+y)(2-x-y)} \le z \right) \, dy \, dx.$$   Since the the integrand is symmetric in $x$ and $y$, as well as symmetric about the line $x + y = 1$ (since replacing $(x,y)$ with $(1-y, 1-x)$ gives the same function), we can restrict our attention to the triangular region with vertices $(0, 0), (1/2, 1/2), (1, 0)$.  The transformation $$u = x-y, \quad v = x+y$$ turns this region into $(0,0), (0,1), (1,1)$ in $(u,v)$-space, and the integral becomes $$\Pr[Z \le z] = 4 \int_{v=0}^1 \int_{u=0}^v \mathbb 1 \left( \frac{u}{v(2-v)} \le z \right) \cdot \frac{1}{2} \, du \, dv, \tag{1}$$ where the $4$ comes from the symmetry argument above, and the $1/2$ is the Jacobian of the transformation from $(x,y)$ to $(u,v)$.  We consider the locus of points satisfying the equation $$u = z v(2-v)$$ as comprising part of the boundary of the region of integration, the other parts being the aforementioned triangle.  As this is clearly a parabola with vertex at $(u,v) = (z,1)$, axis of symmetry $u = 1$, and passing through $(0,0)$, we note that when $1/2 < z \le 1$, there is a nontrivial point of intersection with the line $u = v$ at $$u = v = 2 - \frac{1}{z}.$$  Therefore, we must consider these cases separately.  When $0 \le z \le 1/2$, the integral $(1)$ is simply the area of the parabolic region given by $$2 \int_{v=0}^1 \int_{u=0}^{zv(2-v)} \, du \, dv = 2 \int_{v=0}^1 zv(2-v) \, dv = \frac{4z}{3}.$$  When $1/2 < z \le 1$, we instead evaluate $$\begin{align}
\Pr[Z \le z]
&= 2 \left( \int_{v=0}^{2 - 1/z} \int_{u=0}^v \, du \, dv + \int_{v= 2-1/z}^1 \int_{u=0}^{zv(2-v)} \, du \, dv \right) \\
&= (2-1/z)^2 + 2 \int_{v=2 - 1/z}^1 zv(2-v) \, dv \\
&= (2-1/z)^2 + 2\frac{-1 + 3z + 2z^3}{3z^2} \\
&= \frac{1-6z+12z^2-4z^3}{3z^2}.
\end{align}$$
Therefore, the full CDF is given by $$F_Z(z) = \begin{cases} 0, & z < 0 \\
\frac{4z}{3}, & 0 \le z \le 1/2 \\
\frac{1-6z+12z^2-4z^3}{3z^2}, & 1/2 < z \le 1 \\
1, & 1 < z. \end{cases}$$
Consequently, the density is
$$f_Z(z) = \begin{cases} \frac{4}{3}, & 0 < z < 1/2 \\ \frac{2(-1 + 3z - 2z^3)}{3z^3}, & 1/2 < z < 1. \end{cases}$$
