按顺序计算互质

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上限.内循环在i [2..m]范围内一次一步地前进,因此它在单个外循环迭代期间执行的操作数量由m/i上限.因此,最内线的迭代总数上限为从i = 2到m的m/i之和.m因子可以移到总和之外以获得上限

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

该和是谐波系列中的部分和,并且它由log(m)上限,因此最内循环迭代的总数是O(m log m).

extendedEratosthenes()旨在通过避免所有划分并保持O(m)内存使用来减少常数因素.countCoprimes()实际上所有人都需要知道数字2 <= i <= m是(a)它是否有重复的素数因子,如果没有,(b)它是否具有偶数或奇数个素数因子.为了计算(b)我们可以利用以下事实:Eratosthenes的筛子有效地"击中"任何给定的i,其不同的素因子按递增顺序排列,因此我们可以稍微翻转(parity字段struct entry)以跟踪是否我有一个偶数或奇数的因素.每个数字都以一个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;
}