Description:
神犇YY虐完数论后给傻×kAc出了一题
给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对
kAc这种傻×必然不会了,于是向你来请教……
多组输入
Input:
第一行一个整数T 表述数据组数
接下来T行,每行两个正整数,表示N, M
T = 10000
N, M <= 10000000
Output
T行,每行一个整数表示第i组数据的结果
Sample Input:
2
10 10
100 100
Sample Output:
30
2791
题意为求出 T 组 i=1∑nj=1∑m[gcd(i,j)∈prime]
首先我们令 n<m ,之后对此式进行化简,将 gcd(i,j) 设为变量 p 得
p=1∑ni=1∑nj=1∑m[gcd(i,j)=p](p∈prime)
化简得
p=1∑ni=1∑⌊pn⌋j=1∑⌊pm⌋[gcd(i,j)=1]
利用莫比乌斯函数性质可得
p=1∑ni=1∑⌊pn⌋j=1∑⌊pm⌋d∣gcd(i,j)∑μ(d)
考虑枚举 d 得
p=1∑nd=1∑⌊pn⌋μ(d)i=1∑⌊pn⌋[d∣i]j=1∑⌊pm⌋[d∣j]
其中 i=1∑⌊pn⌋[d∣i] 为 [1,⌊pn⌋] 中 d 的倍数, j=1∑⌊pm⌋[d∣j] 为 [1,⌊pm⌋] 中 d 的倍数,两项可化简为 ⌊pdn⌋ 与 ⌊pdm⌋
这样可以推导出 i=1∑nj=1∑m[gcd(i,j)∈prime]=d=1∑⌊pn⌋μ(d)⌊pdn⌋⌊pdm⌋(p∈prime,n<m)
其中莫比乌斯函数 μ 可以用线性筛筛出,但是这样 O(n) 枚举 d 的计算复杂度在多组输入中还不是正确的复杂度
考虑继续化简该式,令 T=pd ,代入得
T=1∑np∣T∑μ(pT)⌊Tn⌋⌊Tm⌋
整理得
T=1∑n⌊Tn⌋⌊Tm⌋p∣T∑μ(pT)
发现其中 p∣T∑μ(pT) 部分可以做前缀和预处理,另外 ⌊Tn⌋⌊Tm⌋ 部分使用分块优化就可以愉(jian)快(nan)地 AC 掉这道题目啦
AC代码:
#include <bits/stdc++.h>
using namespace std;
const int maxn = (int)(1e7 + 5);
bool IsPrime[maxn];
int Tot;
int Prime[maxn];
int Mu[maxn];
long long F[maxn];
long long Prefix[maxn];
void Moblus() {
for (int i = 0; i < maxn; ++i) IsPrime[i] = true;
Mu[1] = 1;
for (int i = 2; i < maxn; ++i) {
if (IsPrime[i]) {
Prime[Tot++] = i;
Mu[i] = -1;
}
for (int j = 0; j < Tot && Prime[j] * i < maxn; ++j) {
IsPrime[i * Prime[j]] = false;
if (i % Prime[j] == 0) {
Mu[i * Prime[j]] = 0;
break;
}
Mu[i * Prime[j]] = -Mu[i];
}
}
for (int i = 0; i < Tot; ++i) {
for (int j = 1; Prime[i] * j < maxn; ++j) {
F[j * Prime[i]] += Mu[j];
}
}
for (int i = 1; i < maxn; ++i) Prefix[i] = Prefix[i - 1] + F[i];
}
long long Solve(int N, int M) {
if (N > M) swap(N, M);
long long Res = 0;
for (int t = 1, tp; t <= N; t = tp + 1) {
tp = min(N / (N / t), M / (M / t));
Res += (Prefix[tp] - Prefix[t - 1]) * (N / t) * (M / t);
}
return Res;
}
int T;
int N, M;
int main(int argc, char *argv[]) {
Moblus();
scanf("%d", &T);
for (int Case = 1; Case <= T; ++Case) {
scanf("%d%d", &N, &M);
printf("%lld\n", Solve(N, M));
}
return 0;
}