莫比乌斯反演

P3455 [POI2007]ZAP-Queries

题意：求$gcd(x,y)=2$的对数，$x,y$满足$1\le x\le a,1\le y\le b$
思路：有T组数据，T最大是50000，a、b最大也是50000，直接用莫比乌斯反演，会超时。所以采用整除分块的方法，复杂度可以从$O(n)$ 降为 $O(\sqrt{n})$
莫比乌斯反演倍数公式：
设$f(d)$为$gcd(x,y)=d$ 且 $1\le x\le a,1\le y\le b$的对数
则$F(d)$为$d|gcd(x,y)$ 且 $1\le x\le a,1\le y\le b$的对数
$F(d)$可以直接求解：$F(d)=\frac{a}{d}\times \frac{b}{d}$
反演可得：$$f(n)=\sum _{n|d}u(\frac{d}{n})F(d)$$

其中 $f(1)={\sum }_{1|d}^{min(a,b)}u(d)F(d)$就是我们求解的答案

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#define ll long long
using namespace std;
const int maxn=5e4+10,INF=0x3f3f3f3f;
const int mod=1e9+7;

int T;
const int N=5e4+5;
int prime[maxn],visit[maxn];
int mu[maxn],pref[maxn],cnt;

void Mobius()
{
	memset(visit,0,sizeof(visit));
	memset(mu,0,sizeof(mu));
	cnt=0,mu[1]=1,visit[1]=1;
	for(int i=2;i<N;++i)
	{
		if(!visit[i])
			prime[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*prime[j]<N;++j)
		{
			visit[i*prime[j]]=1;
			if(i%prime[j]==0)
				break;
			mu[i*prime[j]]=-mu[i];
		}			
	}	
	
	for(int i=1;i<N;++i)
		pref[i]=pref[i-1]+mu[i];
	
}

ll solve(ll a,ll b,ll k)
{
	a/=k,b/=k;
	ll ans=0,r;
	if(a>b)
		swap(a,b);
	for(int l=1;l<=a;l=r+1)
	{
		r=min(a/(a/l),b/(b/l));
		ans+=1ll*(pref[r]-pref[l-1])*(a/l)*(b/l);
	}
	return ans;
	
}

int main()
{
	Mobius();
	scanf("%d",&T);
	while(T--)
	{
		ll a,b,k;
		scanf("%lld%lld%lld",&a,&b,&k);
		ll ans=solve(a,b,k);
		printf("%lld\n",ans);
	}  
	return 0;
}