algorithm - 计算序列中的互素数

Question

有一个 n <= 10^6 个整数序列，所有整数都不超过 m <= 3*10^6，我想计算其中有多少个互质对。如果两个数的最大公约数是 1，则它们互质。

它可以在 O(n^2 log n) 中轻松完成，但这显然是减慢速度的方式，因为极限表明更接近 O(n log n)。可以快速完成的一件事是分解出所有数字，并在每个数字中抛出多次出现的相同素数，但这并没有带来任何显着的改进。我还考虑过计算相反数 - 具有公约数的对。它可以分组进行 - 首先计算所有对，它们的最小公素因数是 2，然后是 3、5 等等，但在我看来这就像另一个死胡同。

Answer 1

根据您的回答，我想出了一个稍微快一些的替代方案。在我的工作 PC 上，我的 C++ 实现(底部)需要大约 350ms 来解决任何问题实例；在我的旧笔记本电脑上，只需 1 秒多一点。该算法避免了所有除法和取模运算，仅使用 O(m) 空间。

与您的算法一样，基本思想是通过枚举每个不包含重复因子的数字 2 <= i <= m 恰好一次来应用包含 - 排除原则，并且对于每个这样的 i，计算数字的数量可被 i 整除的输入，并从总数中加上或减去它。关键区别在于我们可以“愚蠢地”完成计数部分，只需测试输入中是否出现 i 的每个可能倍数，这仍然只需要 O(m log m) 时间。

最里面的线c += v[j].freq;是多少次在countCoprimes()重复？外循环体对每个不包含重复质因数的数 2 <= i <= m 执行一次；此迭代计数通常以 m 为上限。内循环在 [2..m] 范围内一次前进 i 步，因此它在单个外循环迭代期间执行的操作数上限为 m/i。因此最内层线的总迭代次数上限为从i=2到m的m/i的总和。可以将 m 因子移到总和之外以获得上限

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

那个和是调和级数的部分和，并且it is upper-bounded by log(m) , 所以最内层循环迭代的总数是 O(m log m)。

extendedEratosthenes()旨在通过避免所有除法并保持 O(m) 内存使用来减少常数因子。全部countCoprimes()实际上需要知道一个数字 2 <= i <= m 是 (a) 它是否有重复的质因数，如果没有，(b) 它是否有偶数或奇数个质因数。为了计算 (b)，我们可以利用这样一个事实，即埃拉托色尼筛法有效地“命中”任何给定的 i 及其不同的素因子，因此我们可以稍微翻转一下(parity 中的 struct entry 字段) ) 来跟踪 i 的因子数是偶数还是奇数。每个数字都以 prod 开头字段等于 1；为了记录 (a)，我们通过设置其 prod 来简单地“剔除”任何包含质数平方作为因数的数。字段为 0。此字段有双重用途:如果 v[i].prod == 0 ，说明我被发现有重复因素；否则它包含迄今为止发现的(必然不同的)因素的产物。这个(相当小的)效用是它允许我们在 m 的平方根处停止主筛循环，而不是一直到 m: 到目前为止，对于没有重复因子的任何给定 i，要么v[i].prod == i ，在这种情况下，我们找到了 i 或 v[i].prod < i 的所有因子，在这种情况下，我必须恰好有一个因子 > sqrt(3000000) 我们还没有考虑到。我们可以使用第二个非嵌套循环找到所有这些剩余的“大因素”。

#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;
}