# 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
7 unsigned 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
22 std::vector<unsigned> primes;
23
24 void 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)
43 unsigned 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
60 int 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.

An alterantive approach can be found here (Proposition 3.2). There, we analyzed the probability that integers are coprime (as opposed to as here).