# The Probability That Two Numbers Are Coprime.

Let be a natural number. Denote by the set

of coprime pairs in , and let be their number. It is well-known that ; this was first proven by Ernesto Cesàro.

We are interested in itself. Even though we know that for some , we do not know how this can be chosen.

First note that , where

and that , where is Euler's function. We want to find an explicit lower bound for , which allows us to find an explicit lower bound for itself.

As in the proof of Theorem 330 in Hardy and Wright, we have

where is the Möbius function. Hardy and Wright compare that expression to and show that the error is .

More precisely, it is

First,

Next,

Since

we finally obtain

With a simple computer program, I verified that

for all , and the above bound shows

for . The computer verification took less than one second. We therefore obtain:

Theorem.

We have for all . Equality holds if and only if .

The program for verification of the bound is the following C++ program:

 1#include <iostream>
2#include <iomanip>
3#include <cmath>
4#include <vector>
5
6// Simple implementation of Euclid's algorithm
7unsigned long long gcd(unsigned long long a, unsigned long long b)
8{
9    while (true)
10    {
11        if (a == 0)
12            return b;
13        b %= a;
14        if (b == 0)
15            return a;
16        a %= b;
17    }
18}
19
20// Compute all primes up to a bound
21
22std::vector<unsigned> primes;
23
24void create_primes(unsigned max)
25{
26    // TODO: replace with sieve of Eratosthenes
27    primes.push_back(2);
28    for (unsigned p = 3; p <= max; p += 2)
29    {
30        bool prime = true;
31        for (unsigned i = 0; i < primes.size(); ++i)
32            if ((p % primes[i]) == 0)
33            {
34                prime = false;
35                break;
36            }
37        if (prime)
38            primes.push_back(p);
39    }
40}
41
42// Compute the Euler phi function (for not too large arguments)
43unsigned phi(unsigned n)
44{
45    if (n < 2)
46        return 1;
47    unsigned r = n;
48    bool found_divisor = false;
49    for (unsigned i = 0; i < primes.size(); ++i)
50        if ((n % primes[i]) == 0)
51        {
52            r = r / primes[i] * (primes[i] - 1);
53            found_divisor = true;
54        }
55    if (!found_divisor)
56        --r;
57    return r;
58}
59
60int main()
61{
62    // Compute probabilities up to:
63    unsigned N = 932453;
64    // unsigned N = 10000000;
65
66    // Compute enough primes
67    create_primes(sqrt(N));
68    std::cout << "Prepared " << primes.size() << " primes\n";
69
70    // Compute probabilities
71    unsigned long long C = 1;
72    for (unsigned n = 2; n < N; ++n)
73    {
74        C += phi(n); // number of pairs (x, y) with 1 <= x <= y <= n such that gcd(x, y) == 1
75        // Compute rational number representing probability
76        unsigned long long num = 2 * C - 1;
77        unsigned long long den = n;
78        den *= n;
79        unsigned long long g = gcd(num, den);
80        num /= g;
81        den /= g;
82        // Compute floating point approximation of probability
83        long double P = (long double)num / (long double)den;
84        // Output if probability is <= 0.6079
85        if (num <= (unsigned long long)((long double)0.6079 * den))
86            std::cout << n << ": " << C << ", " << num << "/" << den << " \\approx " << std::setprecision(20) << P << "\n";
87    }
88}


On my laptop, it runs in less than one second. For , it requires less than 22 seconds.