数对统计
题目链接:ybt金牌导航8-7-1
题目大意
给你 n,m,求 gcd(x,y)=1 的数对个数。
1<=x<=n,1<=y<=m
思路
莫比乌斯函数
什么东西
首先我们要知道莫比乌斯函数是个什么鬼东西。
首先,我们先不管莫比乌斯函数,先来看一个函数:(F(n)=sum_{d|n}f(d))
(假设 (f(d)) 是一个给出的函数)
那根据定义,我们可以先弄一下:
(F(1)=f(1))
(F(2)=f(1)+f(2))
(F(3)=f(1)+f(3))
(F(4)=f(1)+f(2)+f(4))
(F(5)=f(1)+f(5))
(F(6)=f(1)+f(2)+f(3)+f(6))
(F(7)=f(1)+f(7))
(F(8)=f(1)+f(2)+f(4)+f(8))
那我们考虑用 (F(n)) 来推 (f(n))
(f(1)=F(1))
(f(2)=F(2)-F(1))
(f(3)=F(3)-F(1))
(f(4)=F(4)-F(2))
(f(5)=F(5)-F(1))
(f(6)=F(6)-F(3)-F(2)+F(1))
(f(7)=F(7)-F(1))
(f(8)=F(8)-F(4))
那我们会发现,它只会由它因数的 (F) 值加减或者不要组成。
那我们可以把它弄成这样的形式:
(f(n)=sum_{d|n}mu(d)F(frac{n}{d}))
那我们的莫比乌斯函数 (mu(d)) 就出现了!
定义
- 如果 (d=1),那 (mu(d)=1)
- 如果 (d=p_1p_2...p_k),(p_i) 是互不相同的素数,那就会有 (mu(d)=(-1)^k)
- 如果不满足上面两个条件,那 (mu(d)=0)
如何求
首先,我们可以很明显的看出用定义法求会比较慢,尤其是要求一个区间的。
我们可以考虑用类似 DP 的方法求。
看到跟素数有关,自然想到先用欧拉筛求。
然后我们可以考虑,在欧拉筛枚举最小质因子的时候,我们可以想到你处理 (i imes prime_j) 这个数。
那我们想 (prime_j) 的 (mu) 值已经求出,那我们可以看这个最小的质因子是否已经是 (prime_j) 的因子(就是能否整除)。
如果能整除,那就说明这个质数有两个,那就直接 (mu_{i imes prime_j}=0),那否则就是多一个素数,就是 (mu_{i imes prime_j}=mu_{i} imes-1)。
你会想,啊,如果原来 (prime_j) 已经有一个素数能分解出两个或以上呢?
那因为这样,它的 (mu) 值就一定是 (0),就算乘了 (-1),也还是 (0),就没有问题了。
一些性质
- (sum_{d|n}mu(d)=left{egin{matrix} 1 n=1\ 0 n>1 end{matrix} ight.)
- (varphi(n)=sum_{d|n}dfrac{mu(d) imes n}{d})
关于这道题
我们考虑利用第一条性质:
(ans=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=1])
(egin{aligned} ans &= sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=1]\ &= sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sumlimits_{d|gcd(i,j)}mu(d)\ &= sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sumlimits_{d|i,d|j}mu(d)\ &= sumlimits_{d}sumlimits_{i=1,d|i}^{n}sumlimits_{j=1,d|j}^{m}mu(d)\ &= sumlimits_{d}mu(d)cdotleft lfloor frac{n}{d} ight floorcdotleft lfloor frac{m}{d} ight floor end{aligned})
然后这就是 (O(n)) 的,但是因为它是多次询问,每个询问都是 (O(n)),就还是过不了。
然后看到向下取整,自然想到整除分块。
那就前缀和 (mu) 函数,然后就整除分块处理就好了。
代码
#include<cstdio>
#include<iostream>
#define ll long long
using namespace std;
ll T;
ll n, m;
ll miu[100001], prime[100001];
ll qz[100001];
bool np[100001];
void get_miu() {//同欧拉筛预处理 μ
miu[1] = 1;
for (int i = 2; i <= 100000; i++) {
if (!np[i]) {
prime[++prime[0]] = i;
miu[i] = -1;
}
for (int j = 1; j <= prime[0] && 1ll * i * prime[j] <= 100000ll; j++) {
np[i * prime[j]] = 1;
if (i % prime[j] == 0) {
miu[i * prime[j]] = 0;
break;
}
else miu[i * prime[j]] = miu[i] * -1;
}
}
}
void get_qz() {//前缀和
for (int i = 1; i <= 100000; i++)
qz[i] = qz[i - 1] + miu[i];
}
void work(int n, int m) {
ll ans = 0;
for (int l = 1; l <= n && l <= m; ) {
int r = min(n / (n / l), m / (m / l));//数论分块加速
ans += (qz[r] - qz[l - 1]) * (n / l) * (m / l);
l = r + 1;
}
printf("%lld
", ans);
}
int main() {
get_miu();
get_qz();
scanf("%lld", &T);
while (T--) {
scanf("%lld %lld", &n, &m);
work(n, m);
}
return 0;
}