计算序列中的互质数

Question

有一个 n

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

Answer 1

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

与您的算法一样，基本思想是通过枚举每个不包含重复因子的数字 2 这仍然只需要 O(m log m) 时间。

countCoprimes() 中最内层的行c += v[j].freq; 重复了多少次？对于每个不包含重复素因数的数字 2

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

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

extendedEratosthenes() 旨在通过避免所有除法并保持 O(m) 内存使用来减少常数因子。所有countCoprimes() 实际上需要知道一个数字 2 struct entry 中的parity 字段) 来跟踪 i 是否有偶数或奇数个因子。每个数字都以等于 1 的 prod 字段开头；记录 (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;
}

Answer 2

进一步利用我在问题中提到的想法，实际上我自己想出了一个解决方案。由于你们中的一些人可能对它感兴趣，我将简要介绍一下。它确实在 O(m log m + n) 中工作，我已经在 C++ 中实现了它并进行了测试 - 在不到 5 秒的时间内解决了最大的情况（10^6 个整数）。

我们有 n 个整数，都不大于 m。我们首先进行 Eratosthenes Sieve 将每个直到 m 的整数映射到它的最小素因子，允许我们在 O(log m) 时间内分解出不大于 m 的任何数字。那么对于所有给定的数字 A[i]，只要有一些素数 p 比 A[i] 的幂大于一，我们就将 A[i] 除以它，因为当询问两个数字是否互质时，我们可以省略指数。这让我们发现所有 A[i] 都是不同素数的乘积。

现在，让我们假设我们能够在合理的时间内构造一个表 T，使得 T[i] 是条目数 A[j] 使得 i 除以 A[j]。这在某种程度上类似于@Brainless 在他的第二个答案中采用的方法。快速构建表 T 是我在问题下方的 cmets 中谈到的技术。

从现在开始，我们将遵循包含-排除原则。有了 T，对于每个 i，我们计算 P[i] - 对 (j,k) 的数量，使得 A[j] 和 A[k] 都可以被 i 整除。然后计算答案，将所有 P[i] 相加，在那些 i 有偶数个素因数的 P[i] 之前取减号。请注意，i 的所有主要除数都是不同的，因为对于所有其他索引 i P[i] 等于 0。通过包含-排除，每对将仅计算一次。为了以不同的方式看待这一点，假设一对 A[i] 和 A[j] 恰好共享 k 个公素因数。然后这对将被计算 k 次，然后打折 kC2 次，计算 kC3 次，打折 kC4 次……对于 nCk，请参见牛顿符号。一些数学运算使我们看到所考虑的对将被计数 1 - (1-1)^k = 1 次，从而得出证明。

到目前为止所做的步骤需要 O(m log log m) 用于筛子和 O(m) 用于计算结果。最后要做的是构造数组 T。我们可以对每个 A[i] 只增加 T[j] 以使所有 j 除以 i。由于 A[i] 最多可以有 O(sqrt(A[i])) 个除数（实际上甚至更少），所以我们可以在 O(n sqrt m) 中构造 T。但我们可以做得更好！

取二维数组 W。在每个时刻，以下不变量成立 - 如果对于每个非零 W[i][j]，我们将表 T 中的计数器增加 W[i][j] 以用于所有数字将 i 相除，并且共享 i 在 i 的 j 个最小素数除数中的确切指数，那么 T 将被正确构造。由于这可能看起来有点令人困惑，让我们看看它的实际效果。在开始时，为了使不变量为真，对于每个 A[i]，我们只需增加 W[A[i]][0]。另请注意，不超过 m 的数字最多可以有 O(log m) 个素数除数，因此 W 的整体大小为 O(m log m)。现在我们看到存储在 W[i][j] 中的信息可以通过以下方式“向前推进”：假设 p 是 i 的第 (j+1) 个素数除数，假设它有一个。然后 i 的某个除数可以使 p 的指数与 i 相同，或者更低。第一种情况是 W[i][j+1] - 我们添加另一个必须被除数“完全取走”的素数。第二种情况是 W[i/p][j] 作为 i 的除数，没有最高指数的 p 也必须除以 i/p。就是这样！我们按降序考虑所有 i，然后按升序考虑 j。我们从 W[i][j]“推进”信息。看看如果 i 正好有 j 个素数除数，那么它的信息就不能被推送，但我们真的不需要那个！如果 i 有 j 个素数除数，则 W[i][j] 基本上说：增加 W[i][j] only 数组 T 中的索引 i。所以当所有信息都被推送到每个 W[i] 中的“最后一行”我们通过这些行并完成构造 T。由于 W[i][j] 的每个单元格都被访问过一次，因此该算法需要 O(m log m) 时间，并且也是 O (n) 在开始时。至此，施工结束。以下是实际实现中的一些 C++ 代码：

