[BZOJ2820]YY的GCD

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


首先假定n<m,p为质数

[sumlimits_psumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=p] ]

[sumlimits_psumlimits_{i=1}^{lfloorfrac{n}{p} floor}sumlimits_{j=1}^{lfloorfrac{m}{p} floor}[gcd(i,j)=1] ]

[sumlimits_psumlimits_{i=1}^{lfloorfrac{n}{p} floor}sumlimits_{j=1}^{lfloorfrac{m}{p} floor}sumlimits_{d|i,d|j}mu(d) ]

[sumlimits_psumlimits_{d=1}^{lfloorfrac{n}{p} floor}lfloordfrac{n}{dp} floorlfloordfrac{m}{dp} floor ]

如果直接这样写的话,复杂度为(O(sumlimits_{p}sqrt{dfrac{n}{p}})),成功TLE,因此我们要做点优化

我们令T=dp,并改一下枚举顺序,得到

[sumlimits_{T=1}^nlfloordfrac{n}{T} floorlfloordfrac{m}{T} floorsumlimits_{p|T}mu(dfrac{T}{p}) ]

我们令(f(T)=sumlimits_{p|T}mu(dfrac{T}{p})),那么我们只要预处理出(f(T))和它的前缀和,就可以愉快的分块了

预处理的话可以线筛,也有个稍微懒点的方法,先预处理出所有的(mu(i)),然后枚举每个质数p,更新他们的倍数即可

for (int i=1;i<=prime[0];i++)
	for (int j=1;j<=n/prime[i];j++)
		f[prime[i]*j]+=mu[prime[i]];

由于(sumlimits_{i=1}^ndfrac{n}{i}approx nln n),每个素数平摊复杂度为(O(ln n)),根据素数定理,1(sim)n中素数个数接近(dfrac{n}{ln n}),所以预处理的时间复杂度近似(O(n))

所以总时间复杂度为(O(n+Tsqrt n))

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline int read(){
	int x=0,f=1;char ch=getchar();
	for (;ch<'0'||ch>'9';ch=getchar())	if (ch=='-')    f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar())	x=(x<<1)+(x<<3)+ch-'0';
	return x*f;
}
inline void print(int x){
	if (x>=10)	print(x/10);
	putchar(x%10+'0');
}
const int N=1e7,M=1e6;
int prime[M+10],miu[N+10],f[N+10],tot;
bool inprime[N+10];
void prepare(){
	miu[1]=1;
	for (int i=2;i<=N;i++){
		if (!inprime[i])	prime[++tot]=i,miu[i]=-1;
		for (int j=1;j<=tot&&i*prime[j]<=N;j++){
			inprime[i*prime[j]]=1;
			if (i%prime[j]==0){
				miu[i*prime[j]]=0;
				break;
			}
			miu[i*prime[j]]=-miu[i];
		}
	}
	for (int i=1;i<=tot;i++)
		for (int j=1;prime[i]*j<=N;j++)
			f[prime[i]*j]+=miu[j];
	for (int i=1;i<=N;i++)	f[i]+=f[i-1];
}
int main(){
	prepare();
	for (int Data=read();Data;Data--){
		int x=read(),y=read(),n=min(x,y),pos;
		ll Ans=0;
		for (int T=1;T<=n;T=pos+1){
			pos=min(x/(x/T),y/(y/T));
			Ans+=1ll*(f[pos]-f[T-1])*(x/T)*(y/T);
		}
		printf("%lld
",Ans);
	}
	return 0;
}
原文地址:https://www.cnblogs.com/Wolfycz/p/9482194.html