CDF of $\frac{X_1X_2}{X_1+X_2+a}$, where $X_1$ and $X_2$ have exp. distributions Let $X_1$ and $X_2$ be two independent exponential random variables with the PDFs $f_{X_1}(x_1)=\lambda_1 \exp(-\lambda_1 x_1)$ and $f_{X_2}(x_2)=\lambda_2 \exp(-\lambda_2 x_2)$. 
Let $X= \frac{X_1X_2}{X_1+X_2+a}$ be a r.v. for which I want to derive the CDF, where $a$ is a positive constant.
In other words, I need to calculate $\mathbb{P}\{X <x \}$.
Any ideas or hints?
Edit (attempt): 
Using a similar approach to that used here, we have $X_1 \in [0, \infty)$ and  $X_2 \in [0, \frac{x(x_1+a)}{x_1-x})$. The latter interval results from: $\frac{x_1x_2}{x_1+x_2+a} < x$, thus $x_2(x_1-x)<x(x_1+a)$ which implies $x_2 < \frac{x(x_1+a)}{x_1-x}$; but I am not sure if this is correct since $x_1-x$ can sometimes be negative (?). Based on the above, we can write: \begin{align}
 \mathbb{P} \{ \frac{X_1X_2}{X_1+X_2+a} < x \}  & = \int_{x_1=0}^\infty \int_{x_2=0}^{ \frac{x(x_1+a)}{x_1-x}} \lambda_1 \exp(-\lambda_1 x_1) \lambda_2 \exp(-\lambda_2 x_2)  dx_2 dx_1  \\ &= \int_{x_1=0}^\infty \left(1-\exp(-\lambda_2\frac{x(x_1+a)}{x_1-x}) \right) \lambda_1 \exp(-\lambda_1 x_1) dx_1  \\ & = 1-\lambda_1 \int_{x_1=x}^\infty \exp(-\lambda_2\frac{x(x_1+a)}{x_1-x} ) \exp(-\lambda_1 x_1)   dx_1
\end{align}
In the last equality: for the integral of $\lambda_1 \exp(-\lambda_1 x_1) dx_1$, I suppose that $x_1 \in [0,\infty)$, whereas for the second term I suppose that $x_1 \ge x$ (since $x_2$ should be positive) (?).
Is my approach correct ?
 A: Thanks to the indications that @Did gave, I was able to derive the CDF of $Y=\frac{X_1X_2}{X_1+X_2+a}$ (note that in the question I use $X$ instead of $Y$) as follows:  
Based on the identity of events 
\begin{align}
\left[ Y <y \right] = \left[ X_1 < y \right] \cup \left[ X_1 \ge y, X_2 < y (X_1+a)(X_1-y)^{-1} \right],
\end{align}
we get the following
\begin{align}
\nonumber \mathbb{P} \{ Y < y \} &= \mathbb{P}\{X_1 <y\}+ \int_y^\infty \mathbb{P}\{X_2 < y (x+a) (x-y)^{-1} \} f_{X_1}(x) \, dx \\ \nonumber &= 1-e^{-\lambda_1y} + \lambda_1 \int_y^\infty (1-e^{-\lambda_2 y(x+a)(x-y)^{-1}}) e^{-\lambda_1 x} dx \\ \nonumber &= 1-e^{-\lambda_1y}+e^{-\lambda_1y}- \lambda_1 \int_y^\infty  e^{-\lambda_2 y(x+a)(x-y)^{-1}} e^{-\lambda_1 x} dx \\ \nonumber & =_{(i)}  1- \lambda_1 \int_{u=0}^\infty e^{ -\lambda_2 y(u+y+a)u^{-1}} e^{-\lambda_1 (u+y)} du \\ \nonumber  &= 1- \lambda_1 e^{-(\lambda_1+\lambda_2)y} \int_{u=0}^\infty e^{-\lambda_2 y (y+a)u^{-1}} e^{-\lambda_1u} du \\ & =_{(ii)} 1- \lambda_1 e^{-(\lambda_1+\lambda_2)y} \, 2 \, \sqrt{ y(y+a) \lambda_2 \lambda_1^{-1} } \, K_1\left(2 \sqrt{  y(y+a) \lambda_2 \lambda_1  } \right),
 \\ & = 1- e^{-(\lambda_1+\lambda_2)y} \, 2 \, \sqrt{ y(y+a) \lambda_2 \lambda_1 } \, K_1\left(2 \sqrt{  y(y+a) \lambda_2 \lambda_1  } \right)
\end{align}
in which equality (i) is due to the change of variable $u=x-y$ and equality (ii) follows from [Table of Integrals, Series and Products, 7th edition - equation 3.471.9].
A: Let $X$ and $Y$ denote two independent  exponential random variables and suppose that $V = \frac{XY}{X+Y+a}$ where $a > 0$. What is the CDF of $V$?
First, note that $V > 0$. Let $v$ denote a positive constant, and let us try to determine the complementary CDF $P\{V > v\}$ by integrating the joint density of $X$ and $Y$ over that part of the first quadrantVwhere $V$ exceeds $v$.  The set in question, call it $A$, is given by
\begin{align}
A  &= \{(x,y)\colon x > 0, y > 0, \frac{xy}{x+y+a} > v\}\\
&= \{(x,y)\colon x > 0, y > 0, xy > v(x+y+a)\}\\
&= \{(x,y)\colon x > 0, y > 0, xy -vx -vy > av\}\\
&= \{(x,y)\colon x > 0, y > 0, (x-v)(y-v) > v^2 + av\}.
\end{align}
Now, the graph of the hyperbola $xy = b$ consists of two curves confined to the first and third quadrants respectively and passing through the points $\left(\sqrt{b}, \sqrt{b}\right)$ and $\left(-\sqrt{b}, -\sqrt{b}\right)$ respectively. Therefore,  the graph of $(x-v)(y-v) = b$ is just these two curves shifted to the right by $v$ and shifted upwards by $v$, and the two curves now pass through $(\sqrt{b}+v, \sqrt{b}+v)$ and $(-\sqrt{b}+v, -\sqrt{b}+v)$ respectively.  Note that the asymptotes of the curves are $x=v, y=v$. Now, when $b$ equals $v^2+av$, $\sqrt{b} > v$ and so the point $(-\sqrt{b}+v, -\sqrt{b}+v)$ is in the third quadrant. Consequently, the lower branch of the hyperbola does not lie in the first quadrant at all.  (It does cross the $x$ and $y$ axes into the second and fourth quadrants but that is immaterial in this problem). It follows that we can express $A$ as
$$A = \{(x,y)\colon x > v, y > v, (x-v)(y-v) > v^2 + av\}.$$
Hence,
\begin{align}
1-F_V(v) &= P((X,Y)\in A)\\
&= \iint_A f_{X,Y}(x,y)\, \mathrm dy \, \mathrm dx\\
&= \int_v^\infty f_X(x) \left[ \int_{y = \frac{v^2+av}{x-v}}^\infty f_Y(y)\,\mathrm dy \right] \, \mathrm dx.
\end{align}
The inner integral is straightforward to evaluate; the outer one is trickier, needing special functions and tables of integrals to
evaluate.
An alternative calculation given in an answer by the OP (with the
help of many suggestions from @Did) directly evaluates
the CDF $P\{V \leq v\}$ by partitioning the set under consideration
into the events $\{X \leq v\}$ and $\left\{X>v, 0 < Y \leq  v\frac{x+a}{x-v}\right\}$.
