215. 破译密码
题目链接
215. 破译密码
达达正在破解一段密码，他需要回答很多类似的问题：
对于给定的整数 \(a,b\) 和 \(d\)，有多少正整数对 \(x,y\)，满足 \(x \le a，y \le b\)，并且 \(gcd(x,y)=d\)。
作为达达的同学，达达希望得到你的帮助。
输入格式
第一行包含一个正整数 \(n\)，表示一共有 \(n\) 组询问。
接下来 \(n\) 行，每行表示一个询问，每行三个正整数，分别为 \(a,b,d\)。
输出格式
对于每组询问，输出一个正整数，表示满足条件的整数对数。
数据范围
\(1 \le n \le 50000\),
\(1 \le d \le a,b \le 50000\)
输入样例：
2
4 5 2
6 4 3
输出样例：
3
2
提示:\(gcd(x,y)\) 返回 \(x，y\) 的最大公约数。
解题思路
容斥原理，莫比乌斯函数
莫比乌斯函数：\(mobius[i]\)：\(mobius[i]\)，对 \(i\) 分解质因数，若 \(i\) 某一项质因子个数大于 \(1\)，则 \(mobius[i]=0\)，如果质因子种数为奇数则 \(mobius[i]=-1\)，否则 \(mobius[i]=1\)。理解：对于 \(n\) 的欧拉函数，由容斥原理，得 $\phi(n)=n-s_2-s_3-\dots + s_{2,3}+s_{3,5}+\dots - s_{2,3,5}-\dots +\dots $，其中 \(s_{a_1,a_2,\dots a_m}\) 表示与 \(n\) 的公因子为 \(a_1,a_2,\dots a_m\) 的个数，其系数正好为莫比乌斯函数
由于 \(gcd(x,y)=d\)，即 \(gcd(x/d,y/d)=1\)，问题转换为满足 \(x\in[1,a/d],y\in[1,b/d]\) 且 \(gcd(x,y)=1\) 的对数，由容斥原理，答案即为 \((a/d)\times (b/d)-s_2-s_3-\dots +s_6+\dots -\dots\)，令 \(a=a/d,b=b/d\)，答案即为 \(a\times b-(a/2)\times (b/2)-(a/3)\times (b/3)+(a/6)\times (b/6)+\dots - \dots\)，令 \(n=min(a,b)\)，即为 \(\sum_{i=1}^{n}(a/i)\times (b/i)\times mobius[i]\)，其中 \(a/i\) 共有 \(O(\sqrt{a})\) 种取值，简略说下，将 \(a\) 分成两段：\([1,\sqrt{a}],[\sqrt{a},a]\)，\(a/i\) 分成的两段的值的范围都为 \(O(\sqrt{a})\)。
另外还有如下定理：
对于 \(a/x\) 该段取值的最大的一个 \(x\) 为 \(g(x)=a/(a/x)\)，\(a,b\) 表示的段都是 \(O(根号)\) 级别，每次跳到 \(a,b\) 表示的段的较小边界即可
- 时间复杂度：\(O(n\sqrt{n})\)
代码
// Problem: 破译密码
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/217/
// Memory Limit: 64 MB
// Time Limit: 2000 ms
//
// Powered by CP Editor (https://cpeditor.org)
// %%%Skyqwq
#include <bits/stdc++.h>
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
template <typename T> void inline read(T &x) {
int f = 1; x = 0; char s = getchar();
while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
x *= f;
}
const int N=5e4+5;
LL res;
int prime[N],u[N],s[N],v[N],m;
void primes(int n)
{
u[1]=1;
for(int i=2;i<=n;i++)
{
if(v[i]==0)
{
u[i]=-1;
v[i]=prime[++m]=i;
}
for(int j=1;prime[j]*i<=n&&j<=m;j++)
{
if(v[i]<prime[j])break;
if(i%prime[j]==0)
u[i*prime[j]]=0;
else
u[i*prime[j]]=-u[i];
v[i*prime[j]]=prime[j];
}
}
for(int i=1;i<=n;i++)s[i]=s[i-1]+u[i];
}
int g(int a,int b)
{
return a/(a/b);
}
int main()
{
primes(N-1);
int t;
cin>>t;
while(t--)
{
int a,b,d;
cin>>a>>b>>d;
a/=d,b/=d;
int n=min(a,b);
res=0;
for(int l=1,r;l<=n;l=r+1)
{
r=min(g(a,l),g(b,l));
res+=1ll*(s[r]-s[l-1])*(a/l)*(b/l);
}
cout<<res<<'\n';
}
return 0;
}
