Evaluation of $\int_{0}^{1}\int_{0}^{1}\{\frac{1}{\,x}\}\{\frac{1}{x\,y}\}dx\,dy\,$ Let $\{\}$ denote the fractional part function, does the following double integral have a closed-form ? 
$$\int_{0}^{1}\int_{0}^{1}\bigg\{\frac{1}{\,x}\bigg\}\bigg\{\frac{1}{x\,y}\bigg\}dx\,dy\,$$
 A: So we are looking for 
$$ \iint_{(1,+\infty)^2}\frac{\{x\}\{xy\}}{x^2 y^2}\,dx\,dy=\int_{1}^{+\infty}\frac{\{x\}}{x}\int_{x}^{+\infty}\frac{\{t\}}{t^2}\,dt\,dx. $$
We may recall that $\int_{1}^{+\infty}\frac{\{x\}}{x^2}\,dx =\sum_{n\geq 0}\int_{0}^{1}\frac{x\,dx}{(x+n+1)^2}=\sum_{n\geq 0}\left[\log(n+2)-\log(n+1)-\frac{1}{n+2}\right]$ equals $1-\gamma$ and
$$\{t\}=\frac{1}{2}-\sum_{n\geq 1}\frac{\sin(2\pi n t)}{\pi n},\qquad \frac{\{x\}}{x}=\frac{1}{2x}-\sum_{m\geq 1}\frac{\sin(2\pi m x)}{\pi m x} $$
hold almost everywhere, such that
$$ \int_{x}^{+\infty}\frac{\{t\}}{t^2}\,dt = \frac{1}{2x}-\sum_{n\geq 1}\frac{\sin(2\pi n x)-2\pi n x\,\text{Ci}(2\pi n x)}{\pi n x} $$
and our integral equals
$$ \frac{1}{4}-\int_{1}^{+\infty}\frac{1}{2x}\sum_{m\geq 1}\frac{\sin(2\pi m x)}{\pi m x}\,dx-\int_{1}^{+\infty}\frac{1}{2x}\sum_{n\geq 1}\frac{\sin(2\pi n x)-2\pi n x\,\text{Ci}(2\pi n x)}{\pi n x}\,dx $$
plus
$$ \int_{1}^{+\infty}\sum_{m\geq 1}\sum_{n\geq 1}\frac{\sin(2\pi m x)}{\pi m x}\cdot\frac{\sin(2\pi n x)-2\pi n x\,\text{Ci}(2\pi n x)}{\pi n x}\,dx. $$
After some simplification we get
$$ \frac{1}{4}+\overbrace{2\sum_{m\geq 1}\text{Ci}(2m\pi)}^{\frac{1}{2}-\gamma}+\sum_{n\geq 1}\int_{1}^{+\infty}\text{Ci}(2\pi n x)\frac{dx}{x}+\sum_{m,n\geq 1}\frac{\pi\min(m,n)+(m-n)\text{Si}(2\pi(m-n))-(m+n)\text{Si}(2\pi(m+n))}{\pi mn}-2\sum_{m,n\geq 1}\int_{1}^{+\infty}\frac{\sin(2\pi m x)}{\pi m x}\text{Ci}(2\pi n x)\,dx $$
and probably these pieces can be further simplified by recalling that the Laplace transform of $\text{Ci}$ (the cosine integral) is essentially a logarithm and the Laplace transform of $\text{sinc}$ is essentially an arctangent. Anyway, just in its current state the previous representation allows to find arbitrarily accurate approximations of the wanted integral, since the behaviour of $\text{Ci}$ and $\text{sinc}$ over $(1,+\infty)$ is extremely regular.
A simple and non-trivial upper bound can be derived from the Cauchy-Schwarz inequality. We have
$$ \int_{0}^{1}\left\{\frac{1}{x}\right\}^2\,dx = \int_{1}^{+\infty}\frac{\{x\}^2}{x^2}\,dx = -1-\gamma+\log(2\pi),$$
$$ \int_{0}^{1}\int_{0}^{1}\left\{\frac{1}{xy}\right\}^2\,dx\,dy = \int_{1}^{+\infty}\int_{1}^{+\infty}\frac{\{xy\}}{x^2 y^2}\,dx \,dy=\\=1-\gamma+2\sum_{n\geq 1}\int_{1}^{+\infty}\frac{\sin(2\pi n x)}{2\pi n}\cdot\frac{1-\log x}{x^2}\,dx=\\
=1-\gamma+\int_{1}^{+\infty}(1-2\{x\})\frac{1-\log x}{x^2}\,dx\leq \frac{1}{2}$$
hence
$$ \int_{0}^{1}\int_{0}^{1}\left\{\frac{1}{x}\right\}\left\{\frac{1}{xy}\right\}\,dx\,dy \leq \sqrt{\frac{\log(2\pi)-\gamma-1}{2}}<\frac{13}{36}. $$
A: Upper bound on the above double integral 
As we have $\forall (x,y)\in (0;1)^2,   0\leq \{1/{x\,y}\} <1$  then it follows the accompanied inequality :
$$\int_{0}^{1}\int_{0}^{1}\bigg\{\frac{1}{x}\bigg\}\bigg\{\frac{1}{x\,y}\bigg\}dx\,dy\,< \int_{0}^{1}\int_{0}^{1}\bigg\{\frac{1}{x}\bigg\}dx\,dy\ =1-\gamma$$
where $\gamma$ represents the Euler-Mascheroni constant.  
A: One possible approach is to write
$$
I
:= \int_{0}^{1}\int_{0}^{1}\left\{\frac{1}{x}\right\}\left\{\frac{1}{xy}\right\}\,dxdy
= \left(\frac{1}{2} - \gamma\right)\log(2\pi) - \gamma_1 - 2 + C_1 - \frac{C_2}{2},
$$
where $\gamma_1$ is the Stieltjes constant and
\begin{align*}
C_1 &:= \lim_{N\to\infty} \bigg[ \sum_{n=1}^{N} \frac{\log(n!)}{n} - \left( N\log N - 2N + \frac{1}{4}\log^2 N + \frac{1+\log(2\pi)}{2}\log N \right) \bigg], \\
C_2 &:= \lim_{N\to\infty} \bigg[ \sum_{n=1}^{N} \log^2 n - \left( N \log^2 N - 2N \log N + 2N + \frac{1+\gamma}{2}\log^2 N \right) \bigg].
\end{align*}
I would be surprised if these constants are expressed in closed forms, much like Stieltjes constants are not known to be so.

