Is there a way to generate this class of bivariate distributions from two uniform random variables? I need to generate a pair of random variables $U,V$ with distribution
$$
f_{U,V}(u,v)\propto(u+v)^{-4}
$$
in $u,v\in[1,L]$ with $L>1$. I managed to invert the CDF of the marginal distribution for one of the RVs but this already requires finding the root of
a fourth order polynomial.
I'm not familiar at all with the theory of generating random variables, so maybe this is the best I can hope for, but to the untrained eye this looks so "simple" that I would have expected a simple solution method to get $U$ and $V$ from a bunch of uniform random variables. I thought I'd post it here in the hope that someone with a good eye for these things might recognize a simple approach.
 A: I will not go too deeply into details and refer to Luc Devroye's book Non-uniform random variate generation for more ideas. 
Sorry for changing your notation: $U$ is always a uniform variable for me. So let $U_1, U_2, \dots$ be uniform $[0,1]$ random variables. We want to generate $X_1$ and $X_2$ with density $$f_X(x_1,x_2) \propto (x_1+x_2)^{-4}\mathbf{1}_{[1,L]}(x_1)\mathbf{1}_{[1,L]}(x_2).$$ It is more convenient to consider $V_i = X_i-1,i=1,2,$ so that $$f_V(x_1,x_2)\propto (2+x_1+x_2)^{-4}\mathbf{1}_{[0,L-1]}(x_1)\mathbf{1}_{[0,L-1]}(x_2).$$
The main idea is that $V_1$ has a uniform distribution given $V_1+V_2$, and the latter has density $\propto x(2+x)^{-4}$ up to some point. 
Generating for large $L$


*

*Generate a random variable $V$ with density $\propto x(x+2)^{-4}$, $x\ge 0$. You can make this using some rejection algorithm, but it is better to observe that $V/(V+2)$ has density $\propto x(1-x)$, which is related to order statistics. Therefore,
i. Generate $U_1$, $U_2$, $U_3$ and take the middle $U_{(2)}$ of them. 
ii. Set $V = 2U_{(2)}/(1-U_{(2)})$.

*If $V>2L-2$, reject and go to 1.

*Generate $U_4$ and set $V_1 = U_4 V$, $V_2 = (1-U_4) V$.

*If $V_1> L - 1$ or $V_2 > L - 1$, reject both and go to 1.

*Add $1$ to both $V_1$ and $V_2$ and output them. 
Why does this work? Just assume for a moment that the rejection step 2 is not there. Then at st step 3 the random variables $V_1$ and $V_2$ have density $\propto (2+x_1+x_2)^{-4}\mathbf{1}_{[0,+\infty)}(x_1)\mathbf{1}_{[0,+\infty)}(x_2)$. 
The rejection is equivalent to conditioning this to $[0,L-1]\times[0,L-1]$, which will get us what we need. 
Note that the step 2 is not necessary at all: if $V>2L-2$, then it is going to be rejected anyway. So this is just for efficiency: we can drop such values before generating $U_4$ and making some operations.
The fact that $L$ is large will ensure that there are not too many iterations. In fact, for $L=4$ the probability of rejection is around $0.02$-$0.05$.
