Find probability of two randomly picked integers (from n integers) to be relatively prime

Question

I came across this problem of finding said probability and my first attempt was to come up with following algorithm: I am counting number of pairs which are relatively prime.

int rel = 0
int total = n * (n - 1) / 2
for i in [1, n)
    for j in [i+1, n)
        if gcd(i, j) == 1
            ++rel;
return rel / total

which is O(n^2).

Here is my attempt to reducing complexity:

Observation (1): 1 is relatively prime to [2, n] so n - 1 pairs are trivial.

Observation (2): 2 is not relatively prime to even numbers in the range [4, n] so remaining odd numbers are relatively prime to 2, so

 #Relatively prime pairs = (n / 2) if n is even 

                         = (n / 2 - 1) if n is odd.

So my improved algorithm would be:

int total = n * (n - 1) / 2
int rel = 0
if (n % 2) // n is odd
    rel = (n - 1) + n / 2 - 1
else // n is even
    rel = (n - 1) + n / 2

for i in [3, n)
    for j in [i+1, n)
        if gcd(i, j) == 1
            ++rel;
return rel / total

With this approach I could reduce two loops, but worst case time complexity is still O(n^2).

Question: My question is can we exploit any mathematical properties other than above to find the desired probability in linear time?

Thanks.

Answer 1

You'll need to calculate the Euler's Totient Function for all integers from 1 to n. Euler's totient or phi function, φ(n), is a arithmetical function that counts the number of positive integers less than or equal to n that are relatively prime to n.

To calculate the function efficiently, you can use a modified version of Sieve of Eratosthenes.

Here is a sample C++ code -

#include <stdio.h>

#define MAXN 10000000

int phi[MAXN+1];
bool isPrime[MAXN+1];

void calculate_phi() {
    int i,j;
    for(i = 1; i <= MAXN; i++) {
        phi[i] = i;
        isPrime[i] = true;
    }

    for(i = 2; i <= MAXN; i++) if(isPrime[i]) {
        for(j = i+i; j <= MAXN; j+=i) {
            isPrime[j] = false;
            phi[j] = (phi[j] / i) * (i-1);
        }
    }

    for(i = 1; i <= MAXN; i++) {
        if(phi[i] == i) phi[i]--;
    }
}

int main() {
    calculate_phi();
    return 0;
}

It uses the Euler's Product Formula described on the Wikipedia page of Totient Function.

Calculating the complexity of this algorithm is a bit tricky, but it is much less than O(n^2). You can get results for n = 10^7 pretty quickly.

Answer 2

The number of integers in the range 0 .. n that are coprime to n is the Euler totient function of n. You are computing the sum of such values, e.g. called summatory totient function. Methods to compute this sum fast are for example described here. You should easily get a method with a better than quadratic complexity, depending on how fast you implement the totient function.

Even better are the references listed in the encyclopedia of integer sequences: http://oeis.org/A002088, though many of the references require some math skills.

Using these formulas you can even get an implementation that is sublinear.

Answer 3

For each prime p, probability of it dividing a randomly picked number between 1 and n is

[n / p] / n

([x] being the biggest integer not greater than x). If n is large, this is approximately 1/p.

The probability of it dividing two such randomly picked numbers is

([n / p] / n)²

Again, this is 1/p² for large n.

Two numbers are coprime if no prime divides both, so the probability in question is the product

Π_{p is prime}(1 - ([n / p] / n)²)

It is enough to calculate it for all primes less than or equal to n. As n goes to infinity, this product approaches 6/π².

I'm not sure you can use the totient function directly, as described in the other answers.