GCD素数数对计算

Formulating the approach

I’m working through the problem of counting pairs (x, y) where the GCD is a prime number. The solution needs to be efficient and handle up to 1e7 elements in C++.

Refining the approach

I'm working on counting pairs (x, y) where the GCD is a prime by examining x = p * a and y = p * b, ensuring gcd(a, b) = 1.

Piecing together

I’m noting that the formula for f(M) counts ordered pairs of (a, b) with gcd(a, b) = 1 up to M, leveraging the symmetry of numbers and known results.

Orchestrating calculations

Considering an efficient O(N log log N) method using prefix sums and the Mobius function, or applying a known identity for counting coprime pairs to streamline calculations.

Revealing patterns

It's interesting to see that the number of pairs (a, b) with gcd(a, b) = d can be expressed as the sum of f(M/d) for divisors of d.

Piecing together

I’m using the formula f(M) = sum_{k=1 to M} mu(k) * floor(M/k)^2, leveraging the Möbius function and inclusion-exclusion principle to efficiently count coprime pairs up to N.

Calculating pairs

I’m summing f(N/p) for each prime up to N. This involves computing coprime pairs up to N/p, which should be more manageable.

Refining the strategy

I’m now counting pairs (a, b) with gcd(a, b) = 1 and 1 ≤ a, b ≤ floor(N/p) for each prime p. The final task is to efficiently compute the sum for N up to 1e7.

Assessing distinct values

I’m analyzing the number of distinct values of floor(N/p) for each prime p, noting it could be up to O(N), which is quite large.

Leveraging computational methods

I’m examining efficient techniques to calculate floor(N/d) for each prime p, considering a sieve of Eratosthenes and the Möbius function, and preferring prefix sums for optimization.

Summing efficiently

I’m thinking through summing mu(d) * floor(M/d)^2 by splitting the domain into intervals where floor(M/d) remains constant, utilizing prefix sums of mu(d) for efficiency.

Calculating efficiently

I’m outlining the steps: determine distinct values of floor(M/d), calculate f(M) for each value using prefix sums of mu(d), and handle up to 664k distinct queries efficiently.

Adjusting my focus

I’m thinking that handling 2.1 billion operations might be tough. By noting that the number of distinct values of floor(N/d) for d up to N is about 2 sqrt(N), it’s clearer how to approach it.

Reducing the scope

By narrowing down operations to around 20 million, it’s possible to sieve primes up to N, find distinct floor(N/p) values for primes, and sum results for efficient calculation.

Calculating distinct values

I’m working on compressing distinct values of floor(N/d), computing f(value) in O(sqrt(value)) using a prefix sum of the Möbius function, and efficiently summing partial sums.

Mapping out the process

I’m detailing an algorithm: precompute the Möbius function, create prefix sums, and efficiently calculate f(M) for each distinct M using partial sums of the Möbius function.

Generating primes

I’m piecing together primes up to N using a sieve, computing M values, and efficiently calculating f(M) for each distinct M with prefix sums of the Möbius function.

Using 64-bit integers

I’m preparing for a large result by considering 64-bit integers to manage counts and ensure precision for pairs (x, y) where gcd(x, y) is a prime.

Mapping pairs

I’m identifying pairs (x, y) with gcd equal to prime numbers. In the example with N=4, four pairs meet this criterion.

Setting up the sieve

I’m initializing and updating the mu array using the sieve method, careful to set values to 0 where necessary and consider time constraints for efficient computation.

Mapping out the sieve

I’m detailing the sieve process to compute the mu array, identifying primes, and using prefix sums for efficient computation, while managing memory constraints for large arrays.

Managing memory efficiently

I'm mapping out ways to handle large arrays and calculations, like using 32-bit integers for prefix sums of mu and 64-bit for final results, to stay within memory limits.

Focusing on coding

I’m planning the code to handle large inputs, aiming to efficiently compute gcd pairs, and considering how to manage computational capacity and token limits to ensure an effective solution.

Thought about GCD and prime pairs，持续 36 秒

