What is $P(X < Y)$ when joint pdf is $f(x,y) = 2ye^{-xy-2y}\mathbb{1}_{[0, \infty) \times [0, \infty)}(x,y)$? 
Joint probability density is given as $$f(x,y) = 2ye^{-xy-2y}\mathbb{1}_{[0, \infty) \times [0, \infty)}(x,y)$$ What is $\mathbb{P}(X < Y)$? Are $X$ and $Y$ independent? What is $\operatorname{Cov}(X, Y)$?

Now, I drew the area and it seems that I have to calculate:
$$
\int_0^{\infty} dx \int_x^{\infty} f(x,y) \,dy
$$
or equivalently
$$
\int_0^{\infty} dy \int_0^{y}  f(x,y)\, dx
$$
Second integral seems easier at glance but:
$$
\int_0^{\infty} dy \int_0^{y}  2ye^{-xy-2y} \,dx = \int_0^{\infty}-2e^{-y(x+2)} \Bigm|_0^y dy =  \int_0^{\infty} e^{-y^2-2y}(2{e^{y^2}-2)\,dy
}$$
which seems atrocious. Am I missing something or calculating the integrals wrong?
To prove $X$ and $Y$ (in)dependence I'd calculate $f_X(x)$ by integrating over $y$, and $f_Y(y)$ analogously, then multiply both to see if I get $f(x,y)$ back. But I am stuck on those integrals! Any help appreciated.
If they were independent, then $\operatorname{Cov}(X,Y)$ would be $0$. If not... well, I'd be stuck again, probably.
 A: To see whether $X$ and $Y$ are independent, you can write the density of $(X,Y)$ as
$$f_{X,Y}(x,y) = 2ye^{-xy-2y}\mathbb{1}_{[0, \infty) \times [0, \infty)}(x,y)=\underbrace{ye^{-yx}1_{[0,\infty)}(x)}_{f_{X\mid Y}(x\mid y)}\cdot\underbrace{2e^{-2y}1_{[0,\infty)}(y)}_{f_Y(y)}$$
In other words, $X$ conditioned on $Y=y$ has an exponential distribution with mean $1/y$ (this indicates dependence between $X$ and $Y$), where $Y$ itself is distributed as exponential with mean $1/2$.
Now you can use law of total expectation to find $\operatorname{Cov}(X,Y)=\operatorname E[XY]-\operatorname E[X]\operatorname E[Y]$ (provided the expectations exist of course). Here $\operatorname E[X]=\operatorname E\left[\operatorname E[X\mid Y]\right]=\operatorname E\left[\frac1Y\right]=2\int_0^\infty \frac{e^{-2y}}{y}\,dy=\infty$, so the covariance does not exist by definition.
A: The integral you have there is doable if you move things around a bit.
\begin{align}
\int_0^{\infty} 2e^{-y^2-2y}(e^{y^2}-1)\textrm{d}y &=2\left(\int_0^{\infty} e^{-2y}\textrm{d}y-\int_0^{\infty} e^{-(y^2+2y)}\textrm{d}y\right)\\
&=2\left(\frac{1}{2}-\int_0^{\infty} e^{-(y+1)^2+1}\textrm{d}y\right)\\
&=1-2e\sqrt{\pi}\frac{1}{\sqrt{\pi}}\int_{1}^{\infty} e^{-u^2}\textrm{d}u \\
&=1-2e\sqrt{\pi}\frac{1}{\sqrt{2 \pi}}\int_{\sqrt{2}}^{\infty} e^{-\frac{t^2}{2}}\textrm{d}t \\
&=1-2e\sqrt{\pi}(1-\Phi(\sqrt{2})) \\
&=1-2e\sqrt{\pi}\Phi(-\sqrt{2}),
\end{align}
where $\Phi$ is the distribution function of the standard Gaussian distribution and we made the changes of variables $u=y+1$ and $t=\sqrt{2}u$.
Edit: For those interested, this probability seems to be roughly $0.24$
