https://darkbzoj.tk/problem/2301

题意：对于给出的 nn 个询问，每次求有多少个数对 (x,y)(x,y)，满足 a≤x≤b，c≤y≤da≤x≤b，c≤y≤d，且 gcd(x,y)=kgcd(x,y) = k。

容斥原理+莫比乌斯反演+整除分块

在 [a,b],[c,d][a,b],[c,d] 上的答案可以由容斥变为 [1,b],[1,d]−[1,a)[1,d]−[1,b][1,c)+[1,a)[1,c)[1,b],[1,d]-[1,a)[1,d]-[1,b][1,c)+[1,a)[1,c)

则现在要求 ∑i=1n∑j=1m[gcd(i,j)=k]\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=k]。

把 kk 除掉，变为 ∑i=1⌊nk⌋∑j=1⌊mk⌋[gcd(i,j)=1]\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[gcd(i,j)=1]

为方便起见，以下把 ⌊nk⌋\lfloor\frac{n}{k}\rfloor 设为 nn

∑i=1n∑j=1m[gcd(i,j)=1]=∑i=1n∑j=1m∑d∣gcd(i,j)μ(d)\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=1]\\ =\sum_{i=1}^n\sum_{j=1}^m\sum_{d|gcd(i,j)}\mu(d)\\

把对 dd 的 ∑\sum 提出到最外层

=∑d=1min(n,m)∑i=1n∑j=1mμ(d)⋅[d∣i⋀d∣j]=∑d=1min(n,m)μ(d)⋅⌊nd⌋⋅⌊md⌋=\sum_{d=1}^{min(n,m)}\sum_{i=1}^n\sum_{j=1}^m\mu(d)\cdot [d|i\bigwedge d|j]\\ =\sum_{d=1}^{min(n,m)}\mu(d)\cdot\lfloor\frac{n}{d}\rfloor\cdot\lfloor\frac{m}{d}\rfloor

对于 ⌊nd⌋\lfloor\frac{n}{d}\rfloor 可以进行整除分块，而在加了个 ⌊md⌋\lfloor\frac{m}{d}\rfloor 仍然可以，只不过在分块时注意右端点要在两个台阶中小的那个。这也说明知道了 ll 后可以直接通过 r=n/(n/l)r=n/(n/l) 得到 rr ，而不管 ll 是否一定恰好在交界处。

对 μ(d)\mu(d) 预处理出前缀和，这里直接用即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 2e5 + 10;
const int INF = 0x3f3f3f3f;
int mu[N], p[N], tot;
bool flg[N];
void getMu() {
	mu[1] = 1;
	for (int i = 2; i < N; ++i) {
		if (!flg[i]) p[++tot] = i, mu[i] = -1;
		for (int j = 1; j <= tot && i * p[j] < N; ++j) {
			flg[i * p[j]] = 1;
			if (i % p[j] == 0) {
				mu[i * p[j]] = 0;
				break;
			}
			mu[i * p[j]] = -mu[i];
		}
	}
	for (int i = 1; i < N; i++)mu[i] += mu[i - 1];
}
int T, a, b, c, d, k;
int solve(int n, int m) {
	int ans = 0;
	if (n > m)swap(n, m);
	n /= k, m /= k;
	for (int l = 1, r; l <= n; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans += (n / l)*(m / l)*(mu[r] - mu[l - 1]);
	}
	return ans;
}
int main() {
	getMu();
	scanf("%d", &T);
	while (T--) {
		scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
		printf("%d\n", solve(b, d) - solve(b, c - 1) - solve(a - 1, d) + solve(a - 1, c - 1));
	}
	return 0;
}