Double integral with strange Change of Variables I am trying to compute the following integral,
$$\iint_{\mathbb{R}^2} \left(\frac{1-e^{-xy}}{xy}\right)^2 e^{-x^2-y^2}dxdy$$
First I tried substituting $x=r\cos{\theta}, y=r\sin{\theta}$ but it didn't really give me anything.
For the second try, I tried $u=x^2+y^2, v=xy$ but after computing the Jacobian, it got really messy, I couldn't really continue further.
Would there be a way to change the variables so that I would be able to compute this integral? 

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
&\bbox[10px,#ffd]{\ds{\iint_{\mathbb{R}^2}
\pars{1 - \expo{-xy} \over xy}^{2}\expo{-x^{2} - y^{2}}\dd x\,\dd y}}
\\[5mm] = &\
\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\
\overbrace{\pars{\int_{0}^{1}\expo{-xya}\dd a}}
^{\ds{\expo{-xy} - 1 \over -xy}}\
\overbrace{\pars{\int_{0}^{1}\expo{-xyb}\dd b}}
^{\ds{\expo{-xy} - 1 \over -xy}}
\expo{-x^{2} - y^{2}}\dd x\,\dd y
\\[5mm] = &\
\int_{0}^{1}\int_{0}^{1}\
\overbrace{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
\exp\pars{-x^{2} - \bracks{a + b}xy - y^{2}}\dd x\,\dd y}
^{\ds{2\pi \over \root{4 - \pars{a + b}^{2}}}}\
\,\dd a\,\dd b
\\[5mm] = &\
2\pi\int_{0}^{1}\int_{0}^{1}{\dd a\,\dd b \over \root{4 - \pars{a + b}^{2}}} =
\bbx{{4 \over 3}\,\pi\pars{3 - 3\root{3} + \pi}} \approx 3.9603
\end{align}