思路概要（请先阅读本节以理解整体算法，再看代码）

题目要求：
给定 $N$，统计所有满足 $1 \le x,y \le N$ 且 $\gcd(x,y)$ 为质数的有序对 $(x,y)$ 的数量。记该数量为 $\text{Ans}(N)$。

1. 数学转换

令质数 $p = \gcd(x,y)$。则可写

并且

从而对于固定的质数 $p$，能与之匹配的 $(x,y)$ 对的数量即是

其中

因此题目所求总数为

2. 计算 $f(M)$ 的方法

记

这是在区间 $[1,M]$ 上互质对 $(a,b)$ 的个数。
一个经典结果是：

其中 $\mu(\cdot)$ 是米宾函数（Möbius 函数），定义如下：

• $\mu(1) = 1$。
• 若 $n$ 含有平方因子（即存在 $p^2 \mid n$），则 $\mu(n) = 0$。
• 若 $n$ 是无平方因子的整数，并且含有 $k$ 个不同的质因子，则 $\mu(n) = (-1)^k$。

因此我们只要能快速地：

1. 预处理出 $\mu(d)$；
2. 对任意 $M \le N$，能在较快时间内求出 $$f(M) = \sum_{d=1}^{M} \mu(d)\,\Bigl\lfloor \tfrac{M}{d}\Bigr\rfloor^2,$$

就能解决问题。

3. 快速求 $\sum_{d=1}^{M} \mu(d) \lfloor \tfrac{M}{d}\rfloor^2$

直接按照定义做需要 $O(M)$ 时间，若对每个质数 $p$（可能多达 $\pi(N)\approx \tfrac{N}{\ln N}$ 个）都要做一次，会非常大，无法在 $N$ 最大到 $10^7$ 时完成。

但我们可以用分块（或称「分段」）技术，在约 $O(\sqrt{M})$ 时间内计算该和：
将 $\lfloor M/d \rfloor$ 的值记为 $v$。当 $v$ 固定时，$d$ 会落在一个连续区间内。具体做法：

1. 预先用大小为 $N$ 的数组存储前缀和 $\text{prefixMu}[i] = \sum_{k=1}^{i} \mu(k)$。
2. 计算 $\sum_{d=1}^{M} \mu(d) \,\lfloor M/d\rfloor^2$ 时，让变量 i 从 1 递增到 $M$，但不是逐一 +1，而是跳到 $\lfloor M/(M/i)\rfloor$ 的下一个，每次区间内 $\lfloor M/d \rfloor$ 值相同，可将该区间 $[i, r]$ 的 $\mu(d)$ 之和用前缀和一次性求出，再乘以相同的 $\lfloor M/i \rfloor^2$。
   这样整个求和只需要大约 $O(\sqrt{M})$ 次区间计算。

4. 处理所有 $\lfloor N/p\rfloor$ 的去重

在最终答案

中，不同质数 $p$ 可能得到相同的 $\lfloor N/p\rfloor$。我们只需去重这些值，然后对每个去重值 $M$ 计算一次 $f(M)$，放到哈希表或数组里，之后在循环中直接取值相加即可。

当 $N\le 10^7$ 时，质数最多约 66 万多个，但 $\lfloor N/p\rfloor$ 的不同取值个数典型情况下约为 $O(\sqrt{N})$ 级别（更严格的估计与分块技巧类似），大约几千到上万级，计算总量仍在合理范围内（配合 C++ 优化后可以通过）。

5. 主要步骤小结

1. 筛法预处理：
   • 最小质因子筛/线性筛计算出 $\mu(i)$（米宾函数），并顺便记录所有质数。
   • 同时计算出 $\text{prefixMu}[i]$ = $\sum_{k=1}^{i} \mu(k)$ 的前缀和。
