Description
对于给出的 (n) 个询问，每次求有多少个数对 ((x,y))，满足 (a le x le b)，(c le y le d)，且 (gcd(x,y) = k)，(gcd(x,y)) 函数为 (x) 和 (y) 的最大公约数。
其中 (1 le n,k le 5 imes 10^4)，(1 le a le b le 5 imes 10^4)，(1 le c le d le 5 imes 10^4)。
Solution
首先肯定是自然而然地想到分成四块容斥，即
[egin{split}
sum limits_{i=a}^{b} sum limits_{j=c}^{d} left[ gcd(i,j) = k
ight] &= sum limits_{i=1}^{b} sum limits_{j=1}^{d} left[ gcd(i,j) =k
ight] \
&- sum limits_{i=1}^{a} sum limits_{j=1}^{d} left[ gcd(i,j) =k
ight]\
&- sum limits_{i=1}^{b} sum limits_{j=1}^{c} left[ gcd(i,j) =k
ight]\
&+ sum limits_{i=1}^{a} sum limits_{j=1}^{c} left[ gcd(i,j) =k
ight]
end{split}
]
这样，每一部分都是 (sum limits_{i=1}^{m} sum limits_{j=1}^{m} left[ gcd(i,j) = k ight]) 的形式。
接下来考虑如何将式子进行变换：
[sum limits_{i=1}^{m} sum limits_{j=1}^{m} left[ gcd(i,j) = k
ight] = sum limits_{i=1}^{lfloor frac{n}{k}
floor} sum limits_{j=1}^{lfloor frac{m}{k}
floor} left[ gcd(i,j) = 1
ight]
]
看到右边的 (left[ gcd(i,j) = 1
ight]) 是不是感觉很熟悉？根据在莫比乌斯反演学习笔记中提到的反演结论：
[left[ gcd(i,j) = 1
ight] iff sum limits_{d mid gcd(i,j)} mu(d)
]
可以继续化简：
[sum limits_{i=1}^{lfloor frac{n}{k}
floor} sum limits_{j=1}^{lfloor frac{m}{k}
floor} sum limits_{d mid gcd(i,j)} mu(d)
]
考虑变换求和顺序：
[sum limits_{d=1}^{min(n,m)} mu(d) sum limits_{i=1}^{lfloor frac{n}{k}
floor} left[ d mid i
ight] sum limits_{j=1}^{lfloor frac{m}{k}
floor} left[ d mid j
ight]
]
易知在 (left[ 1, n ight]) 中约数包含 (d) 的数有 (lfloor frac{n}{d} floor) 个，于是式子可以化简为：
[sum limits_{d=1}^{min(n,m)} mu(d) lfloor frac{lfloor frac{n}{k}
floor}{d}
floor lfloor frac{lfloor frac{m}{k}
floor}{d}
floor
]
再由莫比乌斯反演学习笔记中提到的引理1，可得：
[sum limits_{d=1}^{min(n,m)} mu(d) lfloor frac{n}{kd}
floor lfloor frac{m}{kd}
floor
]
所以最终，就得到了
[sum limits_{i=a}^{b} sum limits_{j=c}^{d} left[ gcd(i,j) = k
ight] = sum limits_{d=1}^{min(n,m)} mu(d) lfloor frac{n}{kd}
floor lfloor frac{m}{kd}
floor
]
这个式子就可以利用数论分块进行求解了。
Code
Talk is cheap. Show me the code.
#include <bits/stdc++.h>
using namespace std;
int ty() {
int x = 0, f = 1; char ch = getchar();
while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
return x * f;
}
const int _ = 5e4 + 10;
int tot, vis[_], p[_], mu[_], sum[_];
void getMu(int lim = 5e4) {
mu[1] = 1;
for (int i = 2; i <= lim; ++i) {
if (!vis[i]) p[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot && p[j] * i <= lim; ++j) {
vis[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 <= lim; ++i) sum[i] = sum[i - 1] + mu[i];
}
int calc(int n, int m, int k) {
n /= k, m /= k;
int ret = 0;
for (int l = 1, r; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
ret += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
}
return ret;
}
int main() {
getMu();
int T = ty(), a, b, c, d, k;
while (T--) {
a = ty(), b = ty(), c = ty(), d = ty(), k = ty();
printf("%d
", calc(b, d, k) - calc(a - 1, d, k) - calc(b, c - 1, k) +
calc(a - 1, c - 1, k));
}
return 0;
}