Counting coprimes in a sequence

Question

Having a sequence of n <= 10^6 integers, all not exceeding m <= 3*10^6, I'd like to count how many coprime pairs are in it. Two numbers are coprime if their greatest common divisor is 1.

It can be done trivially in O(n^2 log n), but this is obviously way to slow, as the limit suggests something closer to O(n log n). One thing than can be done quickly is factoring out all the numbers, and also throwing out multiple occurences of the same prime in each, but that doesn't lead to any significant improvement. I also thought of counting the opposite - pairs that have a common divisor. It could be done in groups - firstly counting all the pairs that their smallest common prime divisor is 2, then 3, 5, and etc., but it seems to me like an other dead end.

Answer 1

I've come up with a slightly faster alternative based on your answer. On my work PC my C++ implementation (bottom) takes about 350ms to solve any problem instance; on my old laptop, it takes just over 1s. This algorithm avoids all division and modulo operations, and uses only O(m) space.

As with your algorithm, the basic idea is to apply the Inclusion-Exclusion Principle by enumerating every number 2 <= i <= m that contains no repeated factors exactly once, and for each such i, counting the number of numbers in the input that are divisible by i and either adding or subtracting this from the total. The key difference is that we can do the counting part "stupidly", simply by testing whether each possible multiple of i appears in the input, and this still takes just O(m log m) time.

How many times does the innermost line c += v[j].freq; in countCoprimes() repeat? The body of the outer loop is executed once for each number 2 <= i <= m that contains no repeated prime factors; this iteration count is trivially upper-bounded by m. The inner loop advances i steps at a time through the range [2..m], so the number of operations it performs during a single outer loop iteration is upper-bounded by m / i. Therefore the total number of iterations of the innermost line is upper-bounded by the sum from i=2 to m of m/i. The m factor can be moved outside the sum to get an upper bound of

m * sum{i=2..m}(1/i)

That sum is a partial sum in a harmonic series, and it is upper-bounded by log(m), so the total number of innermost loop iterations is O(m log m).

extendedEratosthenes() is designed to reduce constant factors by avoiding all divisions and keeping to O(m) memory usage. All countCoprimes() actually needs to know for a number 2 <= i <= m is (a) whether it has repeated prime factors, and if it doesn't, (b) whether it has an even or odd number of prime factors. To calculate (b) we can make use of the fact that the Sieve of Eratosthenes effectively "hits" any given i with its distinct prime factors in increasing order, so we can just flip a bit (the parity field in struct entry) to keep track of whether i has an even or odd number of factors. Each number starts with a prod field equal to 1; to record (a) we simply "knock out" any number that contains the square of a prime number as a factor by setting its prod field to 0. This field serves a dual purpose: if v[i].prod == 0, it indicates that i was discovered to have repeated factors; otherwise it contains the product of the (necessarily distinct) factors discovered so far. The (fairly minor) utility of this is that it allows us to stop the main sieve loop at the square root of m, instead of going all the way up to m: by now, for any given i that has no repeated factors, either v[i].prod == i, in which case we have found all the factors for i, or v[i].prod < i, in which case i must have exactly one factor > sqrt(3000000) that we have not yet accounted for. We can find all such remaining "large factors" with a second, non-nested loop.

#include <iostream>
#include <vector>

using namespace std;

struct entry {
    int freq;       // Frequency that this number occurs in the input list
    int parity;     // 0 for even number of factors, 1 for odd number
    int prod;       // Product of distinct prime factors
};

const int m = 3000000;      // Maximum input value
int n = 0;                  // Will be number of input values
vector<entry> v;

void extendedEratosthenes() {
    int i;
    for (i = 2; i * i <= m; ++i) {
        if (v[i].prod == 1) {
            for (int j = i, k = i; j <= m; j += i) {
                if (--k) {
                    v[j].parity ^= 1;
                    v[j].prod *= i;
                } else {
                    // j has a repeated factor of i: knock it out.
                    v[j].prod = 0;
                    k = i;
                }
            }
        }
    }
    
    // Fix up numbers with a prime factor above their square root.
    for (; i <= m; ++i) {
        if (v[i].prod && v[i].prod != i) {
            v[i].parity ^= 1;
        }
    }
}

void readInput() {
    int i;
    while (cin >> i) {
        ++v[i].freq;
        ++n;
    }
}

void countCoprimes() {
    __int64 total = static_cast<__int64>(n) * (n - 1) / 2;
    for (int i = 2; i <= m; ++i) {
        if (v[i].prod) {
            // i must have no repeated factors.
            
            int c = 0;
            for (int j = i; j <= m; j += i) {
                c += v[j].freq;
            }
            
            total -= (v[i].parity * 2 - 1) * static_cast<__int64>(c) * (c - 1) / 2;
        }
    }
    
    cerr << "Total number of coprime pairs: " << total << "\n";
}

int main(int argc, char **argv) {
    cerr << "Initialising array...\n";
    entry initialElem = { 0, 0, 1 };
    v.assign(m + 1, initialElem);
    
    cerr << "Performing extended Sieve of Eratosthenes...\n";
    extendedEratosthenes();
    
    cerr << "Reading input...\n";
    readInput();
    
    cerr << "Counting coprimes...\n";
    countCoprimes();
    
    return 0;
}