About the limit $\lim_{n \to +\infty} \frac{1}{n^2} \sum_{1 \le a,b \le n} \frac{1}{ \mathrm{gcd} (a,b)} $ This is not homework.
My question is:

Prove or disprove: $$\lim_{n \to +\infty} \frac{1}{n^2} \sum_{a,b=1}^n \frac{1}{ \mathrm{gcd} (a,b)} = \frac{\zeta(3)}{\zeta(2)}$$

This would represent the probability as $n \to +\infty$ that, after having picked randomly $a,b,c \in \{ 1, \dots ,n \}$, the line
$$ax+by+c=0$$
has some integer coordinate point $(x,y \in \Bbb Z^2)$. Indeed, fixed $a,b \in \{ 1, \dots ,n \}$, then $c=-ax-by$ for some $x,y \in \Bbb Z$ is equivalent to $$c \mathrm{\  is  \ a\  multiple\  of\  gcd}(a,b)$$
because of Bezout identity; this event happens with probability $1/\gcd (a,b)$.
What I tried first: clearly
$$0 \le \frac{1}{n^2} \sum_{a,b=1}^n \frac{1}{ \mathrm{gcd} (a,b)} \le \frac{1}{n^2} \sum_{a,b=1}^n 1 = 1$$
so, our limit belongs to $[0,1]$ (if it exists).
Then I made some euristic argument:
I know that the probability, as $n \to +\infty$, that two random numbers $a,b\in \{ 1, \dots ,n \}$ are coprime is  $$\zeta(2)^{-1} = \frac{6}{\pi^2},$$ thus the probability that $\gcd(a,b)=d$ is equal to
$$\frac{\zeta(2)^{-1}}{d^2} = \frac{6}{d^2\pi^2}$$
Thus our probability is (and here is where I have doubts in my argument)
$$\sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2} \to \frac{1}{\zeta(2)} \cdot \zeta(2) =1$$
EDIT: or maybe my last step should be
$$\sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2} \frac{1}{d} \to \frac{\zeta(3)}{\zeta(2)}$$
 A: There is no need to involve probabilities, just a bit of analytic number theory suffices. I don't know how familiar you are with $L$-functions, so I'll leave some details open for now. If any part needs clarification, let me know.
For every positive integer $n$ define 
$$f(n):=\sum_{a,b=1}^n\frac{1}{\gcd(a,b)}
\qquad\text{ and }\qquad
g(n):=\sum_{d\mid n}\frac{\varphi(\tfrac{n}{d})}{d}$$ so that for all $n>1$ we have
\begin{eqnarray*}
f(n)-f(n-1)
&=&-\frac1n+2\sum_{c=1}^n\frac{1}{\gcd(c,n)}\\
&=&-\frac1n+2\sum_{d\mid n}\frac{\varphi(\tfrac{n}{d})}{d}\\
&=&-\frac1n+2g(n).
\end{eqnarray*}
Because $\lim_{n\to\infty}\frac{1}{n^2}\sum_{k=1}^n-\frac1k=0$, this already shows that
$$\lim_{n\to\infty}\frac{1}{n^2}f(n)=2\lim_{n\to\infty}\frac{1}{n^2}\sum_{k=1}^ng(n).$$
Note that $g$ is a multiplicative function, and in fact $g=\iota\ast\varphi$ where $\iota(n):=\tfrac1n$ and $\ast$ denotes the Dirichlet convolution of multiplicative functions. This immediately shows that the $L$-function associated to $g$ converges for all $s\in\Bbb{C}$ with $\operatorname{Re}(s)>2$ and that for all such $s\in\Bbb{C}$ we have
$$L_g(s)=L_{\iota}(s)L_{\varphi}(s)=\zeta(s+1)L_{\varphi}(s).$$
It is easy to see that $L_{\iota}(s)=\zeta(s+1)$ for these $s$, and it is a basic exercise to show that 
$$L_{\varphi}(s)=\frac{\zeta(s-1)}{\zeta(s)},$$
whenever $\operatorname{Re}(s)>2$. Then by a Tauberian theorem, for example Wiener-Ikehara, we have
\begin{eqnarray*}
\lim_{n\to\infty}\frac{1}{n^2}f(n)
&=&2\operatorname{Res}(L_g,2)\\
&=&2\cdot\operatorname{Res}(\zeta(s+1),2)\cdot\operatorname{Res}(L_{\varphi},2)\\
&=&2\cdot\zeta(3)\cdot\frac{1}{2\zeta(2)}\\
&=&\frac{\zeta(3)}{\zeta(2)}.
\end{eqnarray*}
A: Your doubts are justified:  $\sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2}$ is just a sum of limiting probabilities, so it should (and does) yield $1$ in the limit.  The sum you want to compute is of values of $\frac{1}{\gcd(a,b)}$, so you need to weight the possible values of that expression by their probabilities.  This gives
$$
\sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2}\cdot\frac{1}{d}\to\frac{\zeta(3)}{\zeta(2)}\approx 0.730763
$$
A numerical experimentat (with $n=1000$) gives a result close to $0.731$, so this seems likely to be correct.