2. 收集并去重所有 $\lfloor N / p\rfloor$（其中 $p$ 为质数）到一个数组 vals。
3. 对 vals 中每个元素 m，用前缀和 + 分块技术在 $O(\sqrt{m})$ 时间内计算 $$f(m) \;=\; \sum_{d=1}^{m} \mu(d)\,\Bigl\lfloor \tfrac{m}{d}\Bigr\rfloor^2.$$ 存入字典或数组 fval[m]。
4. 求答案：遍历所有质数 $p \le N$，令 m = floor(N/p)，把 fval[m] 加起来即可。

最终输出此和。

6. 示例校验

题目示例 $N=4$：

• 质数集合：$\{2,3\}$（1不算质数，4不是质数）
• 对应 $\lfloor 4/2\rfloor = 2,\; \lfloor 4/3\rfloor = 1$
• 计算：
  • $f(2)$: $\{(1,1),(1,2),(2,1),(2,2)\}$ 中 gcd=1 的对有 3 个，所以 $f(2)=3$。
    也可用公式： $$f(2) = \sum_{d=1}^{2} \mu(d) \left\lfloor \tfrac{2}{d}\right\rfloor^2 = \mu(1)\cdot (2^2) + \mu(2)\cdot (1^2) = 1\cdot 4 + (-1)\cdot 1 = 3.$$
  • $f(1)$: 只有 (1,1) 这一对，gcd=1，所以 $f(1)=1$。
• 最终 $\text{Ans}(4) = f(2) + f(1) = 3 + 1 = 4$，与样例吻合。

( 注：最后我们要数的是 gcd=质数 的对，这个数经过推导就是 $\sum_{p \le 4} f(\lfloor 4/p\rfloor)$。与直接列举得到 4 对完全一致。)

7. 代码实现

以下是一份相对完整的参考实现，使用了 C++ 常见 OI 写法，头文件 <bits/stdc++.h>，命名空间 std，并做了相应的注释。算法核心就在计算 $f(M)$ 函数时使用的「分块 + $\mu$ 前缀和」技巧。

重要提醒：

1. 若实际比赛/测评环境对内存非常苛刻（如只给 256MB），要注意 1e7 大小数组是否可行，可以考虑进行更细的空间优化或使用分段式计算。
2. 若时间紧张，也可能需进一步优化（如更高效的 IO, 更精简的循环等）。

#include <bits/stdc++.h>
using namespace std;

static const int MAXN = 10000000;  // 题目给出的 N 上限

// 全局数组，视情况可改为静态/动态分配。
// mu[i] 存储 Möbius 函数
// prefixMu[i] 存储 mu 的前缀和： prefixMu[i] = mu[1] + ... + mu[i]
static int mu[MAXN+1];
static int prefixMu[MAXN+1];
// 存放最小质因子或用于标记质数的数组，可用来线性筛。
static int spf[MAXN+1];  
// 存储所有质数
static vector<int> primes;

// 线性筛，计算 mu[] 和 primes。
void computeMuAndPrimes(int n) {
    // mu(1) = 1
    mu[1] = 1;
    for(int i=2; i<=n; i++){
        if(spf[i] == 0) {
            // i 是质数
            spf[i] = i;
            mu[i] = -1;     // 质数的 mu 值为 -1
            primes.push_back(i);
        }
        // 给当前的质数表做“筛”:
        for(int &p: primes){
            if(p > spf[i] || 1LL * p * i > n) break;
            spf[p*i] = p;
            if(i % p == 0){
                // 若 i 已经含有 p, 那么 p^2| (p*i), 则 mu=0
                mu[p*i] = 0;
                break;
            } else {
                mu[p*i] = - mu[i];
            }
        }
    }
    // 计算前缀和 prefixMu
    for(int i=1; i<=n; i++){
        prefixMu[i] = prefixMu[i-1] + mu[i];
    }
}

