BZOJ2301 [HAOI2001]Problem b

题目链接

求值:

$\displaystyle\sum_{i=x}^{n}\sum_{j=y}^{m}[\gcd(i,j)=k]\qquad(1\leqslant T,x,y,n,m,k\leqslant 50000)$

分析

由容斥原理可将原式分成4块处理,每一块式子为

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=k]$

化简得

$\displaystyle\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)=1]$

因为当$\gcd(i,j)=1$时对答案有贡献,所以可以替换为$\varepsilon(\gcd(i,j))$,所以原式被化简为

$\displaystyle\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\varepsilon(\gcd(i,j))$

将$\varepsilon(\gcd(i,j))$用$\text{Dirichlet}$卷积展开得

$\displaystyle\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d\mid\gcd(i,j)}\mu(d)$

变换求和顺序得

$\displaystyle\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}d\mid i\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}d\mid j$

容易得知$\displaystyle 1\sim\lfloor\frac{n}{k}\rfloor$中有$d$的倍数$\displaystyle\lfloor\frac{n}{kd}\rfloor$个,所以原式被化简为

$\displaystyle\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor$

很明显,化简后的式子可以用数论分块解决(过程中默认$n\leqslant m$)

时间复杂度为$\displaystyle\Theta(N+T\sqrt{n})$

$Code$

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
#include <bits/stdc++.h>
#define reg register
#define MAXN 50010

using namespace std;

int mu[MAXN], p[MAXN];
bool flag[MAXN];

inline void init() { // 预处理莫比乌斯函数
int tot = 0;
mu[1] = 1;
for(reg int i = 2, j; i <= MAXN; ++i) {
for(!flag[i] ? p[++tot] = i, mu[i] = -1 : 1, j = 1; j <= tot && i * p[j] <= MAXN; ++j) {
flag[i * p[j]] = 1;
if(!(i % p[j])) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for(reg int i = 1; i <= MAXN; ++i) mu[i] += mu[i - 1];
}

int solve(int n, int m) { // 数论分块
int res = 0;
for(reg int i = 1, j; i <= min(n, m); i = j + 1) j = min(n / (n / i), m / (m / i)), res += (mu[j] - mu[i - 1]) * (n / i) * (m / i);
return res;
}

int main() {
int T, a, b, c, d, k;
init();
for(scanf("%d", &T); T; T--) {
scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
printf("%d\n", solve(b / k, d / k) - solve(b / k, (c - 1) / k) - solve((a - 1) / k, d / k) + solve((a - 1) / k, (c - 1) / k));
}
return 0;
}