BZOJ 1101 Luogu P3455 POI 2007 Zap (莫比乌斯反演+数论分块)

手动博客搬家: 本文发表于20171216 13:34:20, 原地址https://blog.csdn.net/suncongbo/article/details/78819470

URL: (Luogu)https://www.luogu.org/problem/show?pid=3455
(BZOJ)http://www.lydsy.com/JudgeOnline/problem.php?id=1101

题目大意:
有t次询问((tle5e4)), 每次给定a,b,d, 询问有多少对(x,y)满足x<=a, y<=b, gcd(a,b)=d. 0<=d<=a,b<=5e4

思路分析:
首先，需要注意的是，要特殊处理(d=0)的情况，答案为0.
对于(dge1), 采用莫比乌斯反演解决:
先将a/=d, b/=d, 因此只需求gcd(x,y)=1的数的对数。
令F[i]表示(1le xle a,1le yle b)且(i|gcd(x,y))的a,b总数, f[i]表示gcd(x,y)=i的数的对数(此处a,b都已经除以d).因此问题转化为求f(1).
根据莫比乌斯反演公式:$$F(n)=sum_{n|x} f(x), f(n)=sum_{n|x} F(x)mu(frac{x}{n})$$
因此，(f(1)=sum_{1|x} F(x)mu(x))
而显然我们有(F(x)=[frac{a}{x}][frac{b}{x}]), 因此可以(O(1))地求出F(x), 也就可以(O(min(a,b)))地求出f(1)了。(莫比乌斯反演函数(mu(x))可在线性筛中求出)
可是这样还不够。算算复杂度，发现会TLE.
注意到一个性质: 对于(xlesqrt{a}), ([frac{a}{x}])的值变化得很快,([frac{a}{x}])的变化速度远高于(x)的变化速度。而对于(xgtsqrt{a}), ([frac{a}{x}])的值变化得很慢, 远低于(x)的变化速度。因此，我们可以求出所有使得([frac{a}{x}])的值变化的点x， 共有(O(sqrt{n}))个(实际上带一个常数2), 然后我们对b做同样的操作。将所有影响([frac{a}{x}])和([frac{b}{x}])的值的点都从小到大排序记录下来，处理莫比乌斯函数的前缀和, 每一个点代表一个区间，这个区间内所有的数([frac{a}{x}])与([frac{b}{x}])的值分别与这个数([frac{a}{x}])和([frac{b}{x}])相等。然后这一段区间对答案的贡献就是区间的(mu())之和乘以([frac{a}{x}][frac{b}{x}]).

代码实现

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;

const int N = 5e4;
const int NN = 317;
int p[N+2];
bool f[N+2];
int mu[N+2];
int s[N+2];
int g[(NN<<2)+2];
int h[(NN<<2)+2];
int a,b,d,m;

void Mobius()
{
	f[1] = true; mu[1] = 1; m = 0;
	for(int i=2; i<=N; i++)
	{
		if(!f[i]) {p[++m] = i; mu[i] = -1;}
		for(int j=1; p[j]*i<=N; j++)
		{
			f[p[j]*i] = true;
			if(i%p[j]==0)
			{
				mu[i*p[j]] = 0;
				break;
			}
			else mu[i*p[j]] = -mu[i];
		}
	}
}

void merge(int aa,int bb)
{
	int i = 1,j = (aa<<1)+1,k = 1;
	while(i<=(aa<<1) && j<=(aa<<1)+(bb<<1))
	{
		if(h[i]<h[j]) g[k++] = h[i++];
		else g[k++] = h[j++];
	}
	while(i<=(aa<<1)) g[k++] = h[i++];
	while(j<=(aa<<1)+(bb<<1)) g[k++] = h[j++];
}

int main()
{
	int t; scanf("%d",&t);
	Mobius(); s[0] = 0;
	for(int i=1; i<=N; i++) s[i] = s[i-1]+mu[i];
	while(t--)
	{
		scanf("%d%d%d",&a,&b,&d);
		if(d==0) {printf("0
"); continue;}
		if(a>b) swap(a,b);
		a /= d; b /= d;
		int aa = (int)sqrt(a),bb = (int)sqrt(b);
		long long ans = 0ll;
		for(int i=1; i<=aa; i++) h[i] = i;
		for(int i=aa; i>=1; i--) h[(aa<<1)-i+1] = a/i;
 //保证h[]在1~(aa<<1)范围内有序
		for(int i=1; i<=bb; i++) h[i+(aa<<1)] = i;
		for(int i=bb; i>=1; i--) h[(aa<<1)+(bb<<1)-i+1] = b/i;
 //保证h[]在1~(bb<<1)范围内有序
		merge(aa,bb);
 //将[1,aa<<1]与[aa<<1+1,aa<<1+bb<<1]归并起来
		for(int i=1; i<=(aa<<1)+(bb<<1); i++)
		{
			ans += (long long)(s[g[i]]-s[g[i-1]])*(a/g[i])*(b/g[i]);
		}
		printf("%lld
",ans);
	}
	return 0;
}