Count of co-prime pairs from two arrays in less than O(n^2) complexity

Question

I came to this problem in a challenge. There are two arrays A and B both of size of N and we need to return the count of pairs (A[i],B[j]) where gcd(A[i],B[j])==1 and A[i] != B[j]. I could only think of brute force approach which exceeded time limit for few test cases.

for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) {
        if(__gcd(a[i],b[j])==1) {
             printf("%d %d\n", a[i], b[j]);
        }
    }
}

Can you advice time efficient algorithm to solve this.

Edit: Not able to share question link as this was from a hiring challenge. Adding the constraints and input/output format as I remember.

Input -

• First line will contain N, the number of elements present in both arrays.
• Second line will contain N space separated integers, elements of array A.
• Third line will contain N space separated integers, elements of array B.

Output -

• The count of pairs A[i],A[j] as per the conditions.

Constraints -

• 1 <= N <= 10^5
• 1 < A[i],B[j] <= 10^9 where i,j < N

Answer 1

The first step is to use Eratosthenes sieve to calculate the prime numbers up to sqrt(10^9). This sieve can then be used to quickly find all prime factors of any number less than 10^9 (see the getPrimeFactors(...) function in the code sample below).

Next, for each A[i] with prime factors p0, p1, ..., pk, we compute all possible sub-products X - p0, p1, p0p1, p2, p0p2, p1p2, p0p1p2, p3, p0p3, ..., p0p1p2...pk and count them in map cntp[X]. Effectively, the map cntp[X] tells us the number of elements A[i] divisible by X, where X is a product of prime numbers to the power of 0 or 1. So for example, for the number A[i] = 12, the prime factors are 2, 3. We will count cntp[2]++, cntp[3]++ and cntp[6]++.

Finally, for each B[j] with prime factors p0, p1, ..., pk, we again compute all possible sub-products X and use the Inclusion-exclusion principle to count all non-coprime pairs C_j (i.e. the number of A[i]s that share at least one prime factor with B[j]). The numbers C_j are then subtracted from the total number of pairs - N*N to get the final answer.

Note: the Inclusion-exclusion principle looks like this:

C_j = (cntp[p0] + cntp[p1] + ... + cntp[pk]) -
      (cntp[p0p1] + cntp[p0p2] + ... + cntp[pk-1pk]) +
      (cntp[p0p1p2] + cntp[p0p1p3] + ... + cntp[pk-2pk-1pk]) -
      ...

and accounts for the fact that in cntp[X] and cntp[Y] we could have counted the same number A[i] twice, given that it is divisible by both X and Y.

Here is a possible C++ implementation of the algorithm, which produces the same results as the naive O(n^2) algorithm by OP:

// get prime factors of a using pre-generated sieve
std::vector<int> getPrimeFactors(int a, const std::vector<int> & primes) {
    std::vector<int> f;
    for (auto p : primes) {
        if (p > a) break;
        if (a % p == 0) {
            f.push_back(p);
            do {
                a /= p;
            } while (a % p == 0);
        }
    }
    if (a > 1) f.push_back(a);

    return f;
}

// find coprime pairs A_i and B_j
// A_i and B_i <= 1e9
void solution(const std::vector<int> & A, const std::vector<int> & B) {
    // generate prime sieve
    std::vector<int> primes;
    primes.push_back(2);

    for (int i = 3; i*i <= 1e9; ++i) {
        bool isPrime = true;
        for (auto p : primes) {
            if (i % p == 0) {
                isPrime = false;
                break;
            }
        }
        if (isPrime) {
            primes.push_back(i);
        }
    }

    int N = A.size();

    struct Entry {
        int n = 0;
        int64_t p = 0;
    };

    // cntp[X] - number of times the product X can be expressed
    // with prime factors of A_i
    std::map<int64_t, int64_t> cntp;

    for (int i = 0; i < N; i++) {
        auto f = getPrimeFactors(A[i], primes);

        // count possible products using non-repeating prime factors of A_i
        std::vector<Entry> x;
        x.push_back({ 0, 1 });

        for (auto p : f) {
            int k = x.size();
            for (int i = 0; i < k; ++i) {
                int nn = x[i].n + 1;
                int64_t pp = x[i].p*p;

                ++cntp[pp];
                x.push_back({ nn, pp });
            }
        }
    }

    // use Inclusion–exclusion principle to count non-coprime pairs
    // and subtract them from the total number of prairs N*N

    int64_t cnt = N; cnt *= N;

    for (int i = 0; i < N; i++) {
        auto f = getPrimeFactors(B[i], primes);

        std::vector<Entry> x;
        x.push_back({ 0, 1 });

        for (auto p : f) {
            int k = x.size();
            for (int i = 0; i < k; ++i) {
                int nn = x[i].n + 1;
                int64_t pp = x[i].p*p;

                x.push_back({ nn, pp });

                if (nn % 2 == 1) {
                    cnt -= cntp[pp];
                } else {
                    cnt += cntp[pp];
                }
            }
        }
    }

    printf("cnt = %d\n", (int) cnt);
}

Live example

I cannot estimate the complexity analytically, but here are some profiling result on my laptop for different N and uniformly random A[i] and B[j]:

For N = 1e2, takes ~0.02 sec
For N = 1e3, takes ~0.05 sec
For N = 1e4, takes ~0.38 sec
For N = 1e5, takes ~3.80 sec

For comparison, the O(n^2) approach takes:

For N = 1e2, takes ~0.00 sec
For N = 1e3, takes ~0.15 sec
For N = 1e4, takes ~15.1 sec
For N = 1e5, takes too long, didn't wait to finish