FORD(i,SIZE(W)-1,2) //i in descending order
{
    int v = i, p;

    FOR(j,0,SIZE(W[i])-2) //exclude last row
    {
        p = S[v]; //j-th divisor; S[v] - smallest prime divisor of v
        while (v%p == 0) v /= p;

        W[i][j+1] += W[i][j];
        W[i/p][j] += W[i][j];
    }

    T[i] = W[i].back();
}

最后我想说我认为数组 T 可以比我展示的更快、更简单地构建。如果有人对如何做到这一点有一些巧妙的想法，我将不胜感激。

Answer 3

这是一个基于完整序列公式1..n 的想法，可在http://oeis.org/A018805 上找到：

a(n) = 2*( Sum phi(j), j=1..n ) - 1, where phi is Euler's totient function

遍历序列S。对于每个术语，S_i：

  对于S_i 的每个主要因数p：
如果p 的哈希不存在：
创建一个索引为p 的哈希，它指向S 的所有索引集，i 除外，
和一个计数器设置为 1，表示到目前为止 S 有多少项可以被 p 整除
其他：
删除现有索引集中的i 并增加计数器

  将S_i 的主要因子的哈希值按其计数器降序排列。从
开始 最大的计数器（这意味着最小的集合），创建一个索引列表，直到 i 也是
下一个最小集合的成员，直到集合用完。添加剩余数量
列表中的索引到累计总数。

例子：

sum phi' [4,7,10,15,21]

S_0: 4
prime-hash [2:1-4], counters [2:1]
0 indexes up to i in the set for prime 2
total 0

S_1: 7
prime hash [2:1-4; 7:0,2-4], counters [2:1, 7:1]
1 index up to i in the set for prime 7
total 1

S_2: 10
prime hash [2:1,3-4; 5:0-1,3-4; 7:0,2-4], counters [2:2, 5:1, 7:1]
1 index up to i in the set for prime 2, which is also a member 
of the set for prime 5
total 2

S_3: 15
prime hash [2:1,3-4; 5:0-1,4; 7:0,2-4; 3:0-2,4], counters [2:2: 5:2, 7:1, 3:1]
2 indexes up to i in the set for prime 5, which are also members 
of the set for prime 3
total 4

S_4: 21
prime hash [2:1,3-4; 5:0-1,4; 7:0,2-3; 3:0-2], counters [2:2: 5:2, 7:2, 3:2]
2 indexes up to i in the set for prime 7, which are also members 
of the set for prime 3
total 6

6 coprime pairs:
(4,7),(4,15),(4,21),(7,10),(7,15),(10,21)

Answer 4

我建议：

1) 使用 Eratosthene 获取 10^6 以下的已排序素数列表。

2) 对于列表中的每个数字 n，获取它的质因数。将它与另一个数 f(n) 关联起来如下：假设 n 的质因数是 3、7 和 17。那么 f(n) 的二进制表示是：

`0 1 0 1 0 0 1`

第一个数字（此处为 0）与质数 2 相关联，第二个（此处为 1）与质数 3 相关联，依此类推......

因此2 个数 n 和 m 互质当且仅当 f(n) & f(m) = 0。

3) 很容易看出有一个 N 使得对于每个 n : f(n)

`1 1 1 1 1 1 1 1 1 1 1 1 1 1 1`

这里N是上面序列中1的个数。获取这个 N 并对数字列表 f(n) 进行排序。我们称这个列表为 L。 如果要优化：在这个列表中，不要对重复项进行排序，而是存储一个包含 f(n) 和 f(n) 重复次数的对。

4) 以这种方式从 1 迭代到 N：初始化 i = 1 0 0 0 0，并在每次迭代时将数字 1 向右移动，所有其他值保持为 0（使用位移位实现）。

在每次迭代中，迭代 L 以获得 L 中元素 l 的数量 d(i)，使得 i & l != 0（如果使用上述优化，请小心）。换句话说，对于每个 i，获取 L 中不与 i 互质的元素的数量，并将此数字命名为 d(i)。加总

D = d(1) + d(2) + ... + d(N)

5) 这个数字 D 是原始列表中不互质的对数。互质数对数为：

M*(M-1)/2 - D

其中 M 是原始列表中的元素数。这种方法的复杂度是O(n log(n))。

祝你好运！

Answer 5

我之前的回答错了，抱歉。我在这里提出修改：

一旦你得到列表中每个数字的素数除数，将列表中以 p 作为除数的数字的数量 l(p) 与每个素数 p 相关联。比如素数5，列表中能被5整除的数是15、100、255，那么l(5)=3。

要在 O(n logn) 中实现它，迭代列表，并且对于列表中的每个数字，迭代它的主要因子；对于每个素数 p，增加它的 l(p)。

那么不互质且能被p整除的对数为

l(p)*(l(p) - 1) / 2

将这个数字与所有素数 p 相加，您将得到列表中不互素的对数（注意 l(p) 可以为 0）。假设这个和是D，那么答案是

M*(M-1)/2 - D

其中 M 是列表的长度。祝你好运！