求值（多组数据）
$_math_inline$ \sum_{i=x}^{n}\sum_{j=y}^{m}[\gcd(i,j)=k]\qquad (1\leqslant T,x,y,n,m,k\leqslant 5\times 10^4) $math_inline_$根据容斥原理，原式可以分成 $_math_inline$4$math_inline_$ 块来处理，每一块的式子都为
$_math_inline$ \sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=k] $math_inline_$考虑化简该式子
$_math_inline$ \sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)=1] $math_inline_$因为 $_math_inline$\gcd(i,j)=1$math_inline_$ 时对答案才用贡献，于是我们可以将其替换为 $_math_inline$\varepsilon(\gcd(i,j))$math_inline_$ （ $_math_inline$\varepsilon(n)$math_inline_$ 当且仅当 $_math_inline$n=1$math_inline_$ 时值为 $_math_inline$1$math_inline_$ 否则为 $_math_inline$0$math_inline_$ ），故原式化为
$_math_inline$ \sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\varepsilon(\gcd(i,j)) $math_inline_$将 $_math_inline$\varepsilon$math_inline_$ 函数展开得到
$_math_inline$ \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) $math_inline_$变换求和顺序，先枚举 $_math_inline$d\mid gcd(i,j)$math_inline_$ 可得
$_math_inline$ \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 $math_inline_$（其中 $_math_inline$d\mid i$math_inline_$ 表示 $_math_inline$i$math_inline_$ 是 $_math_inline$d$math_inline_$ 的倍数时对答案有 $_math_inline$1$math_inline_$ 的贡献） 易知 $_math_inline$1\sim\lfloor\dfrac{n}{k}\rfloor$math_inline_$ 中 $_math_inline$d$math_inline_$ 的倍数有 $_math_inline$\lfloor\dfrac{n}{kd}\rfloor$math_inline_$ 个，故原式化为
$_math_inline$ \displaystyle\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d) \lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor $math_inline_$很显然，式子可以数论分块求解（注意：过程中默认 $_math_inline$n\leqslant m$math_inline_$ ）。
时间复杂度 ： $_math_inline$\Theta(N+T\sqrt{n})$math_inline_$
Code
1#include <bits/stdc++.h>
2using namespace std;
3const int N = 50000;
4int mu[N + 5], p[N + 5];
5bool flg[N + 5];
6
7void init() {
8 int tot = 0;
9 mu[1] = 1;
10 for (int i = 2; i <= N; ++i) {
11 if (!flg[i]) {
12 p[++tot] = i;
13 mu[i] = -1;
14 }
15 for (int j = 1; j <= tot && i * p[j] <= N; ++j) {
16 flg[i * p[j]] = 1;
17 if (i % p[j] == 0) {
18 mu[i * p[j]] = 0;
19 break;
20 }
21 mu[i * p[j]] = -mu[i];
22 }
23 }
24 for (int i = 1; i <= N; ++i) mu[i] += mu[i - 1];
25}
26
27int solve(int n, int m) {
28 int res = 0;
29 for (int i = 1, j; i <= std::min(n, m); i = j + 1) {
30 j = std::min(n / (n / i), m / (m / i));
31 res += (mu[j] - mu[i - 1]) * (n / i) * (m / i);
32 }
33 return res;
34}
35
36int main() {
37 int T, a, b, c, d, k;
38 init();
39 cin >> T;
40 while (T--) {
41 cin >> a >> b >> c >> d >> k;
42 cout << solve(b / k, d / k)
43 - solve(b / k, (c - 1) / k)
44 - solve((a - 1) / k, d / k)
45 + solve((a - 1) / k, (c - 1) / k) << endl;
46 }
47 return 0;
48}
Code 2
1#include <bits/stdc++.h>
2
3using namespace std;
4typedef long long LL;
5
6namespace SOL {
7 const int MAXN = 50011;
8 const int MAXL = 50000;
9 LL a, b, c, d, k, ans;
10 int prime[MAXN], cnt;
11 bool ok[MAXN];
12 int mobius[MAXN], sum[MAXN]; //莫比乌斯函数及其前缀和
13
14 inline void init() { //递推得到莫比乌斯函数
15 //1的莫比乌斯函数是1；质数的莫比乌斯函数为1；含有相同质因子的数莫比乌斯函数为0；
16 //不含有相同质因子的数的莫比乌斯函数函数为（-1）^k，k为质因子个数
17 mobius[1] = 1;
18 for (int i = 2; i <= MAXL; i++) {
19 if (!ok[i]) {
20 mobius[i] = -1;
21 prime[++cnt] = i;
22 }
23 for (int j = 1; j <= cnt && i * prime[j] <= MAXL; j++) {
24 ok[i * prime[j]] = 1;//标记合数
25 if (i % prime[j]) mobius[i * prime[j]] = -mobius[i];//互质的两个数乘起来得到一个不含有相同质因子的数，质因子个数奇偶性改变，莫比乌斯函数变号
26 else {
27 mobius[i * prime[j]] = 0;
28 break;
29 }//留到后面再筛，此处已经可以break
30 }
31 }
32 for (int i = 1; i <= MAXL; i++) sum[i] = sum[i - 1] + mobius[i];//求莫比乌斯函数的前缀和
33 }
34
35 inline LL solve(LL n, LL m) {//计算a在[1,n]且b在[1,m]中的gcd(a,b)==1的数目
36 n /= k;
37 m /= k;
38 if (n > m) swap(n, m);
39 if (n == 0) return 0;
40 LL i, next1, next2, next;//把相等的部分直接分块一起计算
41 LL tot = 0;
42 for (i = 1; i <= n; i = next) {
43 next1 = n / (n / i);
44 next2 = m / (m / i);
45 next = min(next1, next2);
46 tot += (n / i) * (m / i) * (sum[next] - sum[i - 1]);
47 next++;
48 }
49 return tot;
50 }
51
52 inline void work() {
53 int T;
54 cin >> T;
55 init();
56 while (T--) {
57 cin >> a >> b >> c >> d >> k;
58 ans = solve(b, d) - solve(a - 1, d) - solve(b, c - 1) + solve(a - 1, c - 1);//容斥原理
59 cout << ans << endl;
60 }
61 }
62}
63
64int main() {
65 SOL::work();
66 return 0;
67}
除另有声明外，本博客文章均采用 知识共享 (Creative Commons) 署名 4.0 国际许可协议 进行许可。转载请注明原作者与文章出处。