Let's look at a special case and use natural logs to convert this to a sum. Suppose we have $X_1,X_2 \sim U[0,1]$ with $X_1 X_2 = z.$ Let $R= - \mathrm{ln} \left( X_3 \right)$ and $S= - \mathrm{ln} \left( X_4 \right).$

Then we require $R + S = - \mathrm{ln} \left( X_3 \right) + - \mathrm{ln} \left( X_4 \right) = - \mathrm{ln} \left( z \right).$

It can be determined that $R$ and $S$ are gamma random variables with both parameters equal to 1 (exponential in this case).

So now we have a fixed sum of iid gamma random variables. Applying the general logic outlined in the answer here to the gamma case (How to generate two groups of $n$ random numbers in $U(0,1)$ such that sum of these two groups equal?), we find the conditional distribution of $R$ given $z$ is $$R|z \ \sim \ U \left[ 0,-\mathrm{ln} \left( z \right ) \right]$$

Then the algorithm is:

1. Generate independent realizations $x_1$ and $x_2$ as $U[0,1].$
2. Calculate $z=x_1x_2.$
3. Generate a realization $r$ from $R \sim U[0,- \mathrm{ln} \left( z \right) ]$
4. Calculate $x_3 = e^{-r}$
5. Calculate $x_4 =z/x_3$

Let's look at a special case and use natural logs to convert this to a sum. Suppose we have $X_1,X_2 \sim U[0,1]$ with $X_1 X_2 = z.$ Let $R= - \mathrm{ln} \left( X_3 \right)$ and $S= - \mathrm{ln} \left( X_4 \right).$

Then we require $R + S = - \mathrm{ln} \left( X_3 \right) + - \mathrm{ln} \left( X_4 \right) = - \mathrm{ln} \left( z \right).$

It can be determined that $R$ and $S$ are gamma random variables with both parameters equal to 1 (exponential in this case).

So now we have a fixed sum of iid gamma random variables. Applying the general logic outlined in the answer here to the gamma case (How to generate two groups of $n$ random numbers in $U(0,1)$ such that sum of these two groups equal?), we find the conditional distribution of $R$ given $z$ is $$R|z \ \sim \ U \left[ 0,-\mathrm{ln} \left( z \right ) \right]$$

Then the algorithm is:

1. Generate independent realizations $x_1$ and $x_2$ as $U[0,1].$
2. Calculate $z=x_1x_2.$
3. Generate a realization $r$ from $R \sim U[0,- \mathrm{ln} \left( z \right) ]$
4. Calculate $x_3 = e^{-r}$
5. Calculate $x_4 =z/x_3$

On taking a second look, we can write $-r = \left[ {\mathrm{ln}} \left( z \right) \right] U, $ where $U$ is uniform on $[0,1]. $ Then $e^{-r} = z^U.$

Combining some of the above steps, we can write the process as

1. Generate $x_1,x_2$ as independent $U[0,1]$ realizations and find $z=x_1x_2.$
2. Generate $U \sim U[0,1]$ and calculate $x_3 = z^U$ and $x_4=z^{\left( 1 - U \right)}=z/x_3$