您的位置:首页 > 其它

BZOJ2820 YY的GCD 【莫比乌斯反演】

2017-12-17 14:20 218 查看
2820: YY的GCD

Time Limit: 10 Sec Memory Limit: 512 MB

Submit: 2398 Solved: 1269

[Submit][Status][Discuss]

Description

神犇YY虐完数论后给傻×kAc出了一题给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对kAc这种

傻×必然不会了,于是向你来请教……多组输入

Input

第一行一个整数T 表述数据组数接下来T行,每行两个正整数,表示N, M

Output

T行,每行一个整数表示第i组数据的结果

Sample Input

2

10 10

100 100

Sample Output

30

2791

HINT

T = 10000

N, M <= 10000000

题解

还是经典的莫比乌斯反演

我们设f(n)表示gcd等于n的对数

我们设F(n)表示n|gcd的对数

则有

F(n)=⌊Nn⌋⌊Mn⌋

f(n)=∑n|dμ(dn)F(d)

=∑n|dμ(dn)⌊Nn⌋⌊Mn⌋

=∑
1aaeb
Ni=1μ(i)⌊Nn⌋⌊Mn⌋

ans=∑Np∈prime∑Ni=1μ(i)⌊Nn⌋⌊Mn⌋

=∑NT=1⌊Nn⌋⌊Mn⌋∗∑Np|Tμ(Tp)

对于前半部分∑NT=1⌊Nn⌋⌊Mn⌋我们可以分块O(2N−−√)

后半部分我们可以预处理,枚举每个质数以及它的倍数

由于质数的个数大约是nlogn,而每个质数倍数的枚举均摊是logn的,加上线性筛,总的预处理复杂度O(n)

总复杂度O(T∗N−−√+N)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<bitset>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u]; k != -1; k = ed[k].nxt)
using namespace std;
const int maxn = 10000005,maxm = 100005,INF = 1000000000;
inline int RD(){
int out = 0,flag = 1; char c = getchar();
while (c < 48 || c > 57) {if (c == '-') flag = -1; c = getchar();}
while (c >= 48 && c <= 57) {out = (out << 1) + (out << 3) + c - '0'; c = getchar();}
return out * flag;
}
bool isn[maxn];
int prime[maxn],primei = 0,mu[maxn];
LL A[maxn];
void init(){
mu[1] = 1;
for (int i = 2; i <= 10000000; i++){
if (!isn[i]) prime[++primei] = i,mu[i] = -1;
for (int j = 1; j <= primei && i * prime[j] <= 10000000; j++){
isn[i * prime[j]] = true;
if (i % prime[j] == 0) {mu[i * prime[j]] = 0; break;}
mu[i * prime[j]] = -mu[i];
}
}
REP(i,primei) for (int j = 1; prime[i] * j <= 10000000; j++) A[prime[i] * j] += mu[j];
REP(i,10000000) A[i] += A[i - 1];
}
int main(){
init();
int Tim = RD(),N,M;
LL ans;
while (Tim--){
N = RD(); M = RD(); ans = 0; if (N > M) swap(N,M);
for (int i = 1,nxt; i <= N; i = nxt + 1){
nxt = min(N / (N / i),M / (M / i));
ans += (A[nxt] - A[i - 1]) * (N / i) * (M / i);
}
printf("%lld\n",ans);
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: