What is the pdf of $frac{|x-y|}{(x+y)(2-x-y)}$ when $x,y$ are i.i.d uniform on $[0,1]$?

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 1left(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 + 2frac{-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}$$

Never seen a complex example like this so I do not want to go inside all the calculation required....but perhaps this brainstorming can help.

First, the domain of the rv

$$Z=frac{|X-Y|}{(X+Y)[2-(X+Y)]}$$

is $z in [0;1]$

Now, using the CDF method you have to calculate

$$ bbox[5px,border:2px solid black] { F_Z(z)=mathbb{P}[Z leq z]=mathbb{P}[frac{|X-Y|}{(x+y)[2-(x+y)]} leq z] qquad (1) } $$

As the joint density of $(X,Y) $ is $f_{XY}(x,y)=mathbb{1}_{[0;1]}(x)mathbb{1}_{[0;1]}(y)$, to evaluate the CDF(Z) you can "simply" evaluate the area of the curve expressed in (1) in the unit square $[0;1]times[0;1]$

Starting when $y>x$, that is the upper triangle of the unit square, you have to evaluate the area of the curve

$$zy^2-2z(x-1)y<x(2z+1)-zx^2$$

over the line $x=y$

Same brainstorming for the lower triangle....

I do not know if this hint can help because I didn't go on with the calculations