题目链接

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(n,min(g(a,l),g(b,l)));
			res+=1ll*(s[r]-s[l-1])*(a/l)*(b/l);
		}
		cout<<res<<'\n';
    } 
    return 0;
}