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

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

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

思路分析:
首先，需要注意的是，要特殊处理d=0d=0的情况，答案为0.
对于d≥1d\ge1, 采用莫比乌斯反演解决:
先将a/=d, b/=d, 因此只需求gcd(x,y)=1的数的对数。
令F[i]表示1≤x≤a,1≤y≤b1\le x\le a,1\le y\le b且i|gcd(x,y)i|gcd(x,y)的a,b总数, f[i]表示gcd(x,y)=i的数的对数(此处a,b都已经除以d).因此问题转化为求f(1).
根据莫比乌斯反演公式:

F(n)=∑n|xf(x),f(n)=∑n|xF(x)μ(xn)

F(n)=\sum_{n|x} f(x), f(n)=\sum_{n|x} F(x)\mu(\frac{x}{n})
因此， f(1)=∑1|xF(x)μ(x)f(1)=\sum_{1|x} F(x)\mu(x)
而显然我们有 F(x)=[ax][bx]F(x)=[\frac{a}{x}][\frac{b}{x}], 因此可以 O(1)O(1)地求出F(x), 也就可以 O(min(a,b))O(min(a,b))地求出f(1)了。(莫比乌斯反演函数 μ(x)\mu(x)可在线性筛中求出)
可是这样还不够。算算复杂度，发现会TLE.
注意到一个性质: 对于 x≤a√x\le\sqrt{a}, [ax][\frac{a}{x}]的值变化得很快, [ax][\frac{a}{x}]的变化速度远高于 xx的变化速度。而对于x>a√x\gt\sqrt{a}, [ax][\frac{a}{x}]的值变化得很慢, 远低于 xx的变化速度。因此，我们可以求出所有使得[ax][\frac{a}{x}]的值变化的点x，共有 O(n√)O(\sqrt{n})个(实际上带一个常数2), 然后我们对b做同样的操作。将所有影响 [ax][\frac{a}{x}]和 [bx][\frac{b}{x}]的值的点都从小到大排序记录下来，处理莫比乌斯函数的前缀和, 每一个点代表一个区间，这个区间内所有的数 [ax][\frac{a}{x}]与 [bx][\frac{b}{x}]的值分别与这个数 [ax][\frac{a}{x}]和 [bx][\frac{b}{x}]相等。然后这一段区间对答案的贡献就是区间的 μ()\mu()之和乘以 [ax][bx][\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\n"); 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\n",ans);}return 0;
}

BZOJ 1101 Luogu P3455 POI 2007 Zap (莫比乌斯反演+分块)相关推荐

BZOJ 1101 Luogu P3455 POI 2007 Zap (莫比乌斯反演+数论分块)
BZOJ 1101 Luogu P3455 POI 2007 Zap (莫比乌斯反演+数论分块) 手动博客搬家: 本文发表于20171216 13:34:20, 原地址https://blog.csd ...
luogu P3455 [POI2007]ZAP-Queries （莫比乌斯反演 + 整除分块）
整理的算法模板合集: ACM模板题目传送门本题中数据为5e4,我们只需要筛一次5e4就行了. 双倍经验的P4450 双亲数中数据达到了1e6,我们直接筛1e6的莫比乌斯函数有点不可取,因为只有一组 ...
BZOJ 2301 Problem b（莫比乌斯反演+分块优化）
题目链接:http://acm.hust.edu.cn/vjudge/problem/viewProblem.action?id=37166 题意:对于给出的n个询问,每次求有多少个数对(x,y),满 ...
【bzoj 2820】YY的GCD（莫比乌斯反演）
传送门biu~ 参考博客:莫比乌斯反演 #include<bits/stdc++.h> #define N 10000000 using namespace std; int T,n,m, ...
ACdream 1148(莫比乌斯反演+分块)
传送门:GCD SUM 题意:给出N,M 执行如下程序: long long ans = 0,ansx = 0,ansy = 0; for(int i = 1; i <= N; i ++) ...
BZOJ 1101: [POI2007]Zap( 莫比乌斯反演 )
求 answer = ∑ [gcd(x, y) = d] (1 <= x <= a, 1 <= y <= b) . 令a' = a / d, b' = b / d, 化简一下得 ...
Luogu P6055 [RC-02] GCD（莫比乌斯反演，杜教筛）（这题乐死我了，真就图一乐呗）
整理的算法模板合集: ACM模板点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划 Weblink https://www.luogu.com.cn/problem/P6055 Prob ...
P2257 YY的GCD （莫比乌斯反演）
[题目链接] https://www.luogu.org/problemnew/show/P2257 // luogu-judger-enable-o2 /* -------------------- ...
莫比乌斯函数（bzoj 1101: [POI2007]Zap）
莫比乌斯反演: http://blog.csdn.net/acdreamers/article/details/8542292 课件:https://wenku.baidu.com/view/fbec ...

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

BZOJ 1101 Luogu P3455 POI 2007 Zap (莫比乌斯反演+分块)相关推荐

最新文章

热门文章