// 利用分块 + prefixMu[], 计算 f(M) = sum_{d=1..M} mu[d] * (floor(M/d))^2
// 时间复杂度约 O(sqrt(M))。
long long calcF(long long M) {
    // 特判
    if(M <= 0) return 0;
    // 分块求和
    long long ans = 0;
    for(long long i=1; i<=M; ){
        long long v = M / i;                 // floor(M/i)
        long long r = M / v;                // 使得 floor(M/d) = v 对应 d 在 [i, r]
        // 区间 [i..r] 的 mu(d) 之和
        long long sumMu = prefixMu[r] - prefixMu[i-1];
        ans += sumMu * (long long)v * v;
        i = r + 1;
    }
    return ans;
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int N;  
    cin >> N;
    // 预处理 mu 和 primes
    computeMuAndPrimes(N);

    // 收集所有 floor(N/p) (p 为质数)，并去重
    vector<int> vals;
    vals.reserve(primes.size());
    for(int p : primes){
        if(p > N) break;  // 安全判断
        vals.push_back(N/p);
    }
    // 排序并去重
    sort(vals.begin(), vals.end());
    vals.erase(unique(vals.begin(), vals.end()), vals.end());

    // 预先计算并保存 f(m) 的结果
    // 用一个 map<int,long long> 存，或者用数组(但要注意 max 值).
    // 因为可能出现的 m = floor(N/p) <= N. 直接开数组也行 (下标到N)。
    // 不过为了节省时间，这里只对出现的 m 计算。
    unordered_map<int, long long> fMap;
    fMap.reserve(vals.size()*4);  // 预留一些空间，避免频繁 rehash

    for(int m : vals) {
        fMap[m] = calcF(m);
    }

    // 最终答案 = sum_{p prime <= N} f( floor(N/p) )
    long long answer = 0;
    for(int p : primes){
        if(p > N) break;
        int m = N/p;
        answer += fMap[m];
    }

    cout << answer << "\n";

    return 0;
}

主要要点：

1. computeMuAndPrimes 使用了线性筛法，可以在 $O(N)$ 或 $O(N \log \log N)$ 左右的复杂度下同时得到：
   • spf[i]: i 的最小质因子（若 spf[i]==0 则 i 是质数）；
   • mu[i]: Möbius 函数；
   • primes: 从小到大的质数列表。
2. prefixMu[i]: $\sum_{k=1}^{i} \mu(k)$。后面分块时用前缀和快速求 $\sum_{d=i}^{r} \mu(d)$。
3. calcF(M) 函数是核心：它用到分块思路，只需大约 $O(\sqrt{M})$ 次循环即可完成 $\sum_{d=1}^M \mu(d)\lfloor M/d\rfloor^2$ 。
4. 因为 $\lfloor N/p \rfloor$ 可能重复，只需要对不同的取值计算一次 calcF(m) 即可；之后对每个质数查询并累加。
5. 为避免数据过大，所有与计数相关的变量要用 long long（或 64 位）存储；答案可能接近 $N^2$ 级别（如 $N=10^7$，答案量级可达 $10^{14}$，仍在 64 位范围内）。
6. 若担心内存，可改用更紧凑的数据结构或分段计算。上面是演示完整思路的写法，在一些平台下（如给 1~2GB 内存）通常还能运行通过。

时间与空间复杂度大致评估

• 筛法部分：$O(N)$ 或 $O(N \log \log N)$，对 $N=10^7$ 一般尚可（具体实现需注意常数优化）。
• distinct $\lfloor N/p \rfloor$ 的个数典型在 $O(\sqrt{N})$ 量级；每次 calcF(m) 约 $O(\sqrt{m}) \le O(\sqrt{N})$。故这部分大约 $O(\sqrt{N} \times \sqrt{N}) = O(N)$。
• 整体可在合理时间内完成。

如有需要，可以使用更快的读写手段（scanf/printf + __attribute__((optimize)) 等）以进一步加速。

小结

本题的关键在于：

• 把「$\gcd(x,y)$ 是质数」转化为对所有质数 $p$ 累加「$\gcd(a,b)=1$ 且 $a,b \le \lfloor N/p\rfloor$」的个数；
• 迅速统计区间内互质对数量，可利用 Möbius 函数分块计算；
• 小心实现、避免重复计算。

这样即可在 $N$ 较大的情况下仍保持可行的性能。祝一切顺利!