以小于O(n^2)的复杂度计算两个数组的余素对

Question

我是在一个挑战中遇到这个问题的。有两个大小都为N的数组A和B，我们需要返回对的计数(Ai，Bj)，其中gcd(A[i],B[j])==1和A[i] != B[j]。我只能想到蛮力法，它超过了测试用例的时间限制。

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]);
        }
    }
}

你能建议时间高效的算法来解决这个问题吗？

编辑:无法分享问题链接，因为这是一个招聘挑战。添加约束和我记得的输入/输出格式。

投入-

• 第一行将包含N，这是两个数组中的元素数。
• 第二行将包含N个分隔整数，数组A的元素。
• 第三行将包含N个分隔整数，数组B的元素。

产出-

• 按条件计算对Ai，Aj。

限制-

• 1 <= N <= 10^5
• 1< Ai，Bj <= 10^9其中i，j

Answer 1

第一步是使用埃拉托斯提尼筛计算到sqrt(10^9)的素数。然后，可以使用这个筛子快速查找小于10^9的任何数的所有素因子(参见下面代码示例中的getPrimeFactors(...)函数)。

接下来，对于每个具有素数因子的A[i]，我们计算出所有可能的子产品X - p0, p1, p0p1, p2, p0p2, p1p2, p0p1p2, p3, p0p3, ..., p0p1p2...pk，并在map cntp[X]中对它们进行计数。实际上，映射cntp[X]告诉我们可被X整除的元素的数量，其中X是素数的乘积，其幂为0或1。例如，对于数字A[i] = 12，素数因子是2, 3。我们将计算cntp[2]++、cntp[3]++和cntp[6]++。

最后，对于每个具有素数因子的B[j]，我们再次计算出所有可能的子产品X，并使用包容排除原则来计算所有非对等对C_j (即与B[j]共享至少一个素数的A[i]s的数目)。然后将数字C_j从对的总数- N*N中减去，以得到最终答案。

注:包容排除原则如下：

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

并解释了这样一个事实:在cntp[X]和cntp[Y]中，我们可以将相同的数字A[i]计算两次，因为它可以被X和Y整除。

下面是该算法的一个可能的C++实现，它产生与OP的朴素O(n^2)算法相同的结果：

// 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);
}

实例化

我无法分析估计复杂性，但以下是我的笔记本电脑上针对不同N和一致随机A[i]和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

比较而言，O(n^2)方法采用：

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

Answer 2

Python实现：

    import math
    from collections import defaultdict

    def sieve(MAXN):
        spf = [0 for i in range(MAXN)]
        spf[1] = 1
        for i in range(2, MAXN):
            spf[i] = i
        for i in range(4, MAXN, 2):
            spf[i] = 2
        for i in range(3, math.ceil(math.sqrt(MAXN))):
            if (spf[i] == i):
                for j in range(i * i, MAXN, i):
                    if (spf[j] == j):
                        spf[j] = i
        return(spf)

    def getFactorization(x,spf):
        ret = list()
        while (x != 1):
            ret.append(spf[x])
            x = x // spf[x]
        return(list(set(ret)))

    def coprime_pairs(N,A,B):
        MAXN=max(max(A),max(B))+1
        spf=sieve(MAXN)
        cntp=defaultdict(int)
        for i in range(N):
            f=getFactorization(A[i],spf)
            x=[[0,1]]
            for p in f:
                k=len(x)
                for i in range(k):
                    nn=x[i][0]+1
                    pp=x[i][1]*p
                    cntp[pp]+=1
                    x.append([nn,pp])
        cnt=0
        for i in range(N):
            f=getFactorization(B[i],spf)
            x=[[0,1]]
            for p in f:
                k=len(x)
                for i in range(k):
                    nn=x[i][0]+1
                    pp=x[i][1]*p
                    x.append([nn,pp])
                    if(nn%2==1):
                        cnt+=cntp[pp]
                    else:
                        cnt-=cntp[pp]
        return(N*N-cnt)

    import random
    N=10001
    A=[random.randint(1,N) for _ in range(N)]
    B=[random.randint(1,N) for _ in range(N)]
    print(coprime_pairs(N,A,B))