# Two randomly chosen coprime integers

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.

• Heuristically one should expect that pairs of smaller numbers are more likely to be coprime than pairs of larger numbers, for the simple reason that larger numbers tend to have more prime factors, and hence more chances to have a common factor. – Travis Willse Jul 13 '16 at 20:07
• @Travis Yes, this was my intuition as well. – JasonM Jul 13 '16 at 20:09
• It makes sense that you get a number bigger than $0.75$, since your pdf implies you choose a $1$ half the time, so your two numbers are both bigger than $1$ only a quarter of the time. – Barry Cipra Jul 13 '16 at 20:19
• What was the original problem with answer $\frac6{\pi^2}$? – Mark Fischler Jul 13 '16 at 20:27
• You seem to assume that the event of a given variate being divisible by $p_1$ is independent of the event of that variable being divisible by $p_2$. That is not the case with your specified distribution. A consequence of this flaw is that the product shown is not the probability of $a$ and $b$ being coprime. Indeed, considering all non-coprime pairs with neither number greater than 12, you find that the probability of being coprime is less than $0.867973$. – Mark Fischler Jul 13 '16 at 20:55

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.

• Excellent answer, thank you – JasonM Jul 18 '16 at 17:55

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