Two randomly chosen coprime integers

Question

This is a twist on the problem commonly known to have solution $6/\pi^2$. Suppose when choosing from all natural numbers $\mathbb{N}$, the probability of choosing $n \in \mathbb{N}$ is given by $P(n)=\frac{1}{2^n}$. Now when choosing two natural numbers, what is the probability (in closed form) of choosing two coprime numbers?

Notice, the probability of choosing something divisible by $p$ is $$\frac{1}{2^p}+\frac{1}{2^{2p}}+\frac{1}{2^{3p}}+\frac{1}{2^{4p}}+\ldots=\frac{1}{2^p-1}$$

so the probability of choosing two numbers both divisible by $p$ is $$\frac{1}{(2^p-1)^2}$$

Meaning $$P(a,b;p)=1-\frac{1}{(2^p-1)^2}$$ where $P(a, b;p)$ is the probability that either $a$ or $b$ is not divisible by $p$. Then the answer I'm looking for is $$P(a,b)=\prod_{p\text{ prime}}P(a,b;p)=\prod_{p\text{ prime}}\left(1-\frac{1}{(2^p-1)^2}\right)$$ where $P(a,b)$ is the probability that $a$ and $b$ are coprime.

Anyway, I'm curious about a closed form expression for this number, similar to the original problem I mentioned. Any insight would be very helpful.

---

Edit

As Mark Fischler has pointed out below, this product representation assumes the events of $p|a$ and $p|b$ are independent, which should not be the case. If anyone can also explain a way of constructing a more correct probability, it would be very helpful.

Answer 1

By @miracle173 's answer, we are only left with $P(x\leq n, y\leq n, (x,y)=1)$. We can find an asymptotic formula as an application of Mobius function:

$$ \begin{align} P(x\leq n, y\leq n, (x,y)=1)&= \sum_{\substack{{x\leq n, y\leq n}\\{(x,y)=1}}} 2^{-x-y}\\ &= \sum_{x\leq n, y\leq n} 2^{-x-y} \sum_{d|(x,y)} \mu(d)\\ &=\sum_{d\leq n} \mu(d) \sum_{a\leq \frac nd, b\leq \frac nd} 2^{- d(a+b) }\ \ \ (\textrm{substitute }x=da, \ y=db)\\ &=\sum_{d\leq n}\mu(d) \left( \frac{\frac{1}{2^d}}{1-\frac{1}{2^d}}+O(2^{-n})\right)^2\\ &=\sum_{d\leq n}\mu(d) \frac{1}{(2^d-1)^2}+O(n 2^{-n}). \end{align} $$ Thus, the probability has to converge to $$ \sum_{d=1}^{\infty}\frac{\mu(d)}{(2^d-1)^2}\approx 0.867630801985022350790508146212902422392760107477\ldots $$ according to SAGE.

Answer 2

This is an extended comment and not an answer.

An approximation of the value can be found in the following way.

We have $$\begin{array}\\ P((x,y)=1) \\ = P(x \le n, y \le n | (x,y)=1) \\ +P(x \gt n, y \le n | (x,y)=1) +P(x \le n, y \gt n | (x,y)=1) +P(x \gt n, y \gt n | (x,y)=1) \end{array}$$

and $$\begin{array}\\ &P(x \gt n, y \le n | (x,y)=1) +P(x \le n, y \gt n | (x,y)=1) +P(x \gt n, y \gt n | (x,y)=1) \\ \le &P(x \gt n, y \le n ) +P(x \le n, y \gt n ) +P(x \gt n, y \gt n ) \\ \le & 2 P(x \gt n) P( x \le n ) + P( x \gt n)^2 \\ \lt & 2^{-n+1} \end{array}$$

So $$ \left| P((x,y)=1) - P(x \le n, y \le n | (x,y)=1) \right| < \frac{1}{2^{n-1}} $$

With Maxima I calculated the first 16 digits $0.8676308019850214$

(%i1) sum(sum(if gcd(i,j)=1 then 2^(-i-j) else 0,i,1,n),j,1,n),n=61,numer;
(%o1) 0.8676308019850214