Description
对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。
Input
第一行一个整数n，接下来n行每行五个整数，分别表示a、b、c、d、k
Output
共n行，每行一个整数表示满足要求的数对(x,y)的个数
Sample Input
2
2 5 1 5 1
1 5 1 5 2
2 5 1 5 1
1 5 1 5 2
Sample Output
14
3
3
最近在学习莫比乌斯反演，终于看到一道裸题可以让我A掉了
询问 a≤x≤b，c≤y≤d，且gcd(x,y) = k 的个数
即$F(a, b, c, d) = \sum_a^b{x}\sum_c^d{y}[gcd(x, y) == k]$
很明显我们可以用容斥转换成
$F(a, b, c, d)= F(1, b, 1, d) - F(1, b, 1, c - 1) - F(1, a - 1, 1, d) + F(1, a - 1, 1, c - 1)$
所以只需要求出:
$f(k, b, d) = \sum_1^b{x}\sum_1^d{y}[gcd(x, y) == k]$
通过莫比乌斯反演:
$f(k, b, d) = \sum_ {k | t} g(t)μ(t / k)$
所以
$f(k, b, d) = \sum_ {k | x}μ(\frac{x}{k}){[\frac{b}{x}]}[\frac{d}{x}]$
$ = \sum_ {i = 1}^{[\frac{n}{k}]}μ(i){[\frac{\frac{b}{k}}{i}]}[\frac{\frac{d}{k}}{i}]$
然后预处理μ(n)的前缀和，后面可以$\sqrt{n}$ 算出
总复杂度 $n\sqrt{n}$
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 #include <algorithm> 6 7 #define LL long long 8 9 using namespace std; 10 11 const int MAXN = 5e4 + 10; 12 int N; 13 int t1, t2; 14 int a, b, k, c, d; 15 int cnt = 0; 16 bool f[MAXN]; 17 LL sm[MAXN]; 18 int prime[MAXN], mul[MAXN]; 19 inline LL read() 20 { 21 LL x = 0, w = 1; char ch = 0; 22 while(ch < ‘0‘ || ch > ‘9‘) { 23 if(ch == ‘-‘) { 24 w = -1; 25 } 26 ch = getchar(); 27 } 28 while(ch >= ‘0‘ && ch <= ‘9‘) { 29 x = x * 10 + ch - ‘0‘; 30 ch = getchar(); 31 } 32 return x * w; 33 } 34 35 void init() 36 { 37 mul[1] = 1; 38 for(int i = 2; i <= 50000; i++) { 39 if(!f[i]) { 40 prime[++cnt] = i; 41 mul[i] = -1; 42 } 43 for(int j = 1; j <= cnt && prime[j] * i <= 50000; j++) { 44 f[prime[j] * i] = 1; 45 if(i % prime[j] == 0) { 46 mul[prime[j] * i] = 0; 47 break; 48 } 49 mul[prime[j] * i] = mul[i] * (-1); 50 } 51 } 52 for(int i = 1; i <= 50000; i++) { 53 sm[i] = sm[i - 1] + mul[i]; 54 } 55 } 56 int main() 57 { 58 init(); 59 60 N = read(); 61 while(N--) { 62 long long ans = 0; 63 a = read(), b = read(), c = read(), d = read(), k = read(); 64 t1 = b, t2 = d; 65 if(t1 > t2) 66 swap(t1, t2); 67 t1 = t1 / k, t2 = t2 / k; 68 for(int i = 1, j; i <= t1; i = j + 1) { 69 j = min(t1 / (t1 / i), t2 / (t2 / i)); 70 ans += (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i); 71 } 72 t1 = a - 1, t2 = d; 73 if(t1 > t2) 74 swap(t1, t2); 75 t1 = t1 / k, t2 = t2 / k; 76 for(int i = 1, j; i <= t1; i = j + 1) { 77 j = min(t1 / (t1 / i), t2 / (t2 / i)); 78 ans -= (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i); 79 } 80 t1 = b, t2 = c - 1; 81 if(t1 > t2) 82 swap(t1, t2); 83 t1 = t1 / k, t2 = t2 / k; 84 for(int i = 1, j; i <= t1; i = j + 1) { 85 j = min(t1 / (t1 / i), t2 / (t2 / i)); 86 ans -= (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i); 87 } 88 t1 = a - 1, t2 = c - 1; 89 if(t1 > t2) 90 swap(t1, t2); 91 t1 = t1 / k, t2 = t2 / k; 92 for(int i = 1, j; i <= t1; i = j + 1) { 93 j = min(t1 / (t1 / i), t2 / (t2 / i)); 94 ans += (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i); 95 } 96 printf("%lld\n", ans); 97 } 98 return 0; 99 }