Three coprime integers, Geometric Probability

Problem 56.1

Let $a, b, c$ be three integers chosen independently and uniformly from $\{1, 2, \ldots, N\}$ . As $N \to \infty$ , what is the probability that $\gcd(a, b, c) = 1$ ?

Solution. Repeat the Euler-product argument of Chapter 52. For each prime $p$ , let $A_p = \{ p \mid a \;\text{and}\; p \mid b \;\text{and}\; p \mid c \}$ . Since each of $a, b, c$ is divisible by $p$ with probability $1/p + O(1/N)$ , and the three divisibility events are independent, $\mathbb{P}(A_p) \to \frac{1}{p^3} \quad \text{as } N \to \infty.$ The triple $(a, b, c)$ has $\gcd = 1$ iff $\bigcap_p A_p^{c}$ , and the events $\{A_p\}$ across distinct primes are essentially independent in the large- $N$ limit (joint divisibility factorises by the Chinese remainder theorem). Hence $\mathbb{P}(\gcd = 1) \to \prod_p \left(1 - \frac{1}{p^3}\right).$

To justify passing from finitely many primes to all primes, first retain only primes $p\leq M$ . Their joint divisibility probabilities converge as stated. Uniformly in $N$ , the probability that a prime $p>M$ divides all 3 sampled integers is at most $p^{-3}$ , so the union bound makes the omitted probability at most $\sum_{m>M}m^{-3}$ , which tends to $0$ as $M\to\infty$ .

By the Euler product for the Riemann zeta function, $\prod_p (1 - p^{-s})^{-1} = \zeta(s)$ , so $\prod_p \left(1 - \frac{1}{p^3}\right) = \frac{1}{\zeta(3)}.$ Therefore $\boxed{\mathbb{P}(\gcd(a, b, c) = 1) \to \frac{1}{\zeta(3)} \approx 0.8319.} \qedhere$

The constant $\zeta(3) = 1.20205\ldots$ is called Apéry’s constant; Roger Apéry proved in 1978 that it is irrational, a result that surprised number theorists because the analogous $\zeta(2) = \pi^2/6$ is a transcendental number by Lindemann’s theorem, but $\zeta(3)$ has no known closed form in terms of $\pi$ or other classical constants. So unlike the two-variable case of Chapter 52 (where the answer is $6/\pi^2$ ), the three-variable coprime probability is a constant with no known simpler form.

The general formula, for $k$ integers uniformly distributed, is $1/\zeta(k)$ ; the limit as $k \to \infty$ is $1$ (all large- $k$ -tuples are asymptotically coprime, since the probability that some fixed prime divides all of them decays as $p^{-k}$ ).

import numpy as np
Nmax, N = 10**4, 10**6
a, b, c = np.random.randint(1, Nmax+1, (3, N))
g = np.gcd(np.gcd(a, b), c)
p = (g == 1).mean()
zeta3 = 1.2020569031595942  # Apery's constant
print(f"sim: {p:.5f}   exact: {1/zeta3:.5f}")
# sim: 0.83191   exact: 0.83191