Computing product distribution If $X$ and $Y$ have the following joint density function
$$f(x,y)=xe^{-x(y+1)}$$ 
For $0<x,0<y$, compute the density function of $Z=XY$.

I am having some problems with this question because $f_X=e^{-x}$ and $f_Y=\frac{1}{(y+1)^2}$, implying that $X$ and $Y$ are not independent. 

Trying to integrate $$\iint_{xy\le a} xe^{-x(y+1)} \mathrm{d}(x,y)$$
I ran into some difficult integrands that made me doubt whether this was the right approach.

I thought also about finding the distribution of the sum of $\ln X$ and $\ln Y$ but this also gave some difficult integrands.

I am unsure as to how to proceed with this question. What approach could I take for this problem?
 A: While your answer is good, note that it is more straightforward to switch the order of integration.  For $a>0$ we have
$$P(XY \leq a)= \int_A x e^{-x(y+1)}dydx$$ where $A=\{(x,y):x,y>0,xy \leq a\}$.  Writing this as an iterated integral:
$$\begin{align}P(XY \leq a) &=\int_0^{\infty}\int_0^{a/x}xe^{-x(y+1)}\mathrm{d}y\mathrm{d}x \\
  & =\int_0^{\infty} xe^{-x} \int_0^{a/x} e^{-xy} \mathrm{d}y \mathrm{d}x \\
 &=\int_0^\infty xe^{-x} \big[\frac{-1}{x}e^{-yx}\big]_{y=0}^{a/x} \mathrm{d}x \\
 & =\int_0^\infty e^{-x}(1-e^{-a}) \mathrm{d}x \\
 & =(1-e^{-a}) \big[-e^{-x} \big]_{x=0}^\infty \\
 & =1-e^{-a} \end{align}$$
Since if $a \leq 0$ we have $P(XY \leq a)=0$, you may recognize this as the CDF of the exponential distribution with parameter 1. If not, simply differentiate it to get the density of XY:
$$f_{XY}(a)=\frac{d}{da} P(XY \leq a)=\frac{d}{da} \big[ 1-e^{-a}\big] =e^{-a}$$
which holds for $a>0$.  For $a \leq 0$, $f_{XY}(a)=0$.
A: It seems I made a mistake in my calculation of the integral. I am posting this answer so that if anyone else runs into this question in Sheldon Ross' Book, First Course in Probability, they will not have the same problems: 
simply calculate $$\begin{align}\iint_{xy \le a}xe^{-x(y+1)}\mathrm{d}x\mathrm{d}y &=\int_0^{\infty}\int_0^{a/y}xe^{-x(y+1)}\mathrm{d}x\mathrm{d}y \\ I & =\int_0^{\infty} \frac{1-\exp\left[-\frac{a}{y}(y+1)\right]\left(\frac ay(y+1)+1\right)}{(y+1)^2} \mathrm{d}y\end{align}$$
We use a modified version of integration by parts, namely: $$\int \frac{u\mathrm{d}v}{v^2}=\int \frac{\mathrm{d}u}{v}-\frac uv $$
And set $v=y+1$ and $u=1-\exp\left[-\frac{a}{y}(y+1)\right]\left(\frac ay(y+1)+1\right)$ which implies that $\mathrm{d}v=\mathrm{d}y$ and $\mathrm{d}u =-e^{-a}e^{-\frac ay} \frac {a^2}{y^3} (y+1) \mathrm{d}y $
Substituting into the modified integration by parts we get
$$I=-e^{-a}\int e^{-\frac ay} \frac {a^2}{y^3} \mathrm{d}y-\frac {1-\exp\left[-\frac{a}{y}(y+1)\right]\left(\frac ay(y+1)+1\right)}{y+1}$$
After setting $w= \frac ay$, the first integral in this new expression becomes $$e^{-a} \int w e^{-w} \mathrm{d}w$$ 
Integrating by parts and substituting into $I$ we obtain
$$I=-\exp \left[-\frac ay (y+1) \right]\left( \frac ay +1 \right)-\frac {1-\exp\left[-\frac{a}{y}(y+1)\right]\left(\frac ay(y+1)+1\right)}{y+1}$$
Simplifying yields:
$$I=-\frac{1+y\exp\left[-\frac{a}{y}(y+1)\right]}{y+1}$$
Taking limits we get that $$\lim_{y \to 0^+} I = -1$$
$$\lim_{y \to \infty} I = -e^{-a}$$
Therefore, we obtain that $\mathbb{P}(XY \le a)=1-e^{-a}=F_Z(a)$ and after taking derivatives, we get that $f_Z(a)=e^{-a}$.
