YY的GCD[bzoj2820]

题目描述 Description

神犇YY虐完数论后给傻×kAc出了一题给定 N,M ,求 1≤x≤N,1≤y≤M 且 gcd(x,y) 为质数的 (x,y) 有多少对kAc这种傻×必然不会了，于是向你来请教……多组输入

输入描述 Input Description

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

输出描述 Output Description

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

样例输入 Sample Input

2

10 10

100 100

样例输出 Sample Output

30

2791

数据范围及提示 Data Size & Hint

T = 10000

N, M <= 10000000

分析

令 n≤m

ans=∑p∈prime∑i=1n∑j=1m(gcd(i,j)==p)=∑p∈prime∑d=1⌊np⌋μ(d)⌊npd⌋⌊mpd⌋

令 T=pd,我们枚举 T

ans=∑T=1n⌊nT⌋⌊mT⌋∑p|Tμ(TP)

∑p|xμ(xp) 可以枚举得到,最后分块求解即可

代码

#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;

typedef long long LL;

const int Size = 10000001;

int prime[Size],mu[Size],r[Size],tx[Size],s[Size];
bool no_prime[Size];
int t;

void get_prime(int);

int main(){

get_prime(Size-1);

scanf("%d",&t);

LL ans;
while( t-- ){

int n,m,delta;

scanf("%d%d",&n,&m);
if(n > m)
swap(n,m);

ans = 0;
for(int i=1;i<=n;i=delta+1){
delta = min(n/(n/i),m/(m/i));
ans = ans+1LL*(n/i)*(m/i)*(s[delta]-s[i-1]);
}

printf("%lld\n",ans);

}

return 0;

}

void get_prime(int x){

mu[1] = 1;

for(int i=2;i<=x;++i){
if(!no_prime[i]){
prime[++prime[0]] = i;
tx[i] = i;
mu[i] = -1;
}
int t=i<<1,j=1;
while(j<=prime[0] && t<=x){
no_prime[t] = true;
tx[t] = i;
if(i%prime[j] == 0){
mu[t] = 0;
break;
}
mu[t] = -mu[i];
t = i*prime[++j];
}
}
s[1] = r[1] = 0;
for(int i=1;i<=prime[0];++i){
for(int j=1;j<=x/prime[i];++j)
r[prime[i]*j] += mu[j];
}
for(int i=2;i<=x;++i)
s[i] = s[i-1]+r[i];

}