Generating identically-distributed random variables with a constraint Is there a way to generate identically-distributed random variables (eg $x_1,x_2,x_3,x_4$) with the following constraint:
$\frac{x_1*x_2}{x_3*x_4} ≡ 1$  
$x \in (0,1)$
Please note that simply sampling $x_1, x_2, x_3$ from the given distribution and then finding $x_4$ as $\frac{x_1*x_2}{x_3}$ would not work since $x_4$, in that case, would not follow the given distribution.
Edit1: One inefficient way of doing it may be generating a large number of $x_1,x_2,x_3,x_4$ and then simply filtering with some tolerance the ones satisfying the constraint. 
The method may be not as inefficient here because of this specific constraint (ratio = 1) since we are going to get lots of instances where the constraint is approximately satisfied. 
 A: In your problem you wish to generate a set of four random variables $X_1, X_2, X_3, X_4$ with identical marginal distributions that have support on the interval $0<x<1$, and subject to the additional constraint $X_1 \cdot X_2 = X_3 \cdot X_4$.  We will refer to the marginal distribution of each of these values as the desired marginal.  You can generate random variables obeying these constraints as follows.

Preliminary: Start by choosing any bivariate symmetric joint density function $p_{A,B}$ that has the desired marginal densities for both variables.  From this joint density, you can define $D = A \cdot B$ and you obtain the corresponding joint density:
$$p_{A, D}(a, d) = \begin{vmatrix} 1 & 0 \\ d/a & a \end{vmatrix} p_{A,B}(a, d/a) = a \cdot p_{A,B}(a, d/a) \quad \quad \text{for all }0<d<a.$$
You then have the corresponding densities:
$$p_{A| D}(a, d) = \frac{p_{A,D}(a, d)}{p_{D}(d)} \quad \quad \quad p_{D}(d) = \int_d^1 a \cdot p_{A,B}(a, d/a) da.$$
Generating your values: Now, once you have got these density functions, you are ready to generate your sample as follows:


*

*Generate a single value $D \sim p_D$;

*Now use this value to generate $X_1, X_3 | D \sim \text{IID }p_{A|D}$;

*Set $X_2 = D / X_1$ and $X_4 = D / X_3$.


This will give you generated random variables $X_1, X_2, X_3, X_4$ with marginal distribution $p_A$ (which is the desired distribution) and subject to the constraint $X_1 \cdot X_2 = X_3 \cdot X_4$.  Note that the method works for any choice of symmetric bivariate distribution $p_{A,B}$ and so the solution to the problem is non-unique.  The particular choice of $p_{A,B}$ determines the dependency structure.

Special case: A simple special case is to choose a starting bivariate joint density where the values are independent with the specified desired marginal distribution.  This gives you the simplified form $p_{A,B}(a,b) = p_A(a) \cdot p_A(b)$, and the above equations simplify to:
$$p_{A, D}(a, d) = a \cdot p_{A}(a) p_A(d/a) \quad \quad \text{for all }0<d<a,$$
and:
$$p_{A| D}(a, d) = \frac{a \cdot p_{A}(a) p_A(d/a)}{p_{D}(d)} \quad \quad \quad p_{D}(d) = \int_d^1 a \cdot p_{A}(a) p_A(d/a) da.$$
Given a desired distribution $p_A$, it should not be difficult to derive these distributions and perform the generation algorithm.  (Note that even though you start with a bivariate density with independence, you still obtain generated values that are not independent, since they obey your linear constraint.)
A: 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:


*

*Generate independent realizations $x_1$ and $x_2$ as $U[0,1].$  

*Calculate $z=x_1x_2.$

*Generate a realization $r$ from $R \sim U[0,- \mathrm{ln} \left( z \right) ]$

*Calculate $x_3 = e^{-r}$

*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


*

*Generate $x_1,x_2$ as independent $U[0,1]$ realizations and find $z=x_1x_2.$

*Generate $U \sim U[0,1]$ and calculate $x_3 = z^U$ and $x_4=z^{\left( 1 - U \right)}=z/x_3$

A: This is an extended comment asking for clarification of @MarkL.Stone 's comment.  If I following the procedure correctly, I don't get x1, x2, x3, and x4 having the same marginal distributions if I start out with x1, x2, and x3 having uniform distributions on (0,1):
n <- 10000
x1 <- runif(n)
x2 <- runif(n)
x3 <- rep(NA,n)
for (i in 1:n) {
    x3[i] = runif(1)
    while ((x1[i]*x2[i]) >= x3[i]) {
      x3[i] <- runif(1)
    }
}
x4 <- x1*x2/x3
par(mfrow=c(2,2))
hist(x1, freq=FALSE, main="x1", ylim=c(0,2.5), las=1)
hist(x2, freq=FALSE, main="x2", ylim=c(0,2.5), las=1)
hist(x3, freq=FALSE, main="x3", ylim=c(0,2.5), las=1)
hist(x4, freq=FALSE, main="x4", ylim=c(0,2.5), las=1)


What did I interpret incorrectly?
