symmetric double-integral on fractional part Let $\{\}$ denotes the fractional part function, does the following double-integral have a closed-form ?
$$\int_{0}^{1}\int_{0}^{1}\bigg\{\frac{1}{x}+\frac{1}{y}\bigg\}dx\,dy$$
 A: It's not a full answer:
Using identity:
$$\left \{ z \right \}=\frac{1}{2}-\sum _{k=1}^{\infty } \frac{\sin (2 k \pi  z)}{k \pi }$$
We can write (with CAS help):
$$\color{red}{\int_{0}^{1}\int_{0}^{1}\bigg\{\frac{1}{x}+\frac{1}{y}\bigg\}dx\,dy}=\\\int _0^1\int _0^1\left(\frac{1}{2}-\sum _{k=1}^{\infty } \frac{\sin \left(2 k \pi  \left(\frac{1}{x}+\frac{1}{y}\right)\right)}{k \pi
   }\right)dydx=\\\frac{1}{2}-\sum _{k=1}^{\infty } \left(4 k \pi ^2 \text{Ci}(2 k \pi )+4 \cos (2 k \pi ) \text{Ci}(2 k \pi )+2 \pi 
   \sin (2 k \pi )-\frac{\sin (4 k \pi )}{k \pi }+8 k \pi  \text{Ci}(2 k \pi ) \text{Si}(2 k \pi )-4 \sin (2 k \pi ) \text{Si}(2 k \pi
   )\right)=\color{red}{\\\frac{3}{2}-2 \gamma +4 \pi  \sum _{k=1}^{\infty } (2 k \text{Ci}(2 k \pi ) \text{Si}(2 k \pi )-k \pi  \text{Ci}(2 k \pi ))}\approx0.495921$$
where: $\text{Ci}$ and $\text{Si}$ is  cosine(sine) integral function.
A: This is not a solution as it does not provide a closed form but we derive a more compact form of the solution, in which the double integral over a fractional part is reduced to a single integral over a smooth integrand.
This is another illustration of the method described in my solution to Evaluation of $\int_{0}^{1}\int_{0}^{1}\{\frac{1}{\,x}\}\{\frac{1}{x\,y}\}dx\,dy\,$
Result
We show here that
$$i:=\int _0^1\int _0^1\{\frac{1}{x}+\frac{1}{y}\}\,dydx= 1-\gamma +\int_0^1 \psi ^{(1)}(\xi +1) \psi ^{(0)}(2-\xi ) \, d\xi\tag{1}$$
where $\{x\}$ is the fracional part of $x$, $\gamma$ is Euler gamma and $\psi ^{(n)}(x)$ is the polygamma function.
Here the integral over the polygamma functions is a $\simeq 17\%$ correction to $1-\gamma$. Numerically we have, respectively, 
$$N(i)=0.42278433509846713 +0.07313656826103414=0.4959209033595013$$
Derivation
To begin with, we simplify the integral with the aim to reduce the fractional part of two variables to that of a single variable 
Letting $x\to \frac{1}{r}$, $y\to \frac{1}{s}$, followed by $r\to u$, $s \to v-u$ leads to
$$i=\int _1^\infty \int _{1+u}^\infty \frac{1}{u^2 (v-u)^2}\{v\}\,dudv\tag{2}$$
Notice that the Jacobian determinant is $1$ and that, since $s\ge 1$, the v-integral has to start at $1+u$.
Next we split the integrals into separate integrals over integer regions, e.g.
$$\int_1^\infty f(u) \,du = \int_1^2 f(u) \,du +\int_2^3 f(u) \,du +... \\= \int_0^1 f(1+\xi) \,d\xi +\int_0^1 f(2+\xi) \,d\xi +... \\= \sum_{k=1}^\infty \int_0^1 f(k+\xi) \,d\xi = \int_0^1 (\sum_{k=1}^\infty f(k+\xi) )\,d\xi$$
in the last step we have optionally intechanged integration and summation.
In our case we let $u=k+\xi$, $v=m+\eta$ with $\{u\}=\xi$ and $\{v\}=\eta$ to obtain
$$i=\sum_{k=1}^\infty \int_0^1 \frac{\,d\xi}{(k+\xi)^2} \left(\int_{1+k+\xi}^{2+k} \frac{\,dv\{v\}}{(v-k-\xi)^2}\\+\sum_{m=2+k}^\infty \int_0^1 \frac{\,d\eta\; \eta}{(m-k-\xi +\eta)^2}\right)$$
The first term in the bracket gives after letting $v=(1+k+\xi)+\eta$, $\{v\}=\xi+\eta$ in the integral
$$\int_{0}^{1-\xi} \frac{\,d\eta(\xi+\eta)}{(v-k-\xi)^2}= \frac{1}{-2+\xi}+\xi +\log(2-\xi)$$
The second term in the bracket is with $m-k=n$
$$\sum_{n=2}^\infty \int_0^1 \frac{\,d\eta\; \eta}{(n-\xi +\eta)^2}\\=\sum_{n=2}^\infty (-\frac{1}{n-\xi +1}-\log (n-\xi )+\log (n-\xi +1))\\=\psi ^{(0)}(3-\xi )-\log (2-\xi )$$
Hence the bracket becomes
$$(\frac{1}{-2+\xi}+\xi +\log(2-\xi))+(\psi ^{(0)}(3-\xi )-\log (2-\xi ))\\=\frac{1}{-2+\xi}+\xi +\psi ^{(0)}(3-\xi ) $$
Notice that the log-term has dropped out.
Now we observe that the bracket simplifies further due to
$$\psi ^{(0)}(3-\xi )= H_{2-\xi}-\gamma$$
and
$$H_{2-\xi}-\frac{1}{-2+\xi}= H_{1-\xi}=\psi ^{(0)}(2-\xi )+\gamma$$
Since the bracket does not depend on $k$ the overall sum can be done 
$$\sum_{k=1}^\infty \int_0^1 \frac{\,d\xi}{(k+\xi)^2}=\psi ^{(1)}(\xi +1)$$
and we are left with the integral
$$i=\int_0^1 \,d\xi \psi ^{(1)}(\xi +1)(\xi  +\psi ^{(0)}(2-\xi )) $$
Finally, the relation
$$\int_0^1 \xi  \psi ^{(1)}(\xi +1) \, d\xi=1-\gamma$$
completes the derivation.
