[BZOJ 2818]Gcd

Description

给定整数N，求1<=x,y<=N且Gcd(x,y)为素数的
数对(x,y)有多少对.

Input

一个整数N

Output

如题

Sample Input

Sample Output

HINT

对于样例(2,2),(2,4),(3,3),(4,2)
1<=N<=10^7

题解

由题，我们需要求

$$sum_{i=1} ^Nsum_{j=1} ^N[gcd(i,j)为质数]$$

我们不妨令$p=gcd(i,j)$，即$p$为质数，

显然我们不可能去枚举每个$i$,$j$，我们考虑去枚举每个质数$p$，原式化为

$$sum_psum_{i=1}^{lfloor {N over p} floor}sum_{j=1}^{lfloor {N over p} floor}[gcd(i,j)=1]$$

我们取出数对$(i,j)$

1.若$i>j$，我们发现对于式子

$$sum_psum_{i=1}^{lfloor {N over p} floor}sum_{j=1}^{lfloor {N over p} floor}[gcd(i,j)=1]$$

可化为

$$sum_psum_{i=1}^{lfloor {N over p} floor}sum_{j=1}^{i-1}[gcd(i,j)=1]$$

取出

$$sum_{i=1}^{lfloor {N over p} floor}sum_{j=1}^{i-1}[gcd(i,j)=1]$$

发现，这其实就是欧拉函数$φ$的定义式（先不考虑$φ(1)=1$）。

那么原式就可以化为

$$sum_psum_{i=1}^{lfloor {N over p} floor}φ(i)$$

结论。

2.若$i=j$，这种条件下只有$i=j=p$符合条件，显然最后我们只要将答案加上素数的个数就可以了。

3.若$i<j$，实际上只要交换$i$，$j$位置即可，我们只需要将1.得出的结论$×2$即可。

 1 #include<map>
 2 #include<set>
 3 #include<ctime>
 4 #include<cmath>
 5 #include<queue>
 6 #include<stack>
 7 #include<cstdio>
 8 #include<string>
 9 #include<vector>
10 #include<cstdlib>
11 #include<cstring>
12 #include<iostream>
13 #include<algorithm>
14 #define LL long long
15 #define RE register
16 #define IL inline
17 using namespace std;
18 const int N=1e7;
19 
20 int n;
21 
22 bool isprime[N+5];
23 int prime[N+5],top;
24 int phi[N+5];
25 IL void Pre();
26 
27 int main()
28 {
29     scanf("%d",&n);
30     Pre();
31     LL ans=0;
32     for (RE int i=1;i<=top;i++)
33     {
34         int lim=n/prime[i];
35         for (RE int j=1;j<=lim;j++) ans+=phi[j];
36     }
37     ans=ans*2+top;
38     printf("%lld
",ans);
39     return 0;
40 }
41 
42 IL void Pre()
43 {
44     for (RE int i=2;i<=n;i++)
45     {
46         if (!isprime[i]) phi[i]=i-1,prime[++top]=i;
47         for (RE int j=1;j<=top&&i*prime[j]<=n;j++)
48         {
49             isprime[i*prime[j]]=1;
50             if (!(i%prime[j])) {phi[i*prime[j]]=phi[i]*prime[j];break;}
51             else phi[i*prime[j]]=phi[i]*phi[prime[j]];
52         }
53     }
54 }

其实统计答案的第二层循环可以用前缀和优化。

 1 //It is made by Awson on 2018.1.12
 2 #include <set>
 3 #include <map>
 4 #include <cmath>
 5 #include <ctime>
 6 #include <queue>
 7 #include <stack>
 8 #include <cstdio>
 9 #include <string>
10 #include <vector>
11 #include <cstdlib>
12 #include <cstring>
13 #include <iostream>
14 #include <algorithm>
15 #define LL long long
16 #define Max(a, b) ((a) > (b) ? (a) : (b))
17 #define Min(a, b) ((a) < (b) ? (a) : (b))
18 #define Swap(a, b) ((a) ^= (b), (b) ^= (a), (a) ^= (b))
19 using namespace std;
20 const int N = 1e7;
21 void read(int &x) {
22     char ch; bool flag = 0;
23     for (ch = getchar(); !isdigit(ch) && ((flag |= (ch == '-')) || 1); ch = getchar());
24     for (x = 0; isdigit(ch); x = (x<<1)+(x<<3)+ch-48, ch = getchar());
25     x *= 1-2*flag;
26 }
27 void write(LL x) {
28     if (x > 9) write(x/10);
29     putchar(x%10+48);
30 }
31 
32 int n, isprime[N+5], prime[N+5], tot;
33 LL phi[N+5], ans;
34 
35 void get_phi() {
36     memset(isprime, 1, sizeof(isprime)); isprime[1] = 0;//, phi[1] = 1;
37     for (int i = 2; i <= n; i++) {
38     if (isprime[i]) phi[i] = i-1, prime[++tot] = i;
39     for (int j = 1; j <= tot && i*prime[j] <= n; j++) {
40         isprime[i*prime[j]] = 0;
41         if (!(i%prime[j])) {phi[i*prime[j]] = phi[i]*prime[j]; break; }
42         else phi[i*prime[j]] = phi[prime[j]]*phi[i];
43     }
44     }
45 }
46 void work() {
47     read(n); get_phi();
48     for (int i = 1; i <= n; i++) phi[i] += phi[i-1];
49     for (int i = 1; i <= tot; i++) ans += phi[n/prime[i]];
50     write(ans*2+tot);
51 }
52 int main() {
53     work();
54     return 0;
55 }