YY的gcd

YY的GCD

原题链接

题目

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

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

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

多组输入

输入输出格式

输入格式：

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

输出格式：

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

说明

• T=10000，N, M <= 10000000

解题思路

仔细读完题之后可以发现本题要求的答案ANS=(sum_{i=1}^{n}sum_{j=1}^{m}[gcd(x,y)=prim])

那么我们开始设出两个非常套路的函数

设(f(d)=sum_{i=1}^n sum_{j=1}^m[gcd(i,j)=d]),(F(n)=sum_{n|d}f(d)=lfloorfrac{N}{n} floorlfloorfrac{M}{n} floor)

(f(d))表示(gcd(i,j)=d)的方案数，(F(n))表示(gcd(i,j)=n)或(k*n,k∈N^*)的方案数

然后我们就可以开始用莫比乌斯函数愉快的推式子了

(f(n)=sum_{n|d}mu(lfloorfrac{d}{n} floor)F(d))

之后开始化简ANS

(Ans=sum_{pin prime}sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=p]=sum_{pin prime}f(p))

将刚刚反演得到的式子带入一下有(Ans=sum_{pin prim}sum_{p|d}mu(lfloorfrac{d}{p} floor)F(d))

由于此时存在一个(mu(lfloor frac dp floor)) 难以化简所以我们转而枚举 (lfloor frac d p floor)得到

(Ans=sum_{pin prim}sum_{d=1}^{min(lfloorfrac{n}{p} floor,lfloorfrac{m}{p} floor)}mu(d)F(dp)=sum_{pin prim}sum_{d=1}^{min(lfloorfrac{n}{p} floor,lfloorfrac{m}{p} floor)}mu(d)lfloorfrac{n}{dp} floorlfloorfrac{m}{dp} floor)

由于现在需要枚举d和p两个数我们观察最后两项(lfloorfrac{n}{dp} floorlfloorfrac{m}{dp} floor)

可以直接枚举(T=d*p)从而得到一个更简单的式子

(Ans=sum_{T=1}^{min(n,m)}sum_{t|T,tin prime}mu(lfloorfrac{T}{t} floor)lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor)

再次化简将最后两项(lfloorfrac{n}{dp} floorlfloorfrac{m}{dp} floor)提到前边即

(Ans=sum_{T=1}^{min(n,m)}lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor(sum_{t|T,tin prime}mu(lfloorfrac{T}{t} floor)))

然后这样我们就已经化简成功了

此时显然我们可以直接O(n)求解但是对于多组数据我们必须找到一种更优秀的写法来求解这个式子。

首先我们可以对于最后面的(sum_{t|T,tin prime}mu(lfloorfrac{T}{t} floor))求一个前缀和设为SUM

然后我们来看这个式子(Ans=sum_{T=1}^{min(n,m)}lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor SUM)

那么此时我们观察前两项这不就是裸的就是整除分块吗？

然后我们预处理出SUM在对前面的两项做一个整除分块时间复杂度(O(Tsqrt n))

代码

#include<bits/stdc++.h>
#define N 10000100
#define ll long long
using namespace std;
inline int read() {
	int n=1,num=0;char ch=getchar();
	while(!isdigit(ch))	{n=(ch=='-')?-1:1;ch=getchar();}
	while(isdigit(ch))	{num=(num<<1)+(num<<3)+(ch^48);ch=getchar();}
	return n*num;
}
int vis[N],prime[N],mu[N],g[N],tot;//g表示推出式子最后的Σ(t/d) ,sum为其前缀和 
ll sum[N];
void get_mu(int n)
{
	mu[1]=1;
	for(int i=2;i<=n;++i){
		if(!vis[i])	{mu[i]=-1;prime[++tot]=i;}
		for(int j=1;j<=tot&&prime[j]*i<=n;++j) {
			vis[i*prime[j]]=1;
			if(!(i%prime[j]))	break;
			else				mu[i*prime[j]]=-mu[i];
		}
	}
	for(int j=1;j<=tot;++j)
		for(int i=1;i*prime[j]<=n;++i)	g[i*prime[j]]+=mu[i];
	for(int i=1;i<=n;++i)	sum[i]=sum[i-1]+(ll)g[i];
}
int n,m,t;
int main()
{
	t=read();
	get_mu(10000000);
	while(t--) {
		n=read();m=read();
		if(n>m)	swap(n,m);
		ll ans=0;
		for(int l=1,r;l<=n;l=r+1) {
			r=min(n/(n/l),m/(m/l));
			ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
		}
		printf("%lld
",ans);
	}
	return 0;
}