YY的GCD

题目描述

神犇YY虐完数论后给傻×kAc出了一题

给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对

kAc这种傻×必然不会了，于是向你来请教……

多组输入

输入输出格式

输入格式：

第一行一个整数T 表述数据组数

接下来T行，每行两个正整数，表示N, M

输出格式：

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

输入输出样例

输入样例#1：

2
10 10
100 100

输出样例#1：

30
2791

说明

T = 10000

N, M <= 10000000

---

好像是莫比乌斯反演的一道入门题

但是我也不会，然后对着题解抄了好久

题目要求(sum_{i=1}^{n}sum_{j=1}^{m}{gcd(i,j)∈prime})

所以设目标函数(f(n) = sum_{i=1}^{n}sum_{j=1}^{m}{[gcd(i,j)=n]})

然后再设(F(n) = sum_{n|d}{f(d)})

F(n)代表的是gcd是n及n的倍数的数目

所以 (Ans = sum_{p∈prime}sum_{i=1}^{n}sum_{j=1}^{m}{gcd(i,j)=p})

(Ans = sum_{p∈prime}f(p))

然后我们对f(p)套用反演公式

(Ans = sum_{p∈prime}sum_{p|d}{μ(frac{d}{p})F(d)})

设i = (frac{d}{p})

(Ans = sum_{p∈prime}sum_{i=1}^{min(n/p,m/p)}{μ(i)F(i*p)})

再设(k = i * p)

然后我们就只需要枚举k就好辣

(Ans = sum_{k=1}^{min(n,m)}sum_{p|k}{μ(k/p)*F(k)})

然后(F(k) = floor(n/k) * floor(m/k))

所以再对上个式子提公因式

把F(k)提出来

(Ans = sum_{k=1}^{min(n,m)}{floor(n/k) * floor(m/k)}sum_{p|k}{μ(k/p)})

然而这样还是不能通过此题

所以要整除分块

算出每个μ(i)对每个位置的贡献(即μ(i)*prime)

然后记录一下(μ)函数的前缀和

就可以(O(sqrt{n}))做了

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
# define LL long long
const int M = 10000005 ;
using namespace std ;
inline int read() {
	char c = getchar() ; int x = 0 , w = 1 ;
	while(c>'9'||c<'0') { if(c=='-') w = -1 ; c = getchar() ; }
	while(c>='0'&&c<='9') { x = x*10+c-'0' ; c = getchar() ; }
	return x*w ;
}
bool Notp[M] ;
int pnum , p[M] , mu[M] ;
int f[M] ;
LL sum[M] ;
inline void Get_Mu(int n) {
	Notp[1] = true ; mu[1] = 1 ;
	for(int i = 2 ; i <= n ; i ++) {
		if(!Notp[i]) p[++pnum] = i , mu[i] = -1 ;
		for(int j = 1 ; j <= pnum && i * p[j] <= n ; j ++) {
			Notp[i * p[j]] = true ;
			if(i % p[j] == 0) {
				mu[i * p[j]] = 0 ;
				break ;
			}
			mu[i * p[j]] = -mu[i] ;
		}
	}
	for(int i = 1 ; i <= pnum ; i ++)
	  for(int j = 1 ; p[i] * j <= n ; j ++)
	    f[p[i] * j] += mu[j] ;
	for(int i = 1 ; i <= n ; i ++) sum[i] = sum[i - 1] + f[i] ;
}

int main() {
	Get_Mu(10000000) ;
	int T = read() , n , m ;
	LL Ans ;
	while(T--) {
		Ans = 0 ;
		n = read() ; m = read() ;
		if(n > m) swap(n , m) ;
		for(int l = 1 , r ; l <= n ; l = r + 1) {
			r = min( n / (n / l) , m / (m / l) ) ;
			Ans += (LL)(n / l) * (m / l) * (sum[r] - sum[l - 1]) ;
		}
		printf("%lld
",Ans) ;
	}
	return 0 ;
}