Derivation of the above formula is not hard. Indeed, we can begin from
$$ I = \int_{1}^{\infty} \frac{\{x\}}{x} \left( \int_{x}^{\infty} \frac{\{y\}}{y^2} \, dy \right) \, dx. $$
Plugging the identity $\int_{x}^{\infty} \frac{\{y\}}{y^2} \, dy = H_{\lfloor x\rfloor} - \log x - \gamma + \frac{\{x\}}{x}$ and evaluating the integral term-by-term,
\begin{align*}
\int_{1}^{N+1} \frac{\{x\}}{x} H_{\lfloor x \rfloor} \, dx
&= \sum_{n=1}^{N} H_n \left( 1 - n ( \log(n+1) - \log n) \right), \\
-\int_{1}^{N+1} \frac{\{x\}}{x} \log x \, dx
&= N - (N+1)\log(N+1) + \frac{1}{2}N\log^2(N+1) - \frac{1}{2} \sum_{n=1}^{N} \log^2 n, \\
-\gamma \int_{1}^{N+1} \frac{\{x\}}{x} \, dx
&= \gamma \left(-N + N\log(N+1) - \log (N!) \right), \\
\int_{1}^{\infty} \left( \frac{\{x\}}{x} \right)^2 \, dx
&= -1 - \gamma + \log(2\pi)
\end{align*}
and a bit of algebra together with Stirling's formula for $\log (N!)$ gives the desired identity above.
