[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为质数

\[\sum\limits_p\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=p] \]

\[\sum\limits_p\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}[\gcd(i,j)=1] \]

\[\sum\limits_p\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}\sum\limits_{d|i,d|j}\mu(d) \]

\[\sum\limits_p\sum\limits_{d=1}^{\lfloor\frac{n}{p}\rfloor}\lfloor\dfrac{n}{dp}\rfloor\lfloor\dfrac{m}{dp}\rfloor \]

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

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

\[\sum\limits_{T=1}^n\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor\sum\limits_{p|T}\mu(\dfrac{T}{p}) \]

我们令\(f(T)=\sum\limits_{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]];

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

所以总时间复杂度为\(O(n+T\sqrt 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\n",Ans);
	}
	return 0